Skip to contents

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",
  ...,
  formula = NULL,
  method = "lm",
  method.args = list(),
  contrasts = "Tukey",
  p.adjust.method = NULL,
  small.p = FALSE,
  adj.method.tag = 4,
  p.digits = 3,
  label.type = "bars",
  fm.cutoff.p.value = 1,
  mc.cutoff.p.value = 1,
  mc.critical.p.value = 0.05,
  label.y = NULL,
  vstep = NULL,
  output.type = NULL,
  na.rm = FALSE,
  orientation = "x",
  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 to 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.

formula

a formula object. Using aesthetic names x and y instead of original variable names.

method

function or character If character, "lm" (or its equivalent "aov"), "rlm" 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 different to lm(), it must accept as a minimum a model formula through its first parameter, and have formal parameters named data, weights, and method, and return a model fit object accepted by function glht().

method.args

named list with additional arguments.

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

small.p

logical If true, use of lower case p instead of capital P as the symbol for P-value in labels.

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 \(R^2\) and P-value in labels.

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

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.

label.y

numeric vector Values in native data units or if character, one of "top" or "bottom". Recycled if too short and truncated if too long.

vstep

numeric in npc units, the vertical displacement step-size used between labels for different contrasts when label.type = "bars".

output.type

character One of "expression", "LaTeX", "text", "markdown" or "numeric".

na.rm

a logical indicating whether NA values should be stripped before the computation proceeds.

orientation

character Either "x" or "y" controlling the default for formula. Support for orientation is not yet implemented but is planned.

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.

Value

A data frame with one row per comparison for label.type = "bars", or a data frame with one row per factor x level for

label.type = "letters" and for 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.

The explanatory variable mapped to the x aesthetic must be a factor as this creates the required grouping. 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, mainly because they require a new geometry that I need to design, implement and add to package 'ggpp'.

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.

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.

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.

Aesthetics

stat_multcomp() understands x and y, to be referenced in the formula and weight passed as argument to parameter weights. A factor must be mapped to x and numeric variables to y, and, if used, to weight. In addition, the aesthetics understood by the geom ("label_pairwise" is the default for label.type = "bars", "text" is the default for label.type = "letters" and for label.type = "LETTERS") are understood and grouping respected.

Computed variables

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,x.left.tip,x.right.tip

x position, numeric.

y

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.

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 label.type = "letters" or label.type = "LETTERS" the returned tibble contains columns listed below.

x,x.left.tip,x.right.tip

x position, numeric.

y

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.

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 .

See also

This statistic uses the implementation of Tests of General Linear Hypotheses in function glht. See summary.glht and p.adjust for the supported and tests and the references therein for the theory behind them.

Examples


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

## labeleld bars

p1 +
  stat_multcomp()


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


# 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 = "bonferroni")


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


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


# no correction, useful only for comparison
p1 +
  stat_multcomp(p.adjust.method = "none")


# sometimes we need to expand the plotting area
p1 +
  stat_multcomp(geom = "text_pairwise") +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.10)))


# 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(c("Delta", "P")),
                size = 2.75)


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


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


# use capital letters
p1 +
  stat_multcomp(label.type = "LETTERS")


# location
p1 +
  stat_multcomp(label.type = "letters",
                label.y = "top")


p1 +
  stat_multcomp(label.type = "letters",
                label.y = 0)


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


# Inspecting the returned data using geom_debug()
# This provides a quick way of finding out the names of the variables that
# are available for mapping to aesthetics with after_stat().

gginnards.installed <- requireNamespace("gginnards", quietly = TRUE)

if (gginnards.installed)
  library(gginnards)

if (gginnards.installed)
p1 +
  stat_multcomp(label.type = "bars",
                geom = "debug")

#> [1] "Summary of input 'data' to 'draw_panel()':"
#>   xmin xmax                        label  hjust PANEL   x  y x.left.tip
#> 1    1    2 italic(P)["HSD"]~`=`~"1.000" center     1 1.5 48          1
#> 2    1    3 italic(P)["HSD"]~`<`~"0.001" center     1 2.0 52          1
#> 3    1    4 italic(P)["HSD"]~`<`~"0.001" center     1 2.5 56          1
#> 4    2    3 italic(P)["HSD"]~`=`~"0.013" center     1 2.5 60          2
#> 5    2    4 italic(P)["HSD"]~`<`~"0.001" center     1 3.0 64          2
#> 6    3    4 italic(P)["HSD"]~`<`~"0.001" center     1 3.5 68          3
#>   x.right.tip coefficients contrasts        tstat      p.value fm.method
#> 1           2  -0.05246914       2-1  -0.02654434 9.999925e-01     lm:qr
#> 2           3  -5.97968433       3-1  -9.79894539 0.000000e+00     lm:qr
#> 3           4 -11.17389771       4-1 -17.74241622 0.000000e+00     lm:qr
#> 4           3  -5.92721519       3-2  -2.99681688 1.326898e-02     lm:qr
#> 5           4 -11.12142857       4-2  -5.60568576 2.007009e-07     lm:qr
#> 6           4  -5.19421338       4-3  -8.19962959 3.241851e-14     lm:qr
#>   fm.class    fm.formula fm.formula.chr mc.adjusted mc.contrast   n   just
#> 1       lm y ~ factor(x)  y ~ factor(x)         HSD         2-1 234 center
#> 2       lm y ~ factor(x)  y ~ factor(x)         HSD         3-1 234 center
#> 3       lm y ~ factor(x)  y ~ factor(x)         HSD         4-1 234 center
#> 4       lm y ~ factor(x)  y ~ factor(x)         HSD         3-2 234 center
#> 5       lm y ~ factor(x)  y ~ factor(x)         HSD         4-2 234 center
#> 6       lm y ~ factor(x)  y ~ factor(x)         HSD         4-3 234 center
#>   stars.label                p.value.label         delta.label
#> 1      "n.s." italic(P)["HSD"]~`=`~"1.000" Delta~`=`~"-0.0525"
#> 2       "***" italic(P)["HSD"]~`<`~"0.001"   Delta~`=`~"-5.98"
#> 3       "***" italic(P)["HSD"]~`<`~"0.001"   Delta~`=`~"-11.2"
#> 4         "*" italic(P)["HSD"]~`=`~"0.013"   Delta~`=`~"-5.93"
#> 5       "***" italic(P)["HSD"]~`<`~"0.001"   Delta~`=`~"-11.1"
#> 6       "***" italic(P)["HSD"]~`<`~"0.001"   Delta~`=`~"-5.19"
#>             t.value.label                default.label
#> 1 italic(t)~`=`~"-0.0265" italic(P)["HSD"]~`=`~"1.000"
#> 2   italic(t)~`=`~"-9.80" italic(P)["HSD"]~`<`~"0.001"
#> 3   italic(t)~`=`~"-17.7" italic(P)["HSD"]~`<`~"0.001"
#> 4    italic(t)~`=`~"-3.0" italic(P)["HSD"]~`=`~"0.013"
#> 5   italic(t)~`=`~"-5.61" italic(P)["HSD"]~`<`~"0.001"
#> 6   italic(t)~`=`~"-8.20" italic(P)["HSD"]~`<`~"0.001"

if (gginnards.installed)
p1 +
  stat_multcomp(label.type = "letters",
                geom = "debug")

#> [1] "Summary of input 'data' to 'draw_panel()':"
#>   xmin xmax                                   label  hjust PANEL   x    y
#> 1   NA   NA "  "*italic(P)["HSD"]^{"crit"}~`=`~0.05 inward     1 0.3 9.44
#> 2   NA   NA                                       c center     1 1.0 9.44
#> 3   NA   NA                                       c center     1 2.0 9.44
#> 4   NA   NA                                       b center     1 3.0 9.44
#> 5   NA   NA                                       a center     1 4.0 9.44
#>   x.left.tip x.right.tip critical.p.value fm.method fm.class    fm.formula
#> 1         NA          NA             0.05     lm:qr       lm y ~ factor(x)
#> 2         NA          NA             0.05     lm:qr       lm y ~ factor(x)
#> 3         NA          NA             0.05     lm:qr       lm y ~ factor(x)
#> 4         NA          NA             0.05     lm:qr       lm y ~ factor(x)
#> 5         NA          NA             0.05     lm:qr       lm y ~ factor(x)
#>   fm.formula.chr mc.adjusted mc.contrast   n
#> 1  y ~ factor(x)         HSD       Tukey 234
#> 2  y ~ factor(x)         HSD       Tukey 234
#> 3  y ~ factor(x)         HSD       Tukey 234
#> 4  y ~ factor(x)         HSD       Tukey 234
#> 5  y ~ factor(x)         HSD       Tukey 234
#>                             letters.label   just
#> 1 "  "*italic(P)["HSD"]^{"crit"}~`=`~0.05 inward
#> 2                                       c center
#> 3                                       c center
#> 4                                       b center
#> 5                                       a center
#>                             default.label
#> 1 "  "*italic(P)["HSD"]^{"crit"}~`=`~0.05
#> 2                                       c
#> 3                                       c
#> 4                                       b
#> 5                                       a

if (gginnards.installed)
p1 +
  stat_multcomp(label.type = "bars",
                output.type = "numeric",
                geom = "debug")

#> [1] "Summary of input 'data' to 'draw_panel()':"
#>   xmin xmax label  hjust PANEL   x  y x.left.tip x.right.tip coefficients
#> 1    1    2  <NA> center     1 1.5 48          1           2  -0.05246914
#> 2    1    3  <NA> center     1 2.0 52          1           3  -5.97968433
#> 3    1    4  <NA> center     1 2.5 56          1           4 -11.17389771
#> 4    2    3  <NA> center     1 2.5 60          2           3  -5.92721519
#> 5    2    4  <NA> center     1 3.0 64          2           4 -11.12142857
#> 6    3    4  <NA> center     1 3.5 68          3           4  -5.19421338
#>   contrasts        tstat      p.value fm.method fm.class    fm.formula
#> 1       2-1  -0.02654434 9.999925e-01     lm:qr       lm y ~ factor(x)
#> 2       3-1  -9.79894539 0.000000e+00     lm:qr       lm y ~ factor(x)
#> 3       4-1 -17.74241622 0.000000e+00     lm:qr       lm y ~ factor(x)
#> 4       3-2  -2.99681688 1.308100e-02     lm:qr       lm y ~ factor(x)
#> 5       4-2  -5.60568576 2.893068e-07     lm:qr       lm y ~ factor(x)
#> 6       4-3  -8.19962959 3.785861e-14     lm:qr       lm y ~ factor(x)
#>   fm.formula.chr mc.adjusted mc.contrast   n   just default.label
#> 1  y ~ factor(x)         HSD         2-1 234 center          <NA>
#> 2  y ~ factor(x)         HSD         3-1 234 center          <NA>
#> 3  y ~ factor(x)         HSD         4-1 234 center          <NA>
#> 4  y ~ factor(x)         HSD         3-2 234 center          <NA>
#> 5  y ~ factor(x)         HSD         4-2 234 center          <NA>
#> 6  y ~ factor(x)         HSD         4-3 234 center          <NA>