# Brms approach to Stan models for nested versus non-nested multi-level models?

**URL:** <https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278>\
**Category:** brms\
**Created:** [June 15, 2019, 3:54pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278 "2019-06-15T15:54:24Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [June 15, 2019, 3:54pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/1 "2019-06-15T15:54:24Z")

</div>

`brms` underlying `Stan` code for a non-nested multi-level model, e.g.,

`make_stancode(y ~ 1 + (1 | re_1) + (1 | re_2) + (1 | re_3), data = d)`

is identical to the `Stan` code for a nested model,

`make_stancode(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)`

whether or not the nesting is natural. It looks like the difference between `brms`-generated `Stan` models are in how they index the random effects. If naturally nested, it seems each random effect assumes the indexing of the level before it. Is that right? And if not-naturally nested, how does `brms` approach indexing?

Here’s some code:

```
# Simulated data
N <- 100
y <- rnorm(N)

# naturally nested random effects
re_1 <- sample(LETTERS[1:5], N, replace = T)
re_2 <- paste0(re_1, sample(seq(5), N, replace = T))
re_3 <- paste0(re_2, sample(seq(5), N, replace = T))
d <- data.frame(y, re_1, re_2, re_3)

# create stan code and data for nonnested and nested models 
# using naturally nested random effects

library(brms)
nonnested_code <- make_stancode(y ~ 1 + (1 | re_1) + (1 | re_2) + (1 | re_3), data = d)
nonnested_data <- make_standata(y ~ 1 + (1 | re_1) + (1 | re_2) + (1 | re_3), data = d)
nested_code <- make_stancode(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)
nested_data <- make_standata(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)

# test for equality in code and data
identical(nonnested_code, nested_code) # TRUE
identical(nonnested_data, nested_data) # TRUE

# not naturally nested random effects
re_2 <- paste0(sample(seq(5), N, replace = T))
re_3 <- paste0(sample(seq(5), N, replace = T))
d <- data.frame(y, re_1, re_2, re_3)

# create stan code and data for nonnested and nested models 
# using naturally nested random effects

nonnested_code <- make_stancode(y ~ 1 + (1 | re_1) + (1 | re_2) + (1 | re_3), data = d)
nonnested_data <- make_standata(y ~ 1 + (1 | re_1) + (1 | re_2) + (1 | re_3), data = d)
nested_code <- make_stancode(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)
nested_data <- make_standata(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)

# test for equality in code and data
identical(nonnested_code, nested_code) # TRUE
identical(nonnested_data, nested_data) # FALSE

writeLines(nested_code) # look at Stan code
str(nonnested_data)
str(nested_data)
```

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 15, 2019, 6:50pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/2 "2019-06-15T18:50:22Z")

</div>

`(1 | A/B)` translates to `(1 | A) + (1 | A:B)` where `A:B` simply means creating a new grouping factor with the levels of `A` and `B` pasted together. This is explained a little bit longer in [https://journal.r-project.org/archive/2018/RJ-2018-017/index.html](https://journal.r-project.org/archive/2018/RJ-2018-017/index.html)

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [June 15, 2019, 8:27pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/3 "2019-06-15T20:27:18Z")

</div>

Thanks. That’s helpful. In looking at the paper, I do have a follow-up question on the paper description,

> The / operator indicates nested grouping structures and expands one grouping factor into two or more when using multiple / within one term. For instance, if the user were to write (1 | g1/g2), brms will expand this to (1 | g1) + (1 | g1:g2). Instead of |, users may use || in grouping terms to prevent group-level correlations from being modeled.

In reviewing underlying Stan code and data from my example, I want to make sure I understand how the group-level correlations are being modeled, and when I compare:

```
nested_code <- make_stancode(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)
nested_data <- make_standata(y ~ 1 + (1 | re_1 / re_2 / re_3), data = d)

nested_codenc <- make_stancode(y ~ 1 + (1 || re_1 / re_2 / re_3), data = d)
nested_datanc <- make_standata(y ~ 1 + (1 || re_1 / re_2 / re_3), data = d)

# test for equality in code and data
identical(nested_code, nested_codenc) # TRUE
identical(nested_data, nested_datanc) # TRUE

```

The model (stan code and stan data) using `|` seems identical to the model using `||`

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 15, 2019, 8:37pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/4 "2019-06-15T20:37:09Z")

</div>

That’s because there is nothing to correlate in your model. Replace `(1 | ...)` by `(1 + x | ...)` and you will see the difference.

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [June 15, 2019, 8:51pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/5 "2019-06-15T20:51:58Z")

</div>

Ah, that makes sense. I misunderstood the sentence as meaning correlation across groups instead of correlation of predictors within each group.

I’d like to have the group `(1 | g1)` provide/share information for the nested group `(1 | g1:g2)`.

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 15, 2019, 8:55pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/6 "2019-06-15T20:55:42Z")

</div>

Then you need to make a conceptual proposal of how this should look like :-)

---

<div class="post-metadata">

**Author:** ![MegBallard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/megballard/32/5155_2.png) [@MegBallard](https://discourse.mc-stan.org/u/MegBallard)\
**Post date:** [June 19, 2019, 12:48am UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/7 "2019-06-19T00:48:30Z")

</div>

Learning Bayesian methods and Stan coding is difficult and something I don’t know a lot about.

Using and interpreting mulitlevel/mixed/hierarchical models is also hard, much harder than most people realize.

Check out general help on specifying multilevel models.

- Ben Bolker’s [FAQ](http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#model-specification)  
\*This oft cited answer in [this stack exchange thread](https://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet)  
\*Doug Bates article about [mixed models](https://arxiv.org/pdf/1406.5823.pdf) in lme4 (he also has a lot of useful slide shows and lecture notes)  
\*There is a large archive of questions in archived mailing list [r-sig-mixed-models] [archive](https://stat.ethz.ch/pipermail/r-sig-mixed-models/)

Obviously these are not bayesian methods but the formula syntax is largely the same except where explicitly cited in the **brms** manual.

Nested levels should always be explicitly coded so any specification you use is at least not going to be messed up by that one step.

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [June 20, 2019, 10:05pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/8 "2019-06-20T22:05:23Z")

</div>

Would:

> [@Trying to make a "three-level nested linear model" run faster](https://discourse.mc-stan.org/t/trying-to-make-a-three-level-nested-linear-model-run-faster/3405/7):
>
> So, with the help of @wds15, I post here the vectorized code of the model, with the non-centered parametrization. I hope someone else will find it useful in future. Many thanks to Kazuki Yoshida that initially shared his code. data { // Define variables in data // Number of level-1 observations int\<lower=0\> Ni; // Number of level-2 clusters int\<lower=0\> Nj; // Number of level-3 clusters int\<lower=0\> Nk; // Cluster IDs int\<lower=1, upper=Nj\> classid[Ni]; int\<lower=1, upper…

Be a start?

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 20, 2019, 10:09pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/9 "2019-06-20T22:09:28Z")

</div>

I don’t see how this can help since the model you linked to does not contain any correlation or similar measure to share information between random effects of different levels. Or am I overlooking something?

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [June 23, 2019, 3:31am UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/10 "2019-06-23T03:31:24Z")

</div>

As a start, it seems it may be useful to just change the priors from independent groups to using the broader group as the mean for the nested group; same for the standard deviations. So, the brms-generated Stan code in the above example,

```
target += student_t_lpdf(sd_1 | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10);
target += normal_lpdf(z_1[1] | 0, 1);

target += student_t_lpdf(sd_2 | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10);
target += normal_lpdf(z_2[1] | 0, 1);

target += student_t_lpdf(sd_3 | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10);
target += normal_lpdf(z_3[1] | 0, 1);

```

Could be changed to,

```
target += student_t_lpdf(sd_1 | 3, 0, 10)
    - 1 * student_t_lccdf(0 | 3, 0, 10);
target += normal_lpdf(z_1[1] | 0, 1);
  
target += student_t_lpdf(sd_2 | 3, sd_1, 10)
    - 1 * student_t_lccdf(0 | 3, sd_1, 10);
target += normal_lpdf(z_2[1][J_2] | z_1[1][J_1], 1);
  
target += student_t_lpdf(sd_3 | 3, sd_2, 10)
    - 1 * student_t_lccdf(0 | 3, sd_2, 10);
target += normal_lpdf(z_3[1][J_3] | z_2[1][J_2], 1);
```

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [June 23, 2019, 9:55am UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/11 "2019-06-23T09:55:43Z")

</div>

For me to seriously consider a proposal, we cannot just change a bunch of lines of Stan code and hope for the model to be still reasonable. Instead it is important to have a clear theory of what is the target, how the desired model looks like, what properties it has and why those makes sense.

With respect to what you are writing above it seems that, at least party, you want to use a centered parameterization (of the random effects), while brms uses the non-centered parameterization throughout. The latter has the advanatage of also applying to non-nested models and to usually having better sampling properties. Apart from that, if both parameterizations exist, they are equivalent.

---

<div class="post-metadata">

**Author:** ![ssp3nc3r](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ssp3nc3r/32/120_2.png) [@ssp3nc3r](https://discourse.mc-stan.org/u/ssp3nc3r)\
**Post date:** [June 23, 2019, 2:44pm UTC](https://discourse.mc-stan.org/t/brms-approach-to-stan-models-for-nested-versus-non-nested-multi-level-models/9278/12 "2019-06-23T14:44:54Z")

</div>

Thank you for your explanations, I understand more fully what is intended with the brms-generated code. I didn’t mean to imply I’d like you to rework what brms does for me. It’s a really helpful tool, and like that I can even use its code to pick up where it left off on the occasions when I’d like a bit different expression.

I just included the code to show one way I think the information could be minimally shared between the groups, independent of whether a centered parameterization was used.