Hello, I am running a hierarchical model that transforms both individual- and group-level parameters, and it is supposed to return (1) an array of parameters fitted on the individual level, and (2) a mean and standard deviation term for each parameter calculated from the group level parameters. However, currently the model only returns the fitted parameters from the individual level. I am wondering why I cannot access the mean and standard deviation variables (ending in _mu and _sd) from the stanfit object, even though they are declared in the parameters {} block?
Here is my model:
data {
int ntr; // number of total data points (N)
int nsub; // total number of participants
int ParticipantID[ntr]; // which subject each data point belongs to
int nResponseTypes; // K
real var1[ntr];
real var2[ntr];
real proximity[ntr]; //used to for both var3 and var4
int actionTaken[ntr]; //y
}
parameters {
// Group level
real<lower=0> invTemp_mu;
real<lower=0> invTemp_sd;
real<lower=0> b_var1_mu;
real<lower=0> b_var1_sd;
real<lower=0> b_var2_mu;
real<lower=0> b_var2_sd;
real<lower=0> b_var3_sd;
real<lower=0> b_var3_mu;
// Individual participants
real<lower=0> invTemp_pr[nsub];
real b_var1_pr[nsub];
real b_var2_pr[nsub];
real b_var3_pr[nsub];
}
transformed parameters {
// Group level transformed variables
real b_var1_pr2[nsub];
real b_var2_pr2[nsub];
real b_var3_pr2[nsub];
real<lower=0> k_mu;
real<lower=0> b_var1_mu_divk;
real<lower=0> b_var2_mu_divk;
real<lower=0> b_var3_mu_divk;
real<lower=0> b_var4_mu_divk;
// Individual level transformed variables
real k[nsub]; // sum of all absolute regression weights
real b_var1[nsub];
real b_var2[nsub];
real b_var3[nsub];
real b_var4[nsub];
real invTemp[nsub];
// First transform the priors to the group level
for (is in 1:nsub) {
invTemp[is] = exp(invTemp_pr[is] * invTemp_sd + log(invTemp_mu)); // what it does invTemp[is] ~normal(invTemp_mu,invTemp_sd)
b_var1_pr2[is] = b_var1_pr[is] * b_var1_sd + b_var1_mu;
b_var2_pr2[is] = b_var2_pr[is] * b_var2_sd + b_var2_mu;
b_var3_pr2[is] = b_var3_pr[is] * b_var3_sd + b_var3_mu;
}
// Then transform to the individual level using the group level transformations
for (is in 1:nsub) {
k[is] = fabs(1) + fabs(b_var1_pr2[is])+fabs(b_var2_pr2[is])+fabs(b_var3_pr2[is]);
b_var4[is] = 1/k[is];
b_var1[is] = b_var1_pr2[is]/k[is];
b_var2[is] = b_var2_pr2[is]/k[is];
b_var3[is] = b_var3_pr2[is]/k[is];
}
k_mu = fabs(1) + fabs(b_var1_mu) + fabs(b_var2_mu) + fabs(b_var3_mu);
b_var4_mu_divk = 1/k_mu;
b_var1_mu_divk = b_var1_mu/k_mu;
b_var2_mu_divk = b_var2_mu/k_mu;
b_var3_mu_divk = b_var3_mu/k_mu;
}
model{
matrix[ntr,nResponseTypes] invTxUtil;
// Priors - these become the priors for the group level
invTemp_mu ~ cauchy(10,2);
b_var1_mu ~ cauchy(0,1);
b_var2_mu ~ cauchy(0,1);
b_var3_mu ~ cauchy(0,1);
invTemp_sd ~ cauchy(0,3);
b_var1_sd ~ cauchy(0,1);
b_var2_sd ~ cauchy(0,1);
b_var3_sd ~ cauchy(0,1);
// individual subjects
invTemp_pr ~ normal(0,1);
b_var1_pr ~ normal(0,1);
b_var2_pr ~ normal(0,1);
b_var3_pr ~ normal(0,1);
// Iterate through each trial
for (itr in 1:ntr) {
// Utility for action1
invTxUtil[itr,1] = invTemp[ParticipantID[itr]]*(b_var1[ParticipantID[itr]]*var1[itr]);
// Utility for action2
invTxUtil[itr,2] = invTemp[ParticipantID[itr]]*(b_var4[ParticipantID[itr]]*proximity[itr]);
// Utility for action3
invTxUtil[itr,3] = invTemp[ParticipantID[itr]]*(b_var2[ParticipantID[itr]]*var2[itr] + b_var3[ParticipantID[itr]]*proximity[itr]);
actionTaken[itr] ~ categorical_logit(to_vector(invTxUtil[itr])); // link utilities to actual choices
}
}
I appreciate any help, thank you.