# Question about Rhat statistics

**URL:** <https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199>\
**Category:** General\
**Tags:** rstan\
**Created:** [June 25, 2020, 12:29am UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199 "2020-06-25T00:29:02Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gen\_Y](https://avatars.discourse-cdn.com/v4/letter/g/ccd318/32.png) [@Gen\_Y](https://discourse.mc-stan.org/u/Gen_Y)\
**Post date:** [June 25, 2020, 12:29am UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199/1 "2020-06-25T00:29:02Z")

</div>

Hi, I’m new to stan. My stan model code is like:

```stan
data {
  int<lower=1> N; // Number of subjects
  int<lower=1> ms; // Number of measurements
  int subject[ms]; // subject ID
  real Y[ms];
  real X[ms];
}
parameters {
  real alpha;
  real<lower=0> sigma_e;
  real mu;
  real<lower=0> sigma;
  vector u[N];
}
transformed parameters {
  vector[ms] theta;
  // construct the theta
  for (i in 1:ms) {
    theta[i] <- alpha*X[i]+u[subject[i]];
  }
}
model {
  // construct random effect
  u ~ normal(mu, sigma);
  for (i in 1:ms) {
    Y[i] ~ normal(theta[i],sigma_e);
  }
  // construct the priors
  alpha ~ normal(0,10);
  sigma_e ~ inv_gamma(0.01,0.01);
  mu ~ normal(0,10);
  sigma ~ inv_gamma(0.01,0.01);
}

```

Something like linear mixed effects model + measurement error.  
I run for 4 chains and 2,000 iterations. The Rhat of alpha and u are not good, sometimes \>2 (I tried to increase iterations but not helpful). And I checked for all estimated theta, the Rhat is ranging from 1.00 to 1.04, I assume it is acceptable.  
So my question is, if I only need theta for subsequent analysis, is it safe to use theta, even though some other parameters’ Rhat not good?

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [June 25, 2020, 3:06am UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199/2 "2020-06-25T03:06:07Z")

</div>

Welcome to the Stan forums.

> [@Gen\_Y](#):
>
> So my question is, if I only need theta for subsequent analysis, is it safe to use theta, even though some other parameters’ Rhat not good?

It’s not advisable to ignore those high Rhats, but I think we can probably help get those Rhat values down. Are you also getting warnings about divergences when fitting this?

A few tips about your Stan code (the first few aren’t related to Rhat but the ones towards the end might be):

> [@Gen\_Y](#):
>
> ```no-highlight
> real Y[ms];
> real X[ms];
> 
> ```

This is fine but I would change it to

```stan
  vector[ms] Y;
  vector[ms] X;

```

so that you can do vector multiplication later in the code (see below).

> [@Gen\_Y](#):
>
> `vector u[N];`

The syntax for vectors (as opposed to arrays) is `vector[N] u;`. What you have should result in a parser error and it shouldn’t even let to run this model. Are you sure this code runs?

> [@Gen\_Y](#):
>
> ```no-highlight
> vector[ms] theta;
> // construct the theta
> for (i in 1:ms) {
> theta[i] <- alpha*X[i]+u[subject[i]];
> }
> 
> ```

If you change `Y` and `X` to vectors as I recommended above then you can vectorize this and replace the loop with this:

```stan
vector[ms] theta = alpha * X + u[subject];

```

In the model block:

> [@Gen\_Y](#):
>
> ```no-highlight
> for (i in 1:ms) {
> Y[i] ~ normal(theta[i],sigma_e);
> }
> 
> ```

can be vectorized and will be faster as

```stan
Y ~ normal(theta, sigma_e);

```

Also, although these `inv_gamma(0.01,0.01)` used to be popular (maybe still are) we recommend against them (see e.g. section 4.3 of [http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf](http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf)). You could try an exponential distribution or even a normal distribution (Stan will automatically make that a half-normal since you have `<lower=0>`).

> [@Gen\_Y](#):
>
> `u ~ normal(mu, sigma);`

This is known as the _centered_ parameterization of a hierarchical model. This can be fine if your data is highly informative, but that is often not the case and it will help to use the _non-centered_ parameterization (here’s the relevant section of the Stan user’s guide: [Stan User’s Guide](https://mc-stan.org/docs/2_23/stan-users-guide/reparameterization-section.html#hierarchical-models-and-the-non-centered-parameterization)). In your case the non-centered parameterization (which is statistically equivalent but often easier to fit) would look like this (also adding in my recommendations from above):

```stan
data {
  int<lower=1> N; 
  int<lower=1> ms; 
  int subject[ms]; 
  vector[ms] Y;
  vector[ms] X;
}
parameters {
  real alpha;
  real<lower=0> sigma_e;
  real mu;
  real<lower=0> sigma;
  vector[N] u_raw; // this will be N(0,1) and scaled and shifted to N(mu, sigma)
}
transformed parameters {
  vector[ms] u = mu + sigma * u_raw; // construct u by scaling and shifting
  vector[ms] theta = alpha * X + u[subject];
}
model {
  Y ~ normal(theta, sigma_e);
  u_raw ~ normal(0, 1); // implies u ~ normal(mu, sigma)
  alpha ~ normal(0,10);
  mu ~ normal(0,10);

  // these aren't wrong but I recommend changing them (see comment above)
  sigma_e ~ inv_gamma(0.01,0.01); 
  sigma ~ inv_gamma(0.01,0.01);  
}

```

---

<div class="post-metadata">

**Author:** ![Gen\_Y](https://avatars.discourse-cdn.com/v4/letter/g/ccd318/32.png) [@Gen\_Y](https://discourse.mc-stan.org/u/Gen_Y)\
**Post date:** [June 25, 2020, 4:46am UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199/3 "2020-06-25T04:46:44Z")

</div>

Hi, jonah. Thanks so much for your detailed response.  
Your revised code looks more harmonious.

> [@jonah](#):
>
> Are you also getting warnings about divergences when fitting this?

Yes, I also got warnings about BFMI was low, and ESS is too low.

> [@jonah](#):
>
> The syntax for vectors (as opposed to arrays) is `vector[N] u;` . What you have should result in a parser error and it shouldn’t even let to run this model. Are you sure this code runs?

This model could run, with a lot warnings and long computation time.

Your recommended references on inverse gamma and non-centered parameterization are quite helpful.

Acutally my model is a little more complicated than this. If my random effects term u = (u0, u1, u2), should I define u like `matrix[N,3] u_raw`, or like `vector[N] u_raw[3]`? In this case, how should I define theta?  
The formula is like

theta\_i(t) = alpha \* X\_i + u\_{i0} + u\_{i1}\*t + u\_{i2}\*t^2

(there is a `real<lower=0> time[ms]` in data for this model)  
Now my code is like

```no-highlight
for (i in 1:ms) {
  theta[i] <- alpha*X[i]+u[subject[i],1]+u[subject[i],2]*time[i]+u[subject[i],3]*time[i]^2;
}

```

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [June 25, 2020, 1:39pm UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199/4 "2020-06-25T13:39:07Z")

</div>

Assuming each of your `u` vectors has its own `mu` and `sigma` (is that true?), how about this?

```stan
data {
  int<lower=1> N; // Number of subjects
  int<lower=1> ms; // Number of measurements
  int subject[ms]; // subject ID
  vector[ms] Y;
  vector[ms] X;
  vector<lower=0>[ms] time; // declare as a vector so you can do multiplication in transformed params
}
parameters {
  real alpha;
  real<lower=0> sigma_e;
  vector[3] mu;
  vector<lower=0>[3] sigma;
  vector[N] u_raw[3];
}
transformed parameters {
  vector[ms] u[3]; 
  vector[ms] theta;
  
  // each u[j] is a vector, each mu[j] and sigma[j] is a scalar
  for (j in 1:3) {
    u[j] = mu[j] + sigma[j] * u_raw[j];
  }
  
  // using .* for element-wise multiplication
  theta = alpha * X + u[1][subject] + u[2][subject] .* time + u[3][subject] .* square(time);
}
model {
  Y ~ normal(theta, sigma_e); // data model
  for (j in 1:3) u_raw[j] ~ normal(0, 1);   
  
  alpha ~ normal(0,10);
  mu ~ normal(0,10);

  // these aren't wrong but I recommend changing them (see comment above)
  sigma_e ~ inv_gamma(0.01,0.01); 
  sigma ~ inv_gamma(0.01,0.01);  

```

---

<div class="post-metadata">

**Author:** ![Gen\_Y](https://avatars.discourse-cdn.com/v4/letter/g/ccd318/32.png) [@Gen\_Y](https://discourse.mc-stan.org/u/Gen_Y)\
**Post date:** [June 25, 2020, 2:27pm UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199/5 "2020-06-25T14:27:02Z")

</div>

Yes, this is exactly what I want to ask. Thanks again.  
BTW, I modified the earlier simple model as you suggested, by vectorization and non-centered parameterization (and half-normal prior for variance), now the Rhat for all the parameters is below 1.01, looks great. Computation time also reduced. Your post helped a lot.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [June 25, 2020, 2:49pm UTC](https://discourse.mc-stan.org/t/question-about-rhat-statistics/16199/6 "2020-06-25T14:49:48Z")

</div>

Great, glad that helped!