Should Species Distribution Models Account for Spatial Autocorrelation? A Test of Model Projections Across Eight Millennia of Climate Change

ABSTRACT


INTRODUCTION
The last decade has witnessed a marked increase in the application of models that project the potential geographic distributions of species by linking observations of species occurrences to environmental predictor variables.These models, commonly called bioclimatic envelope, ecological niche, or species distribution models (hereafter SDMs), are important tools for forecasting impacts of climatic change on biological diversity and for generating conservation plans and climate-change policy (Guisan & Thuiller, 2005).To project future distributions under different, plausible scenarios of climatic change, SDMs use statistical relationships between present-day distributions of species and climate (Elith et al., 2010).Although generally successful at explaining and predicting current distributions of species (Franklin & Miller, 2009), impact assessments derived from SDMs have been criticized for their reliance on a number of largely untested ecological assumptions, methodological issues, and statistical concerns (e.g., Pearson & Dawson, 2003;Dormann, 2007).
Chief among these issues is the failure of most SDMs to account for spatial dependence of occurrence data (Gelfand et al., 2006;Bahn and McGill, 2007;Dormann, 2007;Elith et al., 2010).Spatial autocorrelation arises in ecological data because nearby points tend to be more similar, in physical characteristics and/or species occurrences or abundances, than are pairs of locations that are farther apart (Legendre, 1993).When model assumptions about independent and identically distributed residuals are violated, there could be a bias in the regression parameter estimates, potentially leading to poor inference.Studies illustrate that failure to account for spatial autocorrelation can lead to misidentification of important driving variables and overly optimistic error rates (e.g., Lichstein, et al., 2002;Segurado et al., 2006;Diez & Pulliam, 2007;Dormann, 2007), especially when small-scale patterns of explanatory variables create instability in broad-scale regression parameter estimates (Hawkins et al., 2007).Further, models based solely on spatial interpolation can provide better fits to species range data than models based on explanatory environmental variables (Bahn & McGill, 2007), suggesting that spatial autocorrelation in unmeasured factors (e.g., population processes such as dispersal or underlying resources such as soil moisture) may account for most of the observed distributional patterns.
Analysis of spatial SDMs primarily has focused on predicting current or simulated species' distributions using a hold-out dataset for model validation (Gelfand et al., 2006;Wilson et al., 2010), but projections of spatial SDMs in changing climates over long time scales remain largely untested.Observed changes in species distributions as a result of past climatic dynamics provide a unique opportunity to compare projections of spatial and non-spatial SDMs parameterized with current conditions (Pearman, et al., 2008a;Nogués-Bravo, 2009;Dobrowski et al., 2011, Veloz et al., 2012).
Projections to environmental conditions different from those used to calibrate SDMs are subject to error (Heikkinen et al., 2006) and may not be ecologically meaningful or statistically valid if there are changes in correlations between variables across time and space (Elith et al., 2010) or if species-environment relationships are not conserved (e.g., Fitzpatrick et. al., 2007, Veloz et al., 2012).It also is not known whether it is desirable to project models with spatial random effects based on the partially observed spatial distribution of a species at one time point into a new temporal domain.
In this study, we developed non-spatial and spatial SDMs for two genera of trees in eastern North America.We calibrated the models with current climate data and Forest Inventory and Analysis (FIA) data collected by the United States Forest Service.We then projected the models back in time using paleoclimate simulations and extensive pollen records as independent validation data.Our approach is similar to that of Pearman et al. (2008a), who used fossil pollen to validate SDMs of European trees projected back to a single time in the mid-Holocene (6,000 years before present).However, the availability of new paleoclimate reconstructions, which provide millennial snapshots of historic climate for the last 21,000 years before present, allowed us to validate the models at a much finer temporal resolution.
To assess the usefulness of adding a spatial term to SDMs we consider the following: 1) a spatially-varying intercept model with no climate variables; 2) a non-spatial model with climate variables; and 3) a spatially-varying intercept model with climate variables.As detailed in the Methods Section and Appendix S3, the spatially-varying intercept was introduced via spatial random effects.The rationale for choosing these candidate models is a follows.If climate variables describe a significant portion of the variability in the observed distribution and if these variables change over time, then projections from models with climatic variables will show a conservative shift away from the observed distribution.For the spatially varying intercept model with climate variables, any projected shifts in distributions are tempered by the spatial random effects.Depending on the amount of spatial autocorrelation, spatial random effects act to draw the projected distribution back toward the observed distribution used to calibrate the model.If climate variables do not describe a significant portion of the variability in the observed distribution, then the spatial random effects will keep projected distributions close to the observed distribution, i.e., the only learning for prediction will come from the observed distribution and hence projected probability of species occurrence will be similar to the observed probability of occurrence.With these three candidate models, we were able to tease apart differences due to the spatial random effects alone, the climate variables alone, and their additive effects.We parameterized and estimated model parameters following a Bayesian framework, which provided full posterior distributions for model parameters and allowed us to estimate the uncertainty in our statistical inferences.We focus on two tree genera, Fagus and Tsuga, whose distributions can be readily inferred from fossil pollen and which possess contrasting life histories.
We address three questions.(1) Do non-spatial SDMs of current distributions of Fagus and Tsuga based on climate variables exhibit residual spatial autocorrelation?(2) Do SDMs with spatial random effects that include or exclude climate variables provide better fits to the observed distributions than non-spatial SDMs with climate variables only?(3) Do hindcasted spatial SDMs better predict historic distributions than non-spatial SDMs?

Study genera
We studied two tree genera, Fagus and Tsuga.In eastern North America, Fagus is represented by only one species, F. grandifolia (Ehrh.)(American Beech), and Tsuga by only two, the widespread T. canadensis (L.) Carr.(Eastern Hemlock), and the narrow endemic T. caroliniana Engelm.)(Carolina Hemlock).For both Fagus and Tsuga, the relationship between local abundance of trees and relative abundance of pollen in sediment cores has already been derived (Davis, 1981).Tsuga is a conifer with passively-dispersed cones, whereas Fagus is deciduous with animal-dispersed seeds.

Occurrence data
We used FIA data to describe the current distribution of Fagus and Tsuga.In every 2428 ha of land in the United States classified "forested", there is one permanent FIA plot, each containing four 7.2 m fixed-radius subplots (Woudenberg et al., 2010).In each subplot, all trees >12.7 cm diameter at breast height have been measured periodically since the 1940s; consistent nationwide annual inventories were initiated in 2001.We used data from the most recent full plot inventory (2003 -2008) to calibrate our models.
Historic distributions of Fagus and Tsuga were derived from fossil pollen data in the Neotoma Paleoecology Database (<www.neotomadb.org>).Paleoclimate data (described below) were available at 1 kiloannum before present (kaBP) intervals from 0-21 kaBP, so we focused on millennial historic distributions of Fagus and Tsuga.Given the variation in temporal scale and spatial resolution across study sites and uncertainties associated with radiocarbon aging of pollen from sediment cores (Blauw et al., 2007), we compiled pollen datasets in which Fagus and Tsuga were counted as present at a site if their pollen percentages reached threshold levels at any time within 500 years centered on each historic millennium (Appendix S1).We chose a 500 year window because cross-validation analyses of biostratigraphic ages from recently revised age models for all pollen sites suggested that 500 years is a conservative estimate of temporal uncertainty for sites in the Neotoma database (Blois et al., 2011).To determine the sensitivity of historic tree distributions to the pollen percentage thresholds used to define a genera's presence or absence at a site, we specified low and high thresholds for each genus (Pearman et al., 2008a): 0.5% or 1% for Fagus and 1% or 2% for Tsuga (Davis, 1981).

Extent and resolution
The extent of the study area was the portion of eastern North America with the highest density of pollen data (Fig. 1).This region contained 75,251 FIA sites and up to 379 Neotoma locations, depending on time period considered.Paciorek & McLachlan (2009) found that spatial patterns relating current and past climates to abundances of pollen and trees were unreliable at resolutions below ~50 km, so the climatic predictors for our model (see below) were downscaled to a resolution of 0.5-degrees (~50-80 km depending on latitude).We upscaled the current tree occurrence data for each grid cell in the climate spatial data layers, keeping track of the number of FIA sites per 0.5-degree cell to be used as weights in the models (Appendix S2).Following this aggregation there were a total of 1,419 FIA observations with presence/absence ratios for Fagus and Tsuga of 706/713 and 380/1,039, respectively.The number of aggregated pollen observations varied for each 1 kaBP time period (Fig. 2).Although both paleoclimatic and pollen data extended back 21 kaBP, the total sample size and the number of pollen grains of each genus declined rapidly beyond 8 kaBP (Fig. 2).Thus, our hindcast projections extend only from 1 to 8 kaBP, which allowed us to validate the models using a minimum of 200 grid cells containing observations, and at least 50 of which contain presences for each genus.

Climate data
Modern climate data came from the observed dataset of the Climate Research Unit (CRU), University of East Anglia (Brohan et al., 2006).Paleoclimate data for this study came from a recent transient simulation of the CCSM3 global circulation model (GCM) (Liu et al., 2009).
The standard change-factor approach was employed to statistically downscale and reduce bias in the climate data (Wilby et al., 2004).For each climate variable at each millennial interval, the difference between modeled paleoclimate and modeled modern climate was calculated and then resampled to a 0.5 × 0.5-degrees grid to match the resolution of the CRU observed climate dataset (Mitchell & Jones, 2005).
Decadal averages of seasonal variables were the highest temporal resolution data available from the archived CCSM3 simulations.To get a 'snapshot' of climatic conditions at each millennial time point, decadal averages of seasonal climate variables from the CRU or CCSM3 simulations were calculated for the first 100 years of each millennium (e.g., 8.0 to 7.9 kaBP).Because summaries of modern observed climate are available at centennial scales, these same centennial summaries of paleoclimate were derived to aid comparisons between paleo and modern SDMs.Bioclimatic variables that captured precipitation and temperature averages and seasonalities were used because response surface analyses for Fagus and Tsuga have shown that climatic annual averages, annual ranges, and seasonality were important factors controlling the Holocene migrations of these genera (Bartlein et al., 1986).Specifically, we calculated six bioclimatic variables (Hijmans et al., 2005): annual mean temperature (BIO1), mean diurnal range (BIO2), temperature seasonality (BIO4), temperature annual range (BIO7), annual precipitation (BIO12), and precipitation seasonality (BIO15).
Two of the six calculated bioclimatic variables, temperature seasonality and temperature annual range, had within-time correlations with the other bioclimatic variables ≥0.7, so they were not included as explanatory variables in the models that included environmental predictors (see Appendix S3).The correlations between mean diurnal range and annual precipitation varied between modern and historic times (see Appendix S3), and such changing correlation structures between times could be problematic when projecting models beyond the present (Elith et al., 2010).To determine if sufficient variance in the current distribution was explained by the two remaining variables with stable correlation structures over time (i.e., annual mean temperature and precipitation seasonality), we compared a model with annual mean temperature, precipitation seasonality, mean diurnal range, and annual precipitation with another that included only annual mean temperature and precipitation seasonality.

Model calibration
We used Bayesian generalized linear models (GLMs) to model genera occurrence.While approaches such as neural networks and genetic algorithms have been used for SDMs and although model projections can be sensitive to the type of statistical model employed (Elith et al., 2010), classical approaches do not provide the statistical inferences we sought.Even though GLMs describe a central tendency and not a limiting effect (e.g., of temperature or precipitation extremes), Bayesian spatial GLMs provide exact inference for the random model parameters, including spatial random effects, by estimating entire posterior distributions at both observed and unobserved geographic locations (Gelfand et al., 2006).Because our goal was to compare consistently SDMs with three different specifications (i.e., spatially-varying intercept only (SVI), climate only, and spatially-varying intercept plus climate), we adopted a Bayesian approach in fitting all of the models.Model structure is detailed in Appendix S2; model code is provided in Appendix S4.
Including the SVI has a potential for overfitting as it allows variable intercepts for every location and thus a very flexible spatial fit to the FIA data.As a null model, we also fit a multilevel B-Spline to the FIA data (Lee et al., 1997) using the 'MBA' package of 'R' statistical software to determine whether our hindcasting test for the inclusion of a SVI in the Bayesian models was sufficient.As an exploratory analysis into the strength of the residual spatial dependence in the FIA data, we calculated Moran's I from the residuals of the non-spatial GLMs.This latter analysis was conducted using the Spatial Analyst Tool in ArcMap10 (ESRI, 2011).

Model fit to calibration data
We fit the Bayesian models to 90% of the FIA data (N = 1,277) and randomly selected a 10% holdout dataset (N = 142) to assess predictive performance.We also used DIC to rank the Bayesian models fit to the calibration data (Spiegelhalter et al., 2002).DIC is the sum of the Bayesian deviance (a measure of model fit) and the effective number of parameters (a penalty for model complexity).Lower DIC values indicate better model fit.Models are compared using ΔDIC: where min(DIC) is the DIC value for the model with the best fit (i.e., lowest DIC value).In general, ΔDIC < 2 indicates weak evidence; 5 < ΔDIC < 10 indicates strong evidence, and ΔDIC >10 indicates very strong evidence that one model is preferred over another (Spiegelhalter et al., 2002).

FIA hold-out dataset and pollen validations
When projecting the spatial models back in time for the pollen validation, the random effects serve to draw the projected distributions for each genus back toward that of the observed distribution used for model calibration (i.e., the FIA data) in the new time period (Appendix S2).
To compare the performance of the models in predicting current and projecting past distributions, three measures were calculated using the 'ROCR' package of 'R' statistical software: the Area Under the Curve (AUC) of a Receiver Operating Curve (ROC), false negative rates (FNR), and false positive rates (FPR).The calculation of FNRs and FPRs requires converting the continuous outputs to a binary form using a threshold, in this case the value that maximizes the sum of sensitivity and specificity (Liu et al., 2005;Lobo et al., 2008).
Differences in AUC, FNR, and FPR between models, genera, pollen percentage thresholds, time, and the model × genus interaction were tested with three GLMs.To normalize residuals and reduce heteroskedasticity, AUC, FNR, and FPR were all arcsin transformed.
Model, genera, pollen percentage threshold, and the model × genus interaction entered the GLM as fixed factors, and time entered as a covariate.The model × genus interaction was of particular interest as it tested whether or not different models performed better or worse in hindcasting the presence-absence of the two genera.The data were analyzed with separate GLMs for AUC, FNR, and FPR to facilitate the interpretation of Tukey's Honest Significant Differences post-hoc comparisons at the expense of increasing Type II error rates.Bonferroni corrections of the Pvalues from the tests did not alter the significance of any of the effects.

Parameter estimates and model fit to calibration data
In non-spatial models with two climatic variables (i.e., annual mean temperature and precipitation seasonality) or four climatic variables (i.e., annual mean temperature, mean diurnal range, annual precipitation, and precipitation seasonality), all climatic variables were significant predictors of presence/absence: none of the 95% credible intervals of the parameter estimates included zero (Tables 1, 2).In contrast, in the spatial models some of the climatic explanatory variables were not significant predictors of presence/absence (e.g., annual mean temperature in the Tsuga models with two climatic variables and mean diurnal range in the Fagus model with four climatic variables; Tables 1 & 2).Changes in the magnitude and sign of parameter estimates between non-spatial and spatial models suggested that non-spatial models violated the assumption of independent identically distributed residuals.The residuals of the non-spatial models for both Fagus and Tsuga also exhibited significant positive spatial autocorrelation (Moran's I = 0.604, P < 1 × 10 -7 for Fagus; Moran's I = 0.761, P < 1 × 10 -7 for Tsuga), supporting the conclusion that non-spatial models were inappropriate for these data.
For Fagus, the SVI plus climate model with annual mean temperature and precipitation seasonality had the lowest DIC value and ∆DIC > 10 relative to all other Fagus models (Table 3, Fig. 3).In contrast, for Tsuga, the SVI model with no bioclimatic predictors had the lowest DIC value and ∆DIC > 10 relative to all other Tsuga models (Table 3, Fig. 4).
The non-spatial SDMs for both Fagus and Tsuga that included only annual mean temperature and precipitation seasonality had ∆DIC values >10 relative to the non-spatial models that included annual mean temperature, precipitation seasonality, mean diurnal range, and annual precipitation (Table 3).Given that the correlative relationship between mean diurnal range and annual precipitation was unstable between modern and historic times (see Appendix S3) and that the inclusion of them did not provide a large decrease in the ΔDIC, these two climatic variables were excluded from the models used for prediction that were validated with the 10% holdout FIA dataset and fossil pollen record.

FIA hold-out dataset and pollen validations
For the contemporary 10% hold-out FIA dataset for both genera, the non-spatial model performed worse than the SVI, SVI plus climate, or multilevel B-Spline models in terms of AUC, FNR, and FPR (Table 4; Appendix S5).However, the same was not true when models were hindcasted.Based on AUC, there were significant main effects of model type (non-spatial, SVI, SVI plus climate, FIA B-Spline; F 3,118 = 32.4,P = 2.4 × 10 -15 ), and a significant genus × model interaction (F 3,118 = 13.8,P = 8.8 × 10 -8 ) (Table 4, Appendix S5) on model performance.
For the Fagus hindcasts, on average the non-spatial model had higher AUC values than the spatial models (i.e., SVI and SVI plus climate) and FIA multilevel B-spline models, but the opposite was true for Tsuga.The FNRs in the hindcasting validation varied by model (F 3,118 = 8.1, P = 6.2 × 10 -5 ).The FIA data multilevel B-spline model had the highest FNR and post-hoc comparisons showed that there were no significant differences between the non-spatial and spatial models in FNRs (Table 4, Appendix S5).Similar to the FNRs, the FPRs also varied by model (F 3,118 = 9.0, P = 1.95 × 10 -5 ) (Table 4, Appendix S5).The FIA data multilevel B-spline and the non-spatial models had higher FPRs than the spatial models.There were no significant genus × model interactions for FNRs (F 3,118 = 2.3, P = 0.08) and FPRs (F 3,118 = 1.7,P = 0.18).

Discussion
A key question regarding the application of SDMs to predicting the response of species to climate change is whether the failure to include ecological and evolutionary processes (e.g., dispersal, biotic interactions, readjustment lags) will prove to be problematic (reviewed by Pearson & Dawson, 2003).Depending on the species and its life history, ecological and evolutionary processes may (or may not) lead to its inability to track changes in climate.While there is evidence that vagile organisms (e.g., butterflies) can track rapid climate change (Warren et al., 2001), sessile organisms (e.g., trees) may not readily disperse to newly suitable habitat resulting in limited niche space filling (Svenning & Skov, 2004;Meier et al., 2012).Species undergoing climate driven range expansions coupled with enemy release are hypothesized to be more capable of realizing their potential niche (Hellman, et al., 2012), whereas species limited by a particular resource (e.g., host availability) can be constrained to the spatial distribution of the resource (Merrill et al., 2007).There is evidence that shorter-lived taxa (e.g., insects and herbaceous plants; Woodward, 1990;Thomas et al., 2001) can evolve in response to rapid climate change, but longer-lived taxa that cannot evolve as quickly may experience readjustment lags (Pearson & Dawson, 2003).
For those taxa whose distributions do not shift over time as a result of ecological and evolutionary processes, the inclusion of spatial random effects in SDMs could improve projections by providing a more conservative prediction of distributional shifts, especially when climatic variables do not explain much variability in their observed distributions.Alternatively, when climatic variables explain most of the variability in a taxon's observed distribution and the taxon is capable of tracking climate, then accounting for spatial autocorrelation in SDMs won't provide better projections.In other words, the spatial random effects keep the projected distribution similar to the data used for model calibration, unless the covariates (e.g., climatic variables) suggest otherwise.Further, if the climate variables do not explain much of the variability in the observed distribution and the genera's distribution shifts far from the observed distribution over time, then none of the models defined here will perform well.The predictive abilities of non-spatial and spatial SDMs have rarely been compared with temporally varying validation datasets to test these assertions (Gelfand et al., 2006).
In this study we tested the predictive abilities of non-spatial and spatial SDMs across eight millennia using data from the pollen record (Appendix S1).We found that spatial SDMs had better fits to the calibration data, higher predictive accuracy for a modern hold-out validation dataset, and greater variance in their outputs than non-spatial SDMs (see also Gelfand et al., 2006;Bahn & McGill, 2007).For Fagus, the SVI plus climate model provided a better fit to the calibration data than the SVI model, but the opposite was true for Tsuga.Also for the two climatic variable models, for Fagus there was no change in the sign of the climatic regression coefficients between the non-spatial and spatial models (Table 1), but with Tsuga there was a sign change in the regression coefficient for annual mean temperature between the non-spatial and SVI plus climate models (Table 2).This result suggests that for Tsuga the spatial random effect could be accounting for dependence in the model's residuals across space as several other studies have found that parameter estimates are affected by spatial autocorrelation (Dormann, 2007;Kühn, 2007;Bini et al., 2009;Hodges & Reich, 2010).
In the hindcasting analyses, the SVI and SVI plus climate models performed similarly.
This suggests that the climatic variables do not contribute much to explaining the variability of occurrence relative to that explained by the spatial random effects.AUC values based on fossil pollen indicated that the non-spatial model performed better for Fagus than either of the two spatial models, but the opposite was true for Tsuga.However, FNR values did not differ among the models for either genus, and FPR values were greater for non-spatial models for both genera.
We have more confidence in FNR and FPR values than in AUC values because the latter describe portions of the ROC curve that are rarely encountered and weights omission and commission errors equally (Lobo et al., 2008).With the pollen record, equal weighting of omission and commission errors may not be ideal; we have much more confidence in the presence of pollen grains than in their absence (Blauww et al., 2007;Blois et al. 2011) and false negatives in the pollen record are more problematic than false positives.The lack of differences in false negative rates between models shows that the non-spatial and spatial models have similar FNRs.
Although we have less confidence in actual absences in the pollen data, the FPRs are interesting when considering the ecological and evolutionary processes leading to conserved spatial structure in the distributions of species.The greater FPRs of non-spatial models for both genera suggest that spatial effects may account for smaller-scale climatic spatial structure that is not otherwise estimated in large-scale or averaged temperature and precipitation values (Gelfand et al., 2006;Hawkins et al., 2007).Evidence from the fossil pollen and paleoclimate records suggests that climatic shifts can result in abrupt ecological changes in vegetation that are driven by internal dynamics, such as site-specific environmental characteristics (e.g., soil moisture) or biotic interactions (e.g., competition) that create geographically localized variation in vegetation composition (Williams et al., 2011).Taxon-specific responses to climate forcing also could explain why the SVI model had the lowest DIC for Tsuga and why the two spatial models performed better in regards to both AUC and FPR for Tsuga, but not for Fagus.Approximately 5.5 kaBP Tsuga experienced a range contraction known as the "hemlock decline" potentially due to an abrupt change in climate, a phytophagous insect infestation, or both (Bhiry & Filion, 1996;Foster et al., 2006).If the hemlock decline was due to an abrupt change in climate, then localized ecological changes could have resulted in stronger spatial structure in its distribution.However, decoupling changes in distributions due to climate and spatial structure due to biotic interactions or site-specific abiotic characteristics is difficult because observed spatial structure is (or was) inherently linked to abrupt climate change.
Alternatively, the spatial random effects may have captured a missing covariate, such as an ecological process that generates spatial structure (Clayton et al., 1993;Paciorek, 2010).Such processes could include dispersal, competitive interactions, land-use history, or underlying features of the terrain.For example, if dispersal limitation prevents distributional shifts, then we might expect that spatial SDMs would perform better for dispersal-limited taxa (e.g., Tsuga) that cannot track changes in climate, but not necessarily for taxa with effective dispersal vectors (e.g.,

Fagus
) that can gain dominance by migrating faster to climatically favorable sites (Pearman et al., 2008b).These taxon-specific differences in dispersal mode and degree of dominance could explain why Tsuga seemed to be less responsive to climate over the past 8 millennia than Fagus as evidenced by the better performance over time of the two spatial models in regards to both AUC and FPR for Tsuga, but not for Fagus.Simulation experiments for European trees with spatially explicit process models accounting for changing macroclimate, competition, and habitat connectivity showed that some of the spatial autocorrelation between two time periods may be due to very slow migration rates resulting in severe time lags that are not accounted for in nondynamic and non-spatial SDMs (Meier et al., 2012).Also, Dobrowski et al. (2011) found that non-spatial SDMs fit to widespread plants with more effective dispersal mechanisms had higher predictive accuracy over 75 years of climate change in California than non-spatial SDMs fit to dispersal-limited plants.
Given the results of this study, should researchers include spatial random effects in SDMs?We found that for two long-lived eastern North American trees, spatial models provided better fits to calibration data and lower FPRs, but not necessarily improvements in AUC or the FNR.The better fits of the spatial SDMs may have resulted from the richness of the FIA data used to calibrate the models.Large samples of evenly-dispersed data likely will capture any spatial structure; consequently a spatial SDM should fit well.However, when sample sizes are small, there is less of a chance that the spatial structure will be represented adequately.
Ultimately, whether to include spatial random effects in SDMs will depend on the taxon being modeled, the cost of false positives, and the quality of the data.

Figure 1 .
Figure 1.Map of the study extent in the eastern United States showing Forest Inventory and

Figure 2 .
Figure 2. Numbers of sites with presences (black fill) or absences (white fill) of Fagus (a and c)

Figure 3 .
Figure 3. Maps of a) a surface approximation of the probability of occurrence of Fagus

Figure 4 .
Figure 4. Maps of a) a surface approximation of the probability of occurrence of Tsuga

Table 1 .
Parameter credible intervals (2.5%, 50.0%, and 97.5% percentiles) for the Fagus spatially-varying intercept (SVI), non-spatial (NS2 and NS4) and SVI plus climate (SVI2 and SVI4) models.The numbers two and four in the acronyms for the non-spatial and SVI plus climate models indicate the number of bioclimatic explanatory variables included in the models.

Table 2 .
Parameter credible intervals (2.5%, 50%, and 97.5% percentiles) for the Tsuga spatiallyvarying intercept (SVI), non-spatial (NS2 and NS4) and SVI plus climate (SVI 2 and SVI 4) models.The numbers two and four in the acronyms for the non-spatial and SVI plus climate models indicate the number of bioclimatic explanatory variables included in the models.The two

Table 3 .
Fits of the spatially-varying intercept (SVI), non-spatial, and SVI plus climate SDMs to the modern Forest Inventory and Analysis (FIA) occurrence data for Fagus and Tsuga.Model fit was evaluated with the Deviance Information Criterion (DIC), which is the sum of P D (the effective number of parameters) and the posterior mean of the deviance.To facilitate model comparison, ΔDIC was also calculated, where the model with the lowest DIC has a value of zero and all other models are compared to it.

Table 4 .
Model performance as measured by Area Under the Receiver Operating Curve (AUC), false negative rates (FNR), and false positive rates (FPR) for the non-spatial model, spatiallyvarying intercept (SVI) model, SVI plus climate, and multilevel B-spline fit to modern Fagus and Tsuga occurrence data from the Forest Inventory and Analysis (FIA) data.Predictions of the models for modern time were validated with a 10% hold-out dataset from the FIA data.Hindcasts were validated with data from the fossil pollen record provided by the Neotoma database using the "high" pollen thresholds for both genera.The numbers behind the AUC, FNR, and FPR values in parentheses for the Bayesian models represent the standard error calculated from 1000 random draws from the post burn-in MCMC iterations.For the FIA multilevel Bspline approximation there is no standard error as there were no MCMC iterations to draw from.