Hi,
Is there a way to optionally include both data and parameters in a stan model? I have a model that I would like to allow random effects to be optionally included at a certain level. So I would need a way to optionally supply e.g. the priors information for the covariance matrix for the random effects and their design matrix, and also optionally define the random effects parameters. I’ve tried a few approaches, but each has issues so far:
Option 1:
If I wrap everything that relates to the optional part (Data, parameters etc.) in if statements checking a boolean (that I’ve supplied as part of the data) that indicates whether or not the data and parameters will be present, the model gives error message
PARSER EXPECTED: <one of the following:
a variable declaration, beginning with type,
(int, real, vector, row_vector, matrix, unit_vector,
simplex, ordered, positive_ordered,
corr_matrix, cov_matrix,
cholesky_corr, cholesky_cov
or '}' to close variable declarations>
Which prints out at the line relating to the first occurence of the boolean in the data section:
data{
...
if(boo_b3_l){
matrix[M,rl] Z3_l;
vector[rl] b3_ones;
real<lower=0> tau_Al_sig;
real<lower=0> tau_Al_nu;
real<lower=1> Omega_Al_eta;
}
...
}
Option 2:
Alternatively, if I try to supply the data linking to these optional random effects, but set the number of these random effects rl to zero, e.g. define and supply the design matrix Z3_l as an M by 0 matrix, and only use the boolean when the effects are being sampled, I get error message:
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: std::bad_alloc
[1] "Error in sampler$call_sampler(args_list[[i]]) : std::bad_alloc"
[1] "error occurred during calling the sampler; sampling not done"
Overall
- is it possible to include optional data and optional parameters?
- what is a reliable appropriate way to do so?
Note: I know this particular model can probably already be fitted using existing code, however this is a learning example for me to build up to coding a complex model that cannot be fitted using current code.
Full code:
functions{
}
data {
//dimensions of data
int<lower=1> M; //number of observations
int<lower=1> N; //number of individuals
int<lower=1> K; //number of individuals
int<lower=1> pl; //number long fixed
int<lower=1> ql; //number long indrand
int<lower=0> rl; //number long studrand
//id variables
int<lower=1> idnum[M];
int<lower=1> studnum[M];
//longitudinal outcome information
vector[M] y1;
//longitudinal fixed information
matrix[M,pl] X_l;
real<lower=0> beta_l_sig;
real<lower=0> sigma_e_sig;
real<lower=0> sigma_e_nu;
//booleans
int<lower = 0, upper = 1> boo_b3_l;
//longitudinal indrand information
matrix[M,ql] Z2_l;
vector[ql] b2_ones;
real<lower=0> tau_Dl_sig;
real<lower=0> tau_Dl_nu;
real<lower=1> Omega_Dl_eta;
//longitudinal studrand information
if(boo_b3_l){
matrix[M,rl] Z3_l;
vector[rl] b3_ones;
real<lower=0> tau_Al_sig;
real<lower=0> tau_Al_nu;
real<lower=1> Omega_Al_eta;
}
}
transformed data{
vector[ql] mu_b2;
if(boo_b3_l){
vector[rl] mu_b3;
}
mu_b2 = rep_vector(0,ql);
if(boo_b3_l){
mu_b3 = rep_vector(0,rl);
}
}
parameters{
vector[pl] beta_l;
real<lower=0> sigma_e;
matrix[N,ql] b2;
corr_matrix[ql] Omega_Dl;
vector<lower=0>[ql] tau_Dl;
if(boo_b3_l){
matrix[K,rl] b3;
corr_matrix[rl] Omega_Al;
vector<lower=0>[rl] tau_Al;
}
}
transformed parameters{
vector[M] linpred_y;
matrix[ql,ql] Dl;
if(boo_b3_l){
matrix[rl,rl] Al;
}
linpred_y = X_l*beta_l + (Z2_l .* to_matrix(b2[idnum,]))*b2_ones;
Dl = quad_form_diag(Omega_Dl, tau_Dl);
if(boo_b3_l){
linpred_y=linpred_y + (Z3_l .* to_matrix(b3[studnum,]))*b3_ones;
Al = quad_form_diag(Omega_Al, tau_Al);
}
}
model{
//Priors
tau_Dl ~ student_t(tau_Dl_nu,0,tau_Dl_sig);
Omega_Dl~ lkj_corr(Omega_Dl_eta);
if(boo_b3_l){
tau_Al ~ student_t(tau_Al_nu,0,tau_Al_sig);
Omega_Al~ lkj_corr(Omega_Al_eta);
}
beta_l ~ normal(0,beta_l_sig);
sigma_e ~ student_t(sigma_e_nu,0,sigma_e_sig);
//Likelihood
y1 ~ normal(linpred_y,sigma_e);
for(i in 1:N){
b2[i,] ~ multi_normal(mu_b2, Dl);
}
if(boo_b3_l){
for(i in 1:K){
b3[i,] ~ multi_normal(mu_b3, Al);
}
}
}
Any comments / suggestions would be greatly appreciated - thanks in advance.