# Mark-recapture models with individual covariates

Hi all,

I am trying to construct some Cormack-Jolly_Seber (CJS) models in Stan, but am running into difficulties incorporating covariates. I am following code associated with the book “Bayesian Population Analysis” and the Stan code version of it supplied here. https://github.com/stan-dev/example-models/tree/master/BPA. I am also tyring to look at the code supplied here https://mc-stan.org/docs/2_18/stan-users-guide/mark-recapture-models.html, that appears to mirror the former link. The Stan users guide says that in the code, “All individuals are modeled as having the same capture probability, but this model could be easily generalized to use a logistic regression here based on individual-level inputs to be used as predictors.” In section 7.5.3 of the BPA book, it states that the individual mode "provides the base for modeling survival as a function of an individual covariate Xi (e.g. size of individual). I am interested in incorporating a flexible design matrix into the CJS model that allows for individual level covariates that are predictors of survivorship, but just not sure where to start. I am guessing it would be in the constraints secion of the code, but just not sure how to do it.

here is (I think) the relevant code from the BPA book.

``````// This models is derived from section 12.3 of "Stan Modeling Language
// User's Guide and Reference Manual"

functions {
int first_capture(int[] y_i) {
for (k in 1:size(y_i))
if (y_i[k])
return k;
return 0;
}

int last_capture(int[] y_i) {
for (k_rev in 0:(size(y_i) - 1)) {
// Compoud declaration was enabled in Stan 2.13
int k = size(y_i) - k_rev;
//      int k;
//      k = size(y_i) - k_rev;
if (y_i[k])
return k;
}
return 0;
}

matrix prob_uncaptured(int nind, int n_occasions,
matrix p, matrix phi) {
matrix[nind, n_occasions] chi;

for (i in 1:nind) {
chi[i, n_occasions] = 1.0;
for (t in 1:(n_occasions - 1)) {
// Compoud declaration was enabled in Stan 2.13
int t_curr = n_occasions - t;
int t_next = t_curr + 1;
/*
int t_curr;
int t_next;

t_curr = n_occasions - t;
t_next = t_curr + 1;
*/
t_curr = n_occasions - t;
t_next = t_curr + 1;
chi[i, t_curr] = (1 - phi[i, t_curr])
+ phi[i, t_curr] * (1 - p[i, t_next - 1]) * chi[i, t_next];
}
}
return chi;
}
}

data {
int<lower=0> nind;            // Number of individuals
int<lower=2> n_occasions;     // Number of capture occasions
int<lower=0,upper=1> y[nind, n_occasions];    // Capture-history
}

transformed data {
int n_occ_minus_1 = n_occasions - 1;
//  int n_occ_minus_1;
int<lower=0,upper=n_occasions> first[nind];
int<lower=0,upper=n_occasions> last[nind];

//  n_occ_minus_1 = n_occasions - 1;
for (i in 1:nind)
first[i] = first_capture(y[i]);
for (i in 1:nind)
last[i] = last_capture(y[i]);
}

parameters {
real<lower=0,upper=1> mean_phi;       // Mean survival
real<lower=0,upper=1> mean_p;         // Mean recapture
vector[nind] epsilon;
real<lower=0,upper=5> sigma;
// In case a weakly informative prior is used
//  real<lower=0> sigma;
}

transformed parameters {
matrix<lower=0,upper=1>[nind, n_occ_minus_1] phi;
matrix<lower=0,upper=1>[nind, n_occ_minus_1] p;
matrix<lower=0,upper=1>[nind, n_occasions] chi;
real mu;

// Constraints
mu = logit(mean_phi);
for (i in 1:nind) {
for (t in 1:(first[i] - 1)) {
phi[i, t] = 0;
p[i, t] = 0;
}
for (t in first[i]:n_occ_minus_1) {
phi[i, t] = inv_logit(mu + epsilon[i]);
p[i, t] = mean_p;
}
}

chi = prob_uncaptured(nind, n_occasions, p, phi);
}

model {
// Priors
// Uniform priors are implicitly defined.
//  mean_phi ~ uniform(0, 1);
//  mean_p ~ uniform(0, 1);
//  sigma ~ uniform(0, 5);
// In case a weaily informative prior is used
//  sigma ~ normal(2.5, 1.25);
epsilon ~ normal(0, sigma);

// Likelihood
for (i in 1:nind) {
if (first[i] > 0) {
for (t in (first[i] + 1):last[i]) {
1 ~ bernoulli(phi[i, t - 1]);
y[i, t] ~ bernoulli(p[i, t - 1]);
}
1 ~ bernoulli(chi[i, last[i]]);
}
}
}

generated quantities {
real<lower=0> sigma2;

sigma2 = square(sigma);
}
``````

I am guessing that the design matrix could replace mu, but just not sure how to do it.

Any help would be appreciated.

I think what it’s saying is that instead of modeling the probabilities as parameters with bounds like:

``````real<lower=0.0, upper=1.0> p;
``````

That you have a regression problem that defines your probabilities using unconstrained variables:

``````real a;
real b;
``````

And a link function (in the model or transformed parameters blocks)

``````vector[N] p = inv_logit(a * x + b);
``````

Assuming x is a vector of covariates.

In this case you have ps and phis at every time point, so if you have T timepoints and N individuals, that’s N * T ps and N * T phis.

So maybe:

``````vector[T] p[N];
for(n in 1:N) {
for(t in 1:T) {
p[n, t] = inv_logit(a * x + b); // or whatever you want inside the inv_logit
}
}
``````

I didn’t look super closely at your code, but it seems like you’re doing this. Be careful here:

``````y[i, t] ~ bernoulli(p[i, t - 1]);
``````

Stan indexes from 1, so make sure that t - 1 is what you want.

You could also just fit a mean at every time point to practice that your inv_logit is working as you expect (so don’t include the covariates and make sure this all matches your previous model).

That help at all?

That is helpful, and that was the direction that I was (slowly!) moving in. I simulated some data with alpha and beta coefficients and am able to recover the parameters, so that is a start. I removed the group survivorship to simplify things. Now I would like to treat the group of the individual as a random effect, but not sure how to integrate back into the code. I guess I should ask this as a separate question?

Here is the code. I am not totally happy with the mixing, but I can work on all that later.

``````library(rstan)

#### parameters for stan ####
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

### MCMC settings
ni <- 5000
nt <- 1
nb <- 1000
nc <- 4

### data generation
N = 500 ## number of marked individuals
n.occasions <- 5                 # Number of capture occasions
marked <- rep(N, n.occasions-1)  # Annual number of newly marked individuals
MarkTot = sum(marked)
v.ind = .25

# covariates
b0 = .2
b1 = 1
X1 = runif(MarkTot,-1,1) ## values for some fake covariate, akin to ndvi
X2 = X1 + rnorm(MarkTot,0,.1)
X3 = X1 + rnorm(MarkTot,0,.1)
X4 = X1 + rnorm(MarkTot,0,.1)
X5 = X1 + rnorm(MarkTot,0,.1)

Y1 = inv.logit(b0 + b1*X1 + rnorm(MarkTot,0,v.ind^.5)) # model with error
Y2 = inv.logit(b0 + b1*X2 + rnorm(MarkTot,0,v.ind^.5)) # model with error
Y3 = inv.logit(b0 + b1*X3 + rnorm(MarkTot,0,v.ind^.5)) # model with error
Y4 = inv.logit(b0 + b1*X4 + rnorm(MarkTot,0,v.ind^.5)) # model with error
Y5 = inv.logit(b0 + b1*X5 + rnorm(MarkTot,0,v.ind^.5)) # model with error

X = data.frame(X2,X3,X4,X5)
PHI = data.frame(Y2,Y3,Y4,Y5)
P = matrix(p,ncol=n.occasions - 1,nrow = MarkTot)

# Simulate capture-histories
CH <- simul.cjs(PHI, P, marked)

# Create vector with occasion of marking
get.first <- function(x) min(which(x!=0))
f <- apply(CH, 1, get.first)

stan_data = list(y=CH,n_occasions = ncol(CH),nind = nrow(CH),
X = X, n_occ_minus_1 = ncol(CH) - 1)

### script
script = 'functions {
int first_capture(int[] y_i) {
for (k in 1:size(y_i))
if (y_i[k])
return k;
return 0;
}

int last_capture(int[] y_i) {
for (k_rev in 0:(size(y_i) - 1)) {
int k = size(y_i) - k_rev;
if (y_i[k])
return k;
}
return 0;
}

matrix prob_uncaptured(int nind, int n_occasions,
matrix p, matrix phi) {
matrix[nind, n_occasions] chi;

for (i in 1:nind) {
chi[i, n_occasions] = 1.0;
for (t in 1:(n_occasions - 1)) {
// Compoud declaration was enabled in Stan 2.13
int t_curr = n_occasions - t;
int t_next = t_curr + 1;

chi[i, t_curr] = (1 - phi[i, t_curr])
+ phi[i, t_curr] * (1 - p[i, t_next - 1]) * chi[i, t_next];
}
}
return chi;
}
}

data {
int<lower=0> nind;							// Number of individuals
int<lower=2> n_occasions;		 			// Number of capture occasions
int<lower=0,upper=1> y[nind, n_occasions];	// Capture-history
int n_occ_minus_1; 							//number of capture occasions minus 1
row_vector[n_occ_minus_1] X[nind];			// covariate
}

transformed data {
int<lower=0,upper=n_occasions> first[nind];
int<lower=0,upper=n_occasions> last[nind];

for (i in 1:nind)
first[i] = first_capture(y[i]);
for (i in 1:nind)
last[i] = last_capture(y[i]);
}

parameters {
real<lower=0,upper=1> mean_phi;		// Mean survival
real<lower=0,upper=1> mean_p;			// Mean recapture
real alpha;
real beta;
}

transformed parameters {
matrix<lower=0,upper=1>[nind, n_occ_minus_1] phi;
matrix<lower=0,upper=1>[nind, n_occ_minus_1] p;
matrix<lower=0,upper=1>[nind, n_occasions] chi;
real mu;

// Constraints
for (i in 1:nind) {
for (t in 1:(first[i] - 1)) {
phi[i, t] = 0;
p[i, t] = 0;
}
for (t in first[i]:n_occ_minus_1) {
mu = inv_logit(alpha + beta*X[i,t]);
phi[i, t] = mu;
p[i, t] = mean_p;
}
}

chi = prob_uncaptured(nind, n_occasions, p, phi);
}

model {
// Priors
alpha ~ normal(0,10);
beta ~ normal(0, 10);

// Likelihood
for (i in 1:nind) {
if (first[i] > 0) {
for (t in (first[i] + 1):last[i]) {
1 ~ bernoulli(phi[i, t - 1]);
y[i, t] ~ bernoulli(p[i, t - 1]);
}
1 ~ bernoulli(chi[i, last[i]]);
}
}
}

generated quantities {

}
'

## Initial values
inits <- function() list(mean_phi = runif(1, 0, 1),
mean_p = runif(1, 0, 1))

## Parameters monitored
params <- c("mean_phi", "mean_p")

## Call Stan from R
model.fit  <- stan(model_code = script,
data = stan_data, init = inits, pars = params,
chains = nc, iter = ni, warmup = nb, thin = nt)

stan_trace(model.fit,pars = c("alpha","beta","mean_phi","mean_p"))

## Summarize posteriors
print(model.fit,pars=c("alpha","beta","mean_phi","mean_p"))``````