Can anyone help me with my issue? I simulated using a state-space model to estiamte recovery(p) and survival (phi) by age class.
My estimates are inaccurate.
I get the following:
beta [1] = 0.331
beta [2] = 0.634
delta [1] = 0.469
delta [2] = 0.529
According to my simulated data (or what I think should happen)
beta [1] = 0.8
beta [2] = 0.8
delta [1] = 0.2
delta [2] = 0.1
Any help would be extremely appreciated! Also let me know if you can't read my code.
### Finding first capture year
get.first <- function(x)min(which(x!=0))
known.state.mr <- function (mr) {
state <- matrix (NA, nrow=dim(mr)[1], ncol= dim(mr)[2])
rec <- which(rowSums(mr)==2)
for (i in 1: length(rec)){
n1<- min(which(mr[rec[i],]==1)) #n1 time of capture
n2<- max(which(mr[rec[i],]==1)) # Recovery time
state[rec[i],n1:n2]<-1
state[rec[i],n1]<-NA
state[rec[i],n2:dim(mr)[2]]<- 0
}
return(state)
}
mr.init.z<-function(mr){
ch<-matrix(NA, nrow=dim(mr)[1],ncol=dim(mr)[2])
rec<- which(rowSums(mr)==1)
for(i in 1:length(rec)){
n1<-which(mr[rec[i],]==1)
ch[rec[i], n1:dim(mr)[2]]<-0
ch[rec[i],n1] <- NA
}
return (ch)
}
# Define function to simulate mark-recovery data
simul.mr <- function(S, R, marked){
n.occasions <- dim(S)[2]
MR <- matrix(NA, ncol = n.occasions+1, nrow = sum(marked))
# Define a vector with the occasion of marking
mark.occ <- rep(1:n.occasions, marked)
# Fill the CH matrix
for (i in 1:sum(marked)){
MR[i, mark.occ[i]] <- 1
for (t in mark.occ[i]:n.occasions){
# Bernoulli trial: has individual survived occasion?
sur <- rbinom(1, 1, S[i,t])
if (sur==1) next
rp <- rbinom(1, 1, R[i,t])
if (rp==0){
MR[i,t+1] <- 0
break
}
if (rp==1){
MR[i,t+1] <- 1
break
}
} #t
} #i
# Replace the NA in the file by 0
MR[which(is.na(MR))] <- 0
return(MR)
}
########### Simulate Data #########################
# Juvenile Matrix
n.occasions <- 14 # Number of release occasions
marked <- rep(50, n.occasions) # Annual number of newly
marked individuals
as <- rep(0.8, n.occasions)
ar <- rep(0.2, n.occasions)
# Define matrices with survival and recovery probabilities
aS <- matrix(as, ncol = n.occasions, nrow = sum(marked))
aR <- matrix(ar, ncol = n.occasions, nrow = sum(marked))
# Execute function
MR.A <- simul.mr(aS, aR, marked)
### Juvenile Matrix #####
js <- rep(0.8, n.occasions)
jr <- rep(0.1, n.occasions)
# Define matrices with survival and recovery probabilities
jS <- matrix(js, ncol = n.occasions, nrow = sum(marked))
jR <- matrix(jr, ncol = n.occasions, nrow = sum(marked))
# Execute function
MR.J <- simul.mr(jS, jR, marked)
get.first <- function(x) min(which(x!=0))
f.j <- apply(MR.J, 1, get.first)
f.a <- apply(MR.A, 1, get.first)
# Create matrices X indicating age classes
x.j <- matrix(NA, ncol = dim(MR.J)[2]-1, nrow = dim(MR.J)[1])
x.a <- matrix(NA, ncol = dim(MR.A)[2]-1, nrow = dim(MR.A)[1])
for (i in 1:dim(MR.J)[1]){
for (t in f.j[i]:(dim(MR.J)[2]-1)){
x.j[i,t] <- 2
x.j[i,f.j[i]] <- 1
} #t
} #i
for (i in 1:dim(MR.A)[1]){
for (t in f.a[i]:(dim(MR.A)[2]-1)){
x.a[i,t] <- 2
} #t
} #i
MR <- rbind(MR.J, MR.A)
f <- c(f.j, f.a)
x <- rbind(x.j, x.a)
# Specify model in BUGS language
sink("Model_MR")
cat("
model {
# Priors and constraints
for (i in 1:nind){
for (t in f[i]:(n.occasions-1)){
phi[i,t] <- beta[x[i,t]]
p[i,t] <- delta[x[i,t]]
} #t
} #i
for (u in 1:2){
beta[u] ~ dunif(0,1) # Priors for age-specific survival
delta[u] ~ dunif(0,1)
}
mean.p ~ dunif(0, 1) # Prior for mean recapture
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,f[i]] <- 1
for (t in (f[i]+1):n.occasions){
# State process
z[i,t] ~ dbern(mu1[i,t])
mu1[i,t] <- phi[i,t-1] * z[i,t-1]
# Observation process
y[i,t] ~ dbern(mu2[i,t])
mu2[i,t] <- p[i,t-1] * (z[i,t-1]-z[i,t])
} #t
} #i
}
",fill = TRUE)
sink()
# Bundle data
bugs.data <- list(y = MR,
f = f,
nind = dim(MR)[1],
n.occasions = dim(MR)[2],
z = known.state.mr(MR),
x = x)
# Initial values
inits <- function(){list(z = mr.init.z(MR),
beta = runif(2, 0, 1),
mean.p = runif(1, 0, 1))}
# Parameters monitored
parameters <- c("beta", "mean.p","delta")
# MCMC settings
ni <- 2000
nt <- 3
nb <- 1000
nc <- 3
# Call WinBUGS from R (BRT 3 min)
cjs.age <- jags(bugs.data, inits, parameters, "Model_MR",
n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
print(cjs.age, digits = 3
I figured out my issue... I was simulating the data wrong. Please don't come after me! Here's the corrected code.
The issue was I didn't simulate transition probability for when juveniles become adults.
Now you know!
# Juvenile Matrix
n.occasions <- 10 # Number of release occasions
marked.j <- rep(50, n.occasions) # Annual number of newly marked individuals
marked.a <- rep(50, n.occasions)
marked.j <- rep(200, n.occasions-1) # Annual number of newly marked juveniles
marked.a <- rep(30, n.occasions-1) # Annual number of newly marked adults
phi <- rep(0.8, n.occasions-1) # Recapture
p.juv <- 0.3 # Juvenile annual survival
p.ad <- 0.65
p.j <- c(p.juv, rep(p.ad, n.occasions-2))
p.a <- rep(p.ad, n.occasions-1)
P.J <- matrix(0, ncol = n.occasions-1, nrow = sum(marked.j))
for (i in 1:length(marked.j)){
P.J[(sum(marked.j[1:i])-marked.j[i]+1):sum(marked.j[1:i]),i:(n.occasions-1)] <- matrix(rep(p.j[1:(n.occasions-i)],marked.j[i]), ncol = n.occasions-i, byrow = TRUE)
}
PHI.J <- matrix(rep(phi, sum(marked.j)), ncol = n.occasions-1, nrow = sum(marked.j), byrow = TRUE)
P.A <- matrix(rep(p.a, sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE)
PHI.A <- matrix(rep(phi, sum(marked.a)), ncol = n.occasions-1, nrow = sum(marked.a), byrow = TRUE)
# Execute function
MR.A <- simul.mr(PHI.A, P.A, marked.a)
# Execute function
MR.J <- simul.mr(PHI.J, P.J, marked.j)
get.first <- function(x) min(which(x!=0))
f.j <- apply(MR.J, 1, get.first)
f.a <- apply(MR.A, 1, get.first)
# Create matrices X indicating age classes
x.j <- matrix(NA, ncol = dim(MR.J)[2]-1, nrow = dim(MR.J)[1])
x.a <- matrix(NA, ncol = dim(MR.A)[2]-1, nrow = dim(MR.A)[1])
for (i in 1:dim(MR.J)[1]){
for (t in f.j[i]:(dim(MR.J)[2]-1)){
x.j[i,t] <- 2
x.j[i,f.j[i]] <- 1
} #t
} #i
for (i in 1:dim(MR.A)[1]){
for (t in f.a[i]:(dim(MR.A)[2]-1)){
x.a[i,t] <- 2
} #t
} #i
MR <- rbind(MR.J, MR.A)
f <- c(f.j, f.a)
x <- rbind(x.j, x.a)