Disclaimer
My code is available for use and alteration, but please give credit with a citation & link to the page.
Goal
Share code for simulating multilevel data, with optional interaction term.
Function to simulate 2 level data with an interaction
This function allows you to specify sample sizes, two main effects (main_x & main_z), an interaction, and other features of the data. If you specify interaction = 0, you can simply specify multilevel data with no interaction term. Of course the code could be extended to include as many effects as you need. It also provides a variable (w) that’s not associated with the outcome (outcome).
show code
simulate_interaction_w_2_levels <-
function(
n, ## number level 1 units
j, ## number level 2 units
intercept_lv1, ## intercept at level 1
main_x, ## main effect of x
main_z, ## main effect of z
interaction, ## interaction of x and z
residual_var_sd_lv1, ## standard deviation of residuals at level 1
random_int_mean_lv2, ## mean of random intercept at level 2
random_int_sd_lv2, ## standard deviation of random intercept at level 2,
start_seed = 123
){
## placeholder for level 2 outcomes
outcomes_j <- vector("list", j)
## for each variable, make list
## as long as level 2
## fill each element with a list of level 1 outcomes
x <- vector("list", j)
z <- vector("list", j)
w <- vector("list", j)
## School distribution (intercept at level 2)
## Standard deviation of random intercept at level 2
set.seed(start_seed)
a_j <- rnorm(j, random_int_mean_lv2, random_int_sd_lv2)
for(i in 1:j) {
## make a level 1 predictor variable:
## set multiple seeds that change each iteration
## prevents identical cases with seed
set.seed(start_seed+i)
x[[i]] <- rnorm(n)
set.seed(start_seed-i)
z[[i]] <- rnorm(n)
set.seed(-start_seed + i)
outcomes_j[[i]] <-
rnorm(
n,
intercept_lv1 +
interaction*z[[i]]*x[[i]] +
main_x*x[[i]] +
main_z*z[[i]] +
a_j[i],
## standard deviation of residual variance
residual_var_sd_lv1
)
set.seed(start_seed*197+197*i)
w[[i]] <- rnorm(n)
}
outcomes_df <-
data.frame(
id = rep(1:j, each = n),
x = unlist(x),
z = unlist(z),
outcome = unlist(outcomes_j),
w = unlist(w)
)
return(outcomes_df)
}
show code
# baseline theme for plots here
project_theme <-
theme(
axis.line = element_line(size = 3, lineend = 'round'),
legend.position = 'none',
plot.title.position = 'plot',
panel.grid.minor = element_blank(),
axis.title.x = element_text(size = 18),
axis.text = element_text(size = 16),
axis.title.y = element_text(size = 20),
plot.title = element_text(size = 25),
strip.text = element_text(size = 25),
)
Simulating data
We have 15 level 1 units, 71 level 2 units. The level 1 intercept of the outcome is 4.25 with a main effect of 1.25 for x, 2.00 for z, and 0.75 for their interaction. The function also makes a nuisance variable w, which does not contribute to outcome at all.
show code
sim_dat <-
simulate_interaction_w_2_levels(
n = 15, ## number level 1 units at each level 2
j = 71, ## number level 2 units
intercept_lv1 = 4.25, ## intercept at level 1
main_x = 1.25, ## main effect of x
main_z = 2.00, ## main effect of z
interaction = 0.75, ## interaction of x and z
residual_var_sd_lv1 = 2.00, ## standard deviation of residuals at level 1
random_int_mean_lv2 = 0, ## mean of random intercept at level 2
random_int_sd_lv2 = 1.00, ## standard deviation of random intercept at level 2,
start_seed = 123 ## ensure you can reproduce
) %>%
tibble()
Put formally, the simulation is approximating:
\[ \begin{aligned} \operatorname{outcome}_{i} &\sim N \left(\mu, \sigma^2 \right) \\ \mu &=\alpha_{j[i]} + \beta_{1}(\operatorname{x}) + \beta_{2}(\operatorname{z}) + \beta_{3}(\operatorname{w}) + \beta_{4}(\operatorname{x} \times \operatorname{z}) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for id j = 1,} \dots \text{,J} \end{aligned} \]
where
\[\mu = 4.25\] \[\sigma = 2.00\] \[\beta_{1} = 1.25\] \[\beta_{2} = 2.00\] \[\beta_{3} = 0\] \[\beta_{4} = 0.75\] \[J = 71\] \[\mu_{\alpha_{j}} = 0\] \[\sigma_{\alpha_{j}}= 1.00\]
Describing data & verifying simulation results
Here are the overall averages for the 1065 cases, with 71 level two units.
show code
# A tibble: 1 x 8
x_mean x_sd z_mean z_sd outcome_mean outcome_sd w_mean w_sd
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.0900 0.946 -0.0197 1.00 4.44 3.38 -0.0208 0.993Distribution of Variables
show code
figure_1 <-
sim_dat %>%
pivot_longer(
cols = -id,
names_to = 'variable',
values_to = 'statistic'
) %>%
group_by(variable) %>%
mutate(average_value = mean(statistic)) %>%
ungroup() %>%
ggplot(
aes(
x = statistic,
fill = variable)
) +
colorblindr::scale_fill_OkabeIto() +
colorblindr::scale_color_OkabeIto() +
geom_vline(aes(xintercept = average_value), size = 1.25, linetype = 2) +
geom_density(alpha = 0.7, color = 'black', size = 2) +
facet_wrap(~variable, scales = 'free') +
labs(
y = element_blank(),
x = 'Distribution of Simulated Values',
title = 'Figure 1. Distribution of Simulated Variables',
caption = 'Note the scales are free on each plot\nDotted Lines = Average Values') +
theme_minimal() +
project_theme +
theme(axis.text.y = element_blank())
figure_1
When you group by the id, the simulated averages are normally distributed around the values (Figure 1).
Distribution of Group Averages
show code
figure_2 <-
sim_dat %>%
group_by(id) %>%
summarize(
x = mean(x),
z = mean(z),
outcome = mean(outcome),
w = mean(w)
) %>%
pivot_longer(
cols = -id,
names_to = 'variable',
values_to = 'statistic'
) %>%
group_by(variable) %>%
mutate(average_value = mean(statistic)) %>%
ungroup() %>%
ggplot(
aes(
x = statistic,
fill = variable)
) +
colorblindr::scale_fill_OkabeIto() +
colorblindr::scale_color_OkabeIto() +
geom_vline(aes(xintercept = average_value), size = 1.25, linetype = 2) +
geom_density(alpha = 0.7, color = 'black', size = 2) +
facet_wrap(~variable, scales = 'free') +
labs(
y = element_blank(),
x = 'Distribution of Simulated Values (group averages)',
title = 'Figure 2. Distribution of Group Averages for\nSimulated Variables',
caption = 'Note the scales are free on each plot') +
theme_minimal() +
project_theme +
theme(axis.text.y = element_blank())
figure_2
show code
figure_3a <-
sim_dat %>%
ggplot(aes(x = x, y = outcome)) +
geom_point() +
geom_smooth(method = 'lm') +
theme_minimal() +
project_theme +
labs(
title = 'Figure 3. Association of Outcome & Each Covariate',
subtitle = '3a. Association of Outcome & x'
) +
theme(plot.subtitle = element_text(size = 15))
figure_3b <-
sim_dat %>%
ggplot(aes(x = z, y = outcome)) +
geom_point() +
geom_smooth(method = 'lm') +
theme_minimal() +
project_theme +
labs(subtitle = '3b. Association of Outcome & z') +
theme(plot.subtitle = element_text(size = 15))
figure_3c <-
sim_dat %>%
ggplot(aes(x = x*z, y = outcome)) +
geom_point() +
geom_smooth(method = 'lm') +
theme_minimal() +
project_theme +
labs(subtitle = '3c. Association of Outcome & Interaction') +
theme(plot.subtitle = element_text(size = 15))
ggpubr::ggarrange(figure_3a, figure_3b, figure_3c, ncol = 1)
Fitting multilevel model to the data
The simulation looks to have worked as planned. With larger samples, the effects becomes closer to what we specified, obviously, but that was omitted for brevity. Let’s see what this data looks like with an {lme4} to be sure.
show code
confirmation_lmer <-
lmer(
data = sim_dat,
formula = outcome ~ x * z + w + # equivalent to x + z + x : z
(1 | id)
)
summary(confirmation_lmer)
Linear mixed model fit by REML ['lmerMod']
Formula: outcome ~ x * z + w + (1 | id)
Data: sim_dat
REML criterion at convergence: 4729.8
Scaled residuals:
Min 1Q Median 3Q Max
-3.4943 -0.6663 -0.0260 0.6799 2.9690
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 0.7374 0.8587
Residual 4.5271 2.1277
Number of obs: 1065, groups: id, 71
Fixed effects:
Estimate Std. Error t value
(Intercept) 4.37559 0.12117 36.111
x 1.09797 0.07135 15.389
z 1.95442 0.06737 29.011
w 0.06962 0.06750 1.031
x:z 0.80677 0.07160 11.268
Correlation of Fixed Effects:
(Intr) x z w
x -0.053
z 0.012 -0.014
w 0.013 -0.024 -0.027
x:z -0.005 -0.049 -0.141 -0.003We see how closely the MLM estimates identify the simulated effects. We see small standard errors and large t-values for these estimates, too. As we’d like to see, w does not have significant associations with the outcome.