Correcting confidence intervals for multiple comparisons using the multcomp package

No, you can't specify type of adjustment when computing confidence intervals with multcomp. As discussed in the comments, you can choose between individual confidence intervals and simultaneous confidence intervals.

• Individual confidence intervals: Control the coverage probability of each confidence interval separately, ie., each confidence interval will contain the true parameter 100(1-α)% times if the experiment is repeated.
• Simultaneous confidence intervals: Control the overall coverage probability, ie, all confidence intervals will contain their true parameter 100(1-α)% times if the experiment is repeated. This corresponds to family-wise error rate FWER = α because the probability of one or more false discoveries (type I errors) is α. [However, simultaneous CIs are not equivalent to Bonferroni correction; see the update below.]

multcomp takes into account the correlations between hypotheses to get the coverage of simultaneous confidence intervals right and that's all the adjustment needed. The documentation is somewhat clear about all this.

To compute p-values, use the summary function on a glht object and specify how to adjust the p-values with the test = adjusted(type = "method") argument. The available methods include "single-step" (the default) and "bonferroni".

To compute confidence intervals, use the confint function of a glht object and specify how to compute the critical value with the calpha argument. There are two options: univariate_calpha() to compute univariate confidence intervals and adjusted_calpha() to compute simultaneous confidence intervals. Basically, confint ignores the type = "method" argument, if provided.

If multcomp doesn't adjust simultaneous confidence intervals, why does it seem that something is going on? It's because there is randomness in the computation of the critical value, calpha.

multcomp uses mvtnorm::qmvt to compute calpha, which in turn uses as a stochastic algorithm to find the quantiles a multivariate normal distribution.

In effect this is what you did: You fix the random seed once and call a stochastic function multiple times. You get slightly different critical values.

set.seed(1234)
calpha <- qmvt(level, ...)
calpha <- qmvt(level, ...)
calpha <- qmvt(level, ...)

Instead fix the seed before calling the stochastic function. You get exactly the same critical value each time and the simultaneous confidence intervals won't change.

set.seed(1234)
calpha <- qmvt(level, ...)
set.seed(1234)
calpha <- qmvt(level, ...)
set.seed(1234)
calpha <- qmvt(level, ...)

Let's demonstrate this on your example.

library("broom")
library("multcomp")
library("tidyverse")

options(
  pillar.sigfig = 5
)

mod <- lm(
  formula = mpg ~ cyl + disp + hp + drat + wt + qsec + vs + gear,
  data = mtcars
)

# Usually we aren't interested in testing that the intercept is 0.
# So let's skip it.
K <- diag(length(coef(mod)))[-1, ]
rownames(K) <- names(coef(mod))[-1]
K
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> cyl     0    1    0    0    0    0    0    0    0
#> disp    0    0    1    0    0    0    0    0    0
#> hp      0    0    0    1    0    0    0    0    0
#> drat    0    0    0    0    1    0    0    0    0
#> wt      0    0    0    0    0    1    0    0    0
#> qsec    0    0    0    0    0    0    1    0    0
#> vs      0    0    0    0    0    0    0    1    0
#> gear    0    0    0    0    0    0    0    0    1

modGLHT <- glht(
  model = mod,
  linfct = K
)

methods <- c("none", "single-step", "Shaffer", "Westfall")
names(methods) <- methods

There is randomness in computing the critical value, calpha. So we get slightly different simultaneous confidence intervals for the same set of hypotheses.

set.seed(1234)

methods %>%
  map_dfr(

    # `confint` silently ignores `adjusted`
    ~ tidy(confint(modGLHT, adjusted(type = .))),
    .id = "method"
  ) %>%
  filter(
    contrast == "cyl"
  )
#> # A tibble: 4 × 5
#>   method      contrast estimate conf.low conf.high
#>   <chr>       <chr>       <dbl>    <dbl>     <dbl>
#> 1 none        cyl      -0.52659  -3.3654    2.3122
#> 2 single-step cyl      -0.52659  -3.3665    2.3133
#> 3 Shaffer     cyl      -0.52659  -3.3670    2.3138
#> 4 Westfall    cyl      -0.52659  -3.3624    2.3093

We fix the random seed before computing the critical value calpha, so the simultaneous confidence intervals are computed in a reproducible way.

methods %>%
  map_dfr(

    # `confint` silently ignores `adjusted`
    ~ {
      set.seed(1234)
      tidy(confint(modGLHT, adjusted(type = .)))
    },
    .id = "method"
  ) %>%
  filter(
    contrast == "cyl"
  )
#> # A tibble: 4 × 5
#>   method      contrast estimate conf.low conf.high
#>   <chr>       <chr>       <dbl>    <dbl>     <dbl>
#> 1 none        cyl      -0.52659  -3.3648    2.3116
#> 2 single-step cyl      -0.52659  -3.3648    2.3116
#> 3 Shaffer     cyl      -0.52659  -3.3648    2.3116
#> 4 Westfall    cyl      -0.52659  -3.3648    2.3116

Updates

Bonferroni controls the familiy-wise error rate (FWER) at the α level and simultaneous confidence intervals are 100(1-α)% family-wise confidence intervals. So Bonferroni and simultaneous CIs have the FWER concept in common. But they are not equivalent. Why?

The Bonferroni method divides the p-values of m simultaneous hypotheses by m. The procedure is conservative: it guarantees that the FWER is ≤ α rather than = α and it tends to over-correct, a lot.

In general, we don't know the dependencies between hypotheses. However, we do with regression because we have an estimate of the covariance matrix Var{(β₁, β₂, ..., β_p)}. So we can estimate correlations between hypotheses which are functions of the model coefficients (such as β₁ = 0, β₂ = 0, ..., β_p = 0).

And this is how multcomp computes simultaneous confidence intervals: it takes into account the correlations between hypotheses to correct the critical value calpha "just the right amount". While Bonferroni ignores this information to apply the most stringent correction.

The simultaneous confidence intervals are determined by the set of hypotheses being tested. Changing the other hypotheses can lead to a different confidence interval for the same individual hypothesis because the overall coverage depends in a complex way on the correlations between all hypotheses.

set.seed(1234)

modGLHT <- glht(
  model = mod,
  linfct = c("cyl = 0", "disp = 0")
)

tidy(confint(modGLHT, adjusted(type = .)))
#> # A tibble: 2 × 4
#>   contrast  estimate  conf.low conf.high
#>   <chr>        <dbl>     <dbl>     <dbl>
#> 1 cyl      -0.52659  -2.8330    1.7798  
#> 2 disp      0.016565 -0.014587  0.047716