Weak identifiability in measurement models

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

There are a couple of things that you might want to try.

First, it seems like you’re trying to overconstrain your data with really specific priors like normal(0,.2). If your data fits the model, it’ll fit the correct values without needing a very specific prior. You should be able to try out normal(0,1) or even normal(0,10).

Second, your sigma_y is probably not an exponential which means the mode is at zero. You should try a cauchy or gamma or something.

Third, make sure you start by fitting fake data, then move onto fitting your real data one parameter at a time. This is the most important thing for bespoke Bayesian modeling, but I tend to find myself skipping this step and having to start over every time.

Hope this helps a little.

2 Likes

Ok, back to the drawing board, thank you!

I think there was a problem with the simulation, I reworked it, and now I can recover parameters with a centered parameterization, and N(0,1) on all coefficients. However, I get warnings with ESS and energy, particularly with sigma_y. If I specify a prior away from 0 (i.e. N(1,.3) or something) I can avoid divergences however the other warnings remain. I’ve attached a pairs plot with normal priors on sigma_I and sigma_P, and a gamma(2,1) on sigma_y.
Do you have any ideas how to make this more efficient?