# Help with (partially?) nested data/model

**URL:** <https://discourse.mc-stan.org/t/help-with-partially-nested-data-model/26771>\
**Category:** brms\
**Tags:** specification, hierarchical-model, brms\
**Created:** [March 12, 2022, 4:49am UTC](https://discourse.mc-stan.org/t/help-with-partially-nested-data-model/26771 "2022-03-12T04:49:13Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Eric\_Green](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric_green/32/13884_2.png) [@Eric\_Green](https://discourse.mc-stan.org/u/Eric_Green)\
**Post date:** [March 12, 2022, 4:49am UTC](https://discourse.mc-stan.org/t/help-with-partially-nested-data-model/26771/1 "2022-03-12T04:49:13Z")

</div>

I’m trying to simulate a nested dataset for a trial of a family intervention where families are assigned to treatment or an inactive control. Data will be collected from several eligible members in each family. Thus, there will be nesting of the form:

`family members (d) < families (c) < groups (b) < group facilitators (a)`

 ![coping together trial diagram (3)](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/ed9aeb12e2bb7bd5c1afd852ecea0ccc174c1941.png)

I’m not sure I have my terminology correct, but I believe this is partial nesting. Each lower level unit is in one and only one higher level unit, but families (c) in the control arm won’t have groups (b) or group facilitators (a).

I’m unsure about two things:

1. How to simulate and set the random intercepts for all random effects.
2. How to specify the `brms()` model to account for this partial nesting.

**1. How to simulate and set the random intercepts for all random effects**

Here is my attempt at setting up the data structure:

```nohighlight
library(faux)
library(tidyverse)
library(brms)

set.seed(1)

b0 <- 2
b1 <- 0.2

df1 <- add_random(a = 3) %>% # 3 facilitators
  add_random(b = 2, .nested_in = "a") %>% # 2 groups per facilitator
  add_random(c = 3, .nested_in = "b") %>% # 3 families per group
  add_random(d = 3, .nested_in = "c") %>% # 3 members per family
  add_within("d", time = c("pre", "post")) %>%
  add_between(.by = "b",
              arm = c("treatment", "control")) %>%
# add dummy
  add_recode("arm", "treatment", control = 0, treatment = 1) %>%
  add_recode("time", "post", pre = 0, post = 1) %>%
  select(-arm) %>%
# add random intercepts
  add_ranef("a", u0a = 0.01) %>%
  add_ranef("b", u0b = 0.05) %>%
  add_ranef("c", u0c = 0.1) %>%
  add_ranef("d", u0d = 0.2) %>%
  add_ranef(sigma = 0.4) %>%
# calculate DV
  mutate(dv = b0 + u0a + u0b + u0c + u0d + 
         (b1 * treatment * post) + sigma) %>%
  select(-post) 

```

This yields:

```nohighlight
df1

# A tibble: 108 × 12
   a b c d time treatment u0a u0b u0c u0d sigma dv
   <chr> <chr> <chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 a1 b1 c01 d01 pre 1 -0.00626 0.0798 -0.0305 -0.294 -0.0541 1.69
 2 a1 b1 c01 d01 post 1 -0.00626 0.0798 -0.0305 -0.294 0.471 2.42
 3 a1 b1 c01 d02 pre 1 -0.00626 0.0798 -0.0305 -0.0956 -0.609 1.34
 4 a1 b1 c01 d02 post 1 -0.00626 0.0798 -0.0305 -0.0956 0.238 2.38
 5 a1 b1 c01 d03 pre 1 -0.00626 0.0798 -0.0305 0.0836 0.133 2.26
 6 a1 b1 c01 d03 post 1 -0.00626 0.0798 -0.0305 0.0836 0.425 2.75
 7 a1 b1 c02 d04 pre 1 -0.00626 0.0798 0.151 0.272 -0.122 2.37
 8 a1 b1 c02 d04 post 1 -0.00626 0.0798 0.151 0.272 0.148 2.84
 9 a1 b1 c02 d05 pre 1 -0.00626 0.0798 0.151 -0.0206 0.107 2.31
10 a1 b1 c02 d05 post 1 -0.00626 0.0798 0.151 -0.0206 -0.217 2.19
# … with 98 more rows

```

(I don’t think I’m simulating the pre data correctly)

Let’s look at family `c04`:

```nohighlight
df1 %>% filter(c=="c04")

# A tibble: 6 × 12
  a b c d time treatment u0a u0b u0c u0d sigma dv
  <chr> <chr> <chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 a1 b2 c04 d10 pre 0 -0.00626 0.0165 -0.0621 -0.0789 -0.189 1.68
2 a1 b2 c04 d10 post 0 -0.00626 0.0165 -0.0621 -0.0789 -0.248 1.62
3 a1 b2 c04 d11 pre 0 -0.00626 0.0165 -0.0621 -0.0119 0.0168 1.95
4 a1 b2 c04 d11 post 0 -0.00626 0.0165 -0.0621 -0.0119 -0.364 1.57
5 a1 b2 c04 d12 pre 0 -0.00626 0.0165 -0.0621 0.220 0.0632 2.23
6 a1 b2 c04 d12 post 0 -0.00626 0.0165 -0.0621 0.220 -0.262 1.91

```

This family is assigned to the control arm, so they are not part of any group (b) and therefore are not led by any group facilitator (a). _ **Would I just set `u0a` and `u0b` to `0` where `treatment==0`? Relatedly, would I do the same for everyone where `time=="pre"` since there is no nesting by group/facilitator for the baseline round of data collection?** _

**2. How to specify the `brms()` model to account for this partial (?) nesting**

Assume for the moment that I simulated the data correctly. We can reshape it for analysis as follows:

```nohighlight
df2 <- df1 %>%
  pivot_wider(id_cols = c(a, b, c, d, treatment),
              names_from = time,
              values_from = dv) %>%
  rename(y_pre = pre,
         y_post = post)

```

Which gives us:

```nohighlight
# A tibble: 54 × 7
   a b c d treatment y_pre y_post
   <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
 1 a1 b1 c01 d01 1 1.69 2.42
 2 a1 b1 c01 d02 1 1.34 2.38
 3 a1 b1 c01 d03 1 2.26 2.75
 4 a1 b1 c02 d04 1 2.37 2.84
 5 a1 b1 c02 d05 1 2.31 2.19
 6 a1 b1 c02 d06 1 2.79 2.97
 7 a1 b1 c03 d07 1 2.38 2.94
 8 a1 b1 c03 d08 1 2.06 1.53
 9 a1 b1 c03 d09 1 1.80 1.74
10 a1 b2 c04 d10 0 1.68 1.62
# … with 44 more rows

```

This is how I’ve specified the nesting, but I’m not certain it’s correct since the families in the control arm are not nested in groups or facilitators.

```nohighlight
brm(y_post ~ 0 + Intercept + treatment + y_pre + (1 | a:b:c),
    data = df2, 
    control = list(adapt_delta = 0.9),
    cores = parallel::detectCores())

```

```nohighlight
Group-Level Effects: 
~a:b:c (Number of levels: 18) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.08 0.06 0.00 0.24 1.00 1527 1297

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.76 0.32 1.14 2.36 1.00 2521 2691
treatment 0.27 0.13 0.02 0.53 1.00 4034 3059
y_pre 0.12 0.15 -0.16 0.40 1.00 2675 2711

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.42 0.04 0.35 0.52 1.00 4477 2410

```

Info:

```nohighlight
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats graphics grDevices utils datasets methods  
[7] base     

other attached packages:
 [1] sjPlot_2.8.9 tidybayes_3.0.2 cmdstanr_0.2.2     
 [4] brms_2.16.3  

```

---

<div class="post-metadata">

**Author:** ![Eric\_Green](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric_green/32/13884_2.png) [@Eric\_Green](https://discourse.mc-stan.org/u/Eric_Green)\
**Post date:** [March 12, 2022, 3:23pm UTC](https://discourse.mc-stan.org/t/help-with-partially-nested-data-model/26771/2 "2022-03-12T15:23:50Z")

</div>

I don’t think I’m simulating pre/post data correctly, so I try again below without this additional layer of complexity to focus on my original questions. I also bumped up the sample size to more closely represent our planning (originally I wanted to create a small dataframe to make it easier to see the structure).

```nohighlight
library(faux)
library(tidyverse)
library(brms)

set.seed(1)

b0 <- 2
b1 <- 0.2

# simulate data
df3 <- add_random(a = 10) %>% # 10 facilitators
  add_random(b = 3, .nested_in = "a") %>% # 3 groups per facilitator
  add_random(c = 3, .nested_in = "b") %>% # 3 families per group
  add_random(d = 4, .nested_in = "c") %>% # 4 members per family
  add_between(.by = "b",
              arm = c("treatment", "control")) %>%
# add dummy
  add_recode("arm", "treatment", control = 0, treatment = 1) %>%
  select(-arm) %>%
# add random intercepts
  add_ranef("a", u0a = 0.01) %>%
  add_ranef("b", u0b = 0.05) %>%
  add_ranef("c", u0c = 0.1) %>%
  add_ranef("d", u0d = 0.2) %>%
  add_ranef(sigma = 0.4) %>%
  # calculate DV
  mutate(dv = b0 + u0a + u0b + u0c + u0d + 
           (b1 * treatment) + sigma) 

# fit model
brm(dv ~ 0 + Intercept + treatment + (1 | a:b:c),
    data = df3, 
    control = list(adapt_delta = 0.9),
    cores = parallel::detectCores())

```

This recovers b0, b1, and sigma pretty well.

```nohighlight
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: dv ~ 0 + Intercept + treatment + (1 | a:b:c) 
   Data: df3 (Number of observations: 360) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~a:b:c (Number of levels: 90) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.13 0.04 0.03 0.21 1.00 890 995

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.97 0.04 1.90 2.05 1.00 2937 2879
treatment 0.20 0.06 0.09 0.31 1.00 2860 2715

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.43 0.02 0.40 0.47 1.00 1976 2111

```

**I’m still unsure about two things:**

1. How to simulate and set the random intercepts for all random effects.
2. How to specify the `brms()` model to account for this partial nesting.

---

<div class="post-metadata">

**Author:** ![Eric\_Green](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric_green/32/13884_2.png) [@Eric\_Green](https://discourse.mc-stan.org/u/Eric_Green)\
**Post date:** [March 13, 2022, 4:51pm UTC](https://discourse.mc-stan.org/t/help-with-partially-nested-data-model/26771/3 "2022-03-13T16:51:06Z")

</div>

Going a bit deeper on question 2:

**2. How to specify the `brms()` model to account for this partial nesting**

Ignoring the repeated measures (dropping pre-treatment measurement) in my previous post, I specified:

```no-highlight
brm(dv ~ 0 + Intercept + treatment + (1 | a:b:c),
    data = df3, 
    control = list(adapt_delta = 0.9),
    cores = parallel::detectCores())

```

However, I came across [this post](https://rpsychologist.com/r-guide-longitudinal-lme-lmer#partially-nested-models) by Kristoffer Magnusson showing how to fit a similar model in lme4 that he calls “partially nested”:

> Partially nesting occurs when we have nesting in one group but not the other. For instance, we might compare a treatment group to a wait-list condition. Subjects in the wait-list will not be nested, but subjects in treatment group will be nested within therapists.

This is similar to my setup where families (c) in the control arm will not be nested in groups (b) or facilitators (a).

 ![diagram2](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/1e7165661b67cb19dc88fdc07685ce85c79667b5.png)

**Am I showing the correct translation to my structure?**

- Should it be `a:b:c` or `a/b/c`?
- Do I need more than just `(0 + treatment | a)` or does this account for the subsequent partial nesting? (I noted that Kristoffer did not include `(0 + tx | therapist:subjects)`.)

---

<div class="post-metadata">

**Author:** ![Eric\_Green](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/eric_green/32/13884_2.png) [@Eric\_Green](https://discourse.mc-stan.org/u/Eric_Green)\
**Post date:** [March 16, 2022, 3:02am UTC](https://discourse.mc-stan.org/t/help-with-partially-nested-data-model/26771/4 "2022-03-16T03:02:50Z")

</div>

I think I solved (1) and made progress on (2). I’m having problems with the `brm()` model in (2) and I would really appreciate feedback on the `lmer()` code I’m applying to `brm()`.

**1. How to simulate and set the random intercepts for all random effects.**

A colleague shared a paper by [Candlish et al. (2018)](https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-018-0559-x) that simulated the analysis of partially nested RCTs. They describe one approach to imposing clustering on the non-nested levels of the control arm that “treats each individual in the control arm as their own cluster of size one (singleton clusters)”

Here’s a partial example showing a leader (a) with 2 groups (b) where each group has 2 families (c) and each family has 2 members (d). Notice where `treatment==0` the values for group (b) and leader (a) are unique. This makes each person in the control arm their own cluster for group and leader since these are not relevant levels for the control arm.

```nohighlight
   d c treatment b a    
   <chr> <chr> <dbl> <chr> <chr>
 1 d001 c01 1 b01 l1   
 2 d002 c01 1 b01 l1   
 3 d003 c02 1 b01 l1   
 4 d004 c02 1 b01 l1   
 5 d005 c03 1 b02 l1   
 6 d006 c03 1 b02 l1   
 7 d007 c04 1 b02 l1   
 8 d008 c04 1 b02 l1   
 9 d009 c05 0 d009 d009   
10 d010 c05 0 d010 d010 

```

Here’s how I coded it:

```nohighlight
library(faux)
library(tidyverse)
library(brms)

set.seed(1)

b0 <- 2
b1 <- 0.2

# simulate data
n_leader <- 4 # number of facilitators 
grp_per_lead <- 3 # groups per facilitator (treated)
n_groups <- (n_leader * grp_per_lead)*2 # x2 to give equal number control

# create initial structure (excluding leaders)
df4 <- add_random(b = n_groups) %>%     
  add_random(c = 3, .nested_in = "b") %>% # 3 families per group
  add_random(d = 4, .nested_in = "c") # 4 members per family
  
# assign families to arm 
t_assignment <- df4 %>%
  distinct(c) %>%
  mutate(treatment = c(rep(1, nrow(.)/2),
                       rep(0, nrow(.)/2)))

# create singleton clusters for group 
df5 <- df4 %>%
  left_join(t_assignment) %>%
  mutate(b = case_when(
    treatment==1 ~ b,
    TRUE ~ d
  ))

# assign groups to leaders
a_assignment <- df5 %>% 
  filter(treatment==1) %>%
  distinct(b) %>%
  mutate(a = rep(row_number(), length.out = n(), 
                 each = grp_per_lead)) %>%
  mutate(a = paste0("l", a))

# join back leader assignments and calculate dv
df6 <- df5 %>%
  left_join(a_assignment) %>%
# create singleton clusters for leader 
  mutate(a = case_when(
    treatment==1 ~ a,
    TRUE ~ d
  )) %>%
  arrange(a, b, c, d) %>%
# add random intercepts
  add_ranef("a", u0a = 0.05) %>%
  add_ranef("b", u0b = 0.1) %>%
  add_ranef("c", u0c = 0.2) %>%
  add_ranef(sigma = 0.4) %>%
# calculate DV
  mutate(dv = b0 + 
              b1*treatment +
              u0a*treatment + # is 0 for control
              u0b*treatment + # is 0 for control
              u0c + # apply family effect to all
              sigma # apply sigma to all (homoscedastic model)
  ) %>% 
# finalize
  select(d, c, treatment, b, a, dv) %>%
  arrange(c)

```

**2. How to specify the `brm()` model to account for this partial nesting.**

Still following Candlish et al. (though I have more levels of nesting), I fit a partially nested homoscedastic mixed effects model.

 ![Screen Shot 2022-03-15 at 10.55.36 PM](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/5/58669427af306b5f80600eb66f72d6e3facb6602.png)

Candlish et al. supplemental file ([docx](https://static-content.springer.com/esm/art%3A10.1186%2Fs12874-018-0559-x/MediaObjects/12874_2018_559_MOESM2_ESM.docx)):

 ![Screen Shot 2022-03-15 at 10.57.20 PM](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/8/8a80a0468cdde850761eadf992f4ff1e784dd7d7.jpeg)

My colleague and I produced equivalent results with `lmer()` and Stata with the following code:

```nohighlight
lmer(dv ~ treatment + (0 + treatment | a/b) + (1 | c), 
     data=df6)

mixed dv treatment || c: || b:treatment, nocons || a:treatment, nocons reml dfmethod(sat) iterate(100) stddev

# A tibble: 6 × 6
  effect group term estimate std.error statistic
  <chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed NA (Intercept) 2.01 0.0406 49.6 
2 fixed NA treatment 0.164 0.0853 1.93
3 ran_pars b:a sd__treatment 0.126 NA NA   
4 ran_pars a sd__treatment 0.103 NA NA   
5 ran_pars c sd__(Intercept) 0.0927 NA NA   
6 ran_pars Residual sd__Observation 0.450 NA NA   

```

The `brm()` model threw some warnings and did not recover random effect a.

```nohighlight
brm(dv ~ 0 + Intercept + treatment + (0 + treatment | a/b) + (1 | c),
    data = df6, 
    control = list(adapt_delta = 0.9),
    cores = parallel::detectCores(),
    backend = "cmdstanr")

```

**Does the `lmer()` code translate directly to `brm()` in this case?**