# Random Effects Model

**URL:** <https://discourse.mc-stan.org/t/random-effects-model/14571>\
**Category:** Modeling\
**Tags:** mixed-model, fitting-issues, specification\
**Created:** [April 24, 2020, 3:46am UTC](https://discourse.mc-stan.org/t/random-effects-model/14571 "2020-04-24T03:46:21Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![the\_stan\_man](https://avatars.discourse-cdn.com/v4/letter/t/76d3ee/32.png) [@the\_stan\_man](https://discourse.mc-stan.org/u/the_stan_man)\
**Post date:** [April 24, 2020, 3:46am UTC](https://discourse.mc-stan.org/t/random-effects-model/14571/1 "2020-04-24T03:46:21Z")

</div>

I am attempting to simulate a Gage R&R study using stan. In this study there are parts and operators, as well as part-operator interactions. For this experiment, I am generating the data with known parameters and attempting to extract those parameters using stan. The random effects model I am fitting is shown below:

y\_{ijk} ~ \mu + \alpha\_i + \beta\_j + \alpha\beta\_{ij} + \epsilon\_{ijk}  
where  
\alpha\_i ~ \mathcal{N}(0,\sigma\_\alpha) … (i= 1, 2, …, nO)  
\beta\_j ~ \mathcal{N}(0,\sigma\_\beta) … (j= 1, 2, …, nP)  
\alpha\beta\_{ij} ~ \mathcal{N}(0,\sigma\_{\alpha\beta}) … (ij= 1, 2, …, nOP)  
\epsilon\_{ijk} ~ \mathcal{N}(0,\sigma\_R) … (ijk= 1, 2, …, N)

My goal is to obtain parameter estimates for \sigma\_\alpha,\sigma\_\beta,\sigma\_{\alpha\beta}, and \sigma\_R.

* * *

I generate data with [\sigma\_\alpha, \sigma\_\beta, \sigma\_{\alpha\beta}, \sigma\_R] = [0.03, 0.8, 0.002, 0.04] where [nO, nP, nOP, N] = [3, 10, 30, 90]. The data is attached here:  
[stan\_data.txt](https://discourse.mc-stan.org/uploads/short-url/qpBqSFt60g08RRxB5P37hOlH5Fg.txt) (1.6 KB)

My stan code for this model is below:

```stan
data {
    int<lower=1> N;
    int<lower=1> nP;
    int<lower=1> nO;
    int<lower=1> nOP;

    int oper[N];
    int part[N];
    int inter[N];

    real y[N];
}

parameters {
    vector[nO] O;
    vector[nP] P;
    vector[nOP] OP;

    real<lower=0> sigma_O;
    real<lower=0> sigma_P;
    real<lower=0> sigma_OP;
    real<lower=0> sigma_R;
}

transformed parameters {
    vector[N] theta = sigma_O*O[oper] + sigma_P*P[part] + sigma_OP*OP[inter];
}

model {
    sigma_O ~ gamma(0.01, 0.01);
    sigma_P ~ gamma(0.01, 0.01);
    sigma_OP ~ gamma(0.01, 0.01);
    sigma_R ~ gamma(0.01, 0.01);

    O ~ normal(0, 1);
    P ~ normal(0, 1);
    OP ~ normal(0, 1);

    y ~ normal( theta, sigma_R );
}

```

* * *

When I pass the data to stan, I set y = y - y.mean().

The results I achieve are below:  
sO: 0.00035  
sP: 1.0639  
sOP: 0.00020  
sR: 0.04932

I also generated a similar dataset with the same values of sigmas, but with [nO, nP, nOP, N] = [16, 16, 256, 1280].

The results of that are shown below:  
sO: 0.61582  
sP: 0.20466  
sOP: 8.8469\bullet10^{-5}  
sR: 0.04155

The current model takes between 25 and 40 minutes to sample (4 chains, 2000 iter, 800 warmup).

How can I get the predictions to be more accurate? And how can I get the model to run faster?

Is there a better way for me to formulate my model?

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [April 24, 2020, 11:19am UTC](https://discourse.mc-stan.org/t/random-effects-model/14571/2 "2020-04-24T11:19:36Z")

</div>

Is there any reason that you omitted an intercept from your model? Without it, you are conveying a belief that your sample mean precisely matches the population intercept and that you are infinitely certain of its value.

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [April 24, 2020, 1:40pm UTC](https://discourse.mc-stan.org/t/random-effects-model/14571/3 "2020-04-24T13:40:11Z")

</div>

Seeing as N is much bigger than the other count variables, you can probably get speed-ups from computing only the unique combinations then indexing into that to define theta (see [example here](https://discourse.mc-stan.org/t/case-study-preview-speeding-up-stan-by-reducing-redundant-computation/14228)), but before you get too worried about performance you should pin down why this model is not recovering the data-generating parameters first. The fact that it takes so long suggests to me that something is awry in the first place, as this isn’t a particularly complicated model nor large data. First add an intercept (and a prior on said intercept) and see how that fares. Then if you’re still not getting good recovery, post your data-generating code here too as possibly the model is recovering the data-generating parameters accurately but you’re not generating the data in the way you think you are.

---

<div class="post-metadata">

**Author:** ![the\_stan\_man](https://avatars.discourse-cdn.com/v4/letter/t/76d3ee/32.png) [@the\_stan\_man](https://discourse.mc-stan.org/u/the_stan_man)\
**Post date:** [April 26, 2020, 12:57am UTC](https://discourse.mc-stan.org/t/random-effects-model/14571/4 "2020-04-26T00:57:34Z")

</div>

Thank you Mike, I will look into this and see what progress I can make.