Labels for pairwise multiple comparisons

stat_multcomp fits a linear model by default with stats::lm() but alternatively using other model fit functions. The model is passed to function glht() from package 'multcomp' to fit Tukey, Dunnet or other pairwise contrasts and generates labels based on adjusted P-values.

Usage

stat_multcomp(
  mapping = NULL,
  data = NULL,
  geom = NULL,
  position = "identity",
  ...,
  orientation = NA,
  formula = y ~ factor(x),
  method = "lm",
  method.args = list(),
  contrasts = "Tukey",
  p.adjust.method = NULL,
  fit.seed = NA,
  fm.cutoff.p.value = 1,
  mc.cutoff.p.value = 1,
  mc.critical.p.value = 0.05,
  small.p = getOption("ggpmisc.small.p", default = FALSE),
  adj.method.tag = 4,
  p.digits = NULL,
  label.type = "bars",
  bjust = 0.5,
  label.y = NULL,
  vstep = NULL,
  output.type = NULL,
  na.rm = FALSE,
  parse = NULL,
  show.legend = FALSE,
  inherit.aes = TRUE
)

Arguments

mapping: The aesthetic mapping, usually constructed with aes(). Only needs to be set at the layer level if you are overriding the plot defaults.
data: A layer specific dataset, only needed if you want to override the plot defaults.
geom: The geometric object to use display the data
position: The position adjustment to use for overlapping points on this layer.
...: other arguments passed on to layer. This can include aesthetics whose values you want to set, not map. See layer for more details.
orientation: character The orientation of the layer. The default (NA) automatically determines the orientation from the aesthetic mapping. In the rare event that this fails it can be given explicitly by setting orientation to either "x" or "y".
formula: a formula object. Using aesthetic names x and y instead of original variable names. The rhs must include a call to factor() even if the variable mapped to the x aesthetic is a factor! In both flipped and not flipped plots, x should be in the rhs.
method: function or character If character, "lm", "rlm", "lmrob", "lts", "gls", "ma", "sma", "segreg", "rq" or the name of a model fit function are accepted, possibly followed by the fit function's method argument separated by a colon (e.g. "rlm:M"). If a function is different to lm(), rlm(), ltsReg(), gls(), ma, sma, it must have formal parameters named formula, data, and weights. See Details.
method.args: named list with additional arguments. Not data or weights which are always passed through aesthetic mappings.
contrasts: character vector of length one or a numeric matrix. If character, one of "Tukey" or "Dunnet". If a matrix, one column per level of the factor mapped to x and one row per pairwise contrast.
p.adjust.method: character As the argument for parameter type of function adjusted() passed as argument to parameter test of summary.glht. Accepted values are "single-step", "Shaffer", "Westfall", "free", "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none".
fit.seed: RNG seed argument passed to set.seed(). Defaults to NA, indicating that set.seed() should not be called.
fm.cutoff.p.value: numeric [0..1] The P-value for the main effect of factor x in the ANOVA test for the fitted model above which no pairwise comparisons are computed or labels generated. Be aware that recent literature tends to recommend to consider which testing approach is relevant to the problem at hand instead of requiring the significance of the main effect before applying multiple comparisons' tests. The default value is 1, imposing no restrictions.
mc.cutoff.p.value: numeric [0..1] The P-value for the individual contrasts above which no labelled bars are generated. Default is 1, labelling all pairwise contrasts tested.
mc.critical.p.value: numeric The critical P-value used for tests when encoded as letters.
small.p: logical Flag to switch use of lower case p for p-value.
adj.method.tag: numeric, character or function If numeric, the length in characters of the abbreviation of the method used to adjust p-values. A value of zero, adds no label and a negative value uses as starting point for the abbreviation the word "adjusted". If character its value is used as subscript. If a function, the value used is the value returned by the function when passed p.adjust.method as its only argument.
p.digits: integer Number of digits after the decimal point to use for P-value in labels. If Inf, use exponential notation with three decimal places.
label.type: character One of "bars", "letters" or "LETTERS", selects how the results of the multiple comparisons are displayed. Only "bars" can be used together with contrasts = "Dunnet".
bjust: numeric usually with range 0..1, the justification of the labels relative to the bars.
label.y: numeric with range 0..1 "normalized parent coordinates" (npc units) or character if using geom_text_npc() or geom_label_npc(). If using geom_text() or geom_label() numeric in native data units. In flipped plots, it refers to x rather than y position.
vstep: numeric in npc units, the vertical displacement step-size used between labels for different contrasts when label.type = "bars". In flipped plots, it increases the x rather than y position.
output.type: character One of "expression", "text", "markdown", "marquee", "latex", "latex.eqn", "latex.deqn" or "numeric".
na.rm: a logical indicating whether NA values should be stripped before the computation proceeds.
parse: logical Passed to the geom. If TRUE, the labels will be parsed into expressions and displayed as described in plotmath. Default is TRUE if output.type = "expression" and FALSE otherwise.
show.legend: logical. Should this layer be included in the legends? NA, the default, includes if any aesthetics are mapped. FALSE never includes, and TRUE always includes.
inherit.aes: If FALSE, overrides the default aesthetics, rather than combining with them. This is most useful for helper functions that define both data and aesthetics and shouldn't inherit behaviour from the default plot specification, e.g. borders.

Value

A data frame with one row per comparison for label.type = "bars", or a data frame with one row per factor level for label.type = "letters" and label.type = "LETTERS". Variables (= columns) as described under Computed variables.

Details

This statistic can be used to automatically annotate a plot with P-values for pairwise multiple comparison tests, based on Tukey contrasts (all pairwise), Dunnet contrasts (other levels against the first one) or a subset of all possible pairwise contrasts. See Meier (2022, Chapter 3) for an accessible explanation of multiple comparisons and contrasts with package 'multcomp', of which stat_multcomp() is mostly a wrapper. It supports most model fit methods supported by function glht(). The requirement is that the necessary query functions for the class of the fitted model object are available. If not supported, an error will be triggered by glht(). The computations are done in two stages. First a model is fitted to all the data in a plot panel, followed by the pairwise contrasts. A threshold $P$-value can the set for the main effect of "treatments" as a condition to apply the pairwise tests.

One variable mapped to either x or y aesthetic must be a factor, and the other variable a continuous one. Mapping a factor creates the required grouping and determines the default orientation, i.e., similarly as in stat_boxplot().

Two ways of displaying the outcomes are implemented, and are selected by `"bars"`, `"letters"` or `"LETTERS"` as argument to parameter `label.type`. `"letters"` and `"LETTERS"` can be used only with Tukey contrasts, as otherwise the encoding is ambiguous. As too many bars clutter a plot, the maximum number of factor levels supported for `"bars"` together with Tukey contrasts is five, while together with Dunnet contrasts or contrasts defined by a numeric matrix, no limit is imposed. label.y determines the location of the letters labels or the position of the first pairwise label. If label.y is numeric and has length > 1, values are taken as positions for successive bar labels, and recycled if the vector is too short, ignoring vstep.

stat_multcomp() by default generates character labels ready to be parsed as R expressions but LaTeX (use TikZ device), markdown (use package 'ggtext') and plain text are also supported, as well as numeric values for user-generated text labels. The value of parse is set automatically based on output.type, but if you assemble labels that need parsing from numeric output, the default needs to be overridden. This statistic only generates annotation labels and segments connecting the compared factor levels, or letter labels that discriminate significantly different groups.

Currently, contrasts that involve more than two levels of a factor, such as the average of two treatment levels against a control level are not supported. Similarly, flipping is not yet functional with bar labels because of the lack of a working geom.

Note

R option OutDec is obeyed based on its value at the time the plot is rendered, i.e., displayed or printed. Set options(OutDec = ",") for languages like Spanish or French.

Computed variables

Computed variables and their names vary depending on the method used to fit a model or the output.type in use. They can also depend for a given method on the contrasts.

If output.type = "numeric" and label.type = "bars" the returned tibble contains columns listed below. In all cases if the model fit function used does not return a value, the label is set to character(0L) and the numeric value to NA.

x,xmin,xmax: x position, numeric.
y,ymin,ymax: y position, numeric.
coefficients: Delta estimate from pairwise contrasts, numeric.
contrasts: Contrasts as two levels' ordinal "numbers" separated by a dash, character.
tstat: t-statistic estimates for the pairwise contrasts, numeric.
p.value: P-value for the pairwise contrasts.
p.signif: logical for the pairwise contrasts significance.
fm.method: Set according method used.
fm.class: Most derived class of the fitted model object.
fm.formula: Formula extracted from the fitted model object if available, or the formula argument.
fm.formula.chr: Formula extracted from the fitted model object if available, or the formula argument, formatted as character.
mc.adjusted: The method used to adjust the P-values.
mc.contrast: The type of contrast used for multiple comparisons.
n: The total number of observations or rows in data.
default.label: text label, always included, but possibly NA.

If output.type is not "numeric" the returned data frame includes in addition the following labels:

stars.label: P-value for the pairwise contrasts encoded as "starts", character.
p.value.label: P-value for the pairwise contrasts, character.
delta.label: The coefficient or estimate for the difference between compared pairs of levels.
t.value.label: t-statistic estimates for the pairwise contrasts, character.

If output.type = "numeric" and label.type = "letters" or label.type = "LETTERS" the returned tibble contains columns listed below.

x,xmin,xmax: x position, numeric.
y,ymin,ymax: y position, numeric.
critical.p.value: P-value used in pairwise tests, numeric.
fm.method: Set according method used.
fm.class: Most derived class of the fitted model object.
fm.formula: Formula extracted from the fitted model object if available, or the formula argument.
fm.formula.chr: Formula extracted from the fitted model object if available, or the formula argument, formatted as character.
mc.adjusted: The method used to adjust the P-values.
mc.contrast: The type of contrast used for multiple comparisons.
n: The total number of observations or rows in data.
default.label: text label, always included, but possibly NA.

If output.type is not "numeric" the returned data frame includes in addition the following labels:

letters.label: Letters that distinguish levels based on significance from multiple comparisons test.

Alternatives

stat_signif() in package 'ggsignif' is an earlier and independent implementation of pairwise tests.

Output types

The formatting of character strings to be displayed in plots are marked as mathematical equations. Depending on the geom used, the mark-up needs to be encoded differently, or in some cases mark-up not applied.

"expression": The labels are encoded as character strings to be parsed into R's plotmath expressions.
"LaTeX", "TeX", "tikz", "latex": The labels are encoded as 'LaTeX' maths equations, without the "fences" for switching in math mode.
"latex.eqn": Same as "latex" but enclosed in single $, i.e., as in-line maths.
"latex.deqn": Same as "latex" but enclosed in double $$, i.e., as display maths.
"markdown": The labels are encoded as character strings using markdown syntax, with some embedded HTML.
"marquee": The labels are encoded as character strings using markdown syntax, with 'marquee' supported spans.
"text": The labels are plain ASCII character strings.
"numeric": No labels are generated. This value is accepted by the statistics, but not by the label formatting functions.
NULL: The value used depends on the argument passed to geom.

If geom = "latex" (package 'xdvir') the output type used is "latex.eqn". If geom = "richtext" (package 'ggtext') or geom = "textbox" (package 'ggtext') the output type used is "markdown". If geom = "marquee" (package 'marquee') the output type used is "marquee". For all other values of geom the default is "expression". Invalid values as argument trigger an error.

Which variables are available for mapping?

Computed variables and their names can vary depending on the method used to fit a model or the output.type in use. They can also depend for a given method on other arguments passed when fitting a model or extracting estimates and other computed values. In many cases, when values are not available, the variables are filled with NA values.

In the statistics returning formatted strings for use as annotations, a message is issued by default in interactive R sessions, listing the short names for available formatted labels as recognized by functions use_label() and f_use_label(), except when output.type = "numeric" is passed, in which case the names of all variables accessible by after_stat() within a call to aes() are listed. This default ("nicknames") can be changed by setting R option "ggpmisc.stat.vars.message" to one of "names", "colnames" or "none".

In the statistics that plot a prediction or more generally mainly return numeric variables, a message is issued by default in interactive R sessions, listing the names of all variables accessible by after_stat() within a call to aes() with at least some non-missing values. This default ("colnames") can be changed by setting R option "ggpmisc.stat.vars.message" to "none".

To explore the whole returned data frame for a given input we suggest the use of geom_debug().

References

Meier, Lukas (2022) ANOVA and Mixed Models: A Short Introduction Using R. Chapter 3 Contrasts and Multiple Testing. The R Series. Boca Raton: Chapman and Hall/CRC. ISBN: 9780367704209, doi:10.1201/9781003146216 .

Aesthetics

stat_multcomp() understands the following aesthetics. Required aesthetics are displayed in bold and defaults are displayed for optional aesthetics:

•	`x`
•	`y`
•	`group`	→ inferred
•	`hjust`	→ `after_stat(just)`
•	`label`	→ `after_stat(default.label)`
•	`size`	→ `2.5`
•	`weight`	→ `1`
•	`xmax`	→ `after_stat(ifelse(flipped_aes, NA_real_, xmax))`
•	`xmin`	→ `after_stat(ifelse(flipped_aes, NA_real_, xmin))`
•	`ymax`	→ `after_stat(ifelse(flipped_aes, ymax, NA_real_))`
•	`ymin`	→ `after_stat(ifelse(flipped_aes, ymin, NA_real_))`

Learn more about setting these aesthetics in vignette("ggplot2-specs").

Examples


p1 <- ggplot(mpg, aes(factor(cyl), hwy)) +
  geom_boxplot(width = 0.33)

## labeleld bars

p1 +
  stat_multcomp()


p1 +
  stat_multcomp(adj.method.tag = FALSE)


# test against a control, with first level being the control
# change order of factor levels in data to set the control group
p1 +
  stat_multcomp(contrasts = "Dunnet")


# arbitrary pairwise contrasts, in arbitrary order
p1 +
  stat_multcomp(contrasts = rbind(c(0, 0, -1, 1),
                                  c(0, -1, 1, 0),
                                  c(-1, 1, 0, 0)))


# different methods to adjust the contrasts
p1 +
  stat_multcomp(p.adjust.method = "holm")


p1 +
  stat_multcomp(p.adjust.method = "fdr")


# position of contrasts' bars (based on scale limits)
p1 +
  stat_multcomp(label.y = "bottom")


p1 +
  stat_multcomp(label.y = 11)


# use different labels: difference and P-value from hypothesis tests
p1 +
  stat_multcomp(use_label("Delta", "P"),
                size = 2.75)


# control smallest P-value displayed and number of digits
p1 +
  stat_multcomp(p.digits = 4)


# highlight significant differences
p1 +
  stat_multcomp(aes(alpha = after_stat(p.signif)),
                alpha.target = c("text", "segment", "box.line")) +
  scale_alpha_manual(values =  c(0.33, 1))


# label only significant differences
# but test and correct for all pairwise contrasts!
p1 +
  stat_multcomp(mc.cutoff.p.value = 0.01)


## letters as labels for test results

p1 +
  stat_multcomp(label.type = "letters")


# stricter critical p-value than default used for test
p1 +
  stat_multcomp(label.type = "letters",
                mc.critical.p.value = 0.01)


# justification of label in bar, and text anchor point
# for a compact display

p1 +
  stat_multcomp(geom = "text_pairwise",
                arrow = grid::arrow(ends = "both",
                                    length = unit(1.5, "mm")),
                bjust = 0,
                hjust = 1.05,
                vjust = 0.5,
                vstep = 0.05,
                label.y = 14,
                adj.method.tag = 0) +
  expand_limits(y = 0)


## Flipping is supported
p2 <- ggplot(mpg, aes(hwy, factor(cyl))) +
  geom_boxplot(width = 0.33)

p2 +
  stat_multcomp(label.type = "letters") +
  scale_y_discrete(expand = expansion(add = c(1.2, 0.5)))


p2 +
  stat_multcomp(label.type = "bars",
                geom = "text_pairwise",
                size = 3)