Hi!
I have repeated behavioural data on multiple indivdiuals which are nested in different groups (colonies). My goal is to test among-individual and within-individual (residual) variation in behaviour and also control for colony effect.
I’m trying to build a double hierarchical mixed model, with the following dataset:
data.list =
list(y = (xdata$totdist-mean(xdata$totdist)) / sd(xdata$totdist),
ind = as.integer(xdata$ind),
colony = as.integer(xdata$colid),
headsize = (xdata$head-mean(xdata$head)) / sd(xdata$head),
colsize = (xdata$size-mean(xdata$size)) / sd(xdata$size),
myrnest = (xdata$myrnest-mean(xdata$myrnest)) / sd(xdata$myrnest),
totnest = (xdata$totnest-mean(xdata$totnest)) / sd(xdata$totnest),
habitat = as.numeric(xdata$habitat),
N = length(xdata$totdist),
N_ind = length(unique(xdata$ind)),
N_col = length(unique(xdata$colid)))
with the following stan code:
data{
int<lower=1> N;
int<lower=1> N_ind;
int<lower=1> N_col;
real y[N];
real headsize[N];
real colsize[N];
real myrnest[N];
real totnest[N];
int habitat[N];
int colony[N];
int ind[N];
}
parameters{
vector[2] re_id[N_ind];
vector[N_col] re_col;
vector[6] b;
vector[6] g;
vector<lower=0,upper=10>[2] sd_id;
real<lower=0,upper=10> sd_col;
corr_matrix[2] id_cor;
}
model{
vector[N] sdy;
vector[N] muy;
id_cor ~ lkj_corr(2);
// Priors for random effect variances commented out as set in parameters.
//sd_id ~ uniform(0, 10); //prior implied in ‘parameters’ block
//sd_col ~ uniform(0, 10);
b ~ normal(0,100); // prior for mean fixed effects
g ~ normal(0,100); // prior for residual fixed effects
//setup random effect term
re_id ~ multi_normal( rep_vector(0,2), quad_form_diag(id_cor,sd_id));
re_col ~ normal(0, sd_col);
///////////// Setup of the likelihood model ///////////////////
for ( i in 1:N ) {
// Mean level model
muy[i] = (b[1] + re_id[ind[i],1] + re_col[colony[i]]) + b[2]*headsize[i] +
b[3]*colsize[i]+ b[4]*myrnest[i] + b[5]*totnest[i] +
b[6]*habitat[i];
// Residual level model
// Log-link moved to right-hand side of equation as the exponential
sdy[i] = exp(g[1] + g[2]*headsize[i] + g[3]*colsize[i] + g[4]*myrnest[i] +
b[5]*totnest[i] + b[6]*habitat[i] + re_id[ind[i],2]);
}
y ~ normal( muy , sdy );
}
Although after running the model, I get the following Error message
here are whatever error messages were returned
[[1]]
Stan model ‘activity’ does not contain samples.
[[2]]
Stan model ‘activity’ does not contain samples.
[[3]]
Stan model ‘activity’ does not contain samples.
Warning message:
In .local(object, …) :
some chains had errors; consider specifying chains = 1 to debug
Could you help me to figure out what the problem could be? I’m a newby with Stan and at this point I have no idea what’s wrong.
Thanks in advance.