mixed_from_lmer.Rmdpamlj::pamlmixed() runs a power analysis for a linear or
generalized linear mixed model. Normally you describe the model by hand,
writing a syntax string in which every fixed effect
carries an expected coefficient and every random effect an expected
variance (see ?pamlmixed). This is convenient when you are
planning a study from scratch, but often you already have a fitted
lme4::lmer()/glmer() model — a pilot study, a
previous experiment, or published data you want to use as the basis for
planning a new, adequately powered study. Retyping that model as a
syntax string is tedious and error prone.
pamlmixed() can instead take the fitted model directly,
via the model = argument, and extract everything it needs
from it.
We use lme4’s built-in sleepstudy data:
reaction time measured on 18 subjects over 10 days of sleep deprivation,
with a random intercept and slope for Days by
Subject.
fit <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
summary(fit)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (Days | Subject)
#> Data: sleepstudy
#>
#> REML criterion at convergence: 1743.6
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.9536 -0.4634 0.0231 0.4634 5.1793
#>
#> Random effects:
#> Groups Name Variance Std.Dev. Corr
#> Subject (Intercept) 612.10 24.741
#> Days 35.07 5.922 0.07
#> Residual 654.94 25.592
#> Number of obs: 180, groups: Subject, 18
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 251.405 6.825 36.838
#> Days 10.467 1.546 6.771
#>
#> Correlation of Fixed Effects:
#> (Intr)
#> Days -0.138Passing the fitted model to pamlmixed() computes the
power of the design actually observed in the fit — the same number of
subjects, the same number of observations per subject, the same
fixed-effect coefficients, residual variance and random-effect
(co)variances:
res <- pamlmixed(model = fit, algo = "mc", mcR = 200, set_seed = TRUE, seed = 1)
#>
#> Simulations replications is set to 200. Please set a number greater than 1000 for more stable results.
#> ℹ pamlj: Finding power for n=10 and k=18
res$powertab$asDF
#> effect k n power df df_error sig.level converged singular
#> "1" Days 18 10 0.98 1 35.32932 0.05 0.99 0.265That’s the whole workflow: fit the model as you normally would, then
hand it to pamlmixed(). (mcR = 200 and
set_seed = TRUE are used here only to keep this vignette
fast and reproducible; for a real analysis use the default Monte Carlo
settings, or more repetitions, e.g. mcR = 2000.)
When model is supplied, pamlmixed() fills
in the arguments it would otherwise need by hand:
syntax — the fixed-effect coefficients
(fixef(fit)) and the random-effect variances (from
VarCorr(fit)), turned into pamlj’s model syntax.clusterpars — the number of levels and
cases per cluster for every grouping factor, read off the fit
(ngrps(fit), nobs(fit)).categorical — any factor predictor,
together with the exact contrast coding it was fitted with (so
the extracted coefficients stay correctly interpretable — see
below).sigma2 and
model_type — the residual variance and
whether the model is linear (lmer) or binomial logistic
(glmer).aim defaults to "power" —
computing the power of the observed design — instead of the usual
default "n".Any of these can still be overridden explicitly; an explicit argument
always wins over what was extracted from the model.
pamlmixed() also reports, via message(),
anything about the extraction worth knowing (e.g. a factor’s contrast
coding, or a random-effect covariance that couldn’t be preserved) —
suppress these with verbose = FALSE, as in the examples
above.
The fitted model tells you the power of the study you already
ran. Often the real question is different: how large would a
follow-up study need to be to reach a target power? Override
aim and power:
res <- pamlmixed(model = fit, aim = "n", power = .90,
algo = "mc", mcR = 200, set_seed = TRUE, seed = 1)
#>
#> With aim=n (find # of cases per cluster), the input number of cases are used as starting point. Clusters: Subject
#>
#> Simulations replications is set to 200. Please set a number greater than 1000 for more stable results.
#> ℹ pamlj: Finding number of cases within clusters
#> ◐ pamlj: INT_SEEK: trying n=10 with boundaries [0,1e+05] ◓ pamlj: INT_SEEK: obtained power =0.98894 target=0.9 next steps=-9 sd=3.18759026685619◑ pamlj: INT_SEEK: trying n=3 with boundaries [0,10] ◒ pamlj: INT_SEEK: obtained power =0.64219 target=0.9 next steps=26 sd=3.43278585904358◐ pamlj: INT_SEEK: trying n=9 with boundaries [3,10] ◓ pamlj: INT_SEEK: obtained power =0.9807 target=0.9 next steps=-8 sd=3.75312621183595◑ pamlj: INT_SEEK: trying n=4 with boundaries [3,9] ◒ pamlj: INT_SEEK: obtained power =0.80005 target=0.9 next steps=10 sd=4.05265992077361◐ pamlj: INT_SEEK: trying n=8 with boundaries [4,9] ◓ pamlj: INT_SEEK: obtained power =0.97124 target=0.9 next steps=-7 sd=4.28525853192959◑ pamlj: INT_SEEK: trying n=5 with boundaries [4,8] ◒ pamlj: INT_SEEK: obtained power =0.89476 target=0.9 next steps=1 sd=4.40142786380503◐ pamlj: INT_SEEK: trying n=6 with boundaries [5,8] ◓ pamlj: INT_SEEK: obtained power =0.92595 target=0.9 next steps=-3 sd=4.33020266231576✔ pamlj: INT_SEEK:closed
#>
#>
#> ✔ pamlj: Quick search found n=6 with exit:closed
#> ℹ pamlj: MC refinement starts from n=6 with step scale=6
#> ◐ pamlj: INT_SEEK: trying n=6 with boundaries [0,1e+05] ◓ pamlj: INT_SEEK: obtained power =0.86 target=0.9 next steps=0 sd=1.9090102147448✔ pamlj: INT_SEEK: nosteps
#>
#>
#> ✔ pamlj: MC search found n=6 with exit:nosteps
res$powertab$asDF
#> effect k n power tested df df_error sig.level converged singular
#> "1" Days 18 6 0.86 Tested 1 30.78603 0.05 0.995 0.34By default this searches for the number of cases per cluster
(find = "n", i.e. how many measurement occasions per
subject); pass find = "k" to search for the number of
clusters (subjects) instead.
A model can have several fixed effects, each with its own power; by
default pamlmixed() reports (and, for
aim = "n", targets) the worst-case one. To focus
the search on one specific effect instead, use focus =,
giving the term’s name exactly as it appears in the fitted model
(attr(terms(fit), "term.labels")):
res <- pamlmixed(model = fit, focus = "Days", aim = "n", power = .90,
algo = "mc", mcR = 200, set_seed = TRUE, seed = 1)
#>
#> With aim=n (find # of cases per cluster), the input number of cases are used as starting point. Clusters: Subject
#>
#> Simulations replications is set to 200. Please set a number greater than 1000 for more stable results.
#> ℹ pamlj: Finding number of cases within clusters
#> ◐ pamlj: INT_SEEK: trying n=10 with boundaries [0,1e+05] ◓ pamlj: INT_SEEK: obtained power =0.98894 target=0.9 next steps=-9 sd=3.18759026685619◑ pamlj: INT_SEEK: trying n=3 with boundaries [0,10] ◒ pamlj: INT_SEEK: obtained power =0.64219 target=0.9 next steps=26 sd=3.43278585904358◐ pamlj: INT_SEEK: trying n=9 with boundaries [3,10] ◓ pamlj: INT_SEEK: obtained power =0.9807 target=0.9 next steps=-8 sd=3.75312621183595◑ pamlj: INT_SEEK: trying n=4 with boundaries [3,9] ◒ pamlj: INT_SEEK: obtained power =0.80005 target=0.9 next steps=10 sd=4.05265992077361◐ pamlj: INT_SEEK: trying n=8 with boundaries [4,9] ◓ pamlj: INT_SEEK: obtained power =0.97124 target=0.9 next steps=-7 sd=4.28525853192959◑ pamlj: INT_SEEK: trying n=5 with boundaries [4,8] ◒ pamlj: INT_SEEK: obtained power =0.89476 target=0.9 next steps=1 sd=4.40142786380503◐ pamlj: INT_SEEK: trying n=6 with boundaries [5,8] ◓ pamlj: INT_SEEK: obtained power =0.92595 target=0.9 next steps=-3 sd=4.33020266231576✔ pamlj: INT_SEEK:closed
#>
#>
#> ✔ pamlj: Quick search found n=6 with exit:closed
#> ℹ pamlj: MC refinement starts from n=6 with step scale=6
#> ◐ pamlj: INT_SEEK: trying n=6 with boundaries [0,1e+05] ◓ pamlj: INT_SEEK: obtained power =0.86 target=0.9 next steps=0 sd=1.9090102147448✔ pamlj: INT_SEEK: nosteps
#>
#>
#> ✔ pamlj: MC search found n=6 with exit:nosteps
res$powertab$asDF
#> effect k n power tested df df_error sig.level converged singular
#> "1" Days 18 6 0.86 Tested 1 30.78603 0.05 0.995 0.34Factor predictors are handled automatically. Here we add a 3-level grouping factor to the fixed effects:
sleepstudy2 <- sleepstudy
sleepstudy2$grp <- factor(rep(c("A", "B", "C"), length.out = nrow(sleepstudy)))
fit2 <- lmer(Reaction ~ Days + grp + (Days | Subject), data = sleepstudy2)
res <- pamlmixed(model = fit2, algo = "mc", mcR = 200, set_seed = TRUE, seed = 1)
#> Categorical predictor(s) grp are simulated with the exact contrast coding read from the fitted model's.
#>
#> Simulations replications is set to 200. Please set a number greater than 1000 for more stable results.
#> ℹ pamlj: Finding power for n=10 and k=18
res$powertab$asDF
#> effect k n power df df_error sig.level converged singular
#> "1" Days 18 10 0.97 1 35.64116 0.05 0.99 0.27
#> 2 grp 18 10 0.09 2 157.69289 0.05 0.99 0.27pamlmixed() reads the exact contrast matrix
grp was fitted with directly off the model’s data, rather
than guessing or re-coding it, so the extracted coefficients keep
meaning the same thing they did in the original fit.
Designs with two crossed grouping factors — for instance subjects and
items in a repeated-measures experiment,
(1 | subject) + (1 | item) — are supported as well.
pamlmixed() reproduces the random-effect covariance between
subject/item terms when the fit has one, and, importantly, checks
whether a fixed-effect predictor is constant within one of the
two grouping factors (e.g. a property of the item that never changes
across subjects). When it is, it is simulated as a
between-cluster effect, so the simulated design does
not need to needlessly balloon in size to represent it. A predictor that
genuinely varies within every combination of subject and item
(a true trial-level, counterbalanced factor) has no such shortcut, and
still needs enough real replication in the original data to be simulated
at all.
set.seed(1)
nsubj <- 20; nitem <- 20
cond_per_subject <- factor(sample(c("control", "treatment"), nsubj, replace = TRUE))
d <- expand.grid(subject = factor(1:nsubj), item = factor(1:nitem))
d$cond <- cond_per_subject[as.integer(d$subject)]
d$y <- rnorm(nrow(d), mean = ifelse(d$cond == "treatment", .3, 0))
fit3 <- lmer(y ~ cond + (1 | subject) + (1 | item), data = d)
#> boundary (singular) fit: see help('isSingular')
res <- pamlmixed(model = fit3, algo = "mc", mcR = 200, set_seed = TRUE, seed = 1)
#> Categorical predictor(s) cond are simulated with the exact contrast coding read from the fitted model's.
#> Predictor(s) cond (constant within subject)
#>
#> Simulations replications is set to 200. Please set a number greater than 1000 for more stable results.
#> ℹ pamlj: Finding power for n=1 and k=20
res$powertab$asDF
#> effect k n power df df_error sig.level converged singular
#> "1" cond 20 1 0.85 1 136.4643 0.05 0.995 0.71lmer) and binomial (glmer)
models are supported.For the full set of options — correlated random effects, custom
contrast coding, etc. — see ?pamlmixed.