Skip to contents

Pairwise Comparisons of Estimated Marginal Means

Source: R/emmeans_test.R
emmeans_test.Rd

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.

Usage

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

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