Pairwise Comparisons of Estimated Marginal Means

Source: R/emmeans_test.R

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)

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 objects such as the result of a call to lm(), from which the overall degrees of freedom are to be calculated.

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 estim…¹    se    df conf.…² conf.…³
#> * <fct> <chr> <chr> <chr>  <chr>       <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl>
#> 1 OJ    dose  len   0.5    1               0   -9.47  1.62    54  -12.7   -6.21 
#> 2 OJ    dose  len   0.5    2               0  -12.8   1.62    54  -16.1   -9.57 
#> 3 OJ    dose  len   1      2               0   -3.36  1.62    54   -6.62  -0.104
#> 4 VC    dose  len   0.5    1               0   -8.79  1.62    54  -12.0   -5.53 
#> 5 VC    dose  len   0.5    2               0  -18.2   1.62    54  -21.4  -14.9  
#> 6 VC    dose  len   1      2               0   -9.37  1.62    54  -12.6   -6.11 
#> # … with 4 more variables: statistic <dbl>, p <dbl>, p.adj <dbl>,
#> #   p.adj.signif <chr>, and abbreviated variable names ¹​estimate, ²​conf.low,
#> #   ³​conf.high