sim-lme4.R

library(tidyverse)
library(lme4)

gen_level2 <- function(n = 30, offsetb0 = 0, scaleb0 = 1, offsetb1 = 0, scaleb1 = 1){
  unit <- tibble::tibble(b0 = rnorm(n, offsetb0, scaleb0),
                         b1 = rnorm(n, offsetb1, scaleb1))
  return(unit)
}
gen_level2()

gen_level1 <- function(level2, mean = 0, sd = 1, noise = 1, m = 100){
  reg <- function(...){
    dots <- list(...)
    tibble::tibble(x = rnorm(m, mean, sd), y = dots$b0 + dots$b1*x + rnorm(m, 0, noise))
  }
  purrr::pmap_dfr(level2, reg, .id = "id")
}
set.seed(345)
private <-
  gen_level1(gen_level2(
    n = 20,
    offsetb0 = 3,
    scaleb0 = 30,
    offsetb1 = 12,
    scaleb1 = 5
  ),
  m = 100,
  mean = 5,
  sd = 2,
  noise = 15) %>% 
  mutate(y = scale(y) * 7 + 16,
         fund = "private")

public <-
  gen_level1(gen_level2(
    n = 20,
    offsetb0 = 3,
    scaleb0 = 30,
    offsetb1 = 7,
    scaleb1 = 5
  ),
  m = 100,
  mean = 5,
  sd = 2,
  noise = 20) %>% 
  mutate(y = scale(y) * 7 + 12,
         fund = "public")
dat <- rbind(private, public) %>% 
  filter(between(x, 0, 10),
         y > 0) %>% 
  group_by(fund) %>% 
  mutate(y = round(as.numeric(y)))

groups <- tibble(id = rep(unique(dat$id), 2),
                 fund = rep(unique(dat$fund), each = 20),
                 group = c(rep(1:5, each = 4),
                           rep(5:1, each = 4)),
                 z = round(runif(40, 5, 50) + (fund == "public")*10))

dat <- left_join(dat, groups, by = c("id", "fund"))

dat <- rename(dat, psf = x, stay = y, beds = z)
  
summary(lmer(stay ~ fund + psf + beds + (1 + psf| id), data = dat))

library(ggplot2)
ggplot(dat, aes(psf, stay, color = id, shape = fund)) +
  geom_point(alpha = .5) +
  facet_grid(group ~ fund) +
  theme_minimal() +
  theme(legend.position = "none",
        strip.text.y = element_blank()) +
  xlab("psycho-social functioning") +
  ylab("stay")