Comparison of Means — compare

Performs one or multiple mean comparisons.

Usage

compare_means(
  formula,
  data,
  method = "wilcox.test",
  paired = FALSE,
  id = NULL,
  group.by = NULL,
  ref.group = NULL,
  symnum.args = list(),
  p.adjust.method = "holm",
  p.format.style = "default",
  p.digits = NULL,
  p.leading.zero = NULL,
  p.min.threshold = NULL,
  p.decimal.mark = NULL,
  signif.cutoffs = NULL,
  signif.symbols = NULL,
  ns.symbol = "ns",
  use.four.stars = FALSE,
  ...
)

Arguments

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.

It's also possible to perform the test for multiple response variables at the same time. For example, formula = c(TP53, PTEN) ~ cancer_group.

data

a data.frame containing the variables in the formula.

method

the type of test. Default is wilcox.test. Allowed values include:

t.test (parametric) and wilcox.test (non-parametric). Perform comparison between two groups of samples. If the grouping variable contains more than two levels, then a pairwise comparison is performed.
anova (parametric) and kruskal.test (non-parametric). Perform one-way ANOVA test comparing multiple groups.

paired

a logical indicating whether you want a paired test. Used only in t.test and in wilcox.test.

id

optional character string naming a column that identifies matched subjects for a paired comparison (method = "t.test" or "wilcox.test"). By default (id = NULL) a paired test pairs observations by row order, so a p-value can be wrong if the data are not sorted so that the compared groups align. Providing id pairs the observations by subject id instead (row-order independent), using only the complete pairs (per-comparison pairwise deletion, via rstatix). It works for a two-group, a pairwise (more than two groups) and a ref.group comparison; it is an error to combine id with anova/kruskal.test or with ref.group = ".all.".

group.by

a character vector containing the name of grouping variables.

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

ref.group can be also ".all.". In this case, each of the grouping variable levels is compared to all (i.e. basemean).

symnum.args

a list of arguments to pass to the function symnum for symbolic number coding of p-values. For example, symnum.args <- list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, Inf), symbols = c("****", "***", "**", "*", "ns")).

In other words, we use the following convention for symbols indicating statistical significance:

ns: p > 0.05
*: p <= 0.05
**: p <= 0.01
***: p <= 0.001
****: p <= 0.0001

Note: If signif.cutoffs is provided, it takes precedence over symnum.args.

p.adjust.method

method for adjusting p values (see p.adjust). Has impact only in a situation, where multiple pairwise tests are performed; or when there are multiple grouping variables. 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".

Note that, when the formula contains multiple variables, the p-value adjustment is done independently for each variable.

p.format.style

character string specifying the p-value formatting style. One of: "default" (backward compatible, uses scientific notation), "apa" (APA style, no leading zero), "nejm" (NEJM style), "lancet" (Lancet style), "ama" (AMA style), "graphpad" (GraphPad style), or "scientific" (scientific notation for GWAS). See list_p_format_styles for details.

p.digits

integer specifying the number of decimal places for p-values. If provided, overrides the style default.

p.leading.zero

logical indicating whether to include leading zero before decimal point (e.g., "0.05" vs ".05"). If provided, overrides the style default.

p.min.threshold

numeric specifying the minimum p-value to display exactly. Values below this threshold are shown as "< threshold". If NULL, the selected style's default threshold is used; styles without a threshold show exact values. If provided, overrides the style default.

p.decimal.mark

character string to use as the decimal mark. If NULL, uses getOption("OutDec").

signif.cutoffs

numeric vector of p-value cutoffs in descending order for assigning significance symbols. For example, c(0.10, 0.05, 0.01) means p < 0.10 gets "*", p < 0.05 gets "**", p < 0.01 gets "***". If use.four.stars = TRUE, can include a fourth level (e.g., c(0.10, 0.05, 0.01, 0.001) for "****" at p < 0.001). Default is NULL, which uses the package defaults (backward compatible).

signif.symbols

character vector of symbols corresponding to signif.cutoffs. If NULL, auto-generated as "*", "**", "***" (and "****" if use.four.stars = TRUE). Must have the same length as signif.cutoffs.

ns.symbol

character string for non-significant results. Default is "ns". Use "" (empty string) to show nothing for non-significant results.

use.four.stars

logical. If TRUE and signif.symbols is NULL, allows four stars (****) for the most significant level when signif.cutoffs has 4 levels. Default is FALSE.

...

Other arguments to be passed to the test function.

Value

a data frame with the following columns:

.y.: the y variable used in the test.
group1,group2: the compared groups in the pairwise tests. Available only when method = "t.test" or method = "wilcox.test".
p: the p-value.
p.adj: the adjusted p-value. Default for p.adjust.method = "holm".
p.format: the formatted p-value.
p.format.signif: the formatted p-value with significance symbols.
p.signif: the significance level.
method: the statistical test used to compare groups.

Significance letters (compact letter display)

To label groups with letters instead of p-values or stars - groups that share a letter are not significantly different - compute the pairwise comparisons and derive the letters with rstatix::add_cld() (available without any extra package, as rstatix is already a dependency), then place them with geom_text():


library(ggpubr)
library(rstatix)
library(dplyr)
df <- ToothGrowth
df$dose <- factor(df$dose)

# all-pairwise comparisons, then compact letters (columns: group, cld)
cld <- df %>% tukey_hsd(len ~ dose) %>% add_cld()

# one letter per group, placed above each box
ypos <- df %>% group_by(dose) %>%
  summarise(y.position = max(len) + 2, .groups = "drop")
cld$y.position <- ypos$y.position[match(cld$group, as.character(ypos$dose))]

ggboxplot(df, "dose", "len") +
  geom_text(data = cld, aes(x = group, y = y.position, label = cld))

add_cld() expects an all-pairwise result (tukey_hsd(), dunn_test(), games_howell_test(), pairwise wilcox_test()/t_test(), or compare_means()); it is not defined for a single ref.group comparison.

To instead flag each treatment by which of several controls it differs from (e.g. "a" for a negative control and "b" for a positive control), run one comparison per control and assemble the letters:


trts <- c("trtA", "trtB", "trtC")
cn <- compare_means(value ~ group, df, ref.group = "neg.ctrl")
cp <- compare_means(value ~ group, df, ref.group = "pos.ctrl")
pv <- function(cc, g) dplyr::filter(cc, group1 == g | group2 == g)$p[1]
lab <- sapply(trts, function(g)
  paste0(if (pv(cn, g) < .05) "a" else "",
         if (pv(cp, g) < .05) "b" else ""))
# then place `lab` above each treatment with geom_text()

Examples

# Load data
# :::::::::::::::::::::::::::::::::::::::
data("ToothGrowth")
df <- ToothGrowth

# One-sample test
# :::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ 1, df, mu = 0)
#> # A tibble: 1 × 9
#>   .y.   group1 group2         p   p.adj p.format p.signif method p.format.signif
#>   <chr>  <dbl> <chr>      <dbl>   <dbl> <chr>    <chr>    <chr>  <chr>          
#> 1 len        1 null m… 1.66e-11 1.7e-11 1.7e-11  ****     Wilco… p = 1.7e-11 **…

# Two-samples unpaired test
# :::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ supp, df)
#> # A tibble: 1 × 9
#>   .y.   group1 group2      p p.adj p.format p.signif method   p.format.signif
#>   <chr> <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>    <chr>          
#> 1 len   OJ     VC     0.0645 0.064 0.064    ns       Wilcoxon p = 0.064 ns   

# Two-samples paired test
# :::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ supp, df, paired = TRUE)
#> # A tibble: 1 × 9
#>   .y.   group1 group2       p  p.adj p.format p.signif method   p.format.signif
#>   <chr> <chr>  <chr>    <dbl>  <dbl> <chr>    <chr>    <chr>    <chr>          
#> 1 len   OJ     VC     0.00431 0.0043 0.0043   **       Wilcoxon p = 0.0043 **  

# Paired test pairing by a subject id column (row-order independent)
# :::::::::::::::::::::::::::::::::::::::::
df$id <- rep(1:30, 2) # pairs the two supp levels by subject
compare_means(len ~ supp, df, paired = TRUE, id = "id")
#> # A tibble: 1 × 9
#>   .y.   group1 group2       p  p.adj p.format p.signif method   p.format.signif
#>   <chr> <chr>  <chr>    <dbl>  <dbl> <chr>    <chr>    <chr>    <chr>          
#> 1 len   OJ     VC     0.00431 0.0043 0.0043   **       Wilcoxon p = 0.0043 **  

# Compare supp levels after grouping the data by "dose"
# ::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ supp, df, group.by = "dose")
#> # A tibble: 3 × 10
#>    dose .y.   group1 group2       p p.adj p.format p.signif method  
#>   <dbl> <chr> <chr>  <chr>    <dbl> <dbl> <chr>    <chr>    <chr>   
#> 1   0.5 len   OJ     VC     0.0232  0.023 0.023    *        Wilcoxon
#> 2   1   len   OJ     VC     0.00403 0.004 0.004    **       Wilcoxon
#> 3   2   len   OJ     VC     1       1     1        ns       Wilcoxon
#> # ℹ 1 more variable: p.format.signif <chr>

# pairwise comparisons
# ::::::::::::::::::::::::::::::::::::::::
# As dose contains more thant two levels ==>
# pairwise test is automatically performed.
compare_means(len ~ dose, df)
#> # A tibble: 3 × 9
#>   .y.   group1 group2          p  p.adj p.format p.signif method p.format.signif
#>   <chr> <chr>  <chr>       <dbl>  <dbl> <chr>    <chr>    <chr>  <chr>          
#> 1 len   0.5    1         7.02e-6 1.4e-5 7.0e-06  ****     Wilco… p = 7.0e-06 **…
#> 2 len   0.5    2         8.41e-8 2.5e-7 8.4e-08  ****     Wilco… p = 8.4e-08 **…
#> 3 len   1      2         1.77e-4 1.8e-4 0.00018  ***      Wilco… p = 0.00018 ***

# Comparison against reference group
# ::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ dose, df, ref.group = "0.5")
#> # A tibble: 2 × 9
#>   .y.   group1 group2          p  p.adj p.format p.signif method p.format.signif
#>   <chr> <chr>  <chr>       <dbl>  <dbl> <chr>    <chr>    <chr>  <chr>          
#> 1 len   0.5    1         7.02e-6 7  e-6 7.0e-06  ****     Wilco… p = 7.0e-06 **…
#> 2 len   0.5    2         8.41e-8 1.7e-7 8.4e-08  ****     Wilco… p = 8.4e-08 **…

# Comparison against all
# ::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ dose, df, ref.group = ".all.")
#> # A tibble: 3 × 9
#>   .y.   group1 group2         p   p.adj p.format p.signif method p.format.signif
#>   <chr> <chr>  <chr>      <dbl>   <dbl> <chr>    <chr>    <chr>  <chr>          
#> 1 len   .all.  0.5    0.0000508 0.00015 5.1e-05  ****     Wilco… p = 5.1e-05 **…
#> 2 len   .all.  1      0.764     0.76    0.76404  ns       Wilco… p = 0.76404 ns 
#> 3 len   .all.  2      0.000179  0.00036 0.00018  ***      Wilco… p = 0.00018 ***

# Anova and kruskal.test
# ::::::::::::::::::::::::::::::::::::::::
compare_means(len ~ dose, df, method = "anova")
#> # A tibble: 1 × 7
#>   .y.          p   p.adj p.format p.signif method p.format.signif 
#>   <chr>    <dbl>   <dbl> <chr>    <chr>    <chr>  <chr>           
#> 1 len   9.53e-16 9.5e-16 9.5e-16  ****     Anova  p = 9.5e-16 ****
compare_means(len ~ dose, df, method = "kruskal.test")
#> # A tibble: 1 × 7
#>   .y.               p        p.adj p.format p.signif method      p.format.signif
#>   <chr>         <dbl>        <dbl> <chr>    <chr>    <chr>       <chr>          
#> 1 len   0.00000000148 0.0000000015 1.5e-09  ****     Kruskal-Wa… p = 1.5e-09 **…