Ecological forecasts reveal limitations of common model selection methods: predicting changes in beaver colony densities.

Over the past two decades, there have been numerous calls to make ecology a more predictive science through direct empirical assessments of ecological models and predictions. While the widespread use of model selection using information criteria has pushed ecology toward placing a higher emphasis on prediction, few attempts have been made to validate the ability of information criteria to correctly identify the most parsimonious model with the greatest predictive accuracy. Here, we used an ecological forecasting framework to test the ability of information criteria to accurately predict the relative contribution of density dependence and density-independent factors (forage availability, harvest, weather, wolf [Canis lupus] density) on inter-annual fluctuations in beaver (Castor canadensis) colony densities. We modeled changes in colony densities using a discrete-time Gompertz model, and assessed the performance of four models using information criteria values: density-independent models with (1) and without (2) environmental covariates; and density-dependent models with (3) and without (4) environmental covariates. We then evaluated the forecasting accuracy of each model by withholding the final one-third of observations from each population and compared observed vs. predicted densities. Information criteria and our forecasting accuracy metrics both provided strong evidence of compensatory density dependence in the annual dynamics of beaver colony densities. However, despite strong within-sample performance by the most complex model (density-dependent with covariates) as determined using information criteria, hindcasts of colony densities revealed that the much simpler density-dependent model without covariates performed nearly as well predicting out-of-sample colony densities. The hindcast results indicated that the complex model over-fit our data, suggesting that parameters identified by information criteria as important predictor variables are only marginally valuable for predicting landscape-scale beaver colony dynamics. Our study demonstrates the importance of evaluating ecological models and predictions with long-term data and revealed how a known limitation of information criteria (over-fitting of complex models) can affect our interpretation of ecological dynamics. While incorporating knowledge of the factors that influence animal population dynamics can improve population forecasts, we suggest that comparing forecast performance metrics can likewise improve our knowledge of the factors driving population dynamics.

Introduction

Creating and validating predictions is an essential component of the scientific process. Numerous publications have advocated for ecologists to adopt a bolder scientific approach by integrating prediction as a fundamental element of ecological research, due in large part to increased societal demand to predict how ecological systems will respond to environmental change (Clark et al. 2001, Evans et al. 2012, Mouquet et al. 2015, Petchey et al. 2015, Pennekamp et al. 2017, Houlahan et al. 2017, Dietze et al. 2018, Maris et al. 2018). Central to the calls for a more predictive ecology is the establishment of a standard framework that can guide research and decision‐making efforts (Dietze 2017). Broadly speaking, this predictive framework should involve explicit descriptions of the underlying assumptions and uncertainties of predictions (Beale and Lennon 2012, Gregr and Chan 2014), recognizing the contextual limitations of predictions (Mouquet et al. 2015), and identifying specific validation criteria for predictions (Rykiel 1996, Dietze 2017).

While prediction is not new to ecology (Mouquet et al. 2015), validating ecological predictions has been practiced sporadically. Instead of direct empirical validation, ecologists have largely relied on alternative tools to assess the explanatory and predictive quality of their models. For instance, goodness‐of‐fit measures, such as the coefficient of determination (R 2), are frequently used to assess the amount of variation in the data that a model explains. Goodness‐of‐fit is a useful diagnostic, but it cannot be used to select predictive models from a set of candidates since the fit of a model always increases with the addition of covariates (although this may be moderately accounted for by using an adjusted R 2). Information criteria (e.g., AIC, BIC, DIC) are statistical tools useful for correcting goodness‐of‐fit by including a term that penalizes increases in model complexity (as measured by the number of parameters and sample size; Burnham and Anderson 2002, Clark et al. 2020). The standard model selection framework using information criteria involves fitting competing models to the same data and then assessing these models with information criteria to determine the most parsimonious model. Information criteria are designed to select models with the greatest predictive accuracy (Aho et al. 2014), though different criterion will use different measures of predictive ability (Gelman et al. 2014). Ecologists often use information criteria because they are considered computationally efficient approaches for making relative comparisons of model accuracy. However, researchers infrequently conduct empirical assessments of the goodness‐of‐fit or predictive accuracy of the ‘best’ model(s) selected using information criteria (Mac Nally et al. 2018). There is growing recognition that whenever researchers are interested in making inferences or projections beyond the data used to fit the model, assessments of model transferability, across space or time, should be performed (Wenger and Olden 2012).

When ecologists do validate predictions they tend to use vague or qualitative assessments on in‐sample data rather than quantitative assessments on out‐of‐sample data (Dietze and Lynch 2019). Proper model validation should instead involve (1) fitting models to within‐sample (i.e., training or parameterization) data, (2) making specific forecasts about out‐of‐sample (i.e., testing) data using the fitted models, and (3) evaluating the model forecasts by quantifying the error between the real and predicted observations (Dietze et al. 2018). Unfortunately, access to independent data is often lacking or unavailable (Yates et al. 2018). Due in part to the relative ease at which researchers can use spatial cross‐validation blocking methods that approximate data independence (Roberts et al. 2017), validations of spatial model forecasts (e.g., species distribution models) have been performed more frequently compared to other model classes (e.g., temporal data sets; Yates et al. 2018). Validations of temporal forecasts are lacking because there is often not enough data to adequately capture the complexity and heterogeneity of temporal dynamics (Harris et al. 2018). This is particularly true for forecasts of animal population dynamics (Chevalier and Knape 2020), which are inherently complex and temporally variable. Long‐term field research studies may therefore play a vital role in assessing the limitations of forecasting animal population dynamics.

Understanding the drivers of animal population dynamics is essential for developing useful population forecasts. It is generally accepted that both density‐dependent (e.g., territoriality, competition, and disease) and density‐independent (e.g., food resources, and environmental and demographic stochasticity) mechanisms influence wildlife population dynamics. Density‐independent factors limit population size by temporarily influencing the carrying capacity or causing a change in a demographic rate that is independent of population density (Lande et al. 2003, White 2008); whereas density‐dependent factors have specific effects that are influenced by the current population density, and which in turn influence the rate that a population tends toward equilibrium (i.e., regulation; Sinclair 1989, Turchin 1995, Sinclair and Pech 1996). Identifying how species are affected by density‐dependent and ‐independent factors will improve predictions of how populations will respond to drivers such as harvest and changes in environmental conditions. For instance, we expect that species that are regulated by density‐dependent mechanisms may be comparatively more buffered against temporary environmental perturbations (Băncilă et al. 2016) and less vulnerable to overharvest due to compensatory demographic responses (Boyce et al. 1999).

The degree to which animal populations are influenced by density‐dependent and ‐independent factors, or some combination of both, has been robustly debated for decades, in part due to technical limitations in our ability to estimate their relative influence from temporal population counts (particularly with estimating density dependence; Sinclair 1989, Dennis and Taper 1994, Koons et al. 2015). Previous large‐scale analyses have demonstrated that many populations exhibit density dependence (Sibly et al. 2005, Brook and Bradshaw 2006, but see Knape and de Valpine 2011). However, the predictive ability of these density‐dependent models is often limited (e.g., Glaser et al. 2014), suggesting that the predictive limitations are often due to more complex population dynamics (and noisy data). In response, statistical methods have been developed to adequately handle the complexity and noise within animal population data, spurring new efforts to disentangle how density‐dependent and ‐independent mechanisms affect populations across space and time (e.g., Rotella et al. 2009, Koons et al. 2015, Ross et al. 2015, Ferguson et al. 2017). These advances allow us to partition process from observation error, accurately estimate demographic parameters, and make forecasts from sophisticated population models. But more formal evaluations are needed to determine how well these new statistical models can predict future population dynamics.

Predicting how wildlife populations respond to density‐dependent and ‐independent drivers may be particularly important for species that are of special concern due to their rarity, economic or cultural importance, or outsized ecological role within local environments. One such ecologically important species is the beaver (Castor canadensis and C. fiber), an iconic ecosystem engineer whose abundance and distribution is increasing in North America and Eurasia following centuries of overexploitation. Beavers can substantially affect local environments (Rosell et al. 2005, Johnston 2017), including mitigating the effects of climate for fish and wildlife species (Hood and Bayley 2008) and increasing the habitat heterogeneity, species diversity, and species richness within beaver‐modified environments (Naiman et al. 1988, Wright et al. 2002, Rosell et al. 2005, Windels 2017, Willby et al. 2018). Indeed, as a direct result of the ecological role and ecosystem services that beavers provide, reintroducing beavers has been touted as a viable method for restoring ecosystems in many parts of Eurasia and western United States (Burchsted et al. 2010, Pollock et al. 2014, Stringer and Gaywood 2016, Law et al. 2017, Law et al. 2019, Willby et al. 2018). On the other hand, beaver dam‐building and foraging habits can be destructive to anthropogenic or economically valuable resources (Bhat et al. 1993, Jensen et al. 2001), and introductions of C. canadensis outside their natural range have caused extensive damage to South American ecosystems (Anderson et al. 2006, Anderson and Rosemond 2007). Despite their global range and extensive management history, surprisingly little is known about the landscape‐scale population dynamics of this large, charismatic rodent.

In the present study, we evaluated how density dependence and numerous extrinsic factors (i.e., external processes that affect the species such as predation and weather) influenced annual finite rates of change in beaver colony densities using data from 15 long‐term survey routes in Minnesota, USA. We also assessed the ability of common model selection methods (i.e., information criteria) to identify which model(s) had the best forecasting accuracy in estimating future colony densities using an ecological forecasting framework. Our specific objectives were to (1) estimate the strength of density dependence along beaver survey routes, (2) evaluate the relative influence that extrinsic covariates (forage availability, harvest, weather, and wolf [Canis lupus] density) had on colony finite rates of change, and (3) test the ability of information criteria to distinguish which model(s) can best predict future colony densities by comparing model predictions of future densities with observed densities. While not all the extrinsic factors evaluated here are always functionally density independent (e.g., predation and harvest), we modeled all covariates as density independent in the present study (i.e., without an interaction term with the previous year’s colony density). Typically, economic drivers are more influential in determining annual changes in harvest rates of furbearers than the density of the furbearer population(s), and there is currently no evidence to suggest that wolf predation of beavers is density‐dependent (Gable et al. 2018). Although beaver colony size can fluctuate spatially and temporally (Fryxell 2001, McTaggart and Nelson 2003), given the spatial and temporal scale of our study, using colony density was the only logistically feasible method for estimating beaver population size. There is also a strong precedent in the literature to suggest that colony density is an appropriate proxy for beaver population size (e.g., Hartman 1994, Busher and Lyons 1999, Jarema et al. 2009, Parker and Rosell 2014, Johnston and Windels 2015, Brommer et al. 2017, Ribic et al. 2017).

Methods

Study area

Our study area encompassed approximately the northern half of Minnesota, USA (Fig. 1) within the Laurentian Mixed Forest Province (Cleland et al. 2007). The study area is within the transition zone between temperate deciduous and boreal forest ecoregions, and the vegetative composition varies considerably (MNDNR 2017). Fire‐dependent oak (Quercus spp.) and jack pine (Pinus banksiana) forests are prevalent in the southern and western portions of the study area, while large swaths of black spruce (Picea mariana) bogs and tamarack (Larix laricina) swamps comprise portions of the western and northern sections. Mesic hardwood forests are common throughout the central and eastern sections of the study area, while coniferous forest communities are prevalent in the northeastern section. Human density varies widely throughout the study area, but most survey routes were in undeveloped areas.

Fig. 1

Map of the study area and location of each survey route. The Minnesota wolf population’s range expanded throughout the study’s time frame, as indicated by the range maps created from wolf surveys conducted in 1978–1979, 1988–1989, 1997–1998, and 2003 (no range expansion was found from 1998 to 2003). Results from the 1978–1979 survey indicated route 2 (Hay‐Kelliher) had established wolf packs and route 11 (Itasca) was undergoing recolonization, but these packs were not included in the official range maps. Wolves were not present for route 14 (Kanabec) surveys, which ceased in 1992.

We obtained relevant temperature and precipitation averages for our study from the PRISM Climate Working Group. Average annual precipitation across our study’s time frame (1972–2002) ranged from 616.2 ± 95.8 mm to 773.0 ± 141.4 mm at each beaver survey route, with an average of 66% of total precipitation falling during the growing season (May–September; PRISM Climate Group 2014). Average winter temperatures (December–March) were similar across all routes, ranging from –11.2° ± 2.3°C to –8.4° ± 2.2°C. Average maximum May temperatures (spring green‐up season) ranged from 18.3° ± 2.6°C to 20.0° ± 2.5°C.

Wolves are the main predator of beavers in our study area, and beavers constitute an important food source for wolves during the ice‐free season (Gable et al. 2017). Minnesota’s wolf population expanded during our study after being listed on the Endangered Species Act in 1974 (MNDNR 2001). After listing, the wolf population grew from an estimated low of 750 to approximately 2,450 by 1997–1998, extending their range by nearly 30,000 km2 (Erb and DonCarlos 2009), including recolonizing areas surrounding four beaver survey routes during our study (Fig. 1). Although black bears (Ursus americanus), coyotes (C. latrans), and bobcats (Lynx rufus) are also present within our study area (Hazard 1982), these species were not included in our assessment because there is no evidence to suggest predation rates from these predators can influence beaver population sizes (excluding insular systems; Smith et al. 1994).

Annual beaver colony surveys

The Minnesota Department of Natural Resources (MNDNR) conducted annual surveys by counting active beaver colonies from fixed‐wing aircraft along 25 standardized routes from 1975 to 2002. Following a data selection process (Appendix S2), we retained data from 15 of 25 routes with an average time series length of 22.3 yr (Appendix S1: Fig. S2, Table S1). Observers classified active colonies by identifying the presence of a visible food cache, which is the colony’s winter food source that consists of piles of semi‐submerged logs and twigs and can be seen in the fall just prior to freeze‐up (Payne 1981, Johnston and Windels 2015). Supplementary observations such as fresh mud on dams and/or lodges were also used to determine whether colonies were active. Surveys were conducted between 09:00 and 16:00 in various two‐ and four‐person fixed‐wing aircraft after leaf‐off, but before ice formed on water features (mid‐September–early November).

We digitized and calculated the length of each survey route in ArcGIS (version 10.5.1; Environmental Systems Research Institute 2017) from hand‐drawn maps used by MNDNR personnel as reference. Route lengths (range 55–336 km) and types were variable; three routes were flown in a series of linear transects, seven routes followed waterways exclusively (e.g., lake shores, rivers, streams), and the remaining five routes used a combination of transect and waterway segments (Fig. 1, Appendix S1: Table S1). We digitized each route by inferring the aircraft’s flight path (Appendix S1: Fig. S1), resulting in density estimates of the number of active colonies per kilometer surveyed by the aircraft.

Variable selection

Weather variables

Due to their unique ability to create and maintain ponds, we expected weather variables would have limited influence on fluctuations in beaver colony densities along survey routes. Nevertheless, previous studies have documented weather influences trends in colony density and juvenile survival rate and recruitment (Campbell et al. 2012, 2013, 2017, Ribic et al. 2017). Based on these studies and our belief that spring dispersal may be affected by water levels and connectivity, we included five weather variables in our analysis: (1) mean maximum temperature during the spring green‐up season (May), (2) growing season (May–September) drought index, (3) fall (August–October) drought index, (4) winter severity preceding fall surveys (mean December–March temperature), and (5) a spring (April–June) drought index. Temperature values were obtained from PRISM (PRISM Climate Group 2014) using the package prism (Hart and Bell 2015) in R (version 3.4.4; R Core Team 2018). We used monthly raster files at a 4‐km scale of resolution, averaged values across the entire route using the Zonal Statistics as Table tool in ArcGIS, and used a Python script to summarize multiple monthly PRISM raster files within ArcGIS. We used the Palmer Drought Severity Index (PDSI) values obtained from the U.S. drought portal (National Integrated Drought Information System 2018) to evaluate drought conditions. PDSI values provide a standardized index (range –7 to 7) for estimating the amount of water that is available for plants (Ribic et al. 2017); values <0 indicate drought conditions. Our study area covered three different PDSI climate divisions: North Central (2102), Northeast (2103), and East Central (2106). Routes that crossed multiple divisions (three routes) were assigned PDSI values corresponding to the division with the longest portion. We averaged (mean) all temperature and drought values across their timeframe of interest (e.g., the fall drought value was the mean of August, September, and October monthly values).

Forage quality index

Previous research has demonstrated that reductions in habitat quality can affect colony persistence (Busher and Lyons 1999, Fryxell 2001), so we expected forage quality to be positively correlated with finite rates of change in beaver colony densities. To evaluate the influence of forage quality, we developed an index of woody resource availability by quantifying the relative abundance of preferred forest types near water features along each route. We first delineated a 1‐km habitat buffer around each route, which corresponded to the 800‐m observer sight distance plus an additional 25% buffer to account for forest type characteristics of ponds that may have straddled or been near the sight distance boundary (Appendix S1: Fig. S1). Because beavers generally restrict their foraging to within 30–50 m of the riparian zone (Donkor and Fryxell 1999, Martell et al. 2006), we applied a second 50‐m forage buffer around all water features within the habitat buffer to isolate only those forest type characteristics that were available for beaver foraging (Appendix S1: Fig. S1). We extracted all stream features from the MNDNR hydrography data set (MNDNR 2014) and all lake and wetland features from the Minnesota National Wetland Inventory (NWI; MNDNR 2009); we selected only NWI features that consisted of unconsolidated bottom (i.e., open water) classes within lacustrine and palustrine systems. We used the 1992 National Land Cover Database (NLCD; Vogelmann et al. 2001) as our forest type input. Beavers are generally considered to be “choosy generalists,” preferring deciduous and early successional forest communities (Jenkins 1979, Curtis and Jensen 2004, Breck et al. 2012), so we defined preferred forage as Deciduous and Mixed forest classes to account for the full availability of deciduous trees. We took the total area of deciduous and mixed classes within the forage buffer, divided by the total area within the habitat buffer, to obtain a final index that approximately equates with the relative abundance of high‐quality forage within each route. Although beavers often deplete woody resources locally, we contend that at the spatial scale evaluated here using forest type data is adequate to provide a general index of woody resource availability.

Beaver harvest index

We expected harvest rates may have been high enough to negatively affect finite rates of change in colony densities, since as many as 170,000 beavers were harvested annually in Minnesota during our study (Appendix S1: Fig. S3). To estimate annual harvests, the MNDNR conducted mail surveys and multiplied the mean number of beavers harvested per respondent by the total number of licenses sold. No spatially explicit harvest data exists for our timeframe, only statewide estimates; we assumed that differences in regional harvest effort were consistent across the period of study. There was a limit of 10 beavers per trapper in 1975 and the harvest season was closed in 1976, but there was no harvest limit for beavers from 1977 to 2002. All routes were available to trappers excluding Kabetogama, where trapping ceased in 1975 when Voyageurs National Park was established (thus, the harvest value for all Kabetogama survey years was 0).

Wolf density estimates

We wanted to account for the influence of beavers’ most significant predator (wolves) in our study, but based on recent research (Gable and Windels 2018) we did not expect wolf predation to affect finite rates of change in colony densities along survey routes. We used wolf density estimates as an index to evaluate the influence of predation (akin to Ross et al. 2015). Because wolf densities generally increase linearly with available ungulate prey biomass (Fuller 1989, Fuller et al. 2003, but see Cariappa et al. 2011), we estimated annual wolf densities for each route by calculating ungulate biomass index (UBI) values (Kuzyk and Hatter 2014, Mech and Barber‐Meyer 2015) from available ungulate data. Although ungulates and beavers can sometimes compete for browse, beavers are able to adapt their forage preferences in areas of high ungulate densities, which can buffer the effects of competition (Hood and Bayley 2008). Our assumption in this index is that ungulate density would influence beaver populations via apparent competition, where ungulate abundance would increase wolf populations, which in turn would increase predation on beavers (Gable et al. 2018).

We obtained white‐tailed deer (Odocoileus virginianus) densities from MNDNR pellet survey estimates (Norton 2018; MNDNR, unpublished data), and moose (Alces alces) densities from MNDNR aerial survey estimates (Karns 1982, Lenarz 1998, 2006, Murray et al. 2006). We used the following regression equation presented in Mech and Barber‐Meyer (2015) to estimate annual wolf densitieswhere UBI was calculated by adding the density of deer (no./km2), plus six times the density of moose (no./km2; the number of white‐tailed deer “relative biomass equivalents” presented in Fuller et al. [2003]).

For the four routes that wolves recolonized during our study time period (Cass, Cass‐Crow, Itasca, Southern Pine; Fig. 1), we estimated wolf densities as a proportion of the UBI‐derived density for each year wolves were actively recolonizing the area. We used wolf population recovery data presented by Hayes and Harestad (2000) to estimate how wolf densities within each route reached their predicted densities within 6 yr of establishment. We then used the population estimates from Hayes and Harestad (2000) to estimate UBI‐derived density proportions for each year of the recolonization as follows: 0.12, year 1; 0.28, year 2; 0.52, year 3; 0.76, year 4; 0.84, year 5; and 1.00, year 6. We determined the likely year of initial wolf recolonization using data from a combination of annual scent‐post surveys (Sargeant et al. 2003; MNDNR, unpublished data) and extensive wolf population surveys from 1978 to 1979, 1988 to 1989, and 1997 to 1998 (MNDNR 2001), using the first year of wolf detection within 50 km of each route as the first year of recolonization. We acknowledge the first wolf detection near a route may have been a dispersing individual rather than an established pack, but we believe this method is adequate for estimating recolonization because we know the approximate year wolves colonized each of these four routes during our study time period based on findings from the population surveys.

Data analysis

Incorporating variable time lags

Natal dispersal is thought to be the primary mechanism for beaver population and range expansion (Baker and Hill 2003), so we selected our extrinsic variables and incorporated time lags into our analysis based on how we hypothesized each variable might affect rates of natal dispersal or recruitment. Although population density and harvest can delay the timing of dispersal up to 1–4 yr (Boyce 1981, Mayer et al. 2017), the majority of beavers disperse from their natal colony by age 2 or 3 (van Deelen and Pletscher 1996, Sun et al. 2000, McNew and Woolf 2005, Mayer et al. 2017). Therefore, we incorporated covariate time lags ranging from 0 to 3 yr into our statistical model (variable time lags are specified in Table 2). We assessed collinearity between covariates by calculating variance inflation factors (VIF). All covariates had VIF values <10, indicating multicollinearity was not an issue in our data set (Dormann et al. 2013).

Table 2

Mean and standard deviation (SD) of parameter estimates from the DDcov model.

Parameter	Interpretation	Mean	SD	BCI
a	density‐independent growth	–0.47†	0.09	(–0.66, –0.29)
b	density dependence	–0.64†	0.07	(–0.77, –0.50)
σa	standard deviation in population‐level density‐independent growth	0.28	0.07	(0.17, 0.46)
β1	annual harvest lag 1	0.02	0.02	(–0.01, 0.05)
β2	habitat quality	0.15	0.08	(0.00, 0.31)
β3	estimated wolf density lag 0	0.00	0.01	(–0.01, 0.01)
β4	mean winter temperature lag 0	–0.04†	0.02	(–0.09, –0.01)
β5	mean winter temperature lag 2	0.02	0.02	(–0.01, 0.05)
β6	maximum May temperature lag 2	0.01	0.02	(–0.04, 0.05)
β7	maximum May temperature lag 3	0.02	0.02	(–0.02, 0.05)
β8	spring PDSI lag 0	–0.05†	0.02	(–0.09, –0.02)
β9	growing season PDSI lag 2	0.01	0.03	(–0.05, 0.07)
β10	growing season PDSI lag 3	–0.01	0.03	(–0.07, 0.05)
β11	fall PDSI lag 2	–0.07†	0.03	(–0.13, –0.01)
β12	fall PDSI lag 3	–0.04	0.03	(−0.10, 0.02)

Negative Palmer Drought Severity Index (PDSI) parameter estimates indicate finite rates of change in beaver colony densities were positively correlated with drier seasons, as PDSI values <0 represent drought conditions. The significant negative winter temperature parameter estimate indicates colder winters were positively correlated with finite rates of change.

Effects where the 95% Bayesian credible interval (BCI) did not overlap 0.

Evaluating the influence of density‐dependent and‐independent factors

Because we were interested in evaluating the ability of information criteria to identify the best forecasting model(s), we sought to model beaver colony dynamics using a standard time‐series model in a Bayesian framework, rather than using a hierarchical model; evaluating hierarchical models (e.g., state‐space models) using information criteria is generally inadequate (Link and Sauer 2016). However, we also wanted to assess whether observation error could have significantly affected colony density estimates, and thus our conclusions about the relative strength of density dependence within and among survey routes (Turchin 1995, Freckleton et al. 2006). To do so, we fit both standard and state‐space models to our data (the latter of which can account for observation error) and compared colony densities from each model to assess whether observation error significantly influenced observed colony densities. Full details of this process can be found in Appendix S2. Our results indicated that observation error had an insignificant effect on colony density estimates and parameter estimates of environmental drivers, so we used the standard model framework in the subsequent analysis.

Based on diagnostic plots of the data, we modeled beaver colony dynamics using a model of contest competition, which describes the increasing utilization of available resources with increasing density (Hassell 1975). Our models described changes in the log density, X , = ln(N , /A), where N is the number of colonies observed along route i in year j, and A is the length (km) of route i. We applied the Gompertz model (Dennis and Taper 1994) using the following system of equations: where a is the log finite rate of population change, and b is a parameter for the strength of density dependence, such that b values further away from 0 are indicative of stronger compensatory (negative b values) or depensatory (positive b values) density dependence (Herrando‐Pérez et al. 2012). We also included a random effect to account for variation among routes in the density‐independent reproductive rate (u) that is not accounted for by covariates, and included the effects of environmental covariates (Z , ) on the density‐independent finite rate of change (see Table 2 for a full list of covariates). The variance term, ε , , accounts for unexplained inter‐annual variation in the density along route i in year j. This term is equivalent to the overdispersion term in the negative binomial distribution (Melbourne and Hastings 2008, Ferguson and Ponciano 2014). We note that, for these high colony counts (all >30), the continuous‐state probability density model that we applied is appropriate; however, at lower counts, a discrete‐state model, such as the negative binomial distribution, would be more appropriate to account for lattice effects (Henson et al. 2001).

Model selection using information criteria

We used the deviance information criterion (DIC) and widely available information criterion (WAIC; Watanabe 2010) to compare (1) the full empirical density‐dependent model described above (Eq. 2; DDcov model), (2) a density‐independent model with covariates (DIcov model), (3) a density‐dependent model without covariates (DD model), and (4) a density‐independent model without covariates (DI model). Models were fit using MCMC implemented by Just Another Gibbs Sampler (Plummer 2003) by making 106 draws from the posterior distribution. We thinned our resulting chain every 102 draw due to strong autocorrelation in some parameters (Christensen et al. 2011).

Evaluating model predictions of colony densities

We evaluated the predictive ability of each of our four models by withholding the final one‐third of observations (96 route‐year observations total, an average of 6.4 per route) for each survey route, fitting the models to this reduced data set, then hindcasting the withheld data. We assessed hindcasting accuracy using the root mean squared prediction error (MSPE) of the predicted density. The MSPE for site i is given by , where is the kth predicted density. We then summed the MSPEs across sites to get a total MSPE. In order to determine whether model inferences were consistent between the full data set and withheld data set, we compared parameter estimates from each data set using Deming regression implemented in the R package deming (Therneau 2018), which allows for errors in both dependent and independent variables.

Results

Model selection and evaluating model forecasts

The DDcov model performed best based on DIC, WAIC, and predictive performance (Figure 2, Table 1). However, the colony density prediction error for the DDcov model did not improve appreciably compared to the DD model. The average improvement in predictability across all routes was only 3% (minimum –52%, maximum 42%), with no qualitative differences among routes in density predictions between the DD and DDcov models (Fig. 3). The improvement in predictability was less than we expected, given that high ΔDIC (12.31) and ΔWAIC (14.14) values indicated strong support for the DDcov model (Table 2). We found a similar pattern in the density‐independent models. Information criteria provided stronger support for the DIcov model, but the DIcov model only had an average improvement in predictive ability of 7% relative to the DI model (Table 2). Both density‐dependent models had substantially lower information criteria and MSPE values than the density‐independent models, providing support for density‐dependent dynamics in beaver populations. We found no systematic differences between the effects of the density‐independent covariates estimated from the full data set and the covariates estimated from the holdout data set. All posterior estimates of the covariates were within 1 standard deviation of the one‐to‐one line that showed equal estimates (Appendix S1: Fig. S4), indicating there were no significant changes in the ecological processes regulating beaver populations between our full and holdout data sets. The slope of the line that best explained the relationship between these estimates was 1.24 (SE = 0.19).

Fig. 2

Composite image of observed data (black dots) and model fits for each route. The blue line is the density‐dependent model with covariates (DDcov model), while the black line is the density‐dependent model without covariates (DD model). The number in the upper right corner of each plot corresponds to the route pictured in Fig. 1. Note that the y‐axis limits are different for each route to highlight route‐specific trends. [Correction added on December 09, 2020, after first online publication: Figure 2 has been updated]

Table 1

Comparison of our four models evaluating the influence of density‐dependent and density‐independent factors on finite rates of change in beaver colony densities.

Model	ΔDIC	ΔWAIC	MSPE
Density dependent with covariates (DDcov)	0.00	0.00	1.72
Density dependent without covariates (DD)	12.31	14.14	1.92
Density independent with covariates (DIcov)	86.51	92.54	7.11
Density independent without covariates (DI)	91.22	108.10	7.57

Deviance information criterion (DIC), widely applicable information criterion (WAIC), and the total root mean squared prediction error (MSPE) values are shown for each model. Our results indicate the DDcov model explained the greatest amount of variation and had the lowest MSPE values for predicting beaver colony densities in the holdout data set; however, MSPE values were similar between the DD and DDcov models, indicating an approximately equal ability to predict future colony densities.

Fig. 3

Comparison of prediction errors between the DD model and DDcov model for each of the 15 survey routes. The gray line denotes no difference between the two models. As evident in the figure, prediction errors between the two models are relatively similar.

Parameter estimates for density‐dependent and ‐independent factors

Average estimated strength of density dependence across all routes was  = −0.64 (SD = 0.07, Bayesian credible interval based on 95% of the highest posterior density [BCI] = −0.77 to −0.50; Table 1). Average log finite rate of population change (density‐independent growth) was  = –0.47 (SD = 0.09, BCI = –0.66 to –0.29), with an average variation in the population‐level rates of change of (SD = 0.07, BCI = 0.17–0.46; Table 1). We also found significant geographical variation in the finite rate of population change throughout our study area; details on the methods and results from this analysis can be found in Appendix S3.

Of the 12 covariates we evaluated in the DDcov model, three had posterior probability densities that did not overlap 0 at the 95% BCI (Table 1). Average winter temperature during the same year (i.e., lag 0) was negatively correlated with finite rates of change in colony densities (β4 = –0.04, SD = 0.02, BCI = –0.09 to –0.01). Both spring (lag 0; β8 = –0.05, SD = 0.02, BCI = –0.09 to –0.02) and fall PDSI values (lag 2; β11 = –0.07, SD = 0.03, BCI = –0.13, –0.01) were negatively correlated with finite rates of change in beaver colony densities, indicating colony densities increased in response to both current and past drought conditions (or, alternatively, decreased in response to high‐water conditions; PDSI values < 0 indicate drought) (Fig. 4). All other covariates did not have a significant influence on finite rates of change in colony densities (Table 1). Based on these results, we also conducted a post‐hoc analysis where we fit our data to an intermediate density‐dependent model with only the significant environmental parameters. This limited covariate model performed better than the DD model but worse than the full DDcov model according to information criteria (ΔDIC = 1.91, ΔWAIC = 4.14), but had a nearly identical predictive performance to the DD model (total MSPE = 1.89).

Fig. 4

Relationship between the log finite rates of population change in beaver colony densities and the three significant environmental covariates, (a) average winter temperature (lag 0), (b) fall Palmer Drought Severity Index (PDSI; lag 2), and (c) spring PDSI (lag 0). All covariates were negatively correlated with finite rates of change, indicating beaver colony densities tended to decrease following warmer winters and wetter seasons.

Discussion

By using an ecological forecasting approach, we were able to evaluate each model’s ability to predict the relative contribution of density‐dependent and density‐independent factors on future beaver colony densities. Overall, our results indicate that density‐dependent mechanisms and, to a lesser extent, density‐independent weather variables (winter temperature [lag 0], spring drought [lag 0], fall drought [lag 2]) significantly influenced inter‐annual fluctuations of beaver colony densities at the landscape scale (Table 1). However, despite information criteria indicating strong support for the density‐dependent model with covariates, the predictive ability of the model was only slightly better than the density‐dependent model without covariates likely due to the model over‐fitting our data. Over‐fitting of the complex model is supported by the perplexing direction of the significant weather variables (i.e., that beaver populations respond positively to drought conditions and colder winters), which are contrary to our current understanding of beaver biology. We speculate that a possible biological explanation for these results could be elevated rates of natal dispersal in response to the stressful environmental conditions; however, our hindcasting results suggest that the density‐independent weather variables are nonetheless relatively unimportant for predicting future colony densities.

The role of density‐dependent and ‐independent factors in beaver colony dynamics

In addition to our hindcasting results, our estimate for the average strength of density dependence across all routes (  = −0.64 [SD = 0.07, 95% BCI: −0.77 to −0.50]) provides compelling evidence that density‐dependent mechanisms significantly influenced finite rates of change in beaver colony densities. Several mechanisms may be responsible for driving landscape‐scale density dependence among beaver colonies; for instance, density has been shown to affect the timing of natal dispersal (Mayer et al. 2017), the age at first breeding (Mayer et al. 2017), and rates of mate change in beavers (Mayer et al. 2017). Beaver scent‐marking behavior and intraspecific aggression also act to regulate territory size and distance between colonies (i.e., colony density; Aleksiuk 1968, Bergerud and Miller 1977, Müller‐Schwarze and Heckman 1980, Rosell and Nolet 1997). Territoriality between beaver colonies may therefore function as a mechanism of site‐dependent regulation (Rodenhouse et al. 1997), causing average fecundity rates to decrease as population density increases (as previously found by Payne 1984, Pietrek et al. 2017). However, a recent study suggests levels of conspecific aggression may actually be inversely density‐dependent as a result of delayed dispersal at high densities (Mayer et al. 2020), which may result in lower rates of change at higher densities.

Identifying density dependence as an important mechanism of beaver colony dynamics indicates beaver populations may be buffered against short‐lived environmental fluctuations (Ives 1995), presumably due to compensation (as previously found by Boyce 1981) and their ability to modify and control their own environment to a degree that few other animals can. It is possible that weather may have negatively affected beaver demography as indicated in previous studies (Campbell et al. 2012, 2013), but those effects did not translate to a net effect on colony density due to changes in mortality rates at other life stages (Burnham and Anderson 1984, Lebreton 2005, Băncilă et al. 2016). For many long‐lived species, environmental variation has a greater effect on reproduction and juvenile survival relative to adult survival, which may translate to greater population stability (i.e., the “demographic buffering hypothesis”; Pfister 1998, Gaillard et al. 2000). Additionally, beaver dams can store and maintain areas of open water even during periods of drought (Hood and Bayley 2008), suggesting colonies are naturally insulated from short‐term environmental changes.

Our finding that wolf densities did not influence rates of change in beaver colony densities is consistent with recent research from northern Minnesota that demonstrated beaver populations can be resilient to intense wolf predation pressure (perhaps also due in part to compensation; Gable and Windels 2018). By using wolf density as an estimate of predation pressure we assumed wolf predation rates on beavers were directly related to wolf density, and there is no evidence to suggest this is a poor assumption (Gable and Windels 2018, Gable et al. 2018). But if wolf predation rates on beavers change in relation to wolf density (e.g., if per‐capita predation rates decrease at high wolf densities when ungulate prey is readily available), then we might not have fully captured the role predation has on beaver colony dynamics. Indeed, predation is a strong density‐dependent factor for many rodent populations, which can affect the periodicity and amplitude of population fluctuations (e.g., Erb et al. 2001, Hanski et al. 2001, Korpimäki et al. 2002). However, unlike many predators that are dependent upon rodents as their primary food source, wolves primarily prey on ungulates in boreal forest ecosystems (Fuller et al. 2003). Though beavers are a seasonally important prey item for wolves (Gable et al. 2017, 2018), our analysis suggests wolves are unable to suppress beaver populations in multiple‐prey systems. Little is known about the functional and numerical relationships of wolf–beaver–ungulate systems (Latham et al. 2013, Gable et al. 2018), and future work is needed to determine how variations in wolf and ungulate densities affect wolf predation rates on beavers.

We suggest that because finite rates of change in beaver colony densities appear resilient to moderate exploitation, regulated trapping seasons can be compatible with management programs that desire increases in beaver population sizes. For example, our results indicate that once beaver populations become established through reintroduction or translocation programs, beavers are likely to persist even if subjected to moderate levels of removal (managed harvest or nuisance control). On the other hand, the resilience of beaver populations may be problematic in areas where beavers are unwanted. Given the ongoing trend of declining trapper participation and average pelt prices, we expect beaver populations may increase in areas where public harvests have previously limited beaver populations (e.g., in easily accessible areas), and will generate more conflicts with anthropogenic and natural resources (e.g., coldwater stream systems; Johnson‐Bice et al. 2018). Sustained management efforts will likely be required in areas where beaver engineering causes conflicts or undesirable ecosystem modifications. This is probably not new information for those accustomed to mitigating beaver conflicts, but our results now quantitatively demonstrate the stability of beaver colony densities at the landscape scale.

Ecological forecast assessment and implications for model selection using information criteria

Our study provides an informative example on the importance of empirically validating ecological models. Failing to validate our models would have led to different biological conclusions about the importance of environmental drivers on beaver colony dynamics. Using a hindcasting framework to validate our models revealed two important aspects of predicting beaver colony densities: (1) models that included a parameter for density dependence predicted future colony densities substantially better than density‐independent models and (2) models that included the effects of covariates were only marginally better at predicting future colony densities than comparable models that did not include covariates (Table 1). Our forecasting results therefore indicated that density dependence was clearly important in determining beaver colony densities but that the extrinsic covariates evaluated were relatively insignificant, the latter of which is inconsistent with a standard interpretation of the information criteria results. While incorporating expert biological knowledge of the factors that influence animal population dynamics can improve population forecasts, our analysis demonstrates that the reverse can also be true: comparing validations of population forecasts can improve our knowledge of the factors driving population dynamics.

One of the key assumptions often over‐looked by ecologists using information criteria for model selection is that, when comparing multiple models, the model selection process is designed to identify which model(s) has the greatest predictive accuracy (Aho et al. 2014). This assumption is based on the fact that numerous information criteria like the AIC, DIC, and WAIC, can be interpreted as measures of prediction error for new data points (Gelman et al. 2014; e.g., the AIC is considered asymptotically equivalent to leave‐one‐out cross‐validation; Stone 1977). However, of the few studies that have empirically tested this assumption using population models, most have found that information criteria are not reliable substitutions for empirical measures of predictive error (Link and Sauer 2016, Link et al. 2017, Clark et al. 2020, Wood et al. 2020). Our analysis is consistent with these previous studies. Other critical evaluations of information criteria have found that the predictive performance of model selection using different criteria can be highly variable, depending on the complexity of the model set (Ferguson et al. 2019) or the degree of heterogeneity in the data (Brewer et al. 2016). Thus, information criteria appear well suited for identifying models that best fit the observed data but may not be reliable for identifying models that best predict out‐of‐sample data (i.e., a lack of generalizability; Clark et al. 2020). The lack of generalizability may often be due to over‐fitting (Wenger and Olden 2012), as we believe occurred here. Over‐fitting likely stems from penalty terms that are too conservative to account for increases in model complexity (Gregr et al. 2018, Clark et al. 2020), the inability of information criteria to account for uncertainty in the measured covariates (Dietze et al. 2018), or uncertainties associated with jointly performing parameter estimation and model selection (Harrell 2015).

Although the most complex model (DDcov) had the lowest forecasting error, the average improvement in forecasting future colony densities was only 3% better than the much simpler DD model (Table 1). Information criteria are often favored by ecologists because they are considered effective measurements of model parsimony (Aho et al. 2014). However, based on results from our study and other recent evaluations (Clark et al. 2020, Wood et al. 2020), we suggest that model selection using information criteria does not always identify the most parsimonious model(s) either. There remains considerable debate about whether ecologists should be using complex or parsimonious models to analyze data (e.g., Evans et al. 2013, Coelho et al. 2019), but the justification for using complex or parsimonious models largely depends on whether the analysis is structured around mechanistic or correlative models (Coelho et al. 2019). If the latter, like in our analysis, then ecologists should strive for model parsimony. Indeed, using a variety of validation approaches, researchers have found simpler models (or models of intermediate complexity) generally perform as well as, or better than, more complex models in predicting out‐of‐sample data (e.g., Ward et al. 2014, Tredennick et al. 2016, Chevalier and Knape 2020, Clark et al. 2020, Wood et al. 2020). But approaches such as borrowing information across sites (when possible) appear to improve the predictive ability of complex models, presumably by reducing over‐fitting (Chevalier and Knape 2020).

Our objective here is not to completely discourage the use of information criteria for prediction but rather to highlight critical limitations with its application. We (and others) have shown that model selection using information criteria does not always lead to the selection of the most parsimonious model, nor the model with the best predictive accuracy. However, information criteria may provide a complementary analytical framework to more formal model evaluations, such as by helping to select important model components that are needed to predict future ecological dynamics (Ferguson et al. 2016) or by reducing the number of competing models to a manageable subset that can be used for testing against out‐of‐sample data (Houlahan et al. 2017). But other emerging forecasting tools may nonetheless provide a more nuanced evaluation of model performance than information criteria. For instance, variance partitioning techniques (i.e., decomposing sources of prediction uncertainty; Dietze 2017) or metrics like the ecological forecast horizon (Petchey et al. 2015) can provide direct assessments of predictive power that can be used to improve predictive models (Pennekamp et al. 2019).

Long‐term studies like ours have been historically important in advancing population ecology and will continue to play a large role going forward (Reinke et al. 2019). Despite the declining abundance of, and funding for, long‐term field studies, they remain essential for fully capturing the temporal heterogeneity of ecological systems and will likely yield unique insights into patterns of model transferability (e.g., by identifying whether there may be temporal variation in the transferability of models; Kleiven et al. 2018, Soininen et al. 2018). Unfortunately, our study is one of many long‐term programs that have been discontinued. Our analysis demonstrates that these historical data sets can nonetheless still provide technical and biological insights into modern ecological issues. Moreover, our study presents one of the first evaluations of forecasting mammalian population dynamics, as most previous efforts to date have been conducted on other taxa (Chevalier and Knape 2020). Broadening evaluations of population forecasts to lesser‐represented species that exhibit variable life‐history strategies may lead to a greater understanding of model forecasting performance specifically (Chevalier and Knape 2020), and the population ecology of species generally.

Conclusion

Standard use and interpretation of information criteria would suggest that our most complex model (density‐dependent with covariates) should have the best ability to predict future beaver colony densities. However, we found that the predictive accuracy of a much simpler model (density‐dependent without covariates) was nearly identical to the complex model. We also showed that, across all performance metrics evaluated, the density‐dependent models performed markedly better than the density‐independent models. Thus, we conclude that while we can reasonably predict the return rate of beaver populations to their carrying capacity through the influence of density‐dependent dynamics, we were unable to adequately predict how environmental factors temporarily drive populations away from carrying capacity, even when the future environmental states could be perfectly predicted. By using an ecological forecasting framework, we revealed how a known, but often‐overlooked limitation of information criteria (over‐fitting of complex models) can affect our interpretations of ecological dynamics. Our analysis clearly demonstrates the importance of validating ecological model predictions using long‐term data sets, highlighting yet another reason why long‐term field research studies are integral to ecology.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.