partial_nesting.R

library(faux)
library(tidyverse)

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)
df1 <- add_random(group = n_groups) %>%     
  add_random(family = 3, .nested_in = "group") %>%    # 3 families per group
  add_random(member = 4, .nested_in = "family")       # 4 members per family

# assign families to arm 
t_assignment <- df1 %>%
  distinct(family) %>%
  mutate(treatment = c(rep(1, nrow(.)/2),
                       rep(0, nrow(.)/2)))

# create singleton clusters for group 
df2 <- df1 %>%
  left_join(t_assignment) %>%
  mutate(group = case_when(
    treatment==1 ~ group,
    TRUE ~ member
  ))

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

# join back leader assignments and calculate dv
df3 <- df2 %>%
  left_join(a_assignment) %>%
# create singleton clusters for leader 
  mutate(leader = case_when(
    treatment==1 ~ leader,
    TRUE ~ member
  )) %>%
  arrange(leader, group, family, member) %>%
  # add random intercepts
  add_ranef("leader", u0l = 0.05) %>%
  add_ranef("group", u0g = 0.1) %>%
  add_ranef("family", u0f = 0.2) %>%
  add_ranef(sigma = 0.4) %>%
# calculate DV
  mutate(dv = b0 + 
              b1*treatment +    # is 0 for control
              u0l*treatment +   # is 0 for control
              u0g*treatment +   # is 0 for control
              u0f +             # apply family effect to all
              sigma             # apply sigma to all (homoscedastic model)
  ) %>% 
  # finalize
  select(member, family, treatment, group, leader, dv) %>%
  arrange(family)