Capture-Recapture with further impact on survivability

Dear Forum,
I have been using Cormack Jolly Seber (CJS) examples in the manual and online to try and estimate posteriors for survivability and prob of recapture, in addition to the influence of a trait which reducing survivability (phi).

p = prob of recapture
phi = prob of survivability
dphi[4] = reduction in survivability based one 4 levels (1=nil, 2=low, 3=med, 4=high)

I created a dataset and for normal p and phi estimations it works using the code in the manual.
Sample of capture recapture data (X)
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20 V21 V22 V23
[1,] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
[3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0
[4,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0
[6,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0

The additional parameter is also included in the generated data as Y. Most animals do NOT exhibit this, but it is present.
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20 V21 V22
[1,] 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1
[2,] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[3,] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3
[4,] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[5,] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
[6,] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

I am getting initialization errors, but have tried to truncate priors to mitigate that as an initial work-around.

Any pearls of wisdom would be greatly appreciated.

Thanks,
Steve

data {
int<lower=1> K;                      // capture events
int<lower=0> I;                      // number of individuals
int<lower=0,upper=1> X[I,K];         // X[i,k]: individual i captured at k (binary matrix)
int<lower=1,upper=4> Y[I,K-1];       // Y[i,k-1]: individual influence on future survivability (1=nil, 2=low, 3=med, 4=high)
}

transformed data {
int<lower=0,upper=K+1> first[I];     // first[i]: ind i first capture
int<lower=0,upper=K+1> last[I];      // last[i]:  ind i last capture

first = rep_array(K+1,I);
last = rep_array(0,I);
for (i in 1:I) {
for (k in 1:K) {
if (X[i,k] == 1) {
if (k < first[i])
first[i] = k;
if (k > last[i])
last[i] = k;
}
}
}

}

parameters {
vector<lower=0,upper=1>[K-1] phi;       // phi[k]: Pr[alive at k + 1 | alive at k]
vector<lower=0,upper=1>[K] p;           // p[k]: Pr[capture at k]
real<lower=0,upper=1> dphi[4];          //upper limit should be < phi... but we I set it as 1?
}

model {
for (i in 1:I) {
if (last[i] > 0) { //implies there are some captures for the animal atleast (put in based on released animals never originally caught)
for (k in (first[i]+1):last[i]) {

//SURVIVABILITY
if (Y[i, k-1] > 1)                                // if individual influencing parameter exists
target += (log(phi[k-1] - dphi[Y[i,k-1]]));     // uses one case, else uses another
else
target += (log(phi[k-1]));

//CAPTUREABILITY
if (X[i,k] == 1)
target += (log(p[k]));       // i captured at k  // increment_log_prob(log(p[k]));
else
target += (log1m(p[k]));     // i not captured at k  // increment_log_prob(log1m(p[k]));
}
}
}

p ~ normal( 0.3 , 0.2 );
phi ~ normal( 0.7 , 0.1 );
dphi[1] ~ uniform( 0.4 , 0.6 ); //This is coded to be redundant above
dphi[2] ~ normal( 0.3 , 0.1 ); // T[0,(min(phi)-0.01)];
dphi[3] ~ normal( 0.6 , 0.1 ); // T[0,(min(phi)-0.01)];
dphi[4] ~ normal( 0.9 , 0.2 ); // T[0,(min(phi)-0.01)];
}


Hi! Can you paste the error message? Thanks for the data and model code.

Hi,

The error results from when trying to initialise.
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.

I suspect it is this:
target += (log(phi[k-1] - dphi[Y[i,k-1]]));

Regards,
Steve

Is it possible that you’re getting negative survival probabilities by subtracting dphi from phi?

This phi[k-1] - dphi[Y[i,k-1]] could give negative values under the current parameterization, if I’m understanding correctly. I wonder whether a different parameterization for your survival effects would help, e.g.,

\phi_{ik} = \text{logit}^{-1}(\alpha + \beta_{y_{i k}}),

where \alpha is an intercept, and \beta is a vector with 4 elements, indexed by y. You could declare alpha and beta in the parameters block, and compute phi as a transformed parameter. The logit transform would ensure that you get survival probabilities between 0 and 1.

1 Like

Also can you post the model call? I want to see if what you settings are and if you specified inits.

Hi,

@mbjoseph Started looking into logit. Seems like a logical way forward, just yet to get it working. Alternatively, perhaps inits can be specified (if I am reading into what @Ara_Winter is saying).

@Ara_Winter
sampling_iterations = 1e5
options(mc.cores = 3)
out = rstan::sampling(
object=model
, data = data
, chains = 3
, iter = sampling_iterations
, warmup = sampling_iterations/2
, thin = 4#2#1
, refresh = sampling_iterations/10
)

Regards,
Steve

Thanks @Stephen_Bradshaw ! Can you try setting your init = 0 in the model call? But I would also take the approach that @mbjoseph outlined.

@Stephen_Bradshaw here’s an example of how you might ensure that the survival probabilities are bounded between 0 and 1 (note the use of log_inv_logit where previously you had log:

data {
int<lower=1> K;                      // capture events
int<lower=0> I;                      // number of individuals
int<lower=0,upper=1> X[I,K];         // X[i,k]: individual i captured at k (binary matrix)
int<lower=1,upper=4> Y[I,K-1];       // Y[i,k-1]: individual influence on future survivability (1=nil, 2=low, 3=med, 4=high)
}

transformed data {
int<lower=0,upper=K+1> first[I];     // first[i]: ind i first capture
int<lower=0,upper=K+1> last[I];      // last[i]:  ind i last capture

first = rep_array(K+1,I);
last = rep_array(0,I);
for (i in 1:I) {
for (k in 1:K) {
if (X[i,k] == 1) {
if (k < first[i])
first[i] = k;
if (k > last[i])
last[i] = k;
}
}
}
}

parameters {
real alpha;
vector[4] beta;
vector<lower=0,upper=1>[K] p;           // p[k]: Pr[capture at k]
}

model {
for (i in 1:I) {
if (last[i] > 0) {
//implies there > 1 captures for the animal
// (put in based on released animals never originally caught)
for (k in (first[i]+1):last[i]) {
// survival
if (Y[i, k-1] > 1)  // if individual influencing parameter exists
target += log_inv_logit(alpha + beta[Y[i, k-1]]);
else
target += log_inv_logit(alpha);

// capture
X[i, k] ~ bernoulli(p[k]);
}
}
}

p ~ beta(1, 1);
alpha ~ normal(0, 1.5);
beta ~ std_normal();
}


1 Like

Thanks a lot for the code snippet! I’ve started implementing it but am yet to get it converge to values which I am expecting (from generated data).