Looking for some advice based the n_eff and the pairs plot for mixed effects models

General
1

Hello everybody,

I am fitting a multiple outcome mixed effects model where the correlation between the outcomes are going to capture by the random effects (random intercept).

The random intercept follows N~(0,sigma_u)
I am having a hard time of selecting an appropriate prior for the sigma_u.
I tried with half-cauchy , inv-gamma priors but so far getting the same warning as follows :

There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low 
3: Examine the pairs() plot to diagnose sampling problems

My summary results are as follows :
This is based on 4 chains ,10000 iterations on each and cauchy(0,4) prior on sigma_u

image
image811×380 30.6 KB

It can be observed that the results are reasonably okay with all parameters except with the sigma_u. For sigma_u the n_eff is very low and the Rhat is close to 1.01.
The pairs plot looks like this for chain1. And it is pretty much same for other chains also.

image
image1216×444 73.6 KB

You can see that the energy parameter is highly correlated with the sigma_u.

Can you suggest anything to improve the results ?
Also based on low n_eff , will this can be improved if I increase the number of iterations.

I am new to Bayesian and Stan. So any advice will be highly appreciate.
Thank you.

##############Update######
My model with following Specifications: two outcomes, separate predictors for each outcome and common random intercepts which will share by each subject across different outcomes

       data {
  int<lower=1> N;
   int<lower=1> K1; 
   int<lower=1> K2; 
  int<lower=0,upper=1> y1[N];
  int<lower=0,upper=1> y2[N];
    
    row_vector[K1] x1[N];
  row_vector[K2] x2[N];
}


parameters {
   real alpha1;
   real alpha2;
    vector[K1] beta1;
    vector[K2] beta2;
    real<lower=0> sigma_u;
    real u[N];
}

model {
  
  u ~ normal(0, sigma_u);
  beta1 ~ normal(0, 100);
  beta2 ~ normal(0, 100);
  alpha1 ~ normal(0, 100);
  alpha2 ~ normal(0, 100);
  sigma_u ~ cauchy(0, 4);
  
  for(i in 1:N){
    
    y1[i] ~ bernoulli(inv_logit(alpha1 + u[i] + x1[i]*beta1 ));
   y2[i] ~ bernoulli(inv_logit(alpha2 + u[i] + x2[i]*beta2 ));
  
  }
}
2

Welcome!

Over here Reparameterizing to avoid low E-BFMI warning this might give you some pointers.

2 Likes
3

@Ara_Winter Thank you for your suggestion. I went through that but it didn’t resolve my issue. I updated the question with the inclusion of the code. Please have a look if possible. Thanks

4

Looking at your pairs plot, there are some strong posterior correlations. alpha2 and beta 2.1 are correlated at -.97. alpha1 and beta 1.1 and beta1.2 are correlated at -.68 and -.78 respectively. This points towards poor identifiability, I believe. Some questions: How much variance is there in your predictors? Are you sure your normal(0,100) priors are reasonable?

3 Likes
5

As Erling said, your priors are likely causing a lot of problems here. With the function inv_logit(x), as soon as x gets larger than ~19, the function will only return 1. So you should look at making those priors significantly smaller.

Additionally, you should try the non-centered parameterisation for u (more information here in the manual), which is as simple as changing the declaration to:

vector<multiplier=sigma_u>[N] u;

Finally, if you change the declaration of x1 and x2 to matrices, you can use the bernoulli_logit_glm function for more efficient and stable sampling.

Putting all of these in place, your model would look like:

data {
  int<lower=1> N;
  int<lower=1> K1; 
  int<lower=1> K2; 
  int<lower=0,upper=1> y1[N];
  int<lower=0,upper=1> y2[N];
  matrix[N,K1] x1;
  matrix[N,K2] x2;
}


parameters {
  real alpha1;
  real alpha2;
  vector[K1] beta1;
  vector[K2] beta2;
  real<lower=0> sigma_u;
  vector<multiplier=sigma_u>[N] u;
}

model {
  u ~ normal(0, sigma_u);
  beta1 ~ normal(0, 5);
  beta2 ~ normal(0, 5);
  alpha1 ~ normal(0, 5);
  alpha2 ~ normal(0, 5);
  sigma_u ~ cauchy(0, 4);

  y1 ~ bernoulli_logit_glm(x1, alpha1 + u, beta1);
  y2 ~ bernoulli_logit_glm(x2, alpha2 + u, beta2);
}

I’ve used normal(0,5) priors here, but you will want to check that those priors are realistic/accurate for your data.

4 Likes
6

@andrjohns Thank you very much for your advice. The results improved dramatically with your suggestions. Now I am not getting the low BFMI warning.