Non‐stationarity, diversity gradients and the metabolic theory of ecology

Tests of mechanisms driving species diversity have usually been conducted using exploratory multiple regression, trying to associate alternative (and usually not mutually exclusive) hypotheses with sets of environmental variables or other surrogates of ecological and historical processes that are correlated with richness patterns (see Willig et al., 2003; Hawkins et al., 2003; Currie et al., 2004). More recently, in the context of the metabolic theory of ecology (MTE hereafter; see Brown et al., 2004), Allen et al. (2002) formally proposed a hypothesis stating that the relationship between ln-transformed species and 1/kT (where k is Boltzmann's constant and T is temperature in Kelvin) has a slope between −0.6 and −0.7 (referring to the model version presented in Brown et al., 2004; see also Allen et al., 2006; Hawkins et al., 2007a). However, Algar et al. (2007) recently showed that for several data sets in North America the relationship between ln-transformed richness and temperature is actually curvilinear, so slopes vary systematically in geographical space. As a consequence, they proposed that Allen et al.'s (2002) model did not give an adequate fit to the data (see also Hawkins et al., 2007a, for a more general test, and Gillooly & Allen, 2007, for a reply) and defined in which parts of the continent the expected slope of the MTE could be found, based on fitted quadratic models of richness against temperature. However, we will show here that even when the linear relationship between ln richness and temperature predicted by Allen et al.'s (2002) model is found, the same problem of spatially variable slopes can appear because of non-stationarity in the data, i.e. the relationship between variables can vary systematically from one locality, or region, to the other. Actually, both non-spatial and spatial regression models assume that results apply equally to the entire study area, but actually they may only provide an ‘average’ of the relationship if non-stationarity occurs in data. Thus, although violation of the stationarity assumption in spatial regression may create difficulties in generalizing their results, application of models that incorporate non-stationarity may be useful in revealing more complex spatial structures in data and relationships between variables in multiple scales. Indeed, geographically weighted regression (GWR) procedures have been developed to deal with this and have recently been applied in macroecological analyses (Fotheringham et al., 2002; Foody, 2004, 2005; Wang et al., 2005; Bickford & Laffan, 2006; see also Jetz et al., 2005 and Foody, 2005 for a discussion about interpretation of GWR in richness data). The advantage of GWR is that instead of estimating a single slope that describes the studied relationship, it generates one regression slope for each locality, or ‘cell’ (called ‘local’ slopes), by using a complex scheme of spatial weights (see Fotheringham et al., 2002 for a detailed description). These slopes (as well as estimates for other regression parameters, such as r2) can then be mapped and used to investigate how the relationship between variables changes across geographical space. Here, we analysed spatial patterns of species richness derived from digital maps for 3019 species of New World amphibians, obtained from the Global Amphibian Assessment (see Stuart et al., 2004), available at http://www.natureserve.org, using a grid with 4187 cells of 1° of latitude and longitude. Our purpose is to expand Algar et al.'s (2007) discussion by using GWR to investigate patterns of non-stationarity in the relationship between ln-richness and temperature expressed as 1/kT. Species richness follows the well-known latitudinal pattern in the Western Hemisphere (e.g. Hawkins et al., 2006 and references therein) (Fig. 1a), and the observed slope (−0.686 ± 0.005) matches MTE predictions quite well (Fig. 2). However, there is some heteroscedasticity caused by a greater variation of richness in warm areas (i.e. areas with low 1/kT values), and a ‘constraint envelope’, commonly found in macroecological studies, may actually better describe the relationship. Although a linear regression model is not entirely statistically adequate to describe the data, the r2 is also very high (= 0.835), and adding a quadratic term did not significantly improve the model fit (r2= 0.837). Although confidence intervals are underestimated because of the strong spatial autocorrelation in data, the relationship remains significant after correcting the degrees of freedom using Dutilleul's (1993) procedure (r=−0.914; d.f. = 9.45; P < 0.01; see also Hawkins et al., 2007a). However, when using a simultaneous autoregressive model (SAR; implemented in sam 2.0; see Rangel et al., 2006), the slope becomes −0.514 ± 0.011 (still close to MTE predictions), but the r2 based on 1/kT drops to 0.127. Although this is still significant at P= 0.05, the differences between the ordinary least squares (OLS) and SAR suggest complex hidden spatial structures at distinct scales and perhaps variation in slopes in different regions (see Diniz-Filho et al., 2003). (a) Patterns of species richness for 3019 amphibian species in the New World, followed by patterns of spatial variation in regression parameters from GWR, including (b) the slope of ln-transformed species richness against 1/kT, (c) the model fit, given by r2, and (d) the GWR residuals. Global relationship between ln-transformed species richness against 1/kT, supporting MTE predictions (model I regression slope =−0.686). We then applied the GWR (implemented in GWR 3.0; see Fotheringham et al., 2002) to detect non-stationarity in the relationship between richness and temperature. The weightings for GWR were defined by a kernel function, with bandwidths interactively established using the Akaike criterion. There was an improvement of GWR with respect to the OLS model (r2 of GWR = 0.914, and F= 25.36 for the improvement of GWR over OLS, taking into account the differences in degrees of freedom between the two models; P << 0.01). Thus, a highly significant non-stationarity was detected in the relationship. Indeed, local slopes of ln-richness against 1/kT vary between −1.33 and +0.97 across the New World, with a median of −0.64 (Fig. 1b). A good match with MTE predictions is found mainly in the middle latitudes of North America, and in small bands of South America. In the Northern Hemisphere, both the explanatory power (r2) and the slopes of the local regression increase almost continuously towards northern latitudes (Fig. 1c), supporting the curvilinearity found by Algar et al. (2007) in this region. Residuals of the GWR do not show spatial autocorrelation (Moran's I in the first distance class equal to 0.073; see Diniz-Filho et al., 2003), although there are higher values clustered along the mountains and coasts (Fig. 1d). The non-stationarity in the relationship between richness and temperature is in agreement with Hawkins et al.'s (2003) conjecture that energy is a more important driver of diversity in higher latitudes (and that water is probably a more limiting factor in the tropics, although this second component was not explicitly tested here; see also Whittaker et al., 2007). In general, spatially structured variations in the effect of multiple predictors may be the underlying cause of the non-stationarity in the relationship between richness and temperature detected here. Another potential explanation of non-stationarity is that effects of temperature on other variables that are assumed to be constant in Allen et al.'s (2002) model (i.e. average body mass and abundance) also vary in space, generating a systematic lack of agreement with the MTE prediction in some regions but not in others. Finally, a nonlinear and spatially structured relationship between energetic drivers of diversity and temperature (that is actually a surrogate of these) (see Huston, 2003; Clarke & Gaston, 2006) may be part of the explanation for non-stationarity. Thus, a GWR model using multiple predictors could be used to further investigate the causes of this non-stationarity in a more complex way. Our GWR analyses of New World amphibians strongly support the recent conclusion by Algar et al. (2007) that MTE models for richness provide variable fits in different regions, with better fits found in mid North America for many groups (for amphibians, trees, blister beetles and tiger beetles). As previously discussed, the curvilinear relationships found by Algar et al. (2007) could be attributed to an absence of important ecological components in Allen et al.'s (2002) model or in part because of violation in the assumptions of constant abundance and body size. However, more importantly, our analysis using GWR shows that this does not occur only when the modelled relationship is curvilinear. Even a linear relationship that is, in principle, consistent with the MTE can hide more complex and non-stationary patterns of spatial variation in the responses of richness to temperature, creating a failure in Allen et al.'s (2002) simple predictions (see Hawkins et al., 2007b). In a recent comparative analysis, Hawkins et al. (2007a) did not find consistent support of the MTE in higher latitudes, but since many data sets were compared, other sources of variation across data, probably related to how different organisms respond to temperature, may make it difficult to detect this better consistency with the MTE at higher latitudes. GWR analyses can be extended to other groups of organisms to check for a correspondence between patterns of non-stationarity which, in turn, could help in elucidating how general the MTE is and in which geographical domains Allen et al.'s (2002) model correctly predicts richness patterns. We thank Brad Hawkins, Luis M. Bini, Andrew Allen and two anonymous reviewers for critical reading of the manuscript and useful discussions about MTE. F.A.S.C. and T.F.L.V.B.R. are supported by a CAPES and CAPES/Fulbright doctoral fellowships, whereas work by J.A.F.D.-F. has been continuously supported by several CNPq and FUNAPE/UFG grants.