Performs pairwise comparisons between groups using the estimated
marginal means. Pipe-friendly wrapper arround the functions emmans() +
contrast() from the emmeans package, which need to be installed
before using this function. This function is useful for performing post-hoc
analyses following ANOVA/ANCOVA tests.
emmeans_test(
data,
formula,
covariate = NULL,
ref.group = NULL,
comparisons = NULL,
p.adjust.method = "bonferroni",
conf.level = 0.95,
model = NULL,
detailed = FALSE
)
get_emmeans(emmeans.test)a data.frame containing the variables in the formula.
a formula of the form x ~ group where x is a
numeric variable giving the data values and group is a factor with
one or multiple levels giving the corresponding groups. For example,
formula = TP53 ~ cancer_group.
(optional) covariate names (for ANCOVA)
a character string specifying the reference group. If specified, for a given grouping variable, each of the group levels will be compared to the reference group (i.e. control group).
If ref.group = "all", pairwise two sample tests are performed for
comparing each grouping variable levels against all (i.e. basemean).
A list of length-2 vectors specifying the groups of
interest to be compared. For example to compare groups "A" vs "B" and "B" vs
"C", the argument is as follow: comparisons = list(c("A", "B"), c("B",
"C"))
method to adjust p values for multiple comparisons. Used when pairwise comparisons are performed. Allowed values include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none". If you don't want to adjust the p value (not recommended), use p.adjust.method = "none".
confidence level of the interval.
a fitted-model object such as the result of a call to
lm(), stats::aov() (including a within-subject Error()
term) or nlme::lme(), on which the estimated marginal means are
computed. Supplying model lets you (i) average over factors that are
in the model but not in formula (e.g. for a factorial design, fit
lm(y ~ a * b) and compare a with formula = y ~ a,
averaging over b), and (ii) run pairwise comparisons for
repeated-measures / mixed designs (pass a within-subject model). When
model = NULL (default), a simple lm() is fitted from
formula (plus the grouping and covariate variables). See examples.
logical value. Default is FALSE. If TRUE, a detailed result is shown.
an object of class emmeans_test.
return a data frame with some the following columns:
.y.: the y variable used in the test.
group1,group2: the
compared groups in the pairwise tests.
statistic: Test
statistic (t.ratio) used to compute the p-value.
df: degrees of
freedom.
p: p-value.
p.adj: the adjusted p-value.
method: the statistical test used to compare groups.
p.signif, p.adj.signif: the significance level of p-values and
adjusted p-values, respectively.
estimate: estimate of the
effect size, that is the difference between the two emmeans (estimated
marginal means).
conf.low,conf.high: Lower and upper bound on a
confidence interval of the estimate.
The returned object has an attribute called args, which is a list holding the test arguments. It has also an attribute named "emmeans", a data frame containing the groups emmeans.
get_emmeans(): returns the estimated marginal means from an object of class emmeans_test
# Data preparation
df <- ToothGrowth
df$dose <- as.factor(df$dose)
# Pairwise comparisons
res <- df %>%
group_by(supp) %>%
emmeans_test(len ~ dose, p.adjust.method = "bonferroni")
res
#> # A tibble: 6 × 10
#> supp term .y. group1 group2 df statistic p p.adj p.adj.signif
#> * <fct> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 OJ dose len 0.5 1 54 -5.83 3.18e- 7 9.53e- 7 ****
#> 2 OJ dose len 0.5 2 54 -7.90 1.43e-10 4.29e-10 ****
#> 3 OJ dose len 1 2 54 -2.07 4.34e- 2 1.30e- 1 ns
#> 4 VC dose len 0.5 1 54 -5.41 1.46e- 6 4.39e- 6 ****
#> 5 VC dose len 0.5 2 54 -11.2 1.13e-15 3.39e-15 ****
#> 6 VC dose len 1 2 54 -5.77 3.98e- 7 1.19e- 6 ****
# Display estimated marginal means
attr(res, "emmeans")
#> # A tibble: 6 × 8
#> supp dose emmean se df conf.low conf.high method
#> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 OJ 0.5 13.2 1.15 54 10.9 15.5 Emmeans test
#> 2 OJ 1 22.7 1.15 54 20.4 25.0 Emmeans test
#> 3 OJ 2 26.1 1.15 54 23.8 28.4 Emmeans test
#> 4 VC 0.5 7.98 1.15 54 5.68 10.3 Emmeans test
#> 5 VC 1 16.8 1.15 54 14.5 19.1 Emmeans test
#> 6 VC 2 26.1 1.15 54 23.8 28.4 Emmeans test
# Show details
df %>%
group_by(supp) %>%
emmeans_test(len ~ dose, p.adjust.method = "bonferroni", detailed = TRUE)
#> # A tibble: 6 × 15
#> supp term .y. group1 group2 null.value estimate se df conf.low
#> * <fct> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 OJ dose len 0.5 1 0 -9.47 1.62 54 -12.7
#> 2 OJ dose len 0.5 2 0 -12.8 1.62 54 -16.1
#> 3 OJ dose len 1 2 0 -3.36 1.62 54 -6.62
#> 4 VC dose len 0.5 1 0 -8.79 1.62 54 -12.0
#> 5 VC dose len 0.5 2 0 -18.2 1.62 54 -21.4
#> 6 VC dose len 1 2 0 -9.37 1.62 54 -12.6
#> # ℹ 5 more variables: conf.high <dbl>, statistic <dbl>, p <dbl>, p.adj <dbl>,
#> # p.adj.signif <chr>
# Marginal means averaged over another factor (e.g. a 2x3 design).
# Fit the full model and pass it with `model =` so that the estimated
# marginal means for `dose` are averaged over `supp` (instead of fitting
# `len ~ dose` alone, which would ignore `supp`):
model <- lm(len ~ supp * dose, data = df)
df %>% emmeans_test(len ~ dose, model = model)
#> NOTE: Results may be misleading due to involvement in interactions
#> # A tibble: 3 × 9
#> term .y. group1 group2 df statistic p p.adj p.adj.signif
#> * <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 dose len 0.5 1 54 -7.95 1.19e-10 3.56e-10 ****
#> 2 dose len 0.5 2 54 -13.5 6.17e-19 1.85e-18 ****
#> 3 dose len 1 2 54 -5.54 9.12e- 7 2.74e- 6 ****
# Repeated-measures / mixed designs: pass a fitted within-subject model
# (e.g. stats::aov() with an Error() term, or nlme::lme()) with `model =`:
# \donttest{
set.seed(123)
d <- data.frame(
id = factor(rep(1:10, 3)),
time = factor(rep(c("t1", "t2", "t3"), each = 10)),
score = rnorm(30)
)
rm_model <- stats::aov(score ~ time + Error(id / time), data = d)
d %>% emmeans_test(score ~ time, model = rm_model)
#> Note: re-fitting model with sum-to-zero contrasts
#> Error in eval(cl): object 'd' not found
# }