M-array CLS model with 3 survival parameters

[Trevor Davies]

Hello All,

I'm still having difficultly getting this model to recapture parameters correctly. 

 I modified the data simulator used in the "integrated populations models"  textbook to ensure I was doing it correctly (I'm not).  I've added the data simulator to the end of this email.

To try to better understand the indexing I've isolated each m-array so as to (hopefully) better understand the issues.  I believe the issue is with the pi.1 line in the likelihood.

Thank you very much for the assistance.

Trevor

##################################
#R CODE 

######################################

################ BUGS code
library(tidyverse)
library(jagsUI)
library(readxl)
library(IPMbook)
library(RMark)
library(secr)
library(readxl)
library(janitor)

## Three survival parameters 

## Simulate data

data.3phi.1p <- data_gen(phi2=c(0.8,0.6,0.4),p=.5,nmarked=1000,nyears=11)
## Put in m-array format

marr <- marrayAge(ch=data.3phi.1p, age=1,mAge=3)

jags.data <- list(marr.1=marr[,,1], marr.2=marr[,,2], marr.a=marr[,,3], nyears=dim(marr)[2], rel.1=rowSums(marr[,,1]), rel.2=rowSums(marr[,,2]),rel.a=rowSums(marr[,,3]))
str(jags.data)

inits <- function(){list(phi1.const=runif(1, 0, 1))}
# Parameters monitored
parameters <- c("phi1.const", "phi2.const", "phia.const", "p.const")
# MCMC settings
ni <- 3000; nb <- 1000; nc <- 3; nt <- 1; na <- 1000
# Call JAGS, check convergence and summarize posteriors
model_3phi_1p <- jags(jags.data, inits, parameters, "model_3phi_1p.stan", n.iter=ni, n.burnin=nb, n.chains=nc,   n.thin=nt, n.adapt=na, parallel=TRUE)
traceplot(model_3phi_1p) # Not shown
print(model_3phi_1p)

############## 3phi; 1p CJS model  CODE

## NOT RECAPTURING PARAMETERS CORRECTLY.

##########################################

model {
# Priors and linear models
phi1.const ~ dunif(0, 1)
phi2.const ~ dunif(0, 1)
phia.const ~ dunif(0, 1)

p.const ~ dunif(0, 1)

for (t in 1:(nyears-1)){
phi1[t] <- phi1.const
phi2[t] <- phi2.const
phia[t] <- phia.const
p[t] <- p.const
}
# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(nyears-1)){
marr.1[t,1:nyears] ~ dmulti(pi.1[t,], rel.1[t])
marr.2[t,1:nyears] ~ dmulti(pi.2[t,], rel.2[t])
marr.a[t,1:nyears] ~ dmulti(pi.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
for (t in 1:(nyears-1)){
# Main diagonal
q[t] <- 1 - p[t] # Probability of non-recapture
pi.1[t,t] <- phi1[t] * p[t]
pi.2[t,t] <- phi2[t] * p[t]
pi.a[t,t] <- phia[t] * p[t]
}

############### Likelihood that I am having a hard time with
# Juvenile age 1: released at t, recaptured at j

for (t in 1:(nyears-1)){
for (j in (t+1):(nyears-1)){
pi.1[t,j] <- phi1[t] * phi2[t+1] * prod(phia[(t+1):j]) * prod(q[t:(j-1)]) * p[j]
} #j
}
# Sub-adult at captured at age two and released at age two
for (t in 1:(nyears-1)){
for (j in (t+1):(nyears-1)){
pi.2[t,j] <- phi2[t] * prod(phia[(t+1):j]) * prod(q[t:(j-1)]) * p[j]
} #j
}

# Adults captured 3+ and released at 3+
for (t in 1:(nyears-1)){
for (j in (t+1):(nyears-1)){
pi.a[t,j] <- prod(phia[t:j]) * prod(q[t:(j-1)]) * p[j]
} #j
}
#############################################################

# Below main diagonal
for (t in 1:(nyears-1)){
for (j in 1:(t-1)){
pi.1[t,j] <- 0
pi.2[t,j] <- 0
pi.a[t,j] <- 0
} #j
} #t

# Last column: probability of non-recapture
for (t in 1:(nyears-1)){
pi.1[t,nyears] <- 1-sum(pi.1[t,1:(nyears-1)])
pi.2[t,nyears] <- 1-sum(pi.2[t,1:(nyears-1)])
pi.a[t,nyears] <- 1-sum(pi.a[t,1:(nyears-1)])
}
}

#######################   

## DATA GENERATION FUNCTION

#######################

data_gen <- function(phi2=c(0.4,.6,.8,.95),p=0.5,nmarked=10,nyears=11){

## Generate generic survival vector
  s.length <- length(phi2)  
  s.vec <- head(phi2,s.length-1)
  const.s <- rep(tail(phi2,1),nyears-s.length)
  phi <- c(s.vec,const.s)
#################
 
f <- rep(1:(nyears-1), each=nmarked)
nind <- length(f) # Total number of marked individuals
# State or ecological process
z <- array(NA, dim=c(nind, nyears)) # Empty alive/dead matrix
# Initial conditions: all individuals alive at f(i)
for (i in 1:nind){
  z[i,f[i]] <- 1
}

set.seed(1) # Initialize the RNGs in R
# Propagate alive/dead process forwards via transition rule:
# Alive individuals survive with probability phi
for (i in 1:nind){
  k <- 1
  for (t in (f[i]+1):nyears){
    z[i,t] <- rbinom(1, 1, z[i,t-1] * phi[k])
    k <- k+1
  } #t
}

# Observation process: simulate observations
y <- array(0, dim=c(nind, nyears))
for (i in 1:nind){
  y[i,f[i]] <- 1
  for(t in (f[i]+1):nyears){
    y[i,t] <- rbinom(1, 1, z[i,t] * p)
  } #t
} #i
return(y)
}

[Trevor Davies]

Seems like I answered my own question.  Putting my code here to close things up and hopefully it'll save someone some time in the future.  

If there's an easier way to code this that I missed please let me know (besides obviously removing the redundant loops that I added for my own clarrity)..

model {
# Priors and linear models
phi1.const ~ dunif(0, 1)
phi2.const ~ dunif(0, 1)
phia.const ~ dunif(0, 1)

p.const ~ dunif(0, 1)

for (t in 1:(nyears-1)){
phi1[t] <- phi1.const
phi2[t] <- phi2.const
phia[t] <- phia.const

p[t] <- p.const # recapture

}
# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(nyears-1)){
marr.1[t,1:nyears] ~ dmulti(pi.1[t,], rel.1[t])
marr.2[t,1:nyears] ~ dmulti(pi.2[t,], rel.2[t])
marr.a[t,1:nyears] ~ dmulti(pi.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
for (t in 1:(nyears-1)){
# Main diagonal
q[t] <- 1 - p[t] # Probability of non-recapture
pi.1[t,t] <- phi1[t] * p[t]
pi.2[t,t] <- phi2[t] * p[t]
pi.a[t,t] <- phia[t] * p[t]
}

# Main diagonal #2 for juveniles (Only 9 rows to fill)
for (t in 2:(nyears-1)){
pi.1[t-1,t] <- phi1[t] * q[t] * phi2[t] * p[t]  

}

# Juvenile age 1: released at t, recaptured at j

# Only 8 rows rows to fill
for (t in 1:(nyears-2)){
for (j in (t+2):(nyears-1)){
pi.1[t,j] <- phi1[t+2] * phi2[t+2] * prod(phia[(t+2):j]) * prod(q[t:(j-1)]) * p[j]

} #j
}


# Sub-adult at captured at age two and released at age two

# Fill entire 10 x 10; nyears == 11

for (t in 1:(nyears-1)){
for (j in (t+1):(nyears-1)){
pi.2[t,j] <- phi2[t] * prod(phia[(t+1):j]) * prod(q[t:(j-1)]) * p[j]
} #j
}

# Adults captured 3+ and released at 3+

# Fill entire 10 x 10; nyears == 11

for (t in 1:(nyears-1)){
for (j in (t+1):(nyears-1)){
pi.a[t,j] <- prod(phia[t:j]) * prod(q[t:(j-1)]) * p[j]
} #j
}
#############################################################


# Below main diagonal
for (t in 1:(nyears-1)){
for (j in 1:(t-1)){
pi.1[t,j] <- 0
pi.2[t,j] <- 0
pi.a[t,j] <- 0
} #j
} #t

# Last column: probability of non-recapture
for (t in 1:(nyears-1)){
pi.1[t,nyears] <- 1-sum(pi.1[t,1:(nyears-1)])
pi.2[t,nyears] <- 1-sum(pi.2[t,1:(nyears-1)])
pi.a[t,nyears] <- 1-sum(pi.a[t,1:(nyears-1)])
}
}

[Michael Schaub]

Hi Trevor,

 

But I think there are still errors in your code – the time indices are not correct. I have added the lines of code below in red, which I hope is correct.

 

Kind regards

 

Michael

pi.1[t-1,t] <- phi1[t] * q[t] * phi2[t+1] * p[t+1]  

}

# Juvenile age 1: released at t, recaptured at j
# Only 8 rows rows to fill
for (t in 1:(nyears-2)){
for (j in (t+2):(nyears-1)){
pi.1[t,j] <- phi1[t+2] * phi2[t+2] * prod(phia[(t+2):j]) * prod(q[t:(j-1)]) * p[j]

pi.1[t,j] <- phi1[t] * phi2[t+1] * prod(phia[(t+2):j]) * prod(q[t:(j-1)]) * p[j]

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/hmecology/8ecf1cab-f5f0-46eb-98d9-be70eb4fe03fn%40googlegroups.com.