My analysis up to this point has used spatial variation in climate normals at 1600-m grid points. Time series analysis on the full 1600-m grid is impractical because each point would have a set of climate observations for each year over the period of the time series. This creates a lot of data, and a lot of redundant variation. An alternate approach is to use the BEC variants to define the ecologically appropriate scale of climate differentiation, and then use the climatic centroid of each BEC variant for time series analysis. The easiest way to do this is to find the 1600-m grid point whose climatic attributes are closest to the actual climatic centroid of the variant. I call this spatial location the climatic “centroid surrogate” for the BEC variant.

Mahalanobis distance vs Standardized Euclidian Distance

To find the “closest” grid point to the climatic centroid, one first has to select a distance measure. Simple Euclidian distance is not appropriate since climate elements have different units (e.g. ^oC vs. mm) and variables with large variance will dominate the distance measurement. The most straightforward approach is to standardize all the variables by subtracting the mean from each observation and dividing by the standard deviation. This converts all variables to units of standard deviation and a mean of zero, and facilitates a balanced measure of Euclidean distance. Another approach is to use Mahalanobis distance which also accounts for the correlation between variables. however, in practice I found mahalanobis distance to be difficult to implement in climate space due to the high collinearity between variables. For the purpose of finding a climatic centroid surrogate, standardized Euclidian distance and mahalanobis distance are likely to be equivalent. For this reason, I used standardized Euclidian distance.

Visual validation in climate space and map space

Figure 1: PCA climate space for the BC coast. BEC variants are labeled at their 2D pseudocentroids in black and full-dimensional (12D) centroid surrogates in blue. The close proximity of the labels provides visual validation that the centroid surrogate is at the climatic centroid.

Figure 2: locations of the climatic centroid surrogate for each coastal BEC variant.

The centroid surrogates are located very close to the BEC variant pseudocentroids of a two-dimensional PCA climate space representing 91% of the variance in the data (Figure 1). This strongly suggests that the calculation was correct and the centroid surrogate is the grid point with minimum distance to the climatic centroid. A map of the centroid surrogates (Figure 2) indicates that they are generally located in representative locations. However, there are several idiosyncrasies:

1. The centroid surrogate for the MHmm2 is located near Terrace in the northern portion of the variant’s geographic range, which extends all the way down to the U.S. border. This is problematic, since historical climate variability is likely to be different in the northern and southern areas of the MHmm2. Climate change forecasts are certainly very different for these two regions. This result suggests that more than one centroid surrogate may be required for variants with large geographic ranges (e.g. CWHvm1, vm2, MHmm1).
2. The centroid surrogate for the CWHwh1 is located in the mapped range of the CWHvh2. This is simply the result of using the BEC7 variant map for classification of the 1600-m grid points and the BEC8 map for mapping. This error will be corrected.
   

   Figure 3: close-up of the centroid surrogate for the CWHdm. It is located on the upper elevation limit of the CWHdm on East Redonda Island. While this location may represent the average climate of the CWHdm, it is not located in a core area of the variant. This result suggests the centroid surrogate selection method needs refinement.

3. The centroid surrogate for the CWHdm is located at the upper elevation limit of this variant’s occurrence on East Redonda Island (Figure 3). Intuitively, centroid surrogates should be located at approximately the geographic core of the BEC variant, rather than in transition zones. I attempted to prevent this problem by restricting centroid surrogate selection to grid points that are within one-half of a standard deviation of the mean elevation of each BEC variant. However, the same grid point was selected for the CWHdm, and resulted in poorer centroid surrogates for for submontane variants. As a result, the original solution shown in Figure 2 will be used for now. the potential for this issue should be kept in mind for other regions of the province as well.

Conclusion

The selection of centroid surrogates highlighted the important point that climatic zones that are considered now to represent relatively uniform climatic conditions may not have divergent historical and future variability. While centroid surrogates are useful for examining climatic differentiation in current conditions, they are not sufficient by themselves for characterizing climate change in BEC variants. Geographically extensive BEC variants likely need to have several centroid surrogates representing different regions of their distribution.

Visualization of historical and projected future climate change at Blind Channel, BC.

Summary

This blog post is a test run of a linear discriminant analysis (LDA) method for climate change visualization. The objective of this analysis is to create a low-dimensional (2D or 3D) climate space that reasonably represents the geographic variation of climate on the BC coast, and use that climate space to contextualize historical climatic variability and projected climate changes at selected locations. The analysis produced following outcomes:

• LDA appears to have produced a meaningful reduced climate space that reasonably represents the relationships between BEC variants. The regional scale appears to be more appropriate than the local (transect) or provincial scale;
• The backdrop of regional climatic variation provided a useful context for historical interannual climate variability and projected climate changes. However, projections into regionally novel climate space were difficult to evaluate due to absence of potential climatic analogs from outside the region;
• There appear to be dramatic sub-regional differences in the scale and novelty of projected climate change. Climate change on the north coast appears to be relatively moderate, and stays within the existing range of regional variation. In contrast, climate change on the south coast appears to be of much greater magnitude and entirely into regionally novel conditions.

The key next steps are to develop and test-run other alternative dimension-reduction methods such as multidimensional scaling, and to develop methods to visualize distant climatic analogs and truly novel climate space.

Introduction

In a previous blog post I discussed the results of a discriminant analysis of BEC variants along an elevational transect in central interior of BC. My conclusion was that only the first discriminant axis was meaningful. The other axes represented non-linearities in the relationships between the climatic variables, and did not have an intuitive climatic interpretation. As a result, discriminant analysis, though useful for classification, is likely of dubious value as a method for creating low-dimensional climate visualizations in a simple climatic gradient. Nevertheless, I was curious how discriminant analysis performs on a larger and more complex climatic landscape.  It is possible that meaningful variation in the second and third dimensions will relegate problematic non-linearities to lesser dimensions.

The purpose of this blog post is to test run an LDA-driven climate change visualization at a regional scale. I do a discriminant analysis on the 1961-1990 climate normals over a 1600-metre grid of coastal British Columbia. The results are evaluated in 2D and 3D climate space. Then I plot historical interannual variability and projected climate change for individual locations to get a glimpse of what climate change visualizations might look like in the reduced climate space.

Biogeoclimatic Climate Classification for the Coast

The coastal region provides a good case study because it has a relatively simple climatic geography, compared to the interior of BC. The climatic variation of the coast can be understood in a simple sense in terms of four climatic gradients:

• dry in the south vs. wet in the north, due to the influence of the Aleutian low;
• wet belts and rainshadows associated with the mountains of Vancouver island, the Olympic Peninsula, Haida Gwaii, and the Coast Range;
• hypermaritime (low seasonality of temperature) on the outer coast vs. submaritime (hot summers, cold winters) inland, due to the moderating influence of the Pacific Ocean; and
• Decreasing annual temperatures with elevation.

Figure 1: An example of the BEC climate classification: BEC variants of Vancouver BC.

The Biogeoclimatic Ecosystem Classification (BEC) includes a detailed classification of the climates of the coast.  The naming conventions of the BEC climate classification of are simple to understand and worth knowing. The vast majority of the coast is classified as the Coastal Western Hemlock Biogeoclimatic Zone. The exceptions are the Coastal Douglas-Fir (CDF) zone in the mediterranean climate of the Georgia Basin, and the Mountain Hemlock (MH) zone in subalpine elevations. Subzones are named with a two-letter code corresponding to their relative precipitation (x=very dry, d=dry, m=moist, w=wet, v=very wet) and continentality (h=hypermaritime, m=maritime, s=submaritime). Each subzone generally has at least two variants. In the CWH, variants are either elevational (1=submontane, 2=montane; CWHmm/vm/wh/ws) or latitudinal (1=southern, 2=central; CWHds/ms/vh). In the MH zone, variants are windward (1) vs. leeward (2). Figure 1 shows a small portion of the BEC variant mapping of the coast to provide a sense of the scale of the BEC climate classification.

Methods

Observations are locations on a 1600-m grid clipped to the geographic range of the CWH, CDF, and MH zones (MHmm1/2 only). I used the following annual climate normals output from ClimateWNA: MAT, MWMT, MCMT, TD, MAP, MSP, AHM, SHM, DD_0, DD5, FFP, PAS, Eref, and CMD. I log-transformed precipitation and heat-moisture variables (MAP, MSP, AHM, SHM, PAS) to normalize their distributions. LDA was done on these data using the “lda” call in the MASS package of R.

Results

Table 1: Correlations between climate variables and the first 5 of 13 linear discriminant functions. Primary correlations are highlighted grey, and the leading correlation is highlighted yellow. “% of trace” is an unrelated index of the percent of between-class variation explained by each discriminant function.

The structure of the LDA is strikingly different from that of the Quesnel-Barkerville transect. The first three discriminant functions (LD1-3) explain a moderate 82% of the between group variance of the coastal BEC variants (bottom of Table 1), which was the same amount explained by LD1 alone in the Quesnel-Barkerville transect. This is no doubt because the coastal region contains many climatic gradients, rather than just one. The correlations between the major discriminant axes (Table 1) are also more easily interpretable than in the transect analysis. LD1 appears to be an index of cold winters, and LD2 appears to be an index of growing season moisture availability (i.e. a negative drought index). LD3 appears to be a negative index of precipitation, especially of winter precipitation.

2D visualization of the coastal climate space

Figure 2: Simplified plot of the coastal grid in the first two discriminant axes, with BEC variants labelled at their climatic centroids.

Figure 3: ClimateBC normals for some reference locations

The relative positions of the BEC variants in the first two dimensions make sense in the context of these climatic indices (Figure 2). The progression from the bottom left to top right (CDFmm to MHmm1) is essentially the sequence of BEC variants moving outwards from the Georgia basin and upwards into the mountains of Vancouver Island and the coast range. The progression from the bottom right to the top left (CWHds1 to wh1) is the sequence of BEC variants from the coast-interior transition zone to the outer coast. The gap between the maritime and submaritime BEC variants in the first dimension is likely due to the presence of the coast mountains; i.e. the narrow glacial valleys connecting these two groups do not contribute enough grid points to fill the transition.  The only pairing of BEC variants that doesn’t make sense in the first two dimensions is the CWHwm (sea level) and MHmm1 (subalpine). Overall, this appears to be a valid low-dimensional representation of the coastal climate space.

Figure 4: Grid points in the third discriminating axis plotted against the first and second axes.

The third dimension of the coastal discriminant climate space is shown two-dimensionally in Figure 4, and below as an interactive three-dimensional visualization. There is a pronounced three-pointed arch shape in the third dimension. However, it is not clear whether this is a cryptic artifact of non-linearities in the data, or a meaningful climatic phenomenon. The correlation structure of the LDA (Table 1) suggests that the third axis is an index of winter dryness. The plot of the second and third dimensions (Figure 4, right) clearly distinguishes the rainshadows of the Olympic Peninsula/Vancouver Island (CDFmm), the Coast Range (CWHds2/ws1), and Haida Gwaii (CWHwh1). Interestingly, the CWHvh1 and vh2 also have intermediate values in the third dimension. This result makes sense, despite the common misconception that the hypermaritime subzones experience the highest rainfall on the coast; large areas of the CWHvh1 and vh2 are outer coastal lowlands that experience more summer fog than adjacent maritime subzones (e.g. CWHvm1) but less frontal precipitation due to the lack of orographic uplift. The CWHwm and MHmm1, which were poorly discriminated in the first two dimensions, are well discriminated in the third dimension. Despite the consistency of these observations with patterns of winter precipitation, the arch-shaped relationship between the first and third dimensions is questionable. It would take further investigation to determine whether this relationship is cryptic or meaningful.  For this reason, I use only the first two dimensions for the remainder of the analysis.

Interactive 3D visualization (drag and zoom to explore the climate space)

Distribution and overlap of BEC variants within the climate space

Table 2: Classification error rates of the discriminant analysis, indicating the number of grid cells that fell within another BEC variant’s domain in full-dimensional climate space. Prior probabilities are simply the proportion of the total number of 1600-m grid cells in each BEC variant.

Figure 6 shows the distribution of grid point normals for three arbitrarily selected BEC variants. The grid points of each BEC variant are disperse widely through the first two dimensions, and overlap with numberous BEC variant centroids. Superficially, this suggests low accuracy of the BEC variant mapping and/or of the ClimateWNA/PRISM data. However, the classification error of the discriminant analysis is 37% (Table 2). This indicates that almost 63% of the area of the BEC variants fall into their own exclusive region of the full-dimensional climate space. The remaining 37% classification error is a combination of errors in the BEC variant mapping, ClimateWNA/PRISM surfaces, and variable selection, in unknown proportions. These results indicate that ClimateWNA data can reasonably differentiate BEC variants in full-dimensional climate space, but that differentiation in low dimensional climate space is poor. The climatic centroid of BEC variants is likely much more meaningful in reduced dimensions than the boundaries between BEC variants.

Figure 5: Spatial variability of ClimateBC output for three BEC variants.

Visualizing historical and projected climate change in the 2D discriminant climate space

Figure 6: locations of the weather stations used for the climate change visualizations. The distribution of the Coastal BEC zones is shown in green.

The locations of three long-term weather stations were selected as case studies for climate change visualization (Figure 6). These sites roughly represent the northern, central, and southern regions of the BC coast, though Blind Channel is in an area that is climatically transitional between the southern and central regions. Although these are weather station locations, it is important to note that all date presented here is from ClimateWNA.

The graph for Stewart (Figure 7) can be used to introduce the elements of the climate change visualizations. The case study location is labeled at the position of the 1961-1990 normals in climate space. The climatic conditions of individual years in the period 1951-2010 at the case study location are shown as empty blue circles. A 30-year moving average of these years, shown as a blue line, shows the drift in climate normals starting with the 1951-1980 normals, and ending with the 1981-2010 normals. Future climate normals projected by two GCMs are shown as red dots and labelled by median year. The red dotted line is a spline interpolation with no climatic meaning except to visually connect the projected normals. Collectively, the intent of all these graphic elements is to contextualize GCM projections in terms of historical climatic variability and regional climatic variation.

Figure 7: Historical climatic variability and GCM-projected climate change at Stewart, BC.

Figure 8: GCM projections for Stewart BC in the third discriminant dimension

Historical interannual climatic variability at Stewart is substantial, though the range of variability generally falls within the 2D climate envelope of the coast. Despite this interannual variability, climate change has been limited over the past 60 years at Stewart, as suggested by the blue moving average of climatic normals. Nevertheless, there is a trend towards droughtier conditions.  Neither of the GCM runs project a continuation of drying trend in the 2011-2040 normal period (labelled “2025”).  Instead, they both project a moderate to rapid shift towards warmer winters initially, then a shift towards warmer winters and droughtier growing seasons. Within the 2D climate space, there are historical climatic analogs for the end point of both projections: conditions in the 2071-2100 normal period (labeled “2085”) are similar to the 1961-1990 climate normals for Agassiz the HadGEM1_A1B(1) run and for Port Renfrew in the cgcm3_A2(4) run. These analogs also apply in the third dimension (Figure 8).

Figure 9: Historical climatic variability and GCM-projected climate change at Blind Channel, BC.

The other two case studies are more than 5 degrees of latitude south of Stewart, and in the neighbouring GCM grid cell. Projected climate changes in these two case studies are more dramatic, especially at Agassiz. The endpoint of both GCM runs is entirely novel to the coast, and well beyond the range of interannual variability over the past 60 years. It is worth noting here that a climate envelope classification based on discriminant analysis would likely classify these novel endpoints as CDFmm due to the lack of available analogs.

Figure 10: Historical climatic variability and GCM-projected climate change at Agassiz, BC.

The most striking pattern of the three case studies is that the apparent scale of climate change is dramatically greater in the south coast. Whether or not this reflects the relative scales of climate change in the GCM runs is yet to be investigated. The most likely reason for this effect is that the second discriminating axis exaggerates changes towards progressively drier conditions. This is entirely possible because CMD was not log-transformed even though it is partially derived from precipitation. To check this hypothesis, I reran the analysis without CMD and Eref. Surprisingly, the results were very similar. Further investigation is required to determine whether the second axis contains distortions that exaggerate climate change in drier areas of the climate space.

Overall evaluation of the method

The purpose of this analysis was to create a low-dimensional (2D or 3D) climate space that reasonably represents the geographic variation of climate on the BC coast, and use that climate space to historical climatic variability and project climate changes. Was the method successful? Here are my thoughts:

• Discriminant analysis of the 1961-1990 normals of the 1600m grid locations produce a 2D climate space that reasonably represents the relationships between BEC variants, with some exceptions. These exceptions were resolved in the 3D climate space. It is likely that this reduced climate space could be improved through careful variable selection and transformation. Notably, CMD is calculated from climate normals in ClimateWNA, which skews its distribution towards zero. For this reason, CMD (and its precursor Eref) should be either dropped from the analysis or recalculated from monthly observations.
• The dimensions of the reduced climate space for the coast are climatically meaningful, and are relatively free of the cryptic artefacts that drove the second and third dimensions of the Quesnel-Barkerville transect. This indicates that there is a “goldilocks”-type sweet spot for the scale of study areas for reduced climate spaces. Transects are too small, because they only represent a single axis of climatic variation. The province is too big, because it contains many important axes of variation that cannot be meaningfully represented in 2D or 3D. The regional scale appears to be “just right”.  Of course, it is unclear whether the “novel conditions” portrayed in the case studies are truly novel or just regionally novel. This is the tradeoff between local vs. global classification that I discussed in the previous blog post.  Nevertheless, it may be possible to resolve this tradeoff with the right method.
• This background of regional variation provides a useful context for interannual variability and climate change. However, the (possibly fatal) flaw of this approach is that it constrains the representation of climate variability and change to the primary dimensions of existing geographic variation (winter temperature and summer drought in this case). Variation in other dimensions is invisible, even though it may be highly biologically relevant. This is a potential problem in situations like the visualization for Stewart, where there appeared to be climatic analogs for the projections in 2D and 3D space. It is much less of a problem where novel conditions are projected in 2D, as in Blind Channel and Agassiz, since it is highly unlikely that analogs would exist in full-dimensional climate space that don’t exist in the reduced space.
• Projections of climate change into novel areas of the climate space are visually dramatic. However, the discriminant axes may conceal non-linearities that distort the representations of these changes. Any axes created through discriminant analysis require extensive validation.
• This analysis was primarily an exercise in visualization. In order for the method to be successful, however, it must also facilitate quantitative analyses of climatic variability, differentiation, and novelty. I’m not yet convinced that discriminant analysis provides reliable axes (metrics) for quantitative purposes.

Next steps

• Develop and test-run some other methods for achieving the same objective. The most obvious alternate approach is multidimensional scaling. Another alternate approach is to use LDA or some other method to assist with supervised variable selection, and use raw climatic variables as visualization axes. For example, it appears from this analysis that MCMT, SHM, and MAP may be viable first, second, and third axes for a more intuitive reduced climate space for the coast region. Ultimately, the method must be compatible with a quantitative analysis of climatic variability, differentiation, and novelty (a la Williams et al. 2007).
• Develop methods to resolve the compromise between local and global classifications. An example of a potential approach is a two-step graphical approach. The first graph would plot grid points for the entire western North America in the discriminant axes of the coastal climate space and portray GCM projections for a focal (BC coastal) location against this background. The axes of the second graph could be the discriminant axes for the climatic analog (Oregon? California?) for the projected conditions of the focal location. This approach would provide a graphical validation of the inferred analog. Localized LDA (described in the previous blog post) may be a viable method for this approach.
• Evaluate the pros and cons of including interannual climatic variability and the GCM projections into the discriminant (or alternate) analysis. The advantage would be that the reduced climate space could portray key axes of temporal variability that don’t exist in the geographic climatic variation. The disadvantage is that it may be difficult to simultaneously represent spatial and temporal variability in a reduced climate space at the regional scale. This approach may only be viable for transects.

References

Williams, J. W., S. T. Jackson, and J. E. Kutzbach. 2007. Projected distributions of novel and disappearing climates by 2100 AD. Proceedings of the National Academy of Sciences of the United States of America 104:5738–42.

To provide a background for historical and projected climate trajectories at barkerville, I am doing a discriminant analysis on the geographical climate variability of the BEC variants along the Quesnel-Barkerville transect. This discriminant analysis creates a distinctive arch-shaped curve in the first two discriminating dimension, and a zig-zag in the third dimension (Figure 1). The arch shape is similar to the results from pilot project discriminant analysis of interannual variation at representative sites along the transect. The purpose of this blog is to answer: What is the meaning of these non-linearities in reduced climate space? Do they indicate a methodological problem?

Figure 1: Biogeoclimatic variants of the Quesnel-Barkerville transect plotted in a 3-dimensional space created using discriminant analysis. Points are the ClimateBC 1961-1990 normals for 1600-m grid cells covering the geographic area of each variant. The points form an arch in the x-y plane (left) and a zig-zag in the x-z plane (right).

Hypothesis

The most likely explanation for these shapes are that they are an inevitable outcome of reducing the dimensionality of a climate space made up of variables with non-linear relationships to each other along a gradient. If all relationships were linear, climatic patterns in the gradient could be reduced to a multidimensional plane. Any two-dimensional representation (e.g. x-y, x-z, or y-z) of this plane would appear as a linear relationship, plus some noise. Indeed, if the linear relationships between the variables are correlated, then the first axis (first rotation) will capture all of these linear relationships, and all of the other axes will simply represent the noise in the rest of the dimensions. However, non-linear relationships between the variables will not be captured in this first rotation. It is likely that second-order non-linearities (arches) will be most common and therefore represented on the y axis. Third-order non-linearities (zigzags) would be captured in the third axis. This sequencing of non-linearities is evident in the Quesnel-Barkerville transect. Notably, the fourth axis contains a fourth order non-linearity.

Non linearities might also be the result of rotations that combine two or more variables that are static along a portion of the transect. Any non-linear variables that have discriminating variability on only one side of the first axis are possibly being combined on the second axis with uncorrelated variables that have discriminating variability on the other side of the first axis. This could also create an arch shape. An example of two variables that might create this effect are PAS and CMD. Both variables can be static (zero) along one extreme of the elevation range (Figure 2). Low elevation variants will be poorly discriminated by PAS, but highly discriminated by CMD. High elevation variants may have the opposite differentiating variables. it seems likely that both could be associated with a secondary discriminating axis, not because they are correlated, but because there is little overlap between the groups (variants) that  they discriminate.

Figure 2: relationship of cumulative moisture deficit to mean annual temperature. CMD is an example of a variable that differentiates only a subset of the groups (BEC variants). It is static (zero) within most of the ESSF. This non-linearity cannot be resolved with a transformation.

Structure of the discriminant functions

Despite the use of 14 climate variables, there are only five linear discriminating (LD) functions because there are only 6 BEC variant groups (Table 1).  The first discriminating function accounts for 72% of the total separation achieved by the analysis, and almost all of the variables are highly correlated with this axis. This dominance of the first LD axis across all variables is expected because the groups are associated with an elevation gradient that strongly affects all variables. Winter variables have lowest correlation, which reflects the poor differentiation of BEC variants during the winter, as observed in previous analyses. The second and third axes both account for enough differentiation to be retained (>5% of trace). The first axis is clearly an index of dry/warm conditions. the second axis has no obvious intuitive climatic interpretation. However, it is notable in a scatterplot matrix (Figure 3) that the variables that are most highly correlated with LD2 (PAS, DD_0, MAP, MAT) have second-order nonlinearities with variables that are highly correlated with LD1 but highly uncorrelated with LD2 (CMD, AHM, SHM). Similarly, the variables that are most highly correlated with LD3 (MCMT, PAS, CMD, TD) have third-order nonlinearities with LD1 variables that are highly uncorrelated with LD3 (MWMT, DD5, SHM, Eref). This pattern supports the main hypothesis stated above, i.e. that the discriminant functions primarily reflect different types (orders) of non-linearities in the relationships between variables.  The relationship of CMD to LD1 and LD2 indicates that the secondary hypothesis (presence of zero values along a portion of the gradient) is not a major factor in creating the arch.

Table 1: structure coefficients of the discriminant functions. These are simply the correlations between the values of each observation on the discriminant functions with their original values on each climate variable.

Figure 3: scatterplot matrix showing the pairwise relationships between individual climate variables.

Are log-linear transformations helpful?

Presumably, some of the non-linearities in the data are due to non-normal distributions in some of the variables. Transformation of these variables should therefore reduce the arch shape in the second dimension.

The discriminant analysis methodology and data used for the Quesnel-Barkerville transect analysis is largely derivative of Hamann and Wang (2006). These authors used a log10-transformation on precipitation and heat-moisture indices. Log10 transformations of these variables appear to be appropriate within the Quesnel-Barkerville BEC variants (Figure 4). The only exception is CMD, which should not be log-transformed due to the large number of zero observations. CMD was not included as an analysis variable in Hamann and Wang (2006). I am not clear on whether it was intentionally excluded, or simply not available in the early versions of ClimateBC used for that study.

Figure 4: histograms of precipitation and heat-moisture indices over the geographic area of the six BEC variants of the Quesnel Barkerville transect.

Discriminant analysis using the transformed precipitation and heat-moisture variables increases the discriminating power (proportion of trace) of the first discriminating function by 7% (Table 2). Nevertheless, the arch relationship between LD1 and LD2 is graphically very similar to the analysis on the untransformed variables (Figure 5). This indicates that the skewed distributions in moisture variables contribute to the arch and zigzag shapes in LD2 and LD3, but that they are not solely responsible. This makes sense, since some temperature variables are correlated with LD2. Nevertheless, it is surprising that the arch shape remains so well preserved in the analysis of transformed variables.

Table 2: Proportion of trace for the discriminant functions (axes) with transformed and untransformed moisture-related variables. This metric describes the amount of the total differentiation achieved by the analysis that is attributable to each discriminant function.

Figure 5: comparison of the first two discriminant axes of the Quesnel-Barkerville BEC variants using untransformed (left) and log-10 transformed (right) precipitation and heat moisture variables (MAP, MSP, PAS, AHM, and SHM).

Is this the arch or horseshoe effect?

Superficially, this situation resembles the arch effect and horseshoe effect that frequently occurs in multivariate analysis of plant communities along environmental gradients. The horseshoe effect is generally regarded as a fatal flaw for use of PCA with species occurrence data in which species have a unimodal response curve to an environmental gradient . The gaussian response curves of each plant species overlap along the gradient, meaning that the species are non-linearly related to each other in species space. This is the cause for the arch effect. In addition, PCA may incorrectly portray sites on the opposite sides of the environmental gradient as being closely related due to a large number of shared species absences at these sites. This causes the ends of the arch to become incurved, and confounds interpretation of the principle components. This situation is called the “horseshoe effect”. Michael Palmer provides a good summary of the horseshoe effect and its causes (http://ordination.okstate.edu/PCA.htm).

Arch effects have been noted in principal components analysis (PCA), correspondence analysis (CoA), principal coordinates analysis (PCoA) and non- metric multidimensional scaling (NMDS) (Podani and Miklos 2002). I couldn’t find any reference to the arch or horseshoe effect in the context of discriminant analysis, but I’m also not aware of a reason that the same principles wouldn’t apply to discriminant analysis. Podani and Miklos (2002) make the important and obvious point that “whereas unimodal species response to long gradients always leads to horseshoe (or arch)-shaped configurations in the first two dimensions, the converse is not true; curvilinear arrangements cannot generally be explained by the Gaussian model.”  The climate variables making up the climate space in the Quesnel-Barkerville transect do not have Gaussian distributions along the elevation/temperature gradient. More generally, they never have zero values at both ends of the gradient. For this reason, there is no reason to believe that the arch shape displayed in the first two discriminant axes is the “arch effect” commonly observed in species space.

Conclusions

This informal investigation appears to have explained the reason for the second-order (arch) and third-order (zig-zag) non-linearities in the second and third discriminating functions for the Quesnel-Barkerville BEC variants. the discriminant functions appear to primarily reflect different types (orders) of non-linearities in the relationships between climate variables, rather than functional differences in individual variables.  this is a logical, and possibly inevitable, result of doing discriminant analysis on groups that are sequentially located along an environmental gradient. This is an important result, because functional differences in climate variables are often sought during interpretation of the meaning of discriminating functions. In this case, functional differences between variables are a secondary consideration. It is unclear how the effect observed on this environmental gradient will play out over larger and more complex sets of BEC variants. Investigating this effect over a broader geographic scale (e.g. all BEC variants of the BC coast) is an important next step.

Climate variables with limits at zero (e.g. CMD) may create non-linearities in the reduced-dimensional climate spaces. However, no major problems were evident from this cursory analysis. The implications of zero-limited variables (e.g. CMD) for discriminant analysis should be investigated.

References

Hamann, A., and T. Wang. 2006. Potential Effects of Climate Change on Ecosystem and Tree Species Distribution in British Columbia. Ecology 87:2773–2786.

Podani, J., and I. Miklos. 2002. Resemblance coefficients and the horseshoe effect in principal coordinates analysis. Ecology 83:3331–3343.

R Code

###R code for Blog 1 (Sept29, 2013) "Non-linearities in reduced-dimension climate spaces along a climatic gradient"
##investigate non-linearities in the linear discriminant analysis of the Quesnel-Barkerville BEC variants

#Read in the ClimateBC 1961-90 Normals on the 1600-m grid provided by Tongli.
mydata <- read.csv("C:\\Users\\Colin\\Documents\\Masters\\Research\\Data\\ClimateBC\\BEC7v      1600_QuesBark_Normal_1961_1990Y.csv", strip.white = TRUE, na.strings = c("NA","") ) 
mydata <- mydata[-(which(mydata$MAT==-9999)),]
colnames(mydata)[1]<-"id400"
colnames(mydata)[2]<-"var"
mydata$var<-gsub(" ","", mydata$var , fixed=TRUE) #strip white space 
VarList<-data.frame(var=as.character(levels(factor(mydata$var))),
      count=as.vector(table(mydata$var)))

#Remove irrelevant and collinear variables
mydata$DD_18<-NULL #remove heating degree days (irrelevant)
mydata$DD18<-NULL #remove cooling degree days (irrelevant)
mydata$NFFD<-NULL #very highly collinear with FFP
mydata$bFFP<-NULL #very highly collinear with FFP
mydata$eFFP<-NULL #very highly collinear with FFP
mydata$EMT<-NULL #not very relevant
mydata$EXT<-NULL #very highly collinear with MCMT and not very relevant

##histograms of untransformed and transformed variables
par(mfrow=c(6,2))
par(mar=c(2,4,4,1))
hist(mydata[,10], main=colnames(mydata)[10])
hist(log10(mydata[,10]), main=paste("log10 (",colnames(mydata)[10],")"))
hist(mydata[,11], main=colnames(mydata)[11])
hist(log10(mydata[,11]), main=paste("log10 (",colnames(mydata)[11],")"))
hist(mydata[,12], main=colnames(mydata)[12])
hist(log10(mydata[,12]), main=paste("log10 (",colnames(mydata)[12],")"))
hist(mydata[,13], main=colnames(mydata)[13])
hist(log10(mydata[,13]), main=paste("log10 (",colnames(mydata)[13],")"))
hist(mydata[,17], main=colnames(mydata)[17])
hist(log10(mydata[,17]), main=paste("log10 (",colnames(mydata)[17],")"))
hist(mydata[,19], main=colnames(mydata)[19])
hist(log10(mydata[,19]), main=paste("log10 (",colnames(mydata)[19],")"))
par(mfrow=c(1,1))
par(mar=c(4,4,2,2))

#Pairs plot of x variables, using an N=500 random subsample to avoid overplotting
mysample <- mydata[sample(1:nrow(mydata), 500, replace=FALSE),]
library("car")
scatterplotMatrix((mysample[,6:length(mysample)]))

#
x <- mysample$MAT
y <- mysample$CMD
xhat<-tapply(x, INDEX = mysample$var, FUN = mean, na.rm=TRUE)
yhat<-tapply(y, INDEX = mysample$var, FUN = mean, na.rm=TRUE)
plot(x,y,col=as.numeric(factor(mysample$var)), pch = as.numeric(factor(mysample$var)),
 xlab = "Mean Annual Temperature", ylab = "Cumulative Moisture Deficit (CMD)", 
 main="")
text(xhat, yhat, levels(factor(mysample$var)), font=2)

library(MASS) #contains lda; also eqscplot(x,y) function, which plots axes of equal scale
library(stats)

#define the matrix of x variables
X<- mydata[,6:length(mydata)]
#transform the 5 precip and heat-moisture variables (optional). 
#X$MAP <- log10(mydata$MAP)
#X$MSP <- log10(mydata$MSP)
#X$AHM <- log10(mydata$AHM)
#X$SHM <- log10(mydata$SHM)
#X$PAS <- log10(mydata$PAS)

##MDA 
z <- lda(X, mydata$var)

##Data frame of discriminant axis scores
M <- as.data.frame(predict(z)$x )

##Plot of BEC variants on first two discriminant axes
x<-M[,1]
y<-M[,2]
z3 <- M[,3]
A<-factor(mydata$var) #points are labelled by mapped bec variant. 
variant<-"ESSFwk1"

# Spinning 3d Scatterplot
library(rgl)
#?plot3d
#?rgl
plot3d(x, y, z3, 
 col = as.numeric(factor(A)), 
 alpha = 0.5,
 ylim=c(0.1,max(y)),
 zlim=c(0.1,max(z3)),
 size=2)
##Labels at the 3D psuedocentroid
xhat<-tapply(x, INDEX = mydata$var, FUN = mean, na.rm=TRUE)
yhat<-tapply(y, INDEX = mydata$var, FUN = mean, na.rm=TRUE)
z3hat<-tapply(z3, INDEX = mydata$var, FUN = mean, na.rm=TRUE)
text3d(xhat, yhat, z3hat, text = levels(factor(A)), font=2)

######################Quantitative Analyses

##Structure coefficients
 #Correlations between predictors and discriminant values 
 #indicate which predictor is most related to discriminant function (not corrected for the            #other variables)
CorXM <- cor(X,M)
write.csv(CorXM, "C:\\Users\\Colin\\Documents\\Masters\\Research\\Data\\CorXM.csv")

##Basic descriptive output for the LDA
#the singular value decomposition gives the ratios of the between- and within-group 
         #standard deviations on the linear discriminant variables. 
z
z$svd