Help with (partially?) nested data/model

[Eric_Green]

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)

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:

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:

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:

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:

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:

# 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.

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

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:

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  

[Eric_Green]

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).

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.

 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.

[Eric_Green]

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:

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 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).

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).)

[Eric_Green]

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) 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.

   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:

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.

Candlish et al. supplemental file (docx):

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

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.

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?