Primary Investigator: Althea ArchMiller
Collaborators: J. Fieberg, B. Dorazio, K. St. Clair
Description This JAGS model takes multiple years of operational survey data and “fixed” sightability model parameters (i.e., \(\Phi^{(s)}=\beta_{0g}^{(s)},\beta_{1g}^{(s)}\)) to estimate abundance for each year using temporally smoothed spline model for mean \(\psi_{h,i,t}\) (i.e., \(\mu_{h,i,t}\)) with stratum-specific intercepts and exchangeable random effects for lambda (\(\lambda_{h,i,j,t}\)).
# Documentation-related
library(devtools)
library(knitr)
library(ezknitr)
Operational survey model
model.TEM <- function(){
# Visual obstruction
for(k in 1:nyrL){ # years
for(i in 1:nstrataL){ # strata
mu.voc[i,k] ~ dunif(0.01,0.99)
lnrho.voc[i,k] ~ dnorm(5,1)%_%T(0.01,10)
rho.voc[i,k] <- exp(lnrho.voc[i,k])
a.voc[i,k] <- mu.voc[i,k]*rho.voc[i,k]
b.voc[i,k] <- rho.voc[i,k]-(mu.voc[i,k]*rho.voc[i,k])
}
}
#### Priors
b0psi ~ dnorm(0,0.04)
b1psi ~ dnorm(0,0.04)
b2psi ~ dnorm(0,0.04)
b3psi ~ dnorm(0,0.04)
b4psi ~ dnorm(0,0.04)
b5psi ~ dnorm(0,0.04)
dummy1h <- c(1,0,0)
dummy2h <- c(0,1,0)
dummy3h <- c(0,0,1)
for(k in 1:nyrL){ # years
for(i in 1:nstrataL){ # strata
logit.mu[i,k] <-
b0psi*dummy1h[i] +
b1psi*splines[k,1] +
b2psi*splines[k,2] +
b3psi*splines[k,3] +
b4psi*dummy2h[i] +
b5psi*dummy3h[i]
mu.psi[i,k] <- exp(logit.mu[i,k])/(1+exp(logit.mu[i,k]))
lnrho.psi[i,k] ~ dnorm(5,1)%_%T(0.01,10)
rho.psi[i,k] <- exp(lnrho.psi[i,k])
a.psi[i,k] <- mu.psi[i,k]*rho.psi[i,k]
b.psi[i,k] <- rho.psi[i,k]-(mu.psi[i,k]*rho.psi[i,k])
}
}
# The probability that each group "belongs to N" (data augmentation)
for(i in 1:Nplots.yr){
psi[i] ~ dbeta(a.psi[h.plots[i],yr.plots[i]],
b.psi[h.plots[i],yr.plots[i]])%_%T(0.0001,0.99)
}
# Prior for lambda, stratum-specific, on which y depends as on 1+pois(lambda[i])
# Exchangeable:
for(i in 1:nstrataL){
sigma.lam[i] ~ dunif(0,3)
prec.lam[i] <- 1/(sigma.lam[i]*sigma.lam[i])
mu.lam[i] ~ dnorm(0,0.1)
}
for(k in 1:nyrL){
for(i in 1:nstrataL){
delta.lam[i,k] ~ dnorm(0, prec.lam[i])%_%T(0.001, 10)
loglam[i,k] <- mu.lam[i] + delta.lam[i,k]
lam[i,k] <- exp(loglam[i,k])
}
}
#### LIKELIHOOD
for(j in 1:Ngroups){ # processed for sampled plots
x[j] ~ dbeta(a.voc[h[j],yr[j]],b.voc[h[j],yr[j]])%_%T(0,0.999)
q[j] ~ dbin(psi[plots[j]],1)
# probability of seeing moose group given voc
g[j] <- exp(b0g + b1g*x[j])/(1+exp(b0g + b1g*x[j]))
# conditional probability: if q = 0, g should be 0
pseen[j] <- q[j]*g[j] # eqn 14
z[j] ~ dbin(pseen[j],1) # eqn 14
# model for number of moose (minus 1)
ym1[j] ~ dpois(lam[h[j],yr[j]]) # eqn 12
y[j] <- ym1[j] + 1
# number of moose (tau) conditional if q not 0
taus[j] <- q[j]*y[j]
}
##### Derived quantities
tau.samp05 <- sum(taus[1:ny05])
tau.samp06 <- sum(taus[(ny05+1):ny06]) # tau for year 2006
tau.samp07 <- sum(taus[(ny06+1):ny07]) # tau for year 2007
tau.samp08 <- sum(taus[(ny07+1):ny08]) # tau for year 2008
tau.samp09 <- sum(taus[(ny08+1):ny09]) # tau for year 2009
tau.samp10 <- sum(taus[(ny09+1):ny10]) # tau for year 2010
tau.samp11 <- sum(taus[(ny10+1):ny11]) # 2011
tau.samp12 <- sum(taus[(ny11+1):ny12]) # 2012
tau.samp13 <- sum(taus[(ny12+1):ny13]) # 2013
tau.samp14 <- sum(taus[(ny13+1):ny14]) # 2014
tau.samp15 <- sum(taus[(ny14+1):ny15]) # 2015
tau.samp16 <- sum(taus[(ny15+1):ny16]) # 2016
}
Session Information
devtools::session_info()
## Session info --------------------------------------------------------------
## setting value
## version R version 3.3.2 (2016-10-31)
## system x86_64, darwin13.4.0
## ui RStudio (1.0.136)
## language (EN)
## collate en_US.UTF-8
## tz America/Chicago
## date 2017-03-29
## Packages ------------------------------------------------------------------
## package * version date source
## devtools * 1.12.0 2016-06-24 CRAN (R 3.3.0)
## digest 0.6.12 2017-01-27 CRAN (R 3.3.2)
## evaluate 0.10 2016-10-11 CRAN (R 3.3.0)
## ezknitr * 0.6 2016-09-16 CRAN (R 3.3.0)
## knitr * 1.15.1 2016-11-22 CRAN (R 3.3.2)
## magrittr 1.5 2014-11-22 CRAN (R 3.3.0)
## markdown 0.7.7 2015-04-22 CRAN (R 3.3.0)
## memoise 1.0.0 2016-01-29 CRAN (R 3.3.0)
## R.methodsS3 1.7.1 2016-02-16 CRAN (R 3.3.0)
## R.oo 1.21.0 2016-11-01 CRAN (R 3.3.0)
## R.utils 2.5.0 2016-11-07 CRAN (R 3.3.0)
## rstudioapi 0.6 2016-06-27 CRAN (R 3.3.0)
## stringi 1.1.2 2016-10-01 CRAN (R 3.3.0)
## stringr 1.1.0 2016-08-19 CRAN (R 3.3.0)
## withr 1.0.2 2016-06-20 CRAN (R 3.3.0)
spun with: ezknitr::ezspin(“programs/ms_programs/a_define_TS_model.R”, out_dir = “output”, fig_dir = “figures”, keep_md = F)