Simulation for random factors

ranfac.R

library(faux)
library(tidyverse)
library(lme4)

# simulate nested data with 12 subjects tested at 6 time points with 100 observations per time point
dat <- add_random(ID = 12) %>%
  add_random(Time = 1:6) %>%
  add_within(rep = 1:100) %>%
  add_ranef(.by = "ID", id_int = 2, id_slope = 1) %>%
  add_ranef(.by = "Time", time_int = 1) %>%
  add_ranef(error = 2) %>%
  mutate(time = as.numeric(gsub("T", "", Time)),
         y = 0 + id_int + (0.4+id_slope)*time + error)

# plot for sense
ggplot(dat, aes(time, y)) +
  geom_point() +
  geom_smooth(method = lm) +
  facet_wrap(~ID)

# only random intercept for ID
m1 <- lmer(y ~ time + (1 | ID), dat)
summary(m1)

# only random intercept for ID:Time
m2 <- lmer(y ~ time + (1 | ID:Time), dat)
summary(m2)

# random intercepts for ID and Time (crossed)
m3 <- lmer(y ~ time + (1 | ID) + (1 | Time), dat)
summary(m3)

# random intercepts for ID and Time (nested)
m4 <- lmer(y ~ time + (1 | ID) + (1 | ID:Time), dat)
summary(m4)

# random intercept plus random slope for ID
m5 <- lmer(y ~ time + (1 + time | ID), dat)
summary(m5)



## Simulation to show inflated fixed effect estimate significant for ID/time model

# sim effect = 0 with some random factor structure
x <- function() {
  dat <- add_random(ID = 12) %>%
    add_random(Time = 1:6) %>%
    add_within(rep = 1:100) %>%
    add_ranef(.by = "ID", id_int = 2, id_slope = 1) %>%
    add_ranef(.by = "Time", time_int = 1) %>%
    add_ranef(error = 2) %>%
    mutate(time = as.numeric(gsub("T", "", Time)),
           y = 0 + id_int + (0+id_slope)*time + error)
  
  # return p-value for effect of time for both models
  list(
    noslope = lmer(y~time + (1 | ID/time), dat) %>% broom.mixed::tidy() %>% pull(p.value) %>% pluck(2),
    slope = lmer(y~time + (1 + time| ID), dat) %>% broom.mixed::tidy() %>% pull(p.value) %>% pluck(2)
  )
}

# run 100 times
df <- map_df(1:100, ~x())

p <- ggplot(df, aes(slope, noslope) ) +
  geom_point() +
  labs(title = "P-values from different models for a null effect",
       x = "y ~ time + (1 + time| ID)",
       y = "y ~ time + (1 | ID/time)")

ggExtra::ggMarginal(p, type = "histogram")