Hi! I am trying to further my understanding of measurement models and trying to write them in raw stan. In this model, there is dependence on V from the outcome and all other variables, I, P, and Y. Currently, I have ‘I’ modeled with hierarchical priors, and ‘P’ is just linked through a linear model. I had also tried hierarchical priors on ‘P’. Although I would welcome feedback on any improvements in the code, the particular problem is with that beta5 and beta2 are not identifiable (no updating from prior). The data are simulated so I know updating should occur. They are coefficients on the same data, so it makes sense to me that they would be difficult to distinguish.
How can I deal with this issue?
Also, the model returns no warnings.
Thank you!
Here is the model
"
data {
int<lower=0> N; // cases
vector[N] I_p; // I observed
vector[N] I_e; // I error
vector[N] V; // measured without error
vector[N] Y; // outcome
vector[N] P_p; // P observed
vector[N] P_e; // P error
}
parameters {
vector[N] I_true;
real <lower=0> sigma_y;
real <lower=0> sigma_I;
real beta1;
real beta2;
real beta3;
real beta4;
real beta5;
real beta6;
real alpha1;
real alpha2;
real alpha3;
}
transformed parameters {
vector[N] P_true;
vector[N] I_t;
P_true = alpha1 + beta1*V; // measurement model for P
I_t = alpha2 + beta2*V + beta3*P_true + I_true*sigma_I; // measurement model for I
}
model {
beta1 ~ normal(0,.3);
alpha1 ~ normal(0, .2);
sigma_I ~ normal(0,.2);
alpha2 ~ normal(0,.3);
beta2 ~ normal(0, .3);
beta3 ~ normal(0, .5);
P_p ~ normal(P_true, P_e);
I_true[N] ~ std_normal();
I_p ~ normal(I_true, I_e);
sigma_y ~ exponential(1);
beta4 ~ normal(0,.5);
beta5 ~ normal(0,.4);
beta6 ~ normal(0,.4);
alpha3 ~ normal(0,.5);
Y ~ normal(alpha3 + beta4*I_t + beta5*V + beta6*P_true, sigma_y);
}
generated quantities {
vector[N] y_pred;
for(i in 1:N) {
y_pred[i] = alpha3 + beta4 * I_t[i] + beta5 * V[i] + beta6 * P_true[i];
}
}
"
I have attached a pairs plot that shows several unhappy pairs