Brittany Mosher Changes in Site Occupancy of Yellow-Rumped Warblers in Response to a Mountain Pine Beetle Outbreak Introduction Background and Justification The pine dominated forests of the Rocky Mountains are of interest to residents of the West for a number of reasons ranging from recreation to aesthetics. Recent outbreaks of mountain pine beetles (Dendroctonus ponderosae) will fundamentally alter this region, impacting management decisions related to fire, logging, and habitat management as well as management and conservation for both game and non-game species (Kaufmann et al. 2008). Appropriate analyses of the impacts of this mountain pine beetle outbreak will be essential in making good management decisions in the face of future outbreaks. The mountain pine beetle is a native and essential part of the natural disturbance cycle of forests in Western North America (Sartwell and Stevens 1975). Fire suppression, climate, and lack of thinning have altered the typical interactions between insects and forests (Caroll 2004, Parker et al. 2006) resulting in larger outbreaks of epidemic proportions. Mountain pine beetles carrying symbiotic blue-stain fungi (Ophiostoma clavigerum and Ophiostoma montium) attack large mature pines, boring holes in trees in order to lay their eggs (Klepzig et al. 2001, Waring 1985). In the process, beetles inoculate the fungi into the sapwood of the host tree, slowing water transport and thereby hastening desiccation and eventual death of infested trees (Waring 1985, Solheim and Krokene 1998). The mountain pine beetle outbreaks have substantial ecological and economic ramifications, and have provoked calls for policy changes of aggressive timber harvesting and fire suppression (Raffa et al. 2008). Mountain pine beetle outbreaks can increase fine woody fuels (< 7.63cm) and litter of an area in the years following a beetle epidemic (Jenkins et al. 2008). Increased fuels and concern for wildfires on public lands has provided a valid reason for many agencies to consider salvage logging as part of their post-beetle management plan. In addition to fuel concerns, timber harvest can be negatively impacted by the mountain pine beetle. There is currently not a market for timber that has been infested by the mountain pine beetle. In British Columbia an estimated 10 million hectares of pine forest and 435 million cubic meters of timber have been lost as a result of mountain pine beetle damage since the outbreak began in 2000 (Walton et al. 2007). In the interest of maintaining revenue from timber harvest in the face of an oncoming beetle epidemic, it may be prudent to log trees prior to the beetle attack. Brittany Mosher The effects of mountain pine beetle outbreaks are beneficial to some native wildlife, including excavators and cavity-nesting birds. Direct benefits are present for species including Hairy Woodpeckers (Picoides villosus) and American Three-toed Woodpeckers (Picoides tridactylus) that consume mountain pine beetle larvae (Jackson et al. 2002 and Leonard 2001). Indirect benefits in the form of more complex forest structure, increased snag density, and increased downed wood and stumps provide habitat and nesting benefits for many other species. Other species that nest and forage in the tree canopy, such as the Yellow-rumped Warbler (Dendroica coronata,) may decline during and after the beetle outbreak. As the pines lose their needles, species like the Yellow-rumped Warbler experience decreased nesting and foraging habitat. The habitat changes resulting from the mountain pine beetle epidemic will also have implications for other important game and non-game wildlife species as well. Occupancy Modeling The objective of this study is to investigate the site occupancy dynamics of the Yellow-rumped Warbler (Picoides villosus) in the Helena National Forest in 2003 (pre-outbreak) and in 2010 (during outbreak). It is hypothesized that occupancy probability for the Yellow-rumped Warbler will be lower in 2010 (during epidemic) than in 2003 (pre-epidemic) since the amount of available nesting and foraging habitat will have decreased. Occupancy probability (ψ) is the probability that a site chosen at random is occupied by the species of interest. The observed occupancy of a site is a binomial response where 1’s indicate presence of the species and 0’s indicate that a) the site is unoccupied or b) the site is occupied but the species went undetected. Since observers detect species imperfectly, we must also estimate detection probability (p). Detection probability is the probability that the species of interest will be detected at a site, given that the site is indeed occupied. Detection probability can be estimated when multiple visits occur at sites within a closed “season”. Covariates for modeling detection can be site specific or survey specific. In a single year study, covariates incorporated to model occupancy are those that may vary by site. Table 1 illustrates an example single-season dataset, where the researchers are studying 4 sites over a season. During the season, they visit each site and look, listen, or search for the species of interest 3 times. In the first 3 sites, the species is detected at least once over the 3 visit season. Since we assume that sites are closed to changes in population, we can infer that the 0’s that occur in sites A, B, and C were instances where we failed to detect the species of interest and are not true absences. In site D, one of two scenarios is occurring: 1) the site is unoccupied or 2) the site is occupied but the observers failed to detect the species of interest in each of 3 surveys. Brittany Mosher Table 1: Example data set for a single season occupancy model. SITE HISTORY A 100 B 101 C 010 D 000 In the maximum likelihood method of calculating occupancy and detection, the first step is to calculate the probability of obtaining each “type” of history, given the data. In this case, the probability of obtaining the history seen for Site C is given by the equation: Pr(hC) = ψ(1-p1) * ψp2*ψ(1-p3). Since we assume the species of interest was at the site throughout the survey season we know that on visit one, the site was occupied but the species was not detected (ψ(1-p1)), on visit two the site was occupied and the species was detected (ψp2) and on visit three the site was occupied and the species was not detected (ψ(1-p3)). The product of these three expressions is the probability of obtaining a history like the one obtained for site C. Probability statements for the histories for site A, B, and C can all be constructed in this way, since we know that the species was present. Site D has a history of 000. In this case we are not sure whether the site was truly unoccupied or whether the site was occupied but the observers failed to detect the species. Both of these possibilities must be accounted for. The equation that represents the probability of obtaining a history like the one for site D is: Pr(hD) = ψ(1-p1)* ψ(1p2)*ψ(1-p3) + (1-ψ). On the left side of the addition sign, we consider the possibility that the site was occupied but the species went undetected on each of the 3 visits. On the right side of the addition sign we illustrate the other possibility; that the site is not occupied in the first place. In a similar way as above, a likelihood can be built that takes into account all observed histories (the data) simultaneously (Equation 1). Equation 1: Model likelihood for a single season occupancy model. N L ( , p | hi ,... h N ) Pr( h i ) [ i 1 K iD k 1 p k (1 p k ) ik i D ik K ][ (1 p k ) (1 )] k 1 ii D Brittany Mosher In this equation the number of sites is i=1,…,N and the number of surveys is k=1,…,K. The number of histories with at least one detection is iD and ik is the number of sites where the species was detected on the kth survey. This maximum likelihood (ML) model is an appropriate means of analyzing occupancy data in many cases, however sometimes it is more advisable to switch to a Bayesian framework. There are 3 main reasons for using a Bayesian method rather than ML to model occupancy. Firstly, the empirical Bayesian methodology is superior to the maximum likelihood method in that it allows us to incorporate informative prior knowledge about the species. Also, the ML model sometimes has difficulty reaching convergence when dealing with rare species. Rare species or species that are not easily detected yield histories with a great deal of 0’s. In these scenarios, the ML method often fails to converge whereas the empirical Bayes method is able to provide parameter estimates. Finally, the Bayesian approach can be used effectively in some more complex occupancy problems. Multi-year and multi-species models have recently been developed using the Bayesian framework so that information can be shared across species, sites, and years (Russell et al. 2009, Ruiz-Gutierrez et al. 2010). Methods 1. Field Methods Study Area The study area is located within the northern Elkhorn Mountains of the Helena National Forest and consists of four study units ranging in size from 130 to 260 ha. The original study was part of an Interior West-wide project (Birds and Burns Network1) designed to evaluate the effectiveness of prescribed fire treatments in reducing fuels and how potential fuel reductions influence habitats and populations of the associated avifuana (e.g., Saab et al. 2006, Fowler et al. 2008, Russell et al. 2009). The data used in this analysis were collected in 2003 and again in 2010. The study units are characteristic of dry mixed coniferous forests, largely dominated by ponderosa pine (Pinus ponderosa) with lesser amounts of Douglas-fir (Pseudotsuga menziesii) and lodgepole pine (Pinus contorta). The study units also contain stands of quaking aspen (Populus tremuloides), found mainly along snow melt drainages and creeks. Point Counts In 2002, we established approximately 20 variable radius point count stations in each of the four study units (76 total). The point count stations were established at locations randomly while also ensuring that points were located 250 m apart and 100 m from the study unit edge. Each station was 1 http://www.rmrs.nau.edu/wildlife/birdsnburns Brittany Mosher visited three times approximately every two weeks between 22 May and 3 July during 2003 and 2010. To ensure a consistent level of bird activity, we began point counts just after the dawn chorus and completed them within five hours. At each point, observers recorded all birds detected during a fiveminute count period, and estimated the distance to each individual observed in distance classes (Russell et al. 2009). Since estimating distances accurately can be difficult (Johnson 2008) we truncated the data for our analysis at 75m to minimize incorrectly categorizing the further distance classes (Alldredge et al. 2007). Vegetation Sampling We measured forest components at 76 random point count stations (18-20 stations per unit) during the summers of 2003 and 2010. Each vegetation plot is 2.5 ha in size and consists of nested plots in a cross plot design of two intersecting 20m x 100m rectangles (Saab et al. 2006). The plot is centered on the point count station. Covariates measured include percent canopy cover using a spherical densitometer, live tree and snag density, and understory vegetation composition. W ild life S n ags (> 23 cm d iam eter at b reast h eigh t) 20 m T h ree 4 m rad iu s su bp lots for can op y, grou n d and sh ru b m easu rem en ts focal 6 m 50 m S ection 2 m Fu el sn ags an d trees (< 23 cm d iam eter at b reast h eigh t) and w il d life logs (> 23 cm large en d diam eter) Figure 1: Vegetation plot diagram. W ild life T rees (> 23 cm D iam eter at b reast h eigh t) Brittany Mosher 2. Model The model developed was a single-season single-species occupancy model based on detection/nondetection point count data. The model, originated by Kery and Royle (2008), assumes that true occupancy is a latent variable which we observe imperfectly. The latent state is denoted z(i) such that z(i) is 1 for sites occupied by the species of interest and 0 for unoccupied sites. The model assumes that if the species of interest was detected at a given site during one of the surveys, the site is occupied. Multiple visits allow us to differentiate non-detection and true absence. The model assumes that observers record no false detections, that populations at a site are closed over a season, and that detections at sites are independent of one another. Models were constructed for occupancy and detection probabilities and were used for the 2003 and 2010 data. The Yellow-rumped Warbler site-specific occurrence model was specified with occupancy (the latent state) as a Bernoulli random variable z(i) ~Bern(ψi), where ψi is the site-specific probability of occupancy. The model assumes that occupancy probability is linearly related to a series of readily available vegetation covariates on the logit scale, as shown below. The covariates were screened for second order correlations, and in situations where correlation exceeded |0.5|, one of the two covariates was eliminated. The covariates were centered before being incorporated into the model. The vegetation covariates are: live tree density(β1—trees that are >23 cm in diameter at breast height [dbh]), snag density (β1—snags that are >23 cm dbh), and live fuel density (β1--where a fuel is >1m tall but <23 cm dbh). logit(ψ[i]) <- β0 + β1*trees.cov[i] + β2*snags.cov[i] + β3*fuels.cov[i] The detection model specifies detection as a site and survey specific Bernoulli random variable x(i,k)~Bern(zi*pi,k), where pi,k is site and/or survey specific detection probability. The relationship shows that if the site is unoccupied (zi=0), there is no chance of detecting the species. The logit link was used to model a linear relationship between detection probability and several covariates. The covariate jdate is day of year, and is the only survey-specific covariate included. logit(p[i,k]) <- α0 + α1*jdate[i,k] + α2*trees.cov[i] + α3*snags.cov[i] + α4*fuels.cov[i] The models were analyzed using a Bayesian approach and vague priors in WinBUGS. The equations were implemented as a single season and single species model for 2003 and then again for 2010. Results were compared. WinBUGS implementation and code are in the appendix. Brittany Mosher Results The main question of interest in this analysis was whether or not the occupancy probability for the canopy specialist Yellow-rumped Warbler changed as a result of the mountain pine beetle outbreak. Table 2 shows the results of the occupancy estimates and 95% credible intervals for 2003 and for 2010. The occupancy estimate for average values of the covariates in 2010 is 0.53 (95%CI: 0.50,0.60) while the estimate for 2003 is 0.69 (95%CI: 0.62,0.73). These non-overlapping credible intervals indicate that there is a difference in occupancy between the two years. Table 3 shows the coefficient estimates for detection and occupancy on the logit scale, along with standard deviations, Monte Carlo error, and 95% credible intervals for the 2003 analysis (Table 4 shows the 2010 analysis). Figure 2 shows the 3 chains over 1500 iterations for 2003 (Figure 4 for the 2010 analysis). Figure 3 shows the 2D marginal plot for the estimate of occupancy at mean values of covariates for 2003 (Figure 5 for the 2010 analysis). Figure 6 shows an example of one of the vegetation covariate coefficients that shows poor convergence. The prior for this particular parameter was U(-10,10) and the marginal plot shows that not very much information was gained from the data to inform this vague prior. Such was the case for several of the vegetation covariates in both years of the study. Widening the prior to U(-20,20) and increasing the number of iterations improved the convergence for this particular parameter (Figure 7). Some more work on tweaking priors and altering chains, thin rate, iterations, and burn in needs to be done before any interpretation of the effect sizes for various covariates occurs. Table 2: Occupancy estimates at mean values of the covariates in 2003 versus 2010. The non-overlapping credible intervals indicate that there may truly be a difference in occupancy between these time periods. YEAR Occupancy Estimate Lower 95% CI Upper 95% CI 2003 0.69 0.62 0.73 2010 0.53 0.50 0.60 Table 3: Model results. a1-a4 denote effect sizes of covariates in the detection model, whereas b1-b3 represent effect sizes for covariates in the occupancy model for the 2003 (pre-outbreak) season. node a1 a2 a3 a4 b1 b2 b3 deviance p0 psi0 mean 0.01211 -0.7144 0.1212 -0.6876 0.6133 0.03636 0.2181 234.9 0.5168 0.6935 sd 0.01491 5.555 5.763 3.21 2.87 2.87 2.768 9.02 0.01532 0.03161 MC error 2.798E-4 0.1015 0.1269 0.07736 0.0548 0.05285 0.06584 0.193 5.008E-4 6.044E-4 2.5% -0.01674 -9.573 -9.512 -6.678 -4.673 -4.682 -4.681 218.2 0.5004 0.6176 median 0.0122 -1.005 0.08708 -0.7574 0.8257 0.02894 0.3447 234.7 0.5125 0.7024 97.5% 0.04189 9.221 9.57 5.93 4.848 4.739 4.701 253.1 0.5568 0.7298 start 501 501 501 501 501 501 501 501 501 501 sample 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 Brittany Mosher psi0 chains 1:3 0.8 0.7 0.6 0.5 501 750 1000 1250 1500 iteration Figure 2: Example diagnostic from WinBUGS, showing 3 well-mixed chains and convergence for 2003. psi0 chains 1:3 sample: 3000 30.0 20.0 10.0 0.0 0.5 0.6 0.7 Figure 3: Marginal plot of occupancy at mean values of the covariates for 2003. Table 4: Model results. a1-a4 denote effect sizes of covariates in the detection model, whereas b1-b3 represent effect sizes for covariates in the occupancy model for the 2010 (during outbreak) season. node a1 a2 a3 a4 b1 b2 b3 deviance p0 psi0 mean -0.05335 -1.166 0.4826 -2.081 -0.05679 -0.09435 1.364 134.1 0.5129 0.53 sd 0.02331 5.702 5.852 4.719 2.883 2.907 2.636 9.865 0.01259 0.02635 MC error 4.754E-4 0.1251 0.1081 0.1288 0.04757 0.05287 0.05176 0.1931 3.699E-4 5.198E-4 2.5% -0.1016 -9.632 -9.581 -9.506 -4.757 -4.767 -4.262 116.1 0.5003 0.5008 median -0.05299 -1.762 0.8614 -2.525 -0.1163 -0.1146 1.861 133.7 0.5088 0.5229 97.5% -0.009719 9.273 9.479 7.954 4.699 4.713 4.849 154.9 0.5469 0.6005 start 501 501 501 501 501 501 501 501 501 501 sample 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 psi0 chains 1:3 0.7 0.65 0.6 0.55 0.5 501 750 1000 1250 1500 iteration Figure 4: Example diagnostic from WinBUGS, showing 3 well-mixed chains and convergence for 2010. Brittany Mosher psi0 chains 1:3 sample: 3000 30.0 20.0 10.0 0.0 0.45 0.5 0.55 0.6 0.65 Figure 5: Marginal plot of occupancy at mean values of the covariates for 2010. a2 chains 1:3 sample: 3000 0.08 0.06 0.04 0.02 0.0 -20.0 -10.0 0.0 10.0 Figure 6: Marginal plot of effect size of tree density on detection with a prior of U(-10,10) and 1500 iterations. Poor convergence is shown. a2 chains 1:3 sample: 3501 0.04 0.03 0.02 0.01 0.0 -40.0 -20.0 0.0 20.0 Figure 7:Marginal plot of effect size of tree density on detection with a prior of U(-20,20) and 5000 iterations. Convergence here is improved from Figure 6. Discussion While the Yellow-rumped Warbler seems to respond negatively to the mountain pine beetle outbreak, other species may respond differently. Differences in nesting or foraging strategy may influence how certain avian species respond to the outbreak. The time scale investigated also may influence the results found. This particular study is on a short time scale. As the outbreak progresses and trees continue to deteriorate and fall, trends in occupancy may change as well. The next step for me is to use a multi-season (6 years), multi-species (69 species) occupancy model (Russell et al. 2009, Ruiz-Gutierrez et al. 2010) to get a clear picture of how various species and guilds in the avian community respond to the mountain pine beetle outbreak. Brittany Mosher Literature Cited Carroll, A. L., S. W. Taylor, J. Regniere, and L. Safranyik. 2004. Effects of climate change on range expansion by the mountain pine beetle in British Columbia. Pages 223–232 in T. L. Shore, J. E. Brooks, and J. E. Stone, editors. Mountain pine beetle symposium: challenges and solutions. Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre, Information Report BC-X-399, Victoria, British Columbia, Canada. Fowler, J. F, C. H Sieg, B. G Dickson, and V. Saab. 2008. Exotic plant species diversity: Influence of roads and prescribed fire in Arizona ponderosa pine forests. Rangeland Ecology & Management 61, no. 3: 284–293. Jackson, Jerome A., Henri R. Ouellet and Bette J. Jackson. 2002. Hairy Woodpecker (Picoides villosus), The Birds of North America Online (A. Poole, Ed.). Ithaca: Cornell Lab of Ornithology; Retrieved from the Birds of North America Online: http://bna.birds.cornell.edu/bna/species/702 Jenkins, M. J, E. Hebertson, W. Page, and C. A Jorgensen. 2008. Bark beetles, fuels, fires and implications for forest management in the Intermountain West. Forest Ecology and Management 254, no. 1: 16–34. Kaufmann M.R., G.H. Aplet, M. Babler, W.L. Baker, B. Bentz, M. Harrington, B.C. Hawkes, L. Stroh Huckaby, M.J. Jenkins, D.M. Kashian, R.E. Keane, D. Kulakowski, C. McHugh, J. Negron, J. Popp, W.H. Romme, T. Schoennagel, W. Shepperd, F.W. Smith, E. Kennedy Sutherland, D. Tinker, and T.T. Veblen. 2008. The status of our scientific understanding of lodgepole pine and mountain pine beetles – a focus on forest ecology and fire behavior. The Nature Conservancy, Arlington, VA. GFI technical report 2008-2.] Kéry, M. and J.A. Royle. 2008. Hierarchical Bayes estimation of species richness and occupancy in spatially replicated surveys. Journal of Applied Ecology, 45: 589–598. Klepzig, K. D., J. C. Moser, F. J. Lombardero, and M. P. Ayres. 2001. Symbiosis and competition: complex interactions among beetles, fungi and mites. Symbiosis 30: 83–96. Leonard, Jr., David L. 2001. American Three-toed Woodpecker (Picoides dorsalis), The Birds of North America Online (A. Poole, Ed.). Ithaca: Cornell Lab of Ornithology; Retrieved from the Birds of North America Online: http://bna.birds.cornell.edu/bna/species/588 Parker, T. J, K. M Clancy, and R. L Mathiasen. 2006. Interactions among fire, insects and pathogens in coniferous forests of the interior western United States and Canada. Agricultural and Forest Entomology 8, no. 3: 167–189. Raffa, K. F, B. H Aukema, B. J Bentz, A. L Carroll, J. A Hicke, M. G Turner, and W. H Romme. 2008. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: the dynamics of bark beetle eruptions. Bioscience 58, no. 6: 501–517. Ruiz-Gutiérrez, V., E. F. Zipkin, and A. A. Dhondt. 2010. Occupancy dynamics in a tropical bird community: unexpectedly high forest use by birds classified as non-forest species. Journal of Applied Ecology 47(3): 621-630. Brittany Mosher Russell, Robin E., J. Andrew Royle, Victoria A. Saab, John F. Lehmkuhl, William M. Block, and John R. Sauer. 2009. Modeling the effects of environmental disturbance on wildlife communities: avian responses to prescribed fire. Ecological Applications 19: 1253-1263. Saab, V., L. Bate, J. Lehmkuhl, B. Dickson, S. Story, S. Jentsch, and W. Block. 2006. Changes in downed wood and forest structure after prescribed fire in ponderosa pine forests. Pages 477-487 In Andrews, PL, Butler, BW (compilers), How to Measure Success: Conference Proceedings. USDA Rocky Mount. Res. Sta. Proceed. RMRS-P-41., Fort Collins, CO. Sartwell, C., and R. E Stevens. 1975. Mountain pine beetle in ponderosa pine–prospects for silvicultural control in second-growth stands. Journal of Forestry 73, no. 3: 136–140. Solheim, H., and P. Krokene. 1998. Growth and virulence of mountain pine beetle associated blue-stain fungi, Ophiostoma clavigerum and Ophiostoma montium. Botany 76, no. 4: 561–566. Walton, A. et al. Provincial-Level Projection of the Current Mountain Pine Beetle Outbreak http://www.for.gov.bc.ca/hre/bcmpb/Year4.htmæ (British Columbia, Ministry of Forest and Range, Victoria, 2007). Waring, R. H., and G. B. Pitman. 1985. Modifying lodgepole pine stands to change susceptibility to mountain pine beetle attack. Ecology 66, no. 3: 889–897. Appendix #Single Season Occupancy with all Covariates yrwa.histories10 <- yrwa.histories[,17:19] library(R2WinBUGS) # WinBUGS occupancy model ################################################################## #Create one survey covariate (jdate) and one static covariate (pcc) tmpobs <- as.matrix(hemt.obscovs[ ,6:9]) jdate <- NULL pcc <- NULL trees.cov snags.cov stems.cov fuels.cov <<<<- NULL NULL NULL NULL numsites <- length(yrwa.histories10[ ,1]) #identify covariates jdate <- as.matrix(hemt.obscovs[ ,6:8]) pcc <- as.vector(hemt.obscovs[ ,9]) trees.cov <- as.vector(treedens) snags.cov <- as.vector(snagdens) Brittany Mosher #stems.cov <- as.vector(stemdens) fuels.cov <- as.vector(fuelsdens) #center and round covariates if necessary pcc <- as.vector(pcc - mean(pcc)) pcc <- round(pcc, digits=3) is.vector(pcc) trees.cov <-as.vector(trees.cov-mean(trees.cov)) trees.cov <- round(trees.cov, digits=3) snags.cov <-as.vector(snags.cov-mean(snags.cov)) snags.cov <- round(snags.cov, digits=3) #stems.cov <-as.vector(stems.cov-mean(stems.cov)) #stems.cov <- round(stems.cov, digits=3) fuels.cov <-as.vector(fuels.cov-mean(fuels.cov)) fuels.cov <- round(fuels.cov, digits=3) jdate <- as.matrix(jdate - mean(jdate)) jdate <- round(jdate, digits=3) jdate <- as.matrix(jdate, 76, 3) #look at correlation between coefs covs <- cbind(trees.cov, snags.cov, fuels.cov, pcc) covs <- as.matrix(covs) dimnames(covs) <- list(NULL, c("trees", "snags", "fuels", "pcc")) pairs(covs) cor(covs) ## tree density and pcc are NOT correlated in 2010, but we'll leave pcc out to be consistent with the 2003 model # y N K Data needed for WinBUGS <- as.matrix(yrwa.histories10) <- length(y[ ,1]) <- length(y[1, ]) data <- list("y", "N", "K", "jdate", "trees.cov", "snags.cov","fuels.cov") # initial values for WinBUGS zst<-rbinom(numsites,1,0.5) pst <- runif(1,.05,.5) psist <- runif(1, .05, .5) #beta.st <- rnorm(1,0,1) inits <- function(){ list(a0 = pst, b0 = psist, z=zst) } # parameters to trace parameters <- c("p0", "psi0", "b1", "a1", "a2", "a3", "b2", "a4", "b3") # Write WinBUGS model to a text file Brittany Mosher filename = "model.txt" cat(" model{ # prior for occupancy b0 ~ dunif(0,1) logit(psi0) <- b0 # prior for p a0 ~ dunif(0,1) logit(p0) <- a0 #priors for regression coefficients a1 ~ dunif(-10,10) a2 ~ dunif(-10,10) a3 ~ dunif(-10,10) a4 ~ dunif(-10,10) a5 ~ dunif(-10,10) b1 b2 b3 b4 ~ ~ ~ ~ dunif(-5,5) dunif(-5,5) dunif(-5,5) dunif(-5,5) # observation model for (i in 1:N){ for (k in 1:K){ logit(p[i,k]) <- a0 + a1*jdate[i,k] + a4*fuels.cov[i] muy[i,k] <- z[i]*p[i,k] y[i,k] ~ dbern(muy[i,k]) } a2*trees.cov[i] + a3*snags.cov[i] + # state model z[i] ~ dbern(psi[i]) # state model logit(psi[i]) <- b0 + b1*trees.cov[i] + b2*snags.cov[i] + b3*fuels.cov[i] } } ", fill=TRUE, file=filename) # Call WinBUGS, supply arguments, store results in an object named “out” nt = 1 nc = 3 nb = 500 ni = 1500 out <- bugs(data, inits, parameters, "model.txt", n.thin=nt, n.chains=nc, n.burnin=nb, n.iter=ni, bugs.directory = "C:/Program Files (x86)/WinBUGS14/", debug=T) # Clean up. rm(data, filename, inits, K, N, nb, nc, ni, nt, parameters, psist, pst, y) Brittany Mosher

1/--страниц