emmeans
analyseslibrary("Superpower")
library("emmeans")
When conducting exact ANOVA power analyses with
Superpower
it is possible to calculate the power for both
the omnibus \(F\)-tests and planned
contrasts or post-hoc comparisons. It is possible to use this approach
to calculate power for standard contrasts, such as pairwise contrasts
between cells. Here I provide a brief overview of how the output of
ANOVA_exact()
can be used to perform power analyses for
tailored planned contrasts and follow-up tests.
All power analyses for emmeans
-objects are based on the
\(F\)- and \(t\)-values from the analyses of the dataset
simulated by ANOVA_exact()
assuming two-sided testing.
Thus, the emmeans_power()
does not honor adjustments of the
testing procedure due to either one-sided testing (including two
one-sided tests) or corrections for multiple comparisons via the
adjust
option in emmeans
. As noted below, for
the Bonferroni-adjustment this limitation can be overcome by manually
adjusting alpha_level
.
First, we will set up a 2 \(\times\)
3 repeated measures design. When calling ANOVA_exact()
pairwise comparisons of expected marginal means are added by setting
emm = TRUE
and contrast_type = "pairwise"
(default).
# Set up a within design with 2 factors, each with 2 and 3 levels
<- ANOVA_design(
design_result design = "2w*3w",
n = 40,
mu = c(0.3, 0, 0.5, 0.3, 0, 0),
sd = 2,
r = 0.8,
label_list = list("condition" = c("cheerful", "sad"),
"voice" = c("human", "robot", "cartoon"))
)
<- ANOVA_exact(
exact_result
design_result,alpha_level = 0.05,
verbose = FALSE,
emm = TRUE,
contrast_type = "pairwise"
)
The result contains the power calculations for both the omnibus \(F\)-tests and pairwise post-hoc comparisons.
$main_results exact_result
power | partial_eta_squared | cohen_f | non_centrality | |
---|---|---|---|---|
condition | 29.08380 | 0.0507099 | 0.2311251 | 2.083333 |
voice | 50.12584 | 0.0621242 | 0.2573700 | 5.166667 |
condition:voice | 41.62968 | 0.0507099 | 0.2311251 | 4.166667 |
head(exact_result$emm_results)
contrast | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|
condition_cheerful voice_cartoon - condition_sad voice_cartoon | 68.37160 | 0.1381215 | 0.4003204 | 6.25 |
condition_cheerful voice_cartoon - condition_cheerful voice_human | 16.41134 | 0.0250000 | 0.1601282 | 1.00 |
condition_cheerful voice_cartoon - condition_sad voice_human | 16.41134 | 0.0250000 | 0.1601282 | 1.00 |
condition_cheerful voice_cartoon - condition_cheerful voice_robot | 68.37160 | 0.1381215 | 0.4003204 | 6.25 |
condition_cheerful voice_cartoon - condition_sad voice_robot | 68.37160 | 0.1381215 | 0.4003204 | 6.25 |
condition_sad voice_cartoon - condition_cheerful voice_human | 30.99652 | 0.0545455 | 0.2401922 | 2.25 |
The output also contains the emmeans
-object on which
these power calculations are based. By manipulating this object it is
possible to tailor the power analyses to the contrasts desired for the
planned study. That is, based on the dataset simulated with
ANOVA_exact()
we can write out the analysis code for
emmeans
-contrasts, just as we would if we were to analyze
the empirical data, and use the output to perform the corresponding
power analysis.
emmeans
contrastsThe emmeans
reference grid and contrasts are included in
the output of ANOVA_exact()
.
::kable(exact_result$emmeans$emmeans) knitr
condition | voice | emmean | SE | df | lower.CL | upper.CL |
---|---|---|---|---|---|---|
condition_cheerful | voice_cartoon | 0.5 | 0.3162278 | 39 | -0.139631 | 1.139631 |
condition_sad | voice_cartoon | 0.0 | 0.3162278 | 39 | -0.639631 | 0.639631 |
condition_cheerful | voice_human | 0.3 | 0.3162278 | 39 | -0.339631 | 0.939631 |
condition_sad | voice_human | 0.3 | 0.3162278 | 39 | -0.339631 | 0.939631 |
condition_cheerful | voice_robot | 0.0 | 0.3162278 | 39 | -0.639631 | 0.639631 |
condition_sad | voice_robot | 0.0 | 0.3162278 | 39 | -0.639631 | 0.639631 |
::kable(exact_result$emmeans$contrasts) knitr
contrast | estimate | SE | df | t.ratio | p.value |
---|---|---|---|---|---|
condition_cheerful voice_cartoon - condition_sad voice_cartoon | 0.5 | 0.2 | 39 | 2.5 | 0.0167338 |
condition_cheerful voice_cartoon - condition_cheerful voice_human | 0.2 | 0.2 | 39 | 1.0 | 0.3234749 |
condition_cheerful voice_cartoon - condition_sad voice_human | 0.2 | 0.2 | 39 | 1.0 | 0.3234749 |
condition_cheerful voice_cartoon - condition_cheerful voice_robot | 0.5 | 0.2 | 39 | 2.5 | 0.0167338 |
condition_cheerful voice_cartoon - condition_sad voice_robot | 0.5 | 0.2 | 39 | 2.5 | 0.0167338 |
condition_sad voice_cartoon - condition_cheerful voice_human | -0.3 | 0.2 | 39 | -1.5 | 0.1416670 |
condition_sad voice_cartoon - condition_sad voice_human | -0.3 | 0.2 | 39 | -1.5 | 0.1416670 |
condition_sad voice_cartoon - condition_cheerful voice_robot | 0.0 | 0.2 | 39 | 0.0 | 1.0000000 |
condition_sad voice_cartoon - condition_sad voice_robot | 0.0 | 0.2 | 39 | 0.0 | 1.0000000 |
condition_cheerful voice_human - condition_sad voice_human | 0.0 | 0.2 | 39 | 0.0 | 1.0000000 |
condition_cheerful voice_human - condition_cheerful voice_robot | 0.3 | 0.2 | 39 | 1.5 | 0.1416670 |
condition_cheerful voice_human - condition_sad voice_robot | 0.3 | 0.2 | 39 | 1.5 | 0.1416670 |
condition_sad voice_human - condition_cheerful voice_robot | 0.3 | 0.2 | 39 | 1.5 | 0.1416670 |
condition_sad voice_human - condition_sad voice_robot | 0.3 | 0.2 | 39 | 1.5 | 0.1416670 |
condition_cheerful voice_robot - condition_sad voice_robot | 0.0 | 0.2 | 39 | 0.0 | 1.0000000 |
By using emmeans_power()
on the contrasts, we can
reproduce the results of the previous power analysis for the pairwise
comparisons.
head(emmeans_power(exact_result$emmeans$contrasts))
contrast | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|
condition_cheerful voice_cartoon - condition_sad voice_cartoon | 68.37160 | 0.1381215 | 0.4003204 | 6.25 |
condition_cheerful voice_cartoon - condition_cheerful voice_human | 16.41134 | 0.0250000 | 0.1601282 | 1.00 |
condition_cheerful voice_cartoon - condition_sad voice_human | 16.41134 | 0.0250000 | 0.1601282 | 1.00 |
condition_cheerful voice_cartoon - condition_cheerful voice_robot | 68.37160 | 0.1381215 | 0.4003204 | 6.25 |
condition_cheerful voice_cartoon - condition_sad voice_robot | 68.37160 | 0.1381215 | 0.4003204 | 6.25 |
condition_sad voice_cartoon - condition_cheerful voice_human | 30.99652 | 0.0545455 | 0.2401922 | 2.25 |
Now, we can manipulate the emmeans
reference grid to
perform additional power analyses. In the following example, we
calculate the power for the contrasts between sad and cheerful condition
for each voice.
<- emmeans(
simple_condition_effects $emmeans$emmeans,
exact_resultspecs = ~ condition | voice
)
emmeans_power(pairs(simple_condition_effects))
contrast | voice | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|---|
condition_cheerful - condition_sad | voice_cartoon | 68.3716 | 0.1381215 | 0.4003204 | 6.25 |
condition_cheerful - condition_sad | voice_human | 5.0000 | 0.0000000 | 0.0000000 | 0.00 |
condition_cheerful - condition_sad | voice_robot | 5.0000 | 0.0000000 | 0.0000000 | 0.00 |
We may also calculate the power for testing all condition means against an arbitrary constant.
emmeans_power(test(simple_condition_effects, null = 0.5))
condition | voice | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|---|
condition_cheerful | voice_cartoon | 5.000000 | 0.0000000 | 0.0000000 | 0.0 |
condition_sad | voice_cartoon | 33.831141 | 0.0602410 | 0.2531848 | 2.5 |
condition_cheerful | voice_human | 9.462673 | 0.0101523 | 0.1012739 | 0.4 |
condition_sad | voice_human | 9.462673 | 0.0101523 | 0.1012739 | 0.4 |
condition_cheerful | voice_robot | 33.831141 | 0.0602410 | 0.2531848 | 2.5 |
condition_sad | voice_robot | 33.831141 | 0.0602410 | 0.2531848 | 2.5 |
Finally, we can calculate the power for custom contrasts between any linear combination of conditions.
<- contrast(
custom_contrast $emmeans$emmeans,
exact_resultlist(robot_vs_sad_human = c(0, 0, 0, 1, -0.5, -0.5))
)
emmeans_power(custom_contrast)
contrast | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|
robot_vs_sad_human | 39.34971 | 0.0714286 | 0.2773501 | 3 |
Although emmeans_power()
currently ignores adjustments
for multiple comparisons, it is possible to calculate the power for
Bonferroni-corrected tests by adjusting alpha_level
.
<- nrow(as.data.frame(simple_condition_effects))
n_contrasts
emmeans_power(
pairs(simple_condition_effects),
alpha_level = 0.05 / n_contrasts
)
contrast | voice | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|---|
condition_cheerful - condition_sad | voice_cartoon | 40.1555499 | 0.1381215 | 0.4003204 | 6.25 |
condition_cheerful - condition_sad | voice_human | 0.8333333 | 0.0000000 | 0.0000000 | 0.00 |
condition_cheerful - condition_sad | voice_robot | 0.8333333 | 0.0000000 | 0.0000000 | 0.00 |
Similarly, if we want to calculate power for a one-sided test, we can
doubling alpha_level
.
emmeans_power(
pairs(simple_condition_effects)[1],
alpha_level = 2 * 0.05
)
contrast | voice | power | partial_eta_squared | cohen_f | non_centrality |
---|---|---|---|---|---|
condition_cheerful - condition_sad | voice_cartoon | 79.14507 | 0.1381215 | 0.4003204 | 6.25 |
Note, that because power is calculated from the squared \(t\)-value, power is only calculated correctly if the alternative hypothesis is true in the simulated dataset. That is, the difference of the condition means is consistent with the tested directional hypothesis.
Because emmeans
can perform equivalence,
non-superiority, and -inferiority tests, emmeans_power()
can calculate the corresponding power for these tests.
emmeans_power(
pairs(simple_condition_effects, side = "equivalence", delta = 0.3)[2]
)
contrast | voice | power | partial_eta_squared | cohen_f | non_centrality | |
---|---|---|---|---|---|---|
2 | condition_cheerful - condition_sad | voice_human | 30.99652 | 0.0545455 | 0.2401922 | 2.25 |
Note, that because power is calculated from the squared \(t\)-value, power is only calculated
correctly if the alternative hypothesis is true in the simulated
dataset. That is, the difference between the condition means is
consistent with the tested directional hypothesis (smaller than
delta
).
Another useful application of emmeans_power()
is to
joint tests. Lets assume we plan to test the main effect of voice for
each of the two conditions separately using joint_tests()
.
We can then calculate the power for each Bonferroni-corrected \(F\)-test as follows.
<- joint_tests(
voice_by_condition $emmeans$emmeans,
exact_resultby = "condition"
)
emmeans_power(voice_by_condition, alpha_level = 0.05 / 2)
model term | condition | power | partial_eta_squared | cohen_f | non_centrality | |
---|---|---|---|---|---|---|
1 | voice | condition_cheerful | 21.80102 | 0.0751061 | 0.2849651 | 3.167 |
3 | voice | condition_sad | 10.39583 | 0.0370370 | 0.1961161 | 1.500 |