28  Multiple Comparisons

When conducting multiple hypothesis tests simultaneously, we need to account for the increased probability of Type I errors (false positives). The marginaleffects package provides several methods for multiple comparison adjustment through the hypotheses() function and its multcomp argument.

First, we fit a model

library(marginaleffects)
mod <- lm(Ozone ~ factor(Month) - 1, data = airquality)
hypotheses(mod)

           Term Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
 factor(Month)5     23.6       5.76  4.10  < 0.001 14.6  12.3   34.9
 factor(Month)6     29.4       9.79  3.01  0.00263  8.6  10.3   48.6
 factor(Month)7     59.1       5.76 10.27  < 0.001 79.7  47.8   70.4
 factor(Month)8     60.0       5.76 10.41  < 0.001 81.9  48.7   71.2
 factor(Month)9     31.4       5.45  5.77  < 0.001 26.9  20.8   42.1

hypotheses(mod, multcomp = "holm")

           Term Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
 factor(Month)5     23.6       5.76  4.10  < 0.001 13.6  8.82   38.4
 factor(Month)6     29.4       9.79  3.01  0.00263  8.6  4.30   54.6
 factor(Month)7     59.1       5.76 10.27  < 0.001  Inf 44.32   73.9
 factor(Month)8     60.0       5.76 10.41  < 0.001  Inf 45.17   74.8
 factor(Month)9     31.4       5.45  5.77  < 0.001 25.3 17.44   45.5

Now, we compute average counterfactual comparisons.

cmp <- avg_comparisons(mod)

Warning: The `Month` variable is treated as a categorical (factor) variable,
but the original data is of class integer. It is safer and faster to convert
such variables to factor before fitting the model and calling a
`marginaleffects` function.This warning appears once per session.

cmp

 Contrast Estimate Std. Error     z Pr(>|z|)    S  2.5 % 97.5 %
    6 - 5     5.83      11.36 0.513    0.608  0.7 -16.43   28.1
    7 - 5    35.50       8.14 4.359   <0.001 16.2  19.54   51.5
    8 - 5    36.35       8.14 4.463   <0.001 16.9  20.38   52.3
    9 - 5     7.83       7.93 0.988    0.323  1.6  -7.71   23.4

Term: Month
Type:  response 

Then, we apply a Holm adjustment to the p-values by calling the hypotheses() function and setting the multcomp argument to "holm". When we use the multcomp argument, the output will include adjusted p values and family-wise confidence intervals.

hypotheses(cmp, multcomp = "holm")

 Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
     5.83      11.36 0.513    0.647  0.6 -22.0   33.7
    35.50       8.14 4.359   <0.001 14.6  15.5   55.5
    36.35       8.14 4.463   <0.001 14.9  16.4   56.3
     7.83       7.93 0.988    0.647  0.6 -11.6   27.3

Term: Month

Notice that the p values and confidence intervals are larger than previously.

We can specify the conf_level argument explicitly. If we do not specify it, the size of confidence intervals will be inherited from the previous call to a marginaleffects function.

hypotheses(cmp, multcomp = "holm", conf_level = 0.95)

 Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
     5.83      11.36 0.513    0.647  0.6 -22.0   33.7
    35.50       8.14 4.359   <0.001 14.6  15.5   55.5
    36.35       8.14 4.463   <0.001 14.9  16.4   56.3
     7.83       7.93 0.988    0.647  0.6 -11.6   27.3

Term: Month

hypotheses(cmp, multcomp = "holm", conf_level = 0.90)

 Estimate Std. Error     z Pr(>|z|)    S  5.0 % 95.0 %
     5.83      11.36 0.513    0.647  0.6 -18.88   30.5
    35.50       8.14 4.359   <0.001 14.6  17.78   53.2
    36.35       8.14 4.463   <0.001 14.9  18.63   54.1
     7.83       7.93 0.988    0.647  0.6  -9.42   25.1

Term: Month

cmp <- avg_comparisons(mod, conf_level = 0.90)
hypotheses(cmp, multcomp = "holm")

 Estimate Std. Error     z Pr(>|z|)    S  5.0 % 95.0 %
     5.83      11.36 0.513    0.647  0.6 -18.88   30.5
    35.50       8.14 4.359   <0.001 14.6  17.78   53.2
    36.35       8.14 4.463   <0.001 14.9  18.63   54.1
     7.83       7.93 0.988    0.647  0.6  -9.42   25.1

Term: Month