Hi all,
I am hoping to get some ideas from folks - I am working on a multistate model to estimate bird mortality rates in different behavioral states (Following this paper: https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/1365-2656.13902). I was able to fit a simple model with my data, but when I make the transition rates time-dependent, it does not initialize. After some troubleshooting with colleagues we got it to work on some simulated data, but not the real data.
There are no differences in general structure between the simulated and real data and we have thoroughly checked to make sure the real data set is correct. In the full data sets, the distributions of each parameter are fairly similar between the simulated and real data.
I am wondering if anyone might have suggestions about what could be going wrong. Thinking maybe some of the parameters need to be constrained differently, or maybe there is a different way to think about incorporating time into the model. The Stan model and R code to simulate data and run the model are below, and I can also provide a subset of the real data. Thank you for your thoughts!
Stan Model:
data{
int<lower = 0> N; //Number of individuals
int<lower = 0> P; // number of time intervals
int<lower = 1> l; // length of each time interval
real<lower = 0> c_p[P]; //midpoint of each time interval
int<lower = 0> max_U; // Maximum number of detections
int<lower = 0> max_U_plus_one; // max number of detections including last detection to TT
//int<lower = 0> TT; // length of study
int<lower = 0> U[N]; // number of detections for each individual
int<lower = 0> d[N, max_U]; // interval within which each obs occurred
int<lower = 0> state[N, max_U];// states
real<lower = 0> delta[N, max_U_plus_one]; // time difference between each observation and the previous one
real<lower = 0> delta_int[N, max_U]; //number of whole intervals between each observation
real<lower = 0> b_r[N, max_U]; // time between detection and start of the interval it's in
real<lower = 0> b_g[N, max_U]; // time between each detection and end of the interval it's in
matrix<lower = 0>[1,4] f;//initial states
matrix<lower = 0>[4,1] ones;// a 1 for each state
}
parameters{
real<lower = 0> h[3];// hazard rate
real<lower = 0> k[3]; //shape parameter
real<lower = 0> g[3]; // scale parameter
real<lower = 0> mu[3];// detection rate
}
transformed parameters{
matrix[4,4] Q[P]; //transition rate matrix
vector<lower = 0>[4] lambda = rep_vector(0,4); // detection rate vector
real<lower = 0> eta12[P]; // Transition rate from state 1 to state 2
real<lower = 0> eta13[P]; // Transition rate from state 1 to state 3
real<lower = 0> eta23[P]; // Transition rate from state 2 to state 3
lambda[1]=mu[1]; // detection rates
lambda[2]=mu[2];
lambda[3]=mu[3];
for (p in 1:P){
eta12[p] = k[1]* c_p[p] ^ (k[1]-1)/g[1]^k[1];
eta13[p] = k[2]* c_p[p] ^ (k[2]-1)/g[2]^k[2];
eta23[p] = k[3]* c_p[p] ^ (k[3]-1)/g[3]^k[3];
// transition rate matrix - rates not listed are 0
Q[p] = rep_matrix(0,4,4);
Q[p,1,1] = -(eta12[p]+eta13[p]+h[1]); //stay in 1
Q[p,1,2] = eta12[p];// 1 to 2
Q[p,1,3] = eta13[p]; //1 to 3
Q[p,1,4] = h[1]; // hazard 1
Q[p,2,2] = -(eta23[p]+h[2]);
Q[p,2,3] = eta23[p]; //2 to 3
Q[p,2,4] = h[2]; // hazard 2
Q[p,3,3] = -(h[3]); // stay in 3
Q[p,3,4] = h[3]; // hazard 3
}
}
model{
matrix[1,4] acc; // temp variable, holds intermediate values for calculating likelihood
matrix[4,4] Gamma; // temp variable, holds likelihood of each detection
matrix[4,4] Gamma_temp;
//priors (following Rushing 2023 appendix 3)
h ~ gamma(1,4); // hazard rate
k ~ gamma(1,1); // shape parameter
g ~ gamma(1,1); // scale parameter
mu ~ gamma(1,4); // detection rates
//likelihood
for (i in 1:N){
matrix[4,4] Omega[U[i]+1];
if (U[i]>0){ // this section is the equation Gamma[i,u] = exp{(Q-Lambda)* Deltat[i,u]} from paper
for (j in 1:U[i]){
if(j==1){
if(delta_int[i,1]==0){ // if detections are in the same interval
Gamma = diag_post_multiply(matrix_exp((Q[1]-diag_matrix(lambda))*delta[i,1]), lambda);
Omega[1] = rep_matrix(0,4,4);
Omega[1,1,state[i,1]] = Gamma[1, state[i,1]];
}else{
Gamma_temp = matrix_exp((Q[1] - diag_matrix(lambda))*l);
if (delta_int[i,1]>1){
for (p in 1:(d[i,1]-1)){
Gamma_temp*= matrix_exp((Q[p]- diag_matrix(lambda))*l);
} //p
} //if
Gamma_temp *= matrix_exp((Q[d[i,1]]- diag_matrix(lambda)) * b_r[i,1]);
Gamma = diag_post_multiply(Gamma_temp, lambda);
Omega[1] = rep_matrix(0,4,4);
Omega[1,1,state[i,1]] = Gamma[1, state[i,1]];
}
}else{
if(delta_int[i,j]==0){
//if detections are within the same interval
Gamma=diag_post_multiply(matrix_exp((Q[d[i, j]] - diag_matrix(lambda))* delta[i,j]), lambda);
Omega[j]=rep_matrix(0,4,4);
Omega[j, state[i, j-1], state[i,j]] = Gamma[state[i,j-1], state[i,j]];
}else{
Gamma_temp = matrix_exp((Q[d[i, j-1]] - diag_matrix(lambda))* b_g[i, j-1]);
if(delta_int[i,j]>1){
for (p in (d[i, j-1]+1):d[i,j]-1){
Gamma_temp *= matrix_exp((Q[p]-diag_matrix(lambda))*l);
}//p
}//if
Gamma = diag_post_multiply(Gamma_temp, lambda);
Omega[j] = rep_matrix(0,4,4);
Omega[j, state[i, j-1], state[i,j]] = Gamma[state[i, j-1], state[i, j]];
} //else
}//else
} //j
// last detection to end of study
Gamma = matrix_exp((Q[d[i, U[i]]]- diag_matrix(lambda))* b_g[i, U[i]]);
if(d[i, U[i]]<P){
for (p in (d[i, U[i]]+1):P){
Gamma *= matrix_exp((Q[p]-diag_matrix(lambda))*l);
}
Omega[U[i]+1] = rep_matrix(0,4,4);
Omega[U[i]+1, state[i, U[i]], state[i, U[i]]:4] = Gamma[state[i, U[i]], state[i, U[i]]:4];
// likelihood for individual i
acc = f * Omega[1];
for (j in 1: U[i]+1){
acc *= Omega[j];
} //j
target += log(acc*ones);
}
}//i
}//N
}//end
R Code to simulate data and run Stan Model
library(rstan)
library(expm)
# Runs with simulated data, but not real data
#Errors when running with real data:
#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.
#Chain 1:
#Chain 1: Initialization between (-2, 2) failed after 100 attempts.
#Chain 1: Try specifying initial values, reducing ranges of constrained values, or #reparameterizing the model.
#[1] "Error : Initialization failed."
# data simulation (adapted from Rushing 2023)
set.seed(4457)
## Set parameters
N <- 20 # Number of individuals
T <- 100 # Length of study
h <- c(0.0025, 0.01, 0.001) # Hazard rates
gamma <- c(6, 4, 8) # Weibull scale parameter
kappa <- c(8, 10, 8) # Weibull shape parameter
lambda <- c(0.2, 0.1, 0.3, 0) # Detection rates
b <- seq(from = 0, to = T, by = 5)
# Interval mid-points, assuming 5-day intervals
c_p <- seq(from = 2.5, to = 97.5, by = 5)/10
# Number of intervals
P <- length(c_p)
# Indicator of what interval each day is in
tau <- rep(1:length(c_p), each = 5)
# Transition rates
eta12 <- kappa[1] * c_p ^ (kappa[1] - 1)/gamma[1] ^ kappa[1]
eta13 <- kappa[2] * c_p ^ (kappa[2] - 1)/gamma[2] ^ kappa[2]
eta23 <- kappa[3] * c_p ^ (kappa[3] - 1)/gamma[3] ^ kappa[3]
## Define transition rate and detection matrices
Q <- array(0, dim = c(T, 4, 4))
for(t in 1:T){
Q[t, 1, 1] <- -(eta12[tau[t]] + eta13[tau[t]] + h[1])
Q[t, 1, 2] <- eta12[tau[t]]
Q[t, 1, 3] <- eta13[tau[t]]
Q[t, 1, 4] <- h[1]
Q[t, 2, 2] <- -(eta23[tau[t]] + h[2])
Q[t, 2, 3] <- eta23[tau[t]]
Q[t, 2, 4] <- h[2]
Q[t, 3, 3] <- -(h[3])
Q[t, 3, 4] <- h[3]
}
## Simulate true states
s <- matrix(NA, nrow = N, ncol = T)
s[,1] <- 1
for(i in 1:N){ for(t in 2:T){
s[i, t] <- which(rmultinom(1, 1, prob = expm(Q[t,,])[s[i, t - 1],]) == 1)
}
}
# Change state 3 to 0 (dead)
s[s == 4] <- 0
## Simulate encounter histories as a list because
## history length differs among individuals
det_list <- state_list <- vector(mode = "list", length = N)
U1 <- U2 <- U3 <- vector(length = N)
# How many days was each individual in each state?
s1 <- apply(s, 1, function(x) sum(x == 1))
s2 <- apply(s, 1, function(x) sum(x == 2))
s3 <- apply(s, 1, function(x) sum(x == 3))
for(i in 1:N){
U1[i] <- rpois(1, lambda[1] * s1[i]) # Number of detections
dets1 <- sort(runif(U1[i], 0, s1[i])) # Time of detections
state1 <- rep(1, U1[i])
U2[i] <- rpois(1, lambda[2] * s2[i]) # Number of detections
dets2 <- sort(runif(U2[i], 0, s2[i])) + s1[i] # Time of detections
state2 <- rep(2, U2[i])
U3[i] <- rpois(1, lambda[3] * s3[i]) # Number of detections
dets3 <- sort(runif(U3[i], 0, s3[i])) + s1[i] + s2[i] # Time of detections
state3 <- rep(3, U3[i])
det_list[[i]] <- c(dets1, dets2, dets3)
state_list[[i]] <- c(state1, state2, state3)
}
# Number of detections for each individual
U <- U1 + U2 + U3
# Convert detections to matrix
det <- state <- matrix(0, nrow = N, ncol = max(U))
for(i in 1:N){ if(U[i] > 0){
det[i, 1:U[i]] <- det_list[[i]]
state[i, 1:U[i]] <- state_list[[i]]
}
}
# Calculate time difference between detections (including last detection to end of study)
delta <- matrix(0, nrow = N, ncol = max(U) + 1)
delta[,1] <- det[,1]
for(i in 1:N){ if(U[i] > 1){
for(j in 2:U[i]){
delta[i, j] <- det[i, j] - det[i, j - 1]
} }
delta[i, U[i] + 1] <- T - max(det[i,])
}
# Number of occasions between detections; upper/lower intervals
d <- delta_int <- b_g <- b_r <- matrix(0, nrow = N, ncol = max(U))
for(i in 1:N){
if(U[i] > 0){
d[i, 1:U[i]] <- as.numeric(cut(det[i, 1:U[i]], breaks = b))
delta_int[i, 1] <- d[i, 1] - 1
b_r[i, 1] <- det[i, 1] - b[d[i, 1]]
b_g[i, 1] <- b[d[i, 1] + 1] - det[i, 1]
if(U[i] > 1){
for(j in 2:U[i]){
delta_int[i, j] <- d[i, j] - d[i, j - 1]
b_r[i, j] <- det[i, j] - b[d[i, j]]
b_g[i, j] <- b[d[i, j] + 1] - det[i, j]
}
}
}
}
stan_data <- list(N = N,
max_U = max(U),
max_U_plus_one = max(U) + 1,
P = length(c_p),
l = 5,
U = U,
delta = delta,
d = d,
state = state,
delta_int = delta_int,
b_r = b_r,
b_g = b_g,
c_p = c_p,
f = matrix(c(1,0,0,0),1,4),
ones = matrix(c(1,1,1,1), 4, 1))
#stan_data <- readRDS('./stan_data_sub_REAL_051424.rds')
#stan_data <- readRDS('./stan_data_sub_SIM_051424.rds')
# Initial values
initf1 <- function() {
list(h = runif(3, 0, 0.015),
g = runif(3, 5, 15),
k = runif(3, 7, 10),
mu = runif(3, 0, 0.25))
}
# run model - few iterations/1 chain for troubleshooting
fit <- stan(file = "./model.stan",
pars = c("h", "g", "k", "mu"),
data = stan_data,
init = initf1,
warmup = 50,
iter = 100,
chains = 1,
cores = 1,
thin = 1)