Pairwise Comparisons of Estimated Marginal Means

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.

See the Datanovia tutorial One-Way ANOVA in R for a worked walkthrough.

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)

Arguments

data

a data.frame containing the variables in the formula.

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.

covariate

(optional) covariate names (for ANCOVA)

ref.group

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).

comparisons

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"))

p.adjust.method

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".

conf.level

confidence level of the interval.

model

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.

detailed

logical value. Default is FALSE. If TRUE, a detailed result is shown.

emmeans.test

an object of class emmeans_test.

Value

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.

Functions

get_emmeans(): returns the estimated marginal means from an object of class emmeans_test

Examples

if (requireNamespace("emmeans", quietly = TRUE)) {

# 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

# Display estimated marginal means
attr(res, "emmeans")

# Show details
df %>%
 group_by(supp) %>%
 emmeans_test(len ~ dose, p.adjust.method = "bonferroni", detailed = TRUE)

# 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)

# 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: Results may be misleading due to involvement in interactions
#> Note: re-fitting model with sum-to-zero contrasts
#> Error in eval(cl): object 'd' not found

Pairwise Comparisons of Estimated Marginal Means

Arguments

Value

Functions

See also

Examples