IndyPM Overlap Models Pub
elaine metz
2025-06-03
IPM 1 (highest complexity)
This model has demographic stochasticity and the survival probability
and expected number of immigrants is variable over the years. The full
model is included at the bottom of this doc.
I have found several issues with this model/the data so take this with a
grain of salt please.
#head(out1$summary, 14)
#knitr::kable(out1$summary, caption = "results summary")
df <- as.data.frame(out1$summary)
paged_table(df, options = list(rows.print = 14, cols.print = 7))Mean novice and experienced breeder suvival is 0.33 and 0.35–very similar. Mean productivity is 1.22, which seems really low to me. More below when looking at estimates through time.
Assess fit of first IPM with DS and annual variation in expected immigrants and survivial
Fits good. P values are not close to 0 or 1.
Looking at IPM 1 through time
## [1] 1.097426## 2.5% 97.5%
## 0.674879 1.529376## [1] 0.675245## 2.5% 97.5%
## 0.2324537 0.9545400
Above Figure shows annual estimates (posterior means and 95% credible
intervals) of novice and experienced breeders in CERW population in IN,
USA. Immigration rate (I(t+1)/N(t)) and total productivity are
derived.
Overall nest success did not vary much over time, but conditional productivity did. Conditional productivity is the productivity of a nest given nest success.
These graph are focusing on the weakest parameters in the model.
Population composition over time.
Population structure estimates under 4 models of reducing
complexity.
Results Visualization
Models
data <- readRDS(here('data/data_indy_ovrlp.RDS'))
str(data)## List of 7
## $ ch : num [1:233, 1:11] 0 0 0 0 0 0 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:11] "2013" "2014" "2015" "2016" ...
## $ age : num [1:233(1d)] 1 2 2 1 2 1 2 1 2 2 ...
## $ count : int [1:11] 125 77 67 75 52 41 35 40 48 45 ...
## $ f : int [1:232] NA 0 NA NA 0 3 3 0 0 0 ...
## $ d : int [1:232] 1 0 1 1 0 0 0 0 0 0 ...
## $ censored: num [1:232] 2 1 3 2 1 4 4 1 1 1 ...
## $ year : int [1:232] 2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...IPM 1 With Demographic Stochasticity AND survival & expected immigrants Variable over years
Below, not actually juvenile and adult stages, instead, novice and adult breeders.
recap <- c(1,2,3,4,5,6,7,8,9,10)
censored <- data$censored
jags.data <- list(n.years=dim(marr)[2], marr.j=marr[,,1], marr.a=marr[,,2], rel.j=rel.j,
rel.a=rel.a, recap=recap, n.recap=max(recap), f=data$f, d=data$d, censored=data$censored, year=data$year-2012,
C=data$count, pNinit=dUnif(1, 20))
str(jags.data)
# ipm 1 ----
cat(file="model1.txt",
model {
# Priors and linear models
for (t in 1:(n.years-1)){
logit.phij[t] ~ dnorm(l.mean.phij, tau.phij)
phij[t] <- ilogit(logit.phij[t])
logit.phia[t] ~ dnorm(l.mean.phia, tau.phia)
phia[t] <- ilogit(logit.phia[t])
pj[t] <- p.j[recap[t]]
pa[t] <- p.a[recap[t]]
}
for (t in 1:n.years){
logit.nu[t] ~ dnorm(l.mean.nu, tau.nu)
nu[t] <- ilogit(logit.nu[t])
log.kappa[t] ~ dnorm(l.mean.kappa, tau.kappa)
kappa[t] <- exp(log.kappa[t])
log.omega[t] ~ dnorm(l.mean.omega, tau.omega)
omega[t] <- exp(log.omega[t])
}
for (t in 1:(n.recap)){
p.j[t] <- mean.p[1]
p.a[t] <- mean.p[2]
}
# Fix at 0 recapture probability in years without capture/resighting
#p.j[n.recap] <- 0 #######!!! not needed here.
#p.a[n.recap] <- 0 ########!!!!
mean.phij ~ dunif(0, 1)
l.mean.phij <- logit(mean.phij)
mean.phia ~ dunif(0, 1)
l.mean.phia <- logit(mean.phia)
mean.nu ~ dunif(0, 1)
l.mean.nu <- logit(mean.nu)
mean.kappa ~ dunif(1, 10)
l.mean.kappa <- log(mean.kappa)
mean.omega ~ dunif(0.01, 30)
l.mean.omega <- log(mean.omega)
for (i in 1:2){
mean.p[i] ~ dunif(0, 1)
}
sigma.phij ~ dunif(0, 5)
tau.phij <- pow(sigma.phij, -2)
sigma.phia ~ dunif(0, 5)
tau.phia <- pow(sigma.phia, -2)
sigma.nu ~ dunif(0, 5)
tau.nu <- pow(sigma.nu, -2)
sigma.kappa ~ dunif(0, 5)
tau.kappa <- pow(sigma.kappa, -2)
sigma.omega ~ dunif(0, 5)
tau.omega <- pow(sigma.omega, -2)
sigma.f ~ dunif(0, 5)
tau.f <- pow(sigma.f, -2)
# Population count data (state-space model)
# Model for the initial population size: uniform priors
R[1] ~ dcat(pNinit) # Local recruits
S[1] ~ dcat(pNinit) # Surviving adults
I[1] ~ dpois(omega[1]) # Immigrants --> could be using more informative
# Process model over time: our model of population dynamics
for (t in 2:n.years){
R[t] ~ dpois(nu[t-1] * kappa[t-1]/2 * phij[t-1] * N[t-1]) # Local recruits
S[t] ~ dbin(phia[t-1], N[t-1]) # Surviving adults
I[t] ~ dpois(omega[t]) # Immigrants
}
# Observation model
for (t in 1:n.years){
N[t] <- S[t] + R[t] + I[t] # counts
C[t] ~ dpois(N[t])
# GoF for population count data: mean absolute percentage error
C.pred[t] ~ dpois(N[t])
disc.C[t] <- pow(((C[t] - N[t]) / C[t]) * ((C[t] - N[t]) / (C[t] +
0.001)), 0.5) # Add a small number to avoid potential division by 0
discN.C[t] <- pow(((C.pred[t] - N[t]) / (C.pred[t] + 0.001)) *
((C.pred[t] - N[t]) / (C.pred[t] + 0.001)), 0.5)
}
fit.C <- 100 / n.years * sum(disc.C)
fitN.C <- 100 / n.years * sum(discN.C)
# Productivity data (zero-inflated normal with right censoring)
for (i in 1:length(f)){
z[i] ~ dbern(nu[year[i]])
f[i] ~ dnorm(z[i] * kappa[year[i]], tau.f)
delta[i] <- step(f[i] - censored[i])
d[i] ~ dbern(z[i] * delta[i])
# GoF for productivity data: deviance
z.pred[i] ~ dbern(nu[year[i]])
f.pred[i] ~ dnorm(z.pred[i] * kappa[year[i]], tau.f)
f.exp[i] <- kappa[year[i]] * nu[year[i]]
dev[i] <- (f[i] - f.exp[i]) * (f[i] - f.exp[i])
devN[i] <- (f.pred[i] - f.exp[i]) * (f.pred[i] - f.exp[i])
}
fit.f <- 2 * sum(dev)
fitN.f <- 2 * sum(devN)
# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1: (n.years-1)){
marr.j[t,1:n.years] ~ dmulti(pr.j[t,], rel.j[t])
marr.a[t,1:n.years] ~ dmulti(pr.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
for (t in 1:(n.years-1)){
# Main diagonal
qj[t] <- 1-pj[t]
qa[t] <- 1-pa[t]
pr.j[t,t] <- phij[t] * pj[t]
pr.a[t,t] <- phia[t] * pa[t]
# Above main diagonal
for (j in (t+1):(n.years-1)){
pr.j[t,j] <- phij[t] * prod(phia[(t+1):j]) * qj[t] * prod(qa[t:(j-1)]) * pa[j] / qa[t]
pr.a[t,j] <-prod(phia[t:j]) * prod(qa[t:(j-1)]) * pa[j]
} #j
# Below main diagonal
for (j in 1:(t-1)){
pr.j[t,j] <- 0
pr.a[t,j] <- 0
} #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.years-1)){
pr.j[t,n.years] <- 1-sum(pr.j[t,1:(n.years-1)])
pr.a[t,n.years] <- 1-sum(pr.a[t,1:(n.years-1)])
}
# GoF for cap.-recap. data: Freeman-Tukey test statistics
for (t in 1:(n.years-1)){
# Simulated m-arrays
marr.j.pred[t,1:n.years] ~ dmulti(pr.j[t,], rel.j[t])
marr.a.pred[t,1:n.years] ~ dmulti(pr.a[t,], rel.a[t])
# Expected values and test statistics
for (j in 1:n.years){
marr.jE[t,j] <- pr.j[t,j] * rel.j[t]
E.org.j[t,j] <- pow((pow(marr.j[t,j], 0.5) - pow(marr.jE[t,j], 0.5)), 2)
E.new.j[t,j] <- pow((pow(marr.j.pred[t,j], 0.5) - pow(marr.jE[t,j], 0.5)), 2)
marr.aE[t,j] <- pr.a[t,j] * rel.a[t]
E.org.a[t,j] <- pow((pow(marr.a[t,j], 0.5) - pow(marr.aE[t,j], 0.5)), 2)
E.new.a[t,j] <- pow((pow(marr.a.pred[t,j], 0.5) - pow(marr.aE[t,j], 0.5)), 2)
} #j
} #t
fit.CJS <- sum(E.org.j) + sum(E.org.a)
fitN.CJS <- sum(E.new.j) + sum(E.new.a)
# Compute derived parameters: total productivity and immigration rate
for (t in 1:n.years){
rho[t] <- nu[t] * kappa[t]
}
mean.rho <- mean.nu * mean.kappa
for (t in 1:(n.years-1)){
omega.rate[t] <- I[t+1] / N[t]
}
}
)
# Initial values
init.Juv <- censored+1
init.Juv[data$d==0] <- NA
s <- rep(1, length(data$f))
s[data$f==0] <- 0
inits <- function() {list(z=s, f=init.Juv)}
# Parameters monitored
parameters <- c("mean.phij", "mean.phia", "mean.nu", "mean.kappa", "mean.rho", "mean.omega",
"sigma.phij", "sigma.phia", "sigma.nu", "sigma.kappa", "sigma.omega", "sigma.f", "mean.p",
"fit.C", "fitN.C", "fit.f", "fitN.f", "fit.CJS", "fitN.CJS", "phij", "phia", "nu", "kappa",
"omega", "rho", "omega.rate", "R", "S", "I", "N")
# MCMC settings
ni <- 30000; nb <- 10000; nc <- 3; nt <- 10; na <- 1000
#ni <- 3000; nb <- 1000; nc <- 3; nt <- 1; na <- 1000 # ~~~ for testing
# Call JAGS from R (ART 3 min) and check convergence
#out1 <- jags(jags.data, inits, parameters, "model1.txt", n.iter=ni, n.burnin=nb, n.chains=nc,
# n.thin=nt, n.adapt=na, parallel=TRUE)IPM 2 Without Demographic Stochasticity AND survival & expected immigrants Variable over years
# ipm2 ----
# 1. Model with temporally variable survival and immigration, but without dem. stoch. (IPM2)
cat(file="model2.txt",
model {
# Priors and linear models
for (t in 1:(n.years-1)){
logit.phij[t] ~ dnorm(l.mean.phij, tau.phij)
phij[t] <- ilogit(logit.phij[t])
logit.phia[t] ~ dnorm(l.mean.phia, tau.phia)
phia[t] <- ilogit(logit.phia[t])
pj[t] <- p.j[recap[t]]
pa[t] <- p.a[recap[t]]
}
for (t in 1:n.years){
logit.nu[t] ~ dnorm(l.mean.nu, tau.nu)
nu[t] <- ilogit(logit.nu[t])
log.kappa[t] ~ dnorm(l.mean.kappa, tau.kappa)
kappa[t] <- exp(log.kappa[t])
log.omega[t] ~ dnorm(l.mean.omega, tau.omega)
omega[t] <- exp(log.omega[t])
}
for (t in 1:(n.recap)){ # Recapture probability !!!(remved -1)
p.j[t] <- mean.p[1]
p.a[t] <- mean.p[2]
}
# p.j[n.recap] <- 0
# p.a[n.recap] <- 0
mean.phij ~ dunif(0, 1)
l.mean.phij <- logit(mean.phij)
mean.phia ~ dunif(0, 1)
l.mean.phia <- logit(mean.phia)
mean.nu ~ dunif(0, 1)
l.mean.nu <- logit(mean.nu)
mean.kappa ~ dunif(1, 10)
l.mean.kappa <- log(mean.kappa)
mean.omega ~ dunif(0.01, 30)
l.mean.omega <- log(mean.omega)
for (i in 1:2){
mean.p[i] ~ dunif(0, 1)
}
sigma.phij ~ dunif(0, 5)
tau.phij <- pow(sigma.phij, -2)
sigma.phia ~ dunif(0, 5)
tau.phia <- pow(sigma.phia, -2)
sigma.nu ~ dunif(0, 5)
tau.nu <- pow(sigma.nu, -2)
sigma.kappa ~ dunif(0, 5)
tau.kappa <- pow(sigma.kappa, -2)
sigma.omega ~ dunif(0, 5)
tau.omega <- pow(sigma.omega, -2)
sigma.f ~ dunif(0, 5)
tau.f <- pow(sigma.f, -2)
# Population count data (state-space model)
# Model for the initial population size: uniform priors
R[1] ~ dcat(pNinit) # Local recruits
S[1] ~ dcat(pNinit) # Surviving adults
I[1] ~ dpois(omega[1]) # Immigrants
# Process model over time: our model of population dynamics
## here, removing demog stochastic in model.
for (t in 2:n.years){
R[t] <- nu[t-1] * kappa[t-1]/2 * phij[t-1] * N[t-1] # Local recruits
S[t] <- phia[t-1] * N[t-1] # Surviving adults
I[t] <- omega[t] # Immigrants
}
# Observation model
for (t in 1:n.years){
N[t] <- S[t] + R[t] + I[t] # Total number of breeding females
C[t] ~ dpois(N[t])
}
# Productivity data (zero-inflated normal with right censoring)
for (i in 1:length(f)){
z[i] ~ dbern(nu[year[i]])
f[i] ~ dnorm(z[i] * kappa[year[i]], tau.f)
delta[i] <- step(f[i] - censored[i])
d[i] ~ dbern(z[i] * delta[i])
}
# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(n.years-1)){
marr.j[t,1:n.years] ~ dmulti(pr.j[t,], rel.j[t])
marr.a[t,1:n.years] ~ dmulti(pr.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
for (t in 1:(n.years-1)){
# Main diagonal
qj[t] <- 1-pj[t]
qa[t] <- 1-pa[t]
pr.j[t,t] <- phij[t] * pj[t]
pr.a[t,t] <- phia[t] * pa[t]
# Above main diagonal
for (j in (t+1):(n.years-1)){
pr.j[t,j] <- phij[t] * prod(phia[(t+1):j]) * qj[t] * prod(qa[t:(j-1)]) * pa[j] / qa[t]
pr.a[t,j] <- prod(phia[t:j]) * prod(qa[t:(j-1)]) * pa[j]
} #j
# Below main diagonal
for (j in 1:(t-1)){
pr.j[t,j] <- 0
pr.a[t,j] <- 0
} #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.years-1)){
pr.j[t,n.years] <- 1-sum(pr.j[t,1:(n.years-1)])
pr.a[t,n.years] <- 1-sum(pr.a[t,1:(n.years-1)])
}
}
)
# Initial values
init.Juv <- censored+1
init.Juv[data$d==0] <- NA
s <- rep(1, length(data$f))
s[data$f==0] <- 0
inits <- function() {list(z=s, f=init.Juv)}
# Parameters monitored
parameters <- c("mean.phij", "mean.phia", "mean.nu", "mean.kappa", "mean.omega",
"sigma.phij", "sigma.phia", "sigma.nu", "sigma.kappa", "sigma.omega",
"sigma.f", "phij", "phia", "nu", "kappa", "omega", "R", "S", "I", "N", "mean.p")
# MCMC settings
ni <- 30000; nb <- 10000; nc <- 3; nt <- 10; na <- 1000
#ni <- 3000; nb <- 1000; nc <- 3; nt <- 1; na <- 1000 # ~~~ for testing
# Call JAGS from R (ART: 9.0 min)
#out2 <- jags(jags.data, inits, parameters, "model2.txt",
# n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=TRUE)IPM 3 With Demographic Stochasticity AND survival & expected immigrants Constant over years
# ipm 3 ----
# 2. Model with constant survival and immigration, but with dem. stoch. (IPM3)
cat(file="model3.txt",
model {
# Priors and linear models
## !!! Survival and expected immig Constant
for (t in 1:(n.years-1)){
phij[t] <- mean.phij
phia[t] <- mean.phia
pj[t] <- p.j[recap[t]]
pa[t] <- p.a[recap[t]]
}
for (t in 1:n.years){
logit.nu[t] ~ dnorm(l.mean.nu, tau.nu)
nu[t] <- ilogit(logit.nu[t])
log.kappa[t] ~ dnorm(l.mean.kappa, tau.kappa)
kappa[t] <- exp(log.kappa[t])
omega[t] <- mean.omega
}
for (t in 1:(n.recap)){ # Recapture probability !!!(remved -1)
p.j[t] <- mean.p[1]
p.a[t] <- mean.p[2]
}
#p.j[n.recap] <- 0
#p.a[n.recap] <- 0
mean.phij ~ dunif(0, 1)
mean.phia ~ dunif(0, 1)
mean.nu ~ dunif(0, 1)
l.mean.nu <- logit(mean.nu)
mean.kappa ~ dunif(1, 10)
l.mean.kappa <- log(mean.kappa)
mean.omega ~ dunif(0.01, 30)
for (i in 1:2){
mean.p[i] ~ dunif(0, 1)
}
sigma.nu ~ dunif(0, 5)
tau.nu <- pow(sigma.nu, -2)
sigma.kappa ~ dunif(0, 5)
tau.kappa <- pow(sigma.kappa, -2)
sigma.f ~ dunif(0, 5)
tau.f <- pow(sigma.f, -2)
# Population count data (state-space model)
# Model for the initial population size: uniform priors
R[1] ~ dcat(pNinit) # Local recruits
S[1] ~ dcat(pNinit) # Surviving adults
I[1] ~ dpois(omega[1]) # Immigrants
# Process model over time: our model of population dynamics
for (t in 2:n.years){
R[t] ~ dpois(nu[t-1] * kappa[t-1]/2 * phij[t-1] * N[t-1]) # Local recruits
S[t] ~ dbin(phia[t-1], N[t-1]) # Surviving adults
I[t] ~ dpois(omega[t]) # Immigrants
}
# Observation model
for (t in 1:n.years){
N[t] <- S[t] + R[t] + I[t] # Total number of breeding females
C[t] ~ dpois(N[t])
}
# Productivity data (zero-inflated normal with right censoring)
for (i in 1:length(f)){
z[i] ~ dbern(nu[year[i]])
f[i] ~ dnorm(z[i] * kappa[year[i]], tau.f)
delta[i] <- step(f[i] - censored[i])
d[i] ~ dbern(z[i] * delta[i])
}
# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(n.years-1)){
marr.j[t,1:n.years] ~ dmulti(pr.j[t,], rel.j[t])
marr.a[t,1:n.years] ~ dmulti(pr.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
for (t in 1:(n.years-1)){
# Main diagonal
qj[t] <- 1-pj[t]
qa[t] <- 1-pa[t]
pr.j[t,t] <- phij[t] * pj[t]
pr.a[t,t] <- phia[t] * pa[t]
# Above main diagonal
for (j in (t+1):(n.years-1)){
pr.j[t,j] <- phij[t] * prod(phia[(t+1):j]) * qj[t] * prod(qa[t:(j-1)]) * pa[j] / qa[t]
pr.a[t,j] <- prod(phia[t:j]) * prod(qa[t:(j-1)]) * pa[j]
} #j
# Below main diagonal
for (j in 1:(t-1)){
pr.j[t,j] <- 0
pr.a[t,j] <- 0
} #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.years-1)){
pr.j[t,n.years] <- 1-sum(pr.j[t,1:(n.years-1)])
pr.a[t,n.years] <- 1-sum(pr.a[t,1:(n.years-1)])
}
}
)
# Initial values
init.Juv <- censored+1
init.Juv[data$d==0] <- NA
s <- rep(1, length(data$f))
s[data$f==0] <- 0
inits <- function() {list(z=s, f=init.Juv)}
# Parameters monitored
parameters <- c("mean.phij", "mean.phia", "mean.nu", "mean.kappa", "mean.omega",
"sigma.nu", "sigma.kappa", "sigma.f", "phij", "phia", "nu", "kappa", "omega",
"R", "S", "I", "N", "mean.p")
# MCMC settings
ni <- 30000; nb <- 10000; nc <- 3; nt <- 10; na <- 1000
#ni <- 3000; nb <- 1000; nc <- 3; nt <- 1; na <- 1000 # ~~~ for testing
# Call JAGS from R (ART: 3.4 min)
#out3 <- jags(jags.data, inits, parameters, "model3.txt",
# n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=TRUE)IPM 4 Without Demographic Stochasticity AND survival & expected immigrants Constant over years
# ipm 4 ----
# 3. Model with constant survival and immigration, without dem. stoch. (IPM4)
cat(file="model4.txt",
model {
# Priors and linear models
## !!! Survival and expected immig CONSTANT
for (t in 1:(n.years-1)){
phij[t] <- mean.phij
phia[t] <- mean.phia
pj[t] <- p.j[recap[t]]
pa[t] <- p.a[recap[t]]
}
for (t in 1:n.years){
logit.nu[t] ~ dnorm(l.mean.nu, tau.nu)
nu[t] <- ilogit(logit.nu[t])
log.kappa[t] ~ dnorm(l.mean.kappa, tau.kappa)
kappa[t] <- exp(log.kappa[t])
omega[t] <- mean.omega
}
for (t in 1:(n.recap)){ # Recapture probability !!!(remved -1)
p.j[t] <- mean.p[1]
p.a[t] <- mean.p[2]
}
# p.j[n.recap] <- 0
# p.a[n.recap] <- 0
mean.phij ~ dunif(0, 1)
mean.phia ~ dunif(0, 1)
mean.nu ~ dunif(0, 1)
l.mean.nu <- logit(mean.nu)
mean.kappa ~ dunif(1, 10)
l.mean.kappa <- log(mean.kappa)
mean.omega ~ dunif(0.01, 30)
for (i in 1:2){
mean.p[i] ~ dunif(0, 1)
}
sigma.nu ~ dunif(0, 5)
tau.nu <- pow(sigma.nu, -2)
sigma.kappa ~ dunif(0, 5)
tau.kappa <- pow(sigma.kappa, -2)
sigma.omega ~ dunif(0, 5)
tau.omega <- pow(sigma.omega, -2)
sigma.f ~ dunif(0, 5)
tau.f <- pow(sigma.f, -2)
# Population count data (state-space model)
# Model for the initial population size: uniform priors
R[1] ~ dcat(pNinit) # Local recruits
S[1] ~ dcat(pNinit) # Surviving adults
I[1] ~ dpois(omega[1]) # Immigrants
# Process model over time: our model of population dynamics
## No deomg Stochastiticty
for (t in 2:n.years){
R[t] <- nu[t-1] * kappa[t-1]/2 * phij[t-1] * N[t-1] # Local recruits
S[t] <- phia[t-1] * N[t-1] # Surviving adults
I[t] <- omega[t] # Immigrants
}
# Observation model
for (t in 1:n.years){
N[t] <- S[t] + R[t] + I[t] # Total number
C[t] ~ dpois(N[t])
}
# Productivity data (zero-inflated normal with right censoring)
for (i in 1:length(f)){
z[i] ~ dbern(nu[year[i]])
f[i] ~ dnorm(z[i] * kappa[year[i]], tau.f)
delta[i] <- step(f[i] - censored[i])
d[i] ~ dbern(z[i] * delta[i])
}
# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(n.years-1)){
marr.j[t,1:n.years] ~ dmulti(pr.j[t,], rel.j[t])
marr.a[t,1:n.years] ~ dmulti(pr.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
for (t in 1:(n.years-1)){
# Main diagonal
qj[t] <- 1-pj[t]
qa[t] <- 1-pa[t]
pr.j[t,t] <- phij[t] * pj[t]
pr.a[t,t] <- phia[t] * pa[t]
# Above main diagonal
for (j in (t+1):(n.years-1)){
pr.j[t,j] <- phij[t] * prod(phia[(t+1):j]) * qj[t] * prod(qa[t:(j-1)]) * pa[j] / qa[t]
pr.a[t,j] <- prod(phia[t:j]) * prod(qa[t:(j-1)]) * pa[j]
} #j
# Below main diagonal
for (j in 1:(t-1)){
pr.j[t,j] <- 0
pr.a[t,j] <- 0
} #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.years-1)){
pr.j[t,n.years] <- 1-sum(pr.j[t,1:(n.years-1)])
pr.a[t,n.years] <- 1-sum(pr.a[t,1:(n.years-1)])
}
}
)
# Initial values
init.Juv <- censored+1
init.Juv[data$d==0] <- NA
s <- rep(1, length(data$f))
s[data$f==0] <- 0
inits <- function() {list(z=s, f=init.Juv)}
# Parameters monitored
parameters <- c("mean.phij", "mean.phia", "mean.nu", "mean.kappa", "mean.omega",
"sigma.nu", "sigma.kappa", "sigma.f", "phij", "phia", "nu", "kappa", "omega",
"R", "S", "I", "N", "mean.p")
# MCMC settings
ni <- 30000; nb <- 10000; nc <- 3; nt <- 10; na <- 1000
#ni <- 3000; nb <- 1000; nc <- 3; nt <- 1; na <- 1000 # ~~~ for testing
# Call JAGS from R (ART: 3.0 min)
#out4 <- jags(jags.data, inits, parameters, "model4.txt",
# n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=TRUE)