We developed a hierarchical model in a Bayesian context to jointly model both the dynamics of beetle activity intensity over time within our plots, as well as the occurrence–accounting for imperfect detection–of black-backed woodpeckers at those same plots. The model largely follows the structure of a single-species occupancy model [17 (link)], where woodpecker observations of detection or non-detection, yjkt, for survey interval k at site j (where sites are individual survey points) in year t, are assumed to be imperfectly observed representations of the true occurrence status, zjt (present or absent), which is constant across all k survey intervals (i.e., closure is assumed within the <17-minute survey period) but can change from year to year. Observed occurrence of black-backed woodpeckers, yjkt, is thus modeled as
where pjkt is the probability of detection for a given survey at a site. Similarly, the true occurrence status of a site in year t, zjt, is modeled as
where ψjt is the probability of occurrence at a site. In this context, occupancy is defined by the space and time in which the survey is conducted and across which closure is assumed [18 ]. Thus, our calculated probabilities represent the probability of at least one black-backed woodpecker occurring within (or alternatively, ’using’ [19 ]) the detection radius of a survey point (approximately 120 m) during a survey period (median = 7 minutes in duration).
The probabilities of woodpecker detection and occurrence are both modeled as logit-linear functions of covariates chosen a priori. Following previous work studying black-backed woodpeckers with this survey methodology [14 (link), 15 (link), 20 (link)], we expected detection, pjkt, to vary as a function of an intercept and the linear additive combination of a categorical covariate representing the survey type (passive = 0, broadcast = 1), giving
The probability of woodpecker occupancy of a survey point was modeled as a function of five covariates: (1) elevation, (2) latitude, (3) snag density, (4) intensity of beetle larvae activity (as indirectly measured by cumulative beetle sign since the fire; modeled as a latent variable, see below), and (5) an interaction between years-since-fire and the intensity of beetle larvae activity (with the hypothesis that cumulative beetle sign becomes less predictive over time). Snag counts were conducted immediately after completing woodpecker surveys and consisted of counting all snags of different size classes (10–30, 30–60, and >60 cm dbh) within 50 m of each survey point. Size-specific snag counts were aggregated in the field into different categories (≤5, 6–15, 16–30, 31–50, 51–100, >100), which were converted to numerical quantities (1, 6, 16, 31, 51, 101, respectively) for analysis [15 (link)]. Counts across all three size classes were summed to calculate a relative index of snag density (snags/ha). The linear additive model for occupancy in the first year of surveys can be described as,
where β represents intercept and slope parameters. To account for pseudoreplication and temporal autocorrelation derived by sampling sites repeatedly in consecutive years, we added a temporal autocorrelation term [21 (link)], ϕ, which was multiplied by the true occurrence status in year t−1, resulting in the following model for additional post-fire years,
As 2018 was the last year of surveys used in this dataset, and also the only year with in situ beetle activity surveys, all surveys conducted in 2018 held the temporal index of t = 10. Surveys in previous years (t = 1,…,9) were treated as missing data if no surveys occurred at a site in that survey-year. Finally, the intercept (β0,j) was modeled as a random effect for each fire (n = 22), drawn from a hierarchical normal informed by a common mean (μβ0) and precision (τβ0).
Previous analyses of black-backed woodpecker occurrence have demonstrated the importance of the number of years since fire in models of the species’ occurrence [14 (link), 17 (link)], yet our model of woodpecker occupancy does not include an independent effect of time since fire (Eqs4 –5 ). By excluding time since fire from our model of woodpecker occupancy, we assume that all temporal changes in occupancy are due to either (1) temporal changes in habitat quality as indicated by intensity of beetle sign (parameters βbeetle and βageXbeetle), or (2) random stochasticity (the frequency of which is governed by the parameter ϕ).
A novel feature of our multi-trophic model is that we treat cumulative beetle larvae sign at a survey point each year (intensityjt) as a latent (i.e., indirectly observed) variable, which for mathematical simplicity we define as continuous. We are then able to model beetle larvae sign as a function of different environmental variables expected to relate to beetle activity and to account for the known dynamic that beetle sign generally accumulates over time even though overall activity may decline. Thus, we hypothesized that the intensity of beetle sign at a site each year (intensityjt) varies as a function of: (1) the number of years since fire; (2) the proportion of sampled trees per point that were of the genus Pinus; (3) and an interaction between the proportion of pines and years since fire. Based on previous work [13 (link)], we hypothesized that beetle sign increases over time (as sign is generally cumulative and lasting), but that pines would have greater activity in early post-fire years and lower activity in later post-fire years (as pine bark generally decomposes faster than bark of other trees in our study areas). We thus modeled beetle sign intensity as,
We fit this model to observed data collected in 2018, by treating the total sum of beetle sign scores across all surveyed trees per point (max = 6) as a binomially distributed variable, as follows,
where the maximum activity score is a product of the number of trees sampled per point (numTreesj) and the maximum potential activity score per tree (i.e., 8).
We fit the model to the data with JAGS [22 ] using the R statistical programming language version 4.0.2 [23 ] and the package ‘R2jags’ [24 ]. We used vague priors (i.e., normal with μ = 0, τ = 0.1). We ran three chains of 50,000 iterations thinned by 50 with a burn-in of 50,000, yielding a posterior sample of 3,000 across all chains. Convergence was checked visually with traceplots and confirmed with a Gelman-Rubin statistic < 1.1 [25 ]. Inference on parameters was made using 95% Bayesian credible intervals (95 CI).
where pjkt is the probability of detection for a given survey at a site. Similarly, the true occurrence status of a site in year t, zjt, is modeled as
where ψjt is the probability of occurrence at a site. In this context, occupancy is defined by the space and time in which the survey is conducted and across which closure is assumed [18 ]. Thus, our calculated probabilities represent the probability of at least one black-backed woodpecker occurring within (or alternatively, ’using’ [19 ]) the detection radius of a survey point (approximately 120 m) during a survey period (median = 7 minutes in duration).
The probabilities of woodpecker detection and occurrence are both modeled as logit-linear functions of covariates chosen a priori. Following previous work studying black-backed woodpeckers with this survey methodology [14 (link), 15 (link), 20 (link)], we expected detection, pjkt, to vary as a function of an intercept and the linear additive combination of a categorical covariate representing the survey type (passive = 0, broadcast = 1), giving
The probability of woodpecker occupancy of a survey point was modeled as a function of five covariates: (1) elevation, (2) latitude, (3) snag density, (4) intensity of beetle larvae activity (as indirectly measured by cumulative beetle sign since the fire; modeled as a latent variable, see below), and (5) an interaction between years-since-fire and the intensity of beetle larvae activity (with the hypothesis that cumulative beetle sign becomes less predictive over time). Snag counts were conducted immediately after completing woodpecker surveys and consisted of counting all snags of different size classes (10–30, 30–60, and >60 cm dbh) within 50 m of each survey point. Size-specific snag counts were aggregated in the field into different categories (≤5, 6–15, 16–30, 31–50, 51–100, >100), which were converted to numerical quantities (1, 6, 16, 31, 51, 101, respectively) for analysis [15 (link)]. Counts across all three size classes were summed to calculate a relative index of snag density (snags/ha). The linear additive model for occupancy in the first year of surveys can be described as,
where β represents intercept and slope parameters. To account for pseudoreplication and temporal autocorrelation derived by sampling sites repeatedly in consecutive years, we added a temporal autocorrelation term [21 (link)], ϕ, which was multiplied by the true occurrence status in year t−1, resulting in the following model for additional post-fire years,
As 2018 was the last year of surveys used in this dataset, and also the only year with in situ beetle activity surveys, all surveys conducted in 2018 held the temporal index of t = 10. Surveys in previous years (t = 1,…,9) were treated as missing data if no surveys occurred at a site in that survey-year. Finally, the intercept (β0,j) was modeled as a random effect for each fire (n = 22), drawn from a hierarchical normal informed by a common mean (μβ0) and precision (τβ0).
Previous analyses of black-backed woodpecker occurrence have demonstrated the importance of the number of years since fire in models of the species’ occurrence [14 (link), 17 (link)], yet our model of woodpecker occupancy does not include an independent effect of time since fire (Eqs
A novel feature of our multi-trophic model is that we treat cumulative beetle larvae sign at a survey point each year (intensityjt) as a latent (i.e., indirectly observed) variable, which for mathematical simplicity we define as continuous. We are then able to model beetle larvae sign as a function of different environmental variables expected to relate to beetle activity and to account for the known dynamic that beetle sign generally accumulates over time even though overall activity may decline. Thus, we hypothesized that the intensity of beetle sign at a site each year (intensityjt) varies as a function of: (1) the number of years since fire; (2) the proportion of sampled trees per point that were of the genus Pinus; (3) and an interaction between the proportion of pines and years since fire. Based on previous work [13 (link)], we hypothesized that beetle sign increases over time (as sign is generally cumulative and lasting), but that pines would have greater activity in early post-fire years and lower activity in later post-fire years (as pine bark generally decomposes faster than bark of other trees in our study areas). We thus modeled beetle sign intensity as,
We fit this model to observed data collected in 2018, by treating the total sum of beetle sign scores across all surveyed trees per point (max = 6) as a binomially distributed variable, as follows,
where the maximum activity score is a product of the number of trees sampled per point (numTreesj) and the maximum potential activity score per tree (i.e., 8).
We fit the model to the data with JAGS [22 ] using the R statistical programming language version 4.0.2 [23 ] and the package ‘R2jags’ [24 ]. We used vague priors (i.e., normal with μ = 0, τ = 0.1). We ran three chains of 50,000 iterations thinned by 50 with a burn-in of 50,000, yielding a posterior sample of 3,000 across all chains. Convergence was checked visually with traceplots and confirmed with a Gelman-Rubin statistic < 1.1 [25 ]. Inference on parameters was made using 95% Bayesian credible intervals (95 CI).
Free full text: Click here
