One of the principal goals of wildlife research and management is to understand and predict relationships between habitat quality, health of individuals and their ability to survive. Infrequent sampling, non-random loss of individuals due to mortality and variation in capture susceptibility create potential biases with conventional analysis methods. To account for such sampling biases, we used a multi-state analytical approach to assess relationships between habitat, health and survival of grizzly bears Ursus arctos horribilis over a 10-year period along the east slopes of the Canadian Rockies in Alberta, Canada. We defined bear health states by body condition estimated from the relationship between weight and body length. We used a sequential model building process to first account for potential sampling biases, and then explored changes in body condition relative to habitat use and survival. Bears that used regenerating forest habitats (mostly due to forest harvesting) containing a diversity of age classes were more likely to see gains in their body condition, whereas bears that used older forests were more likely to see reductions in body condition. Survival rate was reduced most by road densities which in turn were positively correlated with regenerating forest habitat. Human activities which promote young regenerating forests, such as forest harvesting, therefore promotes improved health (increased body condition) in bears, but are offset by reductions in survival rates. Multi-state analyses represents a robust analytical tool when dealing with complex relationships and sampling biases that arise from dynamic environments.
One of the main challenges in wildlife research is the need to better understand how habitat quality affects the health and survival of individual animals. It is often assumed that selected habitats are optimal in meeting an individuals or species' biological needs and thus result in higher levels of fitness and survival (Krebs 1980). Human alteration of a habitat can negatively affect survival, despite being a selected habitat with potentially high rates of fitness (Nielsen et al. 2006, Nielsen et al. 2008). It is difficult, however, to assess the relationships between wildlife health, habitat and survival or to relate the results of disparate analyses between linkages in a single analysis. There are at least three major challenges in assessing health, habitat and survival. First, measuring the change in health often requires live captures, which are often periodically sampled leading to uneven intervals between measurements. Second, correlations between health and survival may result in non-random dropout of individuals. This, in turn, may bias the estimation of model parameters. Finally, differences in capture susceptibility among individual animals (which may also be influenced by health) may result in biased representation of individuals.
We were interested in relating environmental factors to changes in grizzly bear Ursus arctos horribilis body condition and survival. Grizzly bears, like many other wildlife species, have specific habitat requirements but also traverse areas of differential mortality pressure, which is often due to anthropogenic activity (Benn & Herrero 2002, Nielsen et al. 2004b), in the search of seasonally available food items. Multi-state models (Brownie et al. 1993, White et al. 2006) provide a method that is robust to many of the sampling issues that confronted our analysis. In multi-state models, animals are first subdivided into relative states based upon health, disease or other biological attributes. The model is then used to estimate the probability that an animal will change states as well as its survival rate in a given state. The transition probabilities can be constrained to be a function of other health, environment or temporal (climate) variables. Multi-state models have been used to assess costs of reproduction (Nichols et al. 1994, Nichols & Kendall 1995, Schwartz & White 2008), disease transmission (Conn 2008), evolutionary trade-offs (Nichols & Kendall 1995) as well as movement rates of animals between different geographic areas (Hestbeck et al. 1991, Williams et al. 2002).
In our study, we used multi-state modeling to investigate relationships between grizzly bear body condition and habitat quality, including habitat features influenced by human activity, and to explore trade-offs between body condition and survival across a gradient of habitat quality and anthropogenic activity. Of most interest was the testing of whether habitat quality or anthropogenic features most affected changes in bear condition, and whether changes in body condition affected a bear's ability to survive. We also examined how method of sampling affected estimates of body condition and survival. While this study focuses on a single species, i.e. the grizzly bear, the methods are widely applicable to other wildlife species found in complex and dynamic ecosystems.
Study area
The topography of our study area varies from plains and foothills to high alpine areas (Fig. 1). The majority of roads are found in the foothills in association with settlements and forestry-related habitats. A history of forestry, mining and oil and gas development has created a mosaic of forest types and stand ages, as indicated by regenerating forest habitats and an array of permanent road networks. Our study area also included federal and provincial parks and protected areas such as the Jasper National Park and Whitehorse Wildlands Park in which anthropogenic changes in habitat are uncommon and human access features are fewer in number. Historic and contemporary forest fires are also common to the landscape leading to a diverse mosaic of stand ages (Nielsen et al. 2002, 2006).
Methods
Capture, collaring and body-condition measurements
We captured and collared grizzly bears during 1999-2010 using either helicopter aerial darting, leg-hold snares or culvert traps (Cattet et al. 2003a, 2003b). Aerial captures typically occurred in subalpine habitats in the early spring or in forestry cut blocks where helicopter operations were possible. Most snare sites, on the other hand, were at lower elevations and normally within 100 m of a road. Beginning in 2006, capture efforts (2006-2010) were focused on the use of culvert traps and helicopter aerial darting, with the use of foot snares phased out (Cattet et al. 2008). Over the course of the study, some bears were targeted for recapture in order to replace GPS collars. This resulted in multiple sampling events for some bears. Simply by chance, some bears were also recaptured and sampled more than once a year. We fitted Global Positioning Systems (GPS) radio-collars from Televilt Simplex and Tellus (Lindesberg, Sweden; 1999-2010) and Advanced Telemetry Systems (ATS) (Isanti, Minnesota USA; 1999-2001) on captured bears, which were programmed to acquire a location every 1-4 hours. In addition, we fitted very high frequency (VHF) ear-tag transmitters (ATS) on all captured bears (Graham et al. 2010). During capture, we measured straight-line length and body mass for each bear, as well as other measurements, and we collected biological samples. These were converted into a body-condition index which is the standardized residual from the regression of body mass against body length (Cattet et al. 2002). All capture efforts followed guidelines by Canadian Council on Animal Care (2003) and the American Society of Mammalogists (Gannon et al. 2007), and were approved annually both by the University of Saskatchewan's Committee on Animal Care and Supply and the Alberta Department of Sustainable Resource Development Animal Care Committee.
We defined habitat quality and anthropogenic variables from remote sensing-based land cover mapping and databases of anthropogenic footprints (McDermid et al. 2005, 2009). The main habitat covariates considered were variation in crown closure which is a general index of the diversity of habitat types that bears traversed. Presence of regeneration habitat, which was related to forestry and burns, was also used to describe habitats that potentially contained berry resources as well as higher ungulate densities (Nielsen et al. 2008, Latham et al. 2011, Serrouya et al. 2011). For canopy closure and regeneration, we considered mean and standard deviation values. The variation in canopy closure and regeneration habitat reflected both the relative amount of habitat, but also the variation in values therefore providing an index of overall habitat diversity (Nielsen et al. 2004a). We also considered forest age as an index of the relative seral stage of forested areas.
We summarized habitat covariates at a buffer scale equivalent to a bear's daily distance moved for each age, sex class and season (Table 1). This scale captured the availability of forage as well as ungulate species in areas that the bears traversed. For hypophagia and early hyperphagia (dates prior to July 31), we estimated distinct movement rates for adult males, females and subadults and females with cubs, with similar sex-specific rates in the fall (late hyperphagia; see Table 1; John Boulanger, unpubl. data). We preferred this approach rather than an analysis at a home range given that habitat in the immediate area of GPS locations best described the actual encounter rate of bear with resources. The main anthropogenic covariate considered was road density which was the kilometers of road encountered within a 300-m radius of each GPS location. This scale was chosen since it described the actual areas that bears traversed and the relative risk that they were directly exposed (Archibald et al. 1987). Mean elevation and the terrain ruggedness index (Riley et al. 1999) within the same 300-m radius were used to represent general topographic conditions. The average of each covariate for the duration a bear was collared was then summarized as an individual covariate to provide a general index of habitat, topographic and anthropogenic factors that we hypothesized would affect the ‘state’ of the bear's body condition.
Defining multi-state conditions
We classified bears into body-condition states according to the value of their body-condition index (BCI). Only ‘spring’ captures that occurred before July 31 were considered in this analysis since body mass changes rapidly during hyperphagia and the fact that most captures occurred during the spring of the year. We used mass and length measurements from the first capture of a bear within each year to estimate BCI. These were classified into three body-condition states (low, medium and high) using the 33th and 66th percentile thresholds to define low (< 33rd percentile), medium (≥ 33rd percentile and ≤ 66th percentile) and high (> 66th percentile) condition states based on overall BCI measurements. In addition to measures of body condition, we also compiled an annual list of known mortalities of study animals.
We defined multi-state models with live and dead recoveries using program MARK (White & Burnham 1999). Multi-state models estimate survival rate (S), recapture rate (p), transition probabilities between condition states (Ψ) and reporting rate (r), which is the probability that a bear which dies is reported. We were particularly interested in transition probabilities, which are the probability that a bear will change from one state to another state in the interval between yearly captures. In a full multi-state model, transition probabilities are estimated for a bear moving from one state to each of the other states, and for bears from the other states to transition to the given state (Brownie et al. 1993). We used a logit-link for all analyses. Simulated annealing methods (Goffe et al. 1992) were used to check final models for optimal convergence. Program U-CARE (Choquet et al. 2009) was used to test for goodness-of-fit of the live-capture component of the multi-state models (Pradel et al. 2005). Information-theoretic model selection methods (Burnham & Anderson 1992) were used to evaluate relative support of candidate models.
Multi-state model building
We built models sequentially with initial models describing base biological, survival and sampling variation, followed by variation in the change in body-condition state, and finally whether changes between body-condition states influenced survival. This approach was similar to an analysis of covariance (Milliken & Johnson 2002) in which a base model is first formulated to control for dominant sources of variation in the data therefore maximizing the power to detect relationships between changes in body condition and survival.
Accounting for dominant sources of variation
Models were initially built to account for the effect of sampling on bear recapture rates and reporting rates. Bear capture effort was not uniform across our study area. Therefore, we calculated a covariate for capture effort based on the number of days per year in which live-capture snares were set in the home-range area of a given bear. This covariate was then entered for each bear year in the analysis.
Our goal during this time period was to try to distribute GPS collars across watershed units within our study area in order to determine habitat use in a heterogeneous landscape. Once a bear was collared within a watershed unit we would stop capture efforts there and focus on a new watershed. Thus, the number of days of capture effort was not equally distributed, but it was rather influenced by the length of time it took to catch a bear within a specific area. Aerial capture efforts were focused in areas of open habitat where helicopter pursuit was deemed safe and practical. In addition, there were four bears (out of 111) that only had (VHF) radio ear-tags for a proportion of the time they were monitored usually in cases where a bear was too small or too large for a collar. This potentially affected the reporting of bears that died, and therefore covariates that described the proportion of time in which a bear had only a radio ear-tag were developed under the assumption that bears with only radio ear-tags would have a lesser probability of being reported given that their locations were less certain than GPS radio-collared bears.
The base biology and sampling models attempted to explain variation in survival (based on sex and known covariates that influence survival) and variation in recapture rate based upon capture effort in bear home ranges each year. Other studies of bear survival rates (McLellan et al. 1999, Nielsen et al. 2004b) have suggested that age, sex and the presence of roads strongly affect survival, and therefore, we initially considered the effect of each of these covariates on survival. In addition, it was possible that bears captured in early spring one year and later in the spring in subsequent years may have shown an increase in BCI due to the timing of capture relative to den emergence and the amount of time that bears would have been feeding since den emergence (when BCI was likely lower). To account for this temporal factor, we tested the relative change in Julian day of capture for a bear captured in more than one year. Relative change was the difference between the Julian day of capture for a bear in a given year and the mean Julian day of capture across all years. Using this approach, the relative change covariate was set at 0 for years in which bears were not captured or for bears captured only once.
The base biology model building phase also tested whether simplifying assumptions could be made about changes in the body-condition states. In particular, we tested whether change from one state to another had its own unique probability, or whether it could be assumed that some probabilities were similar. For example, it might be possible that the probabilities of changing from one state to the next (e.g. low to medium, medium to high or vice versa) might be similar whereas probabilities of changing across two states (low to high or high to low) might be different. We built a number of models with the aim of finding the most parsimonious base model for change in body condition.
Assessing the effect of covariates on changes in body condition
The next phase of model building considered the effects of habitat quality and anthropogenic variables on changes in body condition. The most supported base model was used as a starting point. A large number of potential relationships between covariates and the state change in body condition were possible, and therefore, an essential part of the analysis was formation of biologically-based models as a means of optimizing the number of candidate models and sharpening the biological applicability of the analysis (Table 2). For example, variables with a ‘+’ were predicted to influence transition to a higher body-condition state, whereas variables with a ‘-’ were predicted to influence transition to a lower body-condition state. Covariates that described factors that influence higher density of ungulates such as variation in canopy closure, regeneration habitat and variation in regeneration habitat (Nielsen et al. 2010) were hypothesized to increase body condition. In contrast, factors that reduced habitat diversity such as forest age and elevation were hypothesized to negatively affect body condition. Significant correlations existed between road density, regeneration habitat, elevation, canopy closure and terrain ruggedness given that most roaded areas (with regeneration habitat) occurred in lower to mid elevation areas. We therefore entered correlated covariates individually for any given parameter to avoid confounding and showed results graphically to further illustrate the relationship between covariates (R Development Core Team 2009). We included the number of captures as Cattet et al. (2008) showed that bears captured more than once had a lower probability of increasing body condition. This effect was modeled by using number of captures as a covariate for state change in body condition with the prediction that number of captures should reduce the probability of changing to a higher BCI state. See Appendix I for more summaries of the covariates used in the analysis.
Relationships between body-condition states and survival
Once the relationships between change in body condition were defined for biological and environmental factors, we assessed how these relationships influenced survival. We initially considered the general associations between covariates and change in body condition. For more supported variables we also estimated probability curves to describe changes in body condition across two states (i.e. low to high or high to low). After we identified the most supported model describing changes in body condition, we tested whether differences in survival were apparent across states and whether these differences could be attributed to our hypothesized factors listed in Table 2. We were also interested in whether the same predictor variables that were influencing change in body condition were also influencing survival rates in an opposite direction. This implies a trade-off between body condition and survival and supports a hypothesis that bears using actively managed landscapes in Alberta are effectively ‘caught’ in an ecological trap/attractive sink where resources are plentiful or highest, but survival is low (Nielsen et al. 2006, Northrup et al. 2012).
Results
In total, 199 captures between 1999 and 2010 from 111 individual bears were used in this analysis (the average number of captures per bear was 3.15; SD = 1.04, minimum = 1, maximum = 6). In 30 capture events, no body-condition measurements were taken and these records were censored from the analysis leading to 169 capture events with a BCI score. Of the 111 bears, 107 received GPS collars and four only received VHF transmitters. Of the 111 bears, 62 were females and 49 were males with the average age of bears at 7.5 years (SD = 4.7, minimum = 2, maximum = 21). During the study, we recorded 37 mortalities. Of these, 27 were human caused, eight were unknown and two were due to natural causes. Of the 37 mortalities, 20 bears had been previously captured with BCI measurements and were therefore included in the analysis. On average, 345 yearly GPS locations (SD = 640.7; range: 7-4,513) were obtained per bear. These were used to estimate habitat use by individual bears. Goodness-of-fit based on global multi-site and subcomponent tests in U-CARE suggested adequate model fit (Global test: χ2 = 11.54, df = 29, P = 0.98) with no evidence of overdispersion. AICc was therefore used for model selection. Further summaries of the mark-recapture data and detailed results of the U-CARE subcomponent tests are given in Appendix I.
Multi-state model development
Accounting for dominant sources of variation
Initial model building efforts focused on defining a base biology model and a sampling variation model. A model using capture effort (recapture rate and proportion of time with an ear-tag radio-transmitter) was substantially more supported than a model without these sampling covariates and other parameters held constant (ΔAICc = 64.1).
Of the base survival and biology models considered, a model with sex and roads as a survival covariate was the most supported (Table 3; Model 16). In addition, the model with roads as a covariate for survival was more supported than elevation (Model 24). For BCI change, we considered the models that had unique transition probabilities between all states (Model 23), equal transition probabilities between all states (Model 18), the models with unique transition probabilities for low to high and high to low change in body condition (Model 22) and models that had a different probability for increasing condition and decreasing condition (Model 16; denoted as I/D). Of these, a model with different probabilities of increases and decreases in condition was most supported. A model with change in Julian day of capture as a covariate was less supported suggesting that the relative timing of ‘spring’ capture across different years did not have a large impact on changes in BCI (Model 19).
Assessing the effect of covariates on body-condition change
We introduced habitat covariates to assess factors influencing increases or decreases in body condition. We started with factors that influence increases in condition and among these factors the variation in regeneration habitat was most supported. We found that the probability of an increase in body condition was similar between one state (i.e. low to medium) and two states (low to high) models and, thus, there was higher support of Model 1 where the two probabilities were equal. We then tested factors influencing decreases in body condition and found weaker support for all models. We also considered a model with the most supported covariates for BCI increase (Model 1: regenvar) and decrease (Model 5: presence of cubs) and the resulting model (Model 2) had weak support (ΔAICc = 0.47, wi = 0.31).
Predictions plotted for the most supported model of changes in body condition (see Table 3; Model 1) illustrate that bears with higher proportions of habitat with variation in regeneration habitat have a higher probability of increasing their body-condition state (Fig. 2). A female with cubs had a 2.26 times (SE = 1.01, CI = 2.2-2.31) greater chance of reducing her body-condition state (Model 2, see Table 3).
Relating body condition to survival
Given our defined baseline model for changes in body condition (see Table 3, Model 1 and Table 4, Model 3), we assessed if state of body condition influenced survival when contrasted against other environmental and biological covariates (see Table 4). Of the models considered, a model with unique survival rates for the low body-condition state (Model 1) without sex as a covariate, and a model with unique survival rates for the low body-condition state with sex included as a covariate (Model 2) were more supported than the baseline body-condition model (Model 3). A model with different survival rates for the high body-condition state (Model 4) was marginally supported as indicated by ΔAICc values of < 2. Models with state-specific survival rates were more supported than models that assumed age as a covariate (Model 8).
The effect of roads on survival was substantial with survival rates reduced at higher road densities for bears of all condition classes (Fig. 3). The precision in survival rates was, however, lower for bears associated with higher road densities since most bears were associated with low to moderate road densities. The majority of bears captured that were in the low condition state had 0 or lower road density with the distribution of road densities shifting toward higher road density with increasing condition state (Fig. 4). Using the mean values for road density from Figure 4, model-averaged estimates of survival suggest that bears with low body condition had higher survival rates (0.98, SE = 0.02, CI = 0.89-1.00) compared to bears in medium (0.94, SE = 0.031, CI = 0.84-0.98) and high (0.93, SE = 0.033, CI = 0.82-0.98) body condition. In this case, the road density of bears in each condition state, and the estimated difference in survival rate for the low strata (see Table 4, Model 1) affected the mean survival rate for bears as a function of condition state. The reporting rate of mortality of a bear with a GPS collar was 0.47 (SE = 0.10, SE = 0.29-0.67), and with a VHF ear tag it was 0.24 (SE = 0.24, CI = 0.02-0.8).
Synthesis of factors influencing changes in body condition and survival
The relationship between elevation, road density and variation in regeneration habitat illustrates demographic and geographic patterns in BCI change and survival (Fig. 5). Bears with lower levels of regeneration habitat have lower probabilities of increasing their body condition but also have higher probabilities of survival. Alternatively, bears with higher regeneration variation have a higher probability of increasing condition, but also have a higher rate of mortality (see Fig. 5). The net effect is the creation of an opposing elevational gradient in survival and bear condition given that most regeneration habitat and roads occurred at lower elevations (see Fig. 1). The differing survival rates and lower habitat value of low road density and low regeneration areas result in low condition bears occurring in areas of lower road density (see Fig 4).
Discussion
The results of our analysis illustrate the utility of using a multi-state model to explore and estimate trade-offs in health, habitat use and survival. Bears living in mountainous areas had higher survival rates but also had less probability of increasing their body condition. Bears in the foothills, on the other hand, have access to more diversity in forest seral stages related to forest management practices. The foothill bears therefore are more likely to increase in body condition, but also face higher rates of mortality. This has resulted in bears with lower body condition, but higher survival rates, a result that is counter-intuitive in terms of natural selection pressures on populations.
Nielsen et al. (2006) conducted a habitat-based analysis of mortality and habitat selection of grizzly bears which also suggested that higher value habitat was found in foothill areas which also had higher probabilities of mortality. Our analysis extends the work by Nielsen et al. (2006) by estimating the individual bear-level survival rates as well as probability of changing condition across the habitat/risk gradient in our study area. Further work will focus on spatially mapping areas where bears have the highest probability of gaining body condition in contrast to the highest risk of mortality. Unlike previous analyses, this method provides a survival estimate associated with various habitat types, and therefore, it provides a more direct estimate of demographic consequences of different habitat types. This type of tool would be useful for assessing the relative population impacts of management scenarios and delineating source and sink habitats. For example, our results demonstrate that anthropogenic activity potentially increases the overall habitat quality of forested areas (through creation of regeneration habitat). However, this is offset by increased road densities, thus creating sink habitats. Strategies that prohibit access to roaded areas would potentially maximize habitat value, while minimizing mortality risk, thus creating source habitats (see also Nielsen et al. 2006, 2008, Northrup et al. 2012).
Limitation of our analysis and alternative approaches
The main constraint of our analysis is the small sample size of bears and subsequent constraints on model complexity. The 111 bears that we monitored over 10 years were used in this analysis, which resulted in an effective sample size of 162. The maximum number of parameters is roughly the sample size of bears divided by 10 which would put the number of estimable parameters in the model at approximately 10-16. This limited our ability to model more complex state-specific relationships. In addition, sparse data limited the ability to test the goodness-of-fit of the multi-state models. It was possible that undetected heterogeneity and non-independence of detections occurred within the data set which could potentially have created overdispersion of multinomial variances. In this case, model selection may be optimistic, however, we note that the key covariates such as road density and regeneration variation showed strong support in terms of comparisons of models with and without these covariates, and therefore, the stronger signal of these relationships would still be meaningful if moderate overdispersion was present (> 1). Another potential issue in the analysis was that many models had the same number of parameters which could cause misleading inference on covariates using information-theoretic methods. An alternative approach would be to use analysis of deviance and permutation hypothesis tests to test for significant covariate associations (Lebreton et al. 2012). We feel that this approach was less ideal for our analysis which used a sequential-model building approach that was more amenable to information-theoretic methods, given inherent dangers of mixing information-theoretic and hypothesis-testing paradigms (Lukacs et al. 2007).
Despite sample size limitations, we suggest that our approach still allowed useful inference on the complex relationships between bears' demography and condition while accounting for sampling biases such as dropout and variable capture effort across time. We addressed sample size limitations by focusing our models on defined biological relationships (see Table 1) and through the use of an analysis of covariance approach that first defined a base biological model that accounted for major sources of variation in parameters followed by assessment of relationships of most interest. For example, previous research has demonstrated that roads are a dominant factor affecting survival rate and therefore this covariate was introduced early in the analysis to help control for heterogeneity in survival rates caused by the proximity of individual bears to roads. In addition, the capture effort covariate efficiently modeled temporal variation in recapture rates with a single parameter. This approach also ensured that the models considered were defined by biological and sampling-based hypotheses in the spirit of information-theoretic model selection methods (Burnham & Anderson 1992).
Survival rates from our analysis were not efficient since the main source of information about bears is either from live-capture or from documentation of mortality. It is not possible to include observations from the GPS collars which would require a Barker joint live-dead-resight (Barker & White 2001) model (that allows incorporation of observations). However, the Barker models do not incorporate condition states. Another approach that could be used is a multi-state model with an unobserved state for bears that were captured but had no body-condition measurement taken or were observed via radio-telemetry (Conn & Cooch 2008). While this approach could potentially offset issues with missing data, sample-size constraints (the number of individual bears and yearly captures) with our data did not allow the more complex model formulation needed to model uncertainty in transitions between observed and unobserved states. We did perform a cursory analysis to assess if the detected relationships between survival and body condition were similar if unobserved state observations were included by using an unobserved state and rerunning the survival models (in Table 4) with additional information from the recaptures of bears. The results suggested the same general ordering of models with bears having a low body-condition state also having the highest survival rates. In summary, the actual mark-recapture model that is optimal for survival analysis depends on the objectives of the analysis. In our case, we were interested in relationships between health states and survival, and therefore, the use of a less efficient survival model was justified.
Implications of our analysis
The results of our analysis suggested that mortality of bears and subsequent sample dropout could influence analysis of bear body condition. For example, our results suggested that bears in low body condition are more likely to remain in the sample compared to bears with medium and high body condition given that lower-condition bears live in lower risk habitats (see Fig. 4). Therefore, the available sample of bears could potentially favour lower condition bears over time which would suggest decreasing condition of bears. To use this approach as part of a monitoring programme it is important to distribute capture effort over a broad landscape with a diversity of habitat and anthropogenic states. These methods allow for an assessment of whether observed trends in body condition are an artifact of differential survival, or an actual biologically-based trend in the population. In addition, heterogeneity in survival rates caused by differential body condition could also contribute to bias in estimates of survival rate if low body-condition bears are more likely to remain in the collared bear samples than other bears (Zens & Peart 2003). The ability to estimate and compare state-specific survival rates allows a test for this form of heterogeneity.
For population management, the removal of higher body-condition bears could potentially result in decreases in reproductive rates of the population given that reproduction is often influenced by the body condition of bears (Elowe & Dodge 1989, Robbins et al. 2012). In terms of population viability, a covariance would be created where directional variation in survival also affects reproductive rates. The potential covariance between survival and reproductive rates could be explored using multi-state models where a female bear's state is determined by its reproductive class. Schwartz & White (2008) demonstrated the multi-state approach as a means to provide unbiased estimates of reproductive rate in bears. Their approach could be expanded by using body condition as a covariate for transitions between reproductive states, as well as by estimating survival rates for each reproductive state.
One further application of our approach is the extension of inference from DNA mark-recapture projects that collect hair from bears to estimate population size (Woods et al. 1999, Boulanger et al. 2002, 2006). Hair samples potentially also provide inference into bear stress levels through cortisol (Macbeth et al. 2010) as well as relative trophic state through stable isotopes (Ben-David et al. 2004). It should be possible to use stress or isotope scores from repeated samples of bears to assess trends in bear stress or other attributes without live capture. The use of the multi-state model, that also estimates recapture rates, would be perfectly suited for this application. In summary, the multi-state approach represents a robust analytical method for dealing with complex sampling biases when assessing relationships within dynamic environments using a variety of data types.
Acknowledgements
Funding for this work was received from the Foothills Research Institute, National Science and Engineering Research Council (NSERC), Canadian Forestry Company (CANFOR), Conoco Phillips, Shell Canada Ltd., Suncor Energy Inc., Sundance Forest Industry Ltd., Talisman Energy Inc., TransCanada Pipelines Ltd., Weyerhaeuser Company Ltd. and other research partners of the Foothills Research Institute Grizzly Bear Project.
References
Figure 1.
The Foothills grizzly bear study area in Alberta, Canada illustrating elevation (in m), terrain, roads and regeneration habitat. In general, regeneration habitat was created by forestry or forest fires (in non-roaded areas).
Figure 2.
Changes in grizzly bear body condition as a function of regeneration habitat variation. Data points are shown as squares in each figure. Confidence limits are given as grey lines.
Figure 3.
Annual survival rate as a function of road density for medium and high body condition bears (blue line) and low body condition bears (red line). Confidence intervals are shown on predictions (as dashed lines). Estimates are from Model 1 in Table 3.
Figure 4.
Violin plots displaying the distribution of road density as a function of condition state. A boxplot that delineates 25th and 75th percentiles is shown in black within each plot. The diameter of the plot indicates the kernel density of road densities for each state (Hintze & Nelson 1998). The white dot indicates the median for each state. Plots were based on samples sizes of 54, 58 and 57 capture events for low, medium and high condition state bears, respectively.
Figure 5.
Plots with bubble size showing the probability of a bear increasing its BCI value (A) and the probability of survival (B) as a function of elevation, regeneration habitat (increasing BCI) and road density (survival). Each bubble represents an individual bear in the analysis. Bears that were known mortalities are denoted as black bubbles.
Table 1.
Mean daily movement rates (km/day) used for defining landscape scales used to assess habitat based covariate in the Foothills grizzly bear study, Alberta, Canada, during 1999-2010. The number of individual bears used for estimates is given in parentheses. Pooled age class movement rates for Fall (late hyperphagia) were 5.68 for males (N = 72) and 7.61 for females (N = 48).
Table 2.
Predicted relationships between survival, body-condition index (BCI) and biological and environmental covariates for grizzly bears in Alberta, Canada. Positive (‘+’), negative (‘-’), and no relationship (‘0’) were hypothesized for each of the covariates.
Table 3.
Model selection results for hypothesized factors affecting change in body-condition index (BCI) for grizzly bears in Alberta, Canada. Covariates with a (+) modeled increases in BCI, whereas covariates with a (-) modeled decreases in BCI. Sample-size adjusted Akaike Information Criteria (AICc), difference in AICc between most supported and given model (ΔAICc), Akaike weight (wi), the number of parameters (K) and Deviance are shown. A capture effort model was used for recapture rate and a model with specific recovery rates for GPS and ear-tagged bears was used for all models. See Table 2 for information on individual covariates. I/D denotes a model that estimate different transition probability of increase and decrease in BCI. A 12 subcript denotes that probabilities of increase across 1 and 2 states were constrained to be equal. A 1&2 subscript denotes that different probabilities were assumed for an increase of 1 and 2 states.
Table 4.
Multi-state model selection results for analysis of factors influencing survival rates for grizzly bears in Alberta, Canada. The most supported model from Table 2 (Model 3) was used as a baseline model for this analysis. Sample-size adjusted Akaike Information Criteria (AICc), difference in AICc between most supported and given model (ΔAICc), Akaike weight (wi), the number of parameters (K) and Deviance are shown. A capture effort model was used for recapture rate and a model with specific recovery rates for GPS and ear-tagged bears was used for all models.
Appendices
Appendix I. Statistical details for the multi-state analysis
Transition events for the multi-state analysis pooled across all individual bears are summarized in Table I-1. The sequence of transitions is considered by previous and current strata. For example, a bear detection history of LL0M0000000 (with L and M indicating a BCI at capture of low and medium states and a 0 indicating non-capture) would contribute 1 event to the previous low-current low cell and 1 event to the previous low-current medium cell in the table. In general, the number of transitions was relatively low due to bears being detected once. We note that frequencies of bears captured once do not necessarily indicate transients given that capture effort was inconsistent across the study area (as accounted for by the capture effort covariate). Bears that were captured once potentially still contributed to the modeling of survival rate (given that some were eventually redetected as mortalities).
Covariates used in this analysis were correlated due to the elevational and anthropogenic gradient within the study area (see Fig. 1). Table I-2 displays spearman rank correlations for covariates.
An assumption of multi-state models is that detection of individuals are independent of each other and independent across time. Program U-CARE provides detailed tests into this assumption (Choquet et al. 2009, Pradel et al. 2005). The global goodness-of-fit test was the summation of a test of whether past capture of individuals affected future captures (Test 3G) and a test of whether individuals of behavioural response to trapping (Test M). Test 3G was further broken down into tests for transient individuals and non-independence or memory of individuals in terms of previous and future states. All tests were non-significant (Table I-3) with the resulting estimate of c-hat (total/df) being < 1. Low sample sizes (see Table I-1) most likely affected the power of tests to detect violation of assumptions and therefore these results should be interpreted cautiously.
Estimates from the most-supported mark-recapture model (see Table 4, Model 1) are shown in Table I-4. Convergence was also tested using simulated annealing. Standard error estimates reveal that some of the beta parameters are imprecise due to the relatively sparse data used in the analysis.
Table I-1.
Summary of transitions between states and related capture events for the Foothills grizzly bear study, Alberta, Canada, during 1999-2010. Non-recaptures between transitions were not considered in the summary. Other captures are based upon the BCI state of the bear when capture occurred. Mortalities are listed by the last state from the last live capture event before the bear was reported as a mortality.
Table I-2.
Spearman rank correlation for covariates used in the analysis based upon 111 comparisons. The column heading correspond to rows but are abbreviated. Significant correlations (at α < 0.05) are written in italics.
Table I-3.
Summary of U-CARE test results for the Foothills Model Forest grizzly bear study during 1999-2010. Tests 3G and M added up to the global goodness-of-fit test as indicated by double lines. Test statistics (χ2), corresponding degrees of freedom (df) and P-values are given. See Choquet et al. (2009) for more details on component tests.
Table I-4:
Parameter estimates for the most supported multi-state model (see Table 4, Model 1). All parameter estimates are on the logit scale.
