Week 7 (R): Power
There have been some updates to the R code since the videos were recorded. The major changes are:
The code now uses the packaged
simr::powerCurve()
function, instead of the functions I previously made, to obtain power estimates across different levels of sample sizesFor some of the examples, I started with a bigger number of clusters and cluster size, which may give more stable results for continuous predictors.
Load Packages
# To install a package, run the following ONCE (and only once on your computer)
# install.packages("simr")
library(tidyverse)
library(lme4)
library(simr)
library(MuMIn) # you can use `performance::r2()` as an alternative
See also an example from the vignette of the simr
package: https://cran.r-project.org/web/packages/simr/vignettes/fromscratch.html
Note: In the examples of this note I used 100 simulation samples only to save time. In practice, you should use at least 1,000 to obtain a reasonable precision on the power estimation.
Binary Predictor at Level 2
\[ \begin{aligned} Y_{ij} & = \beta_{0j} + e_{ij} \\ \beta_{0j} & = \gamma_{00} + {\color{red}{\gamma_{01}}} T_j + u_{0j} \end{aligned} \]
Preparing simulated data
First, generate the predictor matrix and the cluster ID
# Cluster ID (start with 10 clusters)
<- 10
num_clus <- 1:num_clus
id # Binary lv-2 predictor
<- factor(c("c", "t"))
treat # Level-2 data (note: `tibble` is the tidyverse version of `data.frame`)
<- tibble(id, treat = rep_len(treat, length.out = num_clus))
lv2_dat # Cluster size
<- 5
clus_size # Expand each cluster to include more rows
<- lv2_dat %>%
(sim_dat # seq_len(n()) creates a range from 1 to number of rows
slice(rep(seq_len(n()), each = clus_size)))
Then specify the fixed effect coefficients and the random effect variance:
<- c(0, 0.625) # gamma_00, gamma_01
gams <- matrix(0.25) # conditional ICC = 0.25 / (1 + 0.25) = 0.2
taus <- 1 # sigma, usually fixed to be 1.0 for easier scaling sigma
Use simr::makeLmer()
to build an artificial lmer
object:
<- makeLmer(y ~ treat + (1 | id),
sim_m1 fixef = gams, VarCorr = taus, sigma = sigma,
data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m1)
># Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
># R2m R2c
># [1,] 0.07383343 0.2590667
# Cohen's d
2] / sqrt(taus + sigma^2) # treatment effect / sqrt(tau_0^2 + sigma^2) gams[
># [,1]
># [1,] 0.559017
Obtain simulation-based power
Now, check the power:
<- powerSim(sim_m1,
power_m1 # number of simulated samples; use a small number for
# trial run
nsim = 20,
progress = FALSE) # suppress the progress bar
power_m1
># Power for predictor 'treat', (95% confidence interval):
># 25.00% ( 8.66, 49.10)
>#
># Test: Likelihood ratio
>#
># Based on 20 simulations, (0 warnings, 0 errors)
># alpha = 0.05, nrow = 50
>#
># Time elapsed: 0 h 0 m 1 s
Power was obviously not sufficient as it is far from 80%, so we’ll increase it.
Increasing number of clusters
Increase to 50 clusters
(Note that I set the cache = TRUE
option so that I can reuse the results next time.)
<- extend(sim_m1,
sim_m1_50J along = "id", # add more levels of `id` (i.e., clusters)
n = 50) # increase to 50 clusters
<- powerSim(sim_m1_50J, nsim = 20, progress = FALSE)
power_m1_50J power_m1_50J
># Power for predictor 'treat', (95% confidence interval):
># 95.00% (75.13, 99.87)
>#
># Test: Likelihood ratio
>#
># Based on 20 simulations, (0 warnings, 0 errors)
># alpha = 0.05, nrow = 250
>#
># Time elapsed: 0 h 0 m 1 s
This looks close to 80%. Let’s try several numbers of clusters (10, 20, 30, 40, 50) with a larger nsim
. You can use the simr::powerCurve()
function.
# The first argument needs to be the output from running `extend`
<- powerCurve(
powers_m1
sim_m1_50J,along = "id",
breaks = c(10, 20, 30, 40, 50),
nsim = 100,
progress = FALSE # suppress progress bar (for printing purpose)
)# summary; check the numbers of rows are correct powers_m1
># Power for predictor 'treat', (95% confidence interval),
># by largest value of id:
># 10: 37.00% (27.56, 47.24) - 50 rows
># 20: 51.00% (40.80, 61.14) - 100 rows
># 30: 65.00% (54.82, 74.27) - 150 rows
># 40: 80.00% (70.82, 87.33) - 200 rows
># 50: 90.00% (82.38, 95.10) - 250 rows
>#
># Time elapsed: 0 h 0 m 33 s
plot(powers_m1)
This suggests that ~ 50 clusters are needed, if each cluster contains 5 observations.
Increasing cluster size
Instead of increasing the number of cluster, one can also increase the cluster size (with \(J\) = 50):
<- extend(sim_m1,
sim_m1_25n # need to specify the combination of id and lv-2 predictors
within = "id + treat",
n = 25) # cluster size
<- powerSim(sim_m1_25n, nsim = 100, progress = FALSE)
power_m1_25n power_m1_25n
># Power for predictor 'treat', (95% confidence interval):
># 52.00% (41.78, 62.10)
>#
># Test: Likelihood ratio
>#
># Based on 100 simulations, (0 warnings, 0 errors)
># alpha = 0.05, nrow = 250
>#
># Time elapsed: 0 h 0 m 8 s
Note that the power of increasing cluster size by 5 times is much smaller than that of increasing number of clusters by 5 times. Here is the power curve:
<- powerCurve(
powers_m1_n
sim_m1_25n,within = "id + treat",
breaks = c(5, 10, 15, 20, 25),
nsim = 100,
progress = FALSE
) powers_m1_n
># Power for predictor 'treat', (95% confidence interval),
># by number of observations within id + treat:
># 5: 37.00% (27.56, 47.24) - 50 rows
># 10: 44.00% (34.08, 54.28) - 100 rows
># 15: 49.00% (38.86, 59.20) - 150 rows
># 20: 49.00% (38.86, 59.20) - 200 rows
># 25: 48.00% (37.90, 58.22) - 250 rows
>#
># Time elapsed: 0 h 0 m 35 s
plot(powers_m1_n)
In this case, increasing the cluster size only increases the power to close to 50%. There is an upper limit in power improvement with increasing the cluster size, when the number of cluster is kept constant.
Changing effect size
Finally, we can also explore how effect size affects power
# Change gamma_01
<- c(0, 0.9) # gamma_00, gamma_01
gams <- matrix(0.25) # conditional ICC = 0.25 / (1 + 0.25) = 0.2
taus <- 1 # sigma, usually fixed to be 1.0 for easier scaling
sigma # Use `makeLmer()` again
<- makeLmer(y ~ treat + (1 | id),
sim_m1_large_es fixef = gams, VarCorr = taus, sigma = sigma,
data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m1_large_es)
># R2m R2c
># [1,] 0.1418564 0.3134851
# Cohen's d
2] / sqrt(taus + sigma^2) # treatment effect / sqrt(tau_0^2 + sigma^2) gams[
># [,1]
># [1,] 0.8049845
Let’s also look at varying number of clusters
<- extend(sim_m1_large_es, along = "id", n = 40)
sim_m1_large_es_40J <- powerCurve(
power_m1_large_es along = "id", breaks = c(10, 20, 30, 40),
sim_m1_large_es_40J, nsim = 100, progress = FALSE
) power_m1_large_es
># Power for predictor 'treat', (95% confidence interval),
># by largest value of id:
># 10: 67.00% (56.88, 76.08) - 50 rows
># 20: 86.00% (77.63, 92.13) - 100 rows
># 30: 96.00% (90.07, 98.90) - 150 rows
># 40: 99.00% (94.55, 99.97) - 200 rows
>#
># Time elapsed: 0 h 0 m 28 s
plot(power_m1_large_es)
So, to summarize power for the different scenarios, if we want 80% power, we need
- Treatment effect = 0.625, ~ 50 clusters
- Treatment effect = 0.90, ~ 25 clusters
Continuous Predictor at Level 2
\[ \begin{aligned} Y_{ij} & = \beta_{0j} + e_{ij} \\ \beta_{0j} & = \gamma_{00} + {\color{red}{\gamma_{01}}} W_j + u_{0j} \end{aligned} \]
Preparing Simulated Data
In the previous section, we started with a small simulated data. When we have continuous predictor(s), it is usually best to start with a larger number of clusters.
First, generate the predictor matrix and the cluster ID
# Cluster ID (use 100 clusters)
<- 100
num_clus <- 1:num_clus
id # Continuous lv-2 predictor
<- rnorm(num_clus)
w # Convert to z score so that SD = 1
<- (w - mean(w)) / sd(w) * (num_clus - 1) / num_clus
w # Level-2 data
<- tibble(id, w = w)
lv2_dat # Cluster size
<- 25
clus_size # Expand each cluster to include more rows
<- lv2_dat %>%
(sim_dat # seq_len(n()) creates a range from 1 to number of rows
slice(rep(seq_len(n()), each = clus_size)))
Then specify the fixed effect coefficients and the random effect variance:
<- c(1, 0.3) # gamma_00, gamma_01
gams <- matrix(0.5) # conditional ICC = 0.5 / (1 + 0.5) = 0.33
taus <- 1 # sigma, usually fixed to be 1.0 for easier scaling sigma
Use simr::makeLmer()
to build an artificial lmer
object:
<- makeLmer(y ~ w + (1 | id),
sim_m2 fixef = gams, VarCorr = taus, sigma = sigma,
data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m2)
># R2m R2c
># [1,] 0.05503588 0.3700239
Obtain simulation-based power
Use various numbers of clusters
<- powerCurve(
powers_m2 along = "id", breaks = c(10, 20, 50, 100),
sim_m2, nsim = 100, progress = FALSE
) powers_m2
># Power for predictor 'w', (95% confidence interval),
># by largest value of id:
># 10: 26.00% (17.74, 35.73) - 250 rows
># 20: 37.00% (27.56, 47.24) - 500 rows
># 50: 74.00% (64.27, 82.26) - 1250 rows
># 100: 99.00% (94.55, 99.97) - 2500 rows
>#
># Time elapsed: 0 h 4 m 54 s
plot(powers_m2)
Use various cluster sizes (negligible change in power)
<- powerCurve(
powers_m2_n within = "id + w", breaks = c(5, 10, 25, 50),
sim_m2, nsim = 100, progress = FALSE
) powers_m2_n
># Power for predictor 'w', (95% confidence interval),
># by number of observations within id + w:
># 5: 94.00% (87.40, 97.77) - 500 rows
># 10: 98.00% (92.96, 99.76) - 1000 rows
># 25: 99.00% (94.55, 99.97) - 2500 rows
># NA: 99.00% (94.55, 99.97) - 2500 rows
>#
># Time elapsed: 0 h 6 m 19 s
plot(powers_m2_n)
Binary Predictor at Level 2 With a Continuous Level-1 Covariate (and Random Slope)
Cluster-mean centering needs to be considerd
\[ \begin{aligned} Y_{ij} & = \beta_{0j} + \beta_{1j} X\text{\_cmc}_{ij} + e_{ij} \\ \beta_{0j} & = \gamma_{00} + \gamma_{01} W_j + \gamma_{02} X\text{\_cm}_{j} + u_{0j} \\ \beta_{1j} & = \gamma_{10} + {\color{red}{\gamma_{11}}} W_j + u_{1j} \end{aligned} \]
We consider a covariate with ICC = .2
Preparing Simulated Data
First, generate the simulated data. I use a large \(J\) and a moderate \(n\) here.
# Step 1: Lv-2 data (basically the same as before)
# Cluster ID (start with 10 clusters)
<- 100
num_clus <- 1:num_clus
id # Binary lv-2 predictor
<- factor(c("c", "t"))
treat # Cluster means of X, with ICC(X) = 0.2
<- rnorm(num_clus, mean = 0, sd = sqrt(0.2))
x_cm # Force SD = sqrt(0.2)
<- (x_cm - mean(x_cm)) / sd(x_cm) * sqrt(0.2) * (num_clus - 1) / num_clus
x_cm # Level-2 data
<- tibble(id, treat = rep_len(treat, length.out = num_clus), x_cm)
lv2_dat # Expand each cluster to include more rows
<- 40 # Cluster size
clus_size <- lv2_dat %>%
sim_dat slice(rep(seq_len(n()), each = clus_size))
# Step 2: Add level-1 covariate
# Within-cluster component of X, with sigma^2(X) = 0.8
<- num_clus * clus_size
num_obs <- rnorm(num_obs, mean = 0, sd = 1)
x_cmc <- x_cmc - ave(x_cmc, lv2_dat$id, FUN = mean)
x_cmc # Force SD = sqrt(0.8)
<- x_cmc / sd(x_cmc) * sqrt(0.8) * (num_obs - 1) / num_obs
x_cmc # Add `x_cmc` to `sim_dat`
$x_cmc <- x_cmc sim_dat
Then specify the fixed effect coefficients and the random effect variance:
# gamma_00, gamma_01, gamma_02, gamma_10, gamma_11
<- c(0, 0.3, 0.2, 0.1, 0.2)
gams <- matrix(c(0.25, 0,
taus 0, 0.10), nrow = 2) # tau^2_0 = 0.25, tau^2_1 = 0.1
<- sqrt(.96) # sigma sigma
Use simr::makeLmer()
to build an artificial lmer
object:
<- makeLmer(y ~ treat + x_cm + x_cmc + treat:x_cmc + (x_cmc | id),
sim_m4 fixef = gams, VarCorr = taus, sigma = sigma,
data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m4)
># R2m R2c
># [1,] 0.0528709 0.295127
# Compare to a model without the cross-level interaction
<- makeLmer(y ~ treat + x_cm + x_cmc + (x_cmc | id),
sim_m4null fixef = gams[1:4], VarCorr = taus, sigma = sigma,
data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m4null)
># R2m R2c
># [1,] 0.02919444 0.2775064
# So the cross-level interaction corresponds to an R^2 change of about 2.5%
Obtain simulation-based power
Use various numbers of clusters
(Note: by default the first term in the fixed effect [other than the intercept] will be tested. We’ll specify the interaction instead.)
<- powerCurve(
powers_m4 along = "id", breaks = c(10, 20, 50, 100),
sim_m4, # `method = "anova"` means power for LRT
test = fixed("treat:x_cmc", method = "anova"),
nsim = 100, progress = FALSE
) powers_m4
># Power for predictor 'treat:x_cmc', (95% confidence interval),
># by largest value of id:
># 10: 11.00% ( 5.62, 18.83) - 400 rows
># 20: 20.00% (12.67, 29.18) - 800 rows
># 50: 50.00% (39.83, 60.17) - 2000 rows
># 100: 80.00% (70.82, 87.33) - 4000 rows
>#
># Time elapsed: 0 h 1 m 11 s
plot(powers_m4)
Use various cluster sizes (results not shown)
# Increase to n = 200
<- extend(sim_m4, within = "id + treat + x_cm", n = 200)
sim_m4_200n # Decrease to J = 50
<- extend(sim_m4_200n, along = "id", n = 50)
sim_m4_200n <- powerCurve(
powers_m4_n within = "id + treat + x_cm", breaks = c(40, 80, 200),
sim_m4_200n, # `method = "anova"` means power for LRT
test = fixed("treat:x_cmc", method = "anova"),
nsim = 100, progress = FALSE
)
powers_m4_nplot(powers_m4_n)
Bonus: Growth Modeling
With a binary person-level predictor
\[ \begin{aligned} Y_{ti} & = \beta_{0j} + \beta_{1i} \text{time}_{ti} + e_{ti} \\ \beta_{0i} & = \gamma_{00} + {\color{red}{\gamma_{01}}} W_i + u_{0i} \\ \beta_{1i} & = \gamma_{10} + \gamma_{11} W_i + u_{1i} \end{aligned} \]
Preparing Simulated Data
First, generate the predictor matrix and the cluster ID
# Person ID (start with 50 clusters)
<- 50
num_clus <- 1:num_clus
id # Binary lv-2 predictor
<- factor(c("c", "t"))
treat # Level-2 data
<- tibble(id, treat = rep_len(treat, length.out = num_clus))
lv2_dat # Expand each cluster to include more rows
<- 4 # Number of time points
clus_size # Within-cluster component of X, with sigma^2(X) = 0.8
<- 1:clus_size
time # Expand each cluster to include more rows
<- expand_grid(lv2_dat, time)) (sim_dat
Then specify the fixed effect coefficients and the random effect variance:
# gamma_00, gamma_01, gamma_10, gamma_11
<- c(1, 0, -0.2, 0.3)
gams <- matrix(c(0.5, 0,
taus 0, 0.10), nrow = 2) # tau^2_0 = 0.5, tau^2_1 = 0.1
<- sqrt(1) # sigma sigma
Use simr::makeLmer()
to build an artificial lmer
object:
# I use `I(time - 1)` so that the intercept is for the first time point
<- makeLmer(
sim_m5 ~ treat + I(time - 1) + treat:I(time - 1) + (I(time - 1) | id),
y fixef = gams, VarCorr = taus, sigma = sigma,
data = sim_dat
)# Check R^2
r.squaredGLMM(sim_m5)
># R2m R2c
># [1,] 0.04258501 0.4824784
# Compare to a model without the binary treatment predictor
<- makeLmer(y ~ I(time - 1) + (I(time - 1) | id),
sim_m5null fixef = gams[c(1, 3)], VarCorr = taus, sigma = sigma,
data = sim_dat)
# Check R^2
r.squaredGLMM(sim_m5null)
># R2m R2c
># [1,] 0.02644453 0.4737538
# So treat and its interaction corresponds to an R^2 change of about 2%
Obtain simulation-based power
Use various numbers of clusters. The default is to test the first fixed effect (i.e., treat
, group difference at the initial time point) using a likelihood ratio test. The true value of this coefficient, so the estimated power is actually Type I error in this case.
# Increase to 100 participants
<- extend(sim_m5, along = "id", n = 100)
sim_m5_100J <- powerCurve(
powers_m5 along = "id", breaks = c(10, 20, 50, 100),
sim_m5_100J, nsim = 100, progress = FALSE
) powers_m5
># Power for predictor 'treat', (95% confidence interval),
># by largest value of id:
># 10: 3.00% ( 0.62, 8.52) - 40 rows
># 20: 5.00% ( 1.64, 11.28) - 80 rows
># 50: 4.00% ( 1.10, 9.93) - 200 rows
># 100: 6.00% ( 2.23, 12.60) - 400 rows
>#
># Time elapsed: 0 h 0 m 39 s
plot(powers_m5)
Now test the slope difference (i.e., treat:time
interaction)
<- powerCurve(
powers_m5_int along = "id", breaks = c(10, 20, 50, 100),
sim_m5_100J, test = fixed("treat:I(time - 1)", method = "anova"),
nsim = 100, progress = FALSE
)plot(powers_m5_int)
Use more time points
# Increase to 8 time points
<- extend(sim_m5, along = "time", values = 1:8)
sim_m5_50n <- powerCurve(
powers_m5_n along = "time", breaks = c(4, 6, 8),
sim_m5_50n, nsim = 100, progress = FALSE,
test = fixed("treat:I(time - 1)", method = "anova")
) powers_m5_n
># Power for predictor 'treat:I(time - 1)', (95% confidence interval),
># by largest value of time:
># 4: 49.00% (38.86, 59.20) - 200 rows
># 6: 82.00% (73.05, 88.97) - 300 rows
># 8: 85.00% (76.47, 91.35) - 400 rows
>#
># Time elapsed: 0 h 0 m 25 s
plot(powers_m5_n)
Increasing number of time points helps a little bit for testing the slope difference.