GEOG 566

         Advanced spatial statistics and GIScience

April 26, 2018

Using multivariate statistics to test spatial variability of model performance.

Filed under: Exercise/Tutorial 1 2018 @ 2:42 pm

Question that you asked

Understanding where my model performs well compared to in-situ observations is important for assessing the appropriateness of snow and climate metrics derived from it in comparison to Dall sheep demographic data. For instance, is the model adequately capturing the paths of wind-blown snow that are so critical for winter grazing? Are there arrangements of landscape whose snow evolution is consistently difficult to represent accurately? Answering these questions produces dual benefit; identification of space to improve the model and the possibility to characterise and handle uncertainties in subsequent analyses using model data.


The data I’m using in this analysis are measurements of snow depth taken at a field site in the Wrangell St Elias National Park between the 18thand 24thMarch 2017. These snow depth measurements were taken at relatively high spatial resolution using a GPS enabled Magnaprobe and have been therefore aggregated into 30m by 30m grid cells of mean snow depth, along with other simple statistics (e.g. standard deviation). These grid cells match the output rasters of my spatially distributed snow evolution model enabling comparison of model performance by comparison of observed mean snow depth to modelled snow depth. For simplicity’s sake, I have chosen the middle date of the observation period, the 21stMarch, for the modelled observations. This is done in the knowledge that there was no snowfall or high winds in the entire observation period, both in reality and model space, and that snowpack evolution governed by other processes is relatively minor under such conditions and time-period.  Important to note also is that the model output compared to is the product of >250 iterations of model calibration. The selected model after this process produces the lowest Root Mean Squared Error (RMSE) between model and observed snow depth and water equivalent. This analysis is hence designed to examine the spatial distribution of this error. To do this I have categorised the result from subtracting the mean observed snow depth and the modelled observed snow depth as either being weak or strong, over or under prediction based on whether it is within or over the RMSE, positively or negatively.

The landcover dataset, derived from the National Land Cover Database 2011 Alaska (Homer et al. 2015), alongside a DEM for the model domain provide landscape data.

Fig 1; Spiral of Magnaprobe measurements on Jaeger Mesa

Name of the tool or approach that you used.

To produce landscape metrics I used the following data processing tools in QGIS;

  • GRASS GIS 7 r.slope and r.aspect,
  • Raster Calculator to produce a ‘northerness’ layer – cosine(aspect)

I also used QGIS to aggregate the point snow depth measurements into 30m resolution rasters for mean depth, standard deviation, max/min depth.

To perform the multivariate analysis I used the R package FactoMineR (Lê, S., Josse, J. & Husson, F. 2008). FactoMinerR is a Multivariate Exploratory Data Analysis and Data Mining Tool that enables Multiple Factor Analysis to deal with discrete data such as the nominal landcover data in my analysis. The particular package used is FAMD – Factor Analysis for Mixed Data.

Brief description of steps you followed to complete the analysis.

Data preparation;

For the aggregation of the in-situ measurements I used QGIS tools to create a 30m polygon grid that matched the extent and aligned to my SnowModel domain. I then used the ‘Join Attributes by Location’ tool to intersect the point data to these polygons, selecting mean, min, max, and standard deviation statistics for summary. The Rasterise tool then could be used to convert the polygon layer into ArcInfo Ascii format.

To compare model results vs observations, I first subtracted the model result raster for the 21stMarch 2017 from the mean observed snow depth raster. The resulting raster showing the difference could then be mapped on top of the landscape variables such as elevation, slope etc. (see figs 2 to 3 below) in QGIS.

Fig. 2; Model performance of snow depth mapped against elevation

Fig. 3; Model performance of snow depth mapped against landcover class

Matlab was employed to build the table for Multiple Factor Analysis in R. Using the arcgridread function it was possible to open the arrays of observations, model results and landcover variables, and then build a table holding each sampled cell’s corresponding variables. This table was then converted to .csv format for use in R.

Within R the table was imported as a data frame and the land cover vector converted into factor type.

Data Analysis;

The FAMD package of FactoMineR was used to perform the multivariate analysis, producing tables and plots of the results.

Brief description of results you obtained.

Fig. 4; Screeplot of explained variance by dimension

Table 1; Eigenvalues, percent variance and cumulative variance by dimension

The FAMD analysis allows insight into which dimensions describe the greatest variability, effectively ‘flattening’ the data, reducing a large set of variables to a selection that contains most of the information. In Fig 4 and Table 1 above, we can see that for my data 24 % of the variance is explained by the 1stdimension, 16 % by the 2nddimension. Cumulatively, by the 5thdimension 70% of the variance is described.

Fig. 5; Variable contribution to dimension 1

It is then possible to see which variables most contribute to each variable. Figure 5 shows that elevation and landcover both contribute >20 % of the variance seen in the 1stdimension. Above the 95% significance level, the Observed Mean Snow Depth (obsMeanSNOD) is also seen as a contributor.

Fig. 6; Variable contribution to dimension 2

The second dimension has the greatest contribution from the categorical values of the model performance (diffCategory), and significant contributions by both modelMeanSNOD and landcover, see figure 6.

Fig. 7; Quantitive variables on dimensions 1 and 2 coloured by the cosine of their contribution.

It is possible to then plot both the quantitive, figure 7, and qualitative, figure 8, variables on both the 1stand 2nddimension axis. By doing so some insight into how the variables interact is possible. The high positive standing on the 1stdimension and 2nddimensions by Observed Mean Snow Depth (fig. 7, obsMeanSNOD) is matched by the qualitative variables of coniferous forest and strong under predict. This suggests that the model strongly under predicts when the land cover is coniferous forest and the observed snow depth is high. This is confirmed by the greatest influence of elevation on the 1stdimension, higher elevations in my study domain have lower snow depth and bare or prostrate shrub tundra land cover. Indeed, my model is specifically calibrated to best reproduce snow in these areas as they’re prime Dall Sheep habitat. Looking further at figure 8 it is possible to see that weak underpredict matches the bare land cover class quite well, and similarly weak overpredict lines up nicely with prostrate shrub tundra. On the 2nddimension it is possible to see that erect shrub tundra is in the direction of the Model Mean Snow Depth and strong over-predict. This would suggest that for this land cover the model is producing too much snow.

Fig. 8. Qualitative variables contribution to dimensions 1 and 2

Figures 9 and 10 further describe the patterns explored above by plotting the individual rows, i.e. the pixels where I’ve sampled, on the first 2 dimensions and colouring them by their category for both the categorical description of model performance and land cover. Here however, we see that the inferences aren’t necessarily as clear cut. For example, the erect shrub tundra land cover class does not entirely comprise of strong model overprediction.

Fig. 9; Individuals by model performance category

Fig. 10; Individuals by land cover category

Critique of the method – what was useful, what was not?

This task was useful for quite a number of reasons, not least because it forced me into baby steps use of R. It neatly confirmed my suspicions about where my model is over predicting, high elevations with prostrate shrub tundra, but also identified other areas of poor representation, namely strong under prediction coniferous forest. It too gave me a certain satisfaction that for the majority of pixels in high elevation, low vegetation Dall Sheep terrain, there’s only weak over or under prediction, suggesting my calibration efforts haven’t been without cause.

It is also interesting to observe that the bare land cover class is subject to weak under prediction, whereas the prostrate shrub tundra class tends to be over predicted. This is quite likely due to with how the model treats the roughness of the land cover compared to how rough that land cover is in reality (as well as whether the 30m NLCD land cover classes are actually a good representation of ground conditions). In my study area, the bare patches were more frequently pretty rough scree patches that had the capacity to intercept a fair amount of snow. Conversely the prostrate shrub tundra snow holding height in the model parameters is set at 10 cm, whereas the reality of the land cover is these areas was a patchwork of sparse alpine grasses and sedges, probably not over 5 cm in height. This is useful information that allows me to play with adjusting these parameters.

A surprise is that the other landscape variables, slope, aspect, northerness, have little influence on the variability. I had expected that slope may have had more of an influence given how important it is to wind redistribution of snow

Where I have queries with the method in respect to my particular problem is the selection, type (qualitative vs quantitative) and quantity of the variables included. For instance, I originally left the description of model performance as a continuous variable but found the results not nearly as easily interpretable as when I categorised them. There is also an element of variation within variation that is masked by attempting to flatten the data. Dominant variables take precedence, such as land cover and elevation in my instance, so that the causes of variation within land covers, for example, are hard to see. A further exercise useful to my specific problem would be to just run the analysis on pixels from the bare and prostrate shrub tundra classes and seeing whether landscape variables are greater contributors.



Homer, C.G., Dewitz, J.A., Yang, L., Jin, S., Danielson, P., Xian, G., Coulston, J., Herold, N.D., Wickham, J.D., and Megown, K., 2015, Completion of the 2011 National Land Cover Database for the conterminous United States-Representing a decade of land cover change information. Photogrammetric Engineering and Remote Sensing, v. 81, no. 5, p. 345-354

Lê, S., Josse, J. & Husson, F. (2008). FactoMineR: An R Package for Multivariate AnalysisJournal of Statistical Software25(1). pp. 1-18.

Print Friendly, PDF & Email

No Comments

No comments yet.

RSS feed for comments on this post.

Sorry, the comment form is closed at this time.

© 2019 GEOG 566   Powered by WordPress MU    Hosted by