6
$\begingroup$

I often get reviewers asking me to correct my results for multiple comparisons. It's quite straightforward to adjust both p-values and confidence intervals within models, i.e. via adjusting the familywise error rate, and there are several packages within R that do this (e.g. emmeans, glht). I was aided in my quest by a helpful CV user who cracked open the software I was using and helped me understand what was going on under the hood.

What I don't know is how to adjust confidence intervals for the experimentwise error rate, which I believe is what reviewers are asking for: i.e. adjusting p-values and CIs across multiple, potentially different statistical models. With p-values this is easy enough to do using functions like the p.adjust() function in R, where you can simply enter the p-values for all the tests you want to cover and the function adjusts those p-values. Crucially the p.adjust function has a range of options that provide protection from type-1 error while not being overly draconian and controlling out potentially legitimate effects (e.g. the Benjamini-Hochberg procedure). The draconian procedure I am talking about is of course the Bonferroni method, which has simplicity to recommend it but little else.

Here is an example from the mtcars dataset in R. Four separate simple regressions, one coefficient of interest in each, estimating the effect of cyl, disp, hp, and gear on mpg. Apologies to non-R users. Please don't take offence.

(sum1 <- summary(test1 <- lm(mpg ~ cyl, 
                             data = mtcars)))
# Output
#               Estimate Std. Error t value Pr(>|t|)    
#   (Intercept)  37.8846     2.0738   18.27  < 2e-16 ***
#   cyl          -2.8758     0.3224   -8.92 6.11e-10 ***

(sum2 <- summary(test2 <- lm(mpg ~ disp,
                             data = mtcars)))

# Output
#                Estimate Std. Error t value Pr(>|t|)    
#   (Intercept) 29.599855   1.229720  24.070  < 2e-16 ***
#   disp        -0.041215   0.004712  -8.747 9.38e-10 ***

(sum3 <- summary(test3 <- lm(mpg ~ hp,
                             data = mtcars)))

# Output
#               Estimate Std. Error t value Pr(>|t|)    
#   (Intercept) 30.09886    1.63392  18.421  < 2e-16 ***
#   hp          -0.06823    0.01012  -6.742 1.79e-07 ***


(sum4 <- summary(test4 <- lm(mpg ~ gear,
                             data = mtcars)))

# Output
#             Estimate Std. Error t value Pr(>|t|)   
# (Intercept)    5.623      4.916   1.144   0.2618   
# gear           3.923      1.308   2.999   0.0054 **

Now adjusting for p-values is pretty simple with the p.adjust() function: just extract the p-values for each estimate from the models...

unadjusted_p <- c(sum1$coefficients[2,4],
                  sum2$coefficients[2,4],
                  sum3$coefficients[2,4],
                  sum4$coefficients[2,4])

...and run them through the p_adjust() function, applying the false discovery rate procedure (also known as the Benjamini-Hochberg procedure).

p.adjust(p = unadjusted_p,
         method = "fdr") -> adjusted_p

Now if we compare of the unadjusted and adjusted p-values we can see that adjustments have been made

data.frame(test = 1:4, unadjusted_p, adjusted_p)

# Output
#   test unadjusted_p   adjusted_p
# 1    1 6.112687e-10 1.876065e-09
# 2    2 9.380327e-10 1.876065e-09
# 3    3 1.787835e-07 2.383780e-07
# 4    4 5.400948e-03 5.400948e-03

Now let's leave p-values aside and move on to CIs. Unfortunately if you adjust the p-values but not the CIs you can get confusing situations where a p-value is less than the alpha level (e.g. 0.05) but the confidence interval contains 0. All journals want confidence intervals these days so there's no option but to control CIs for the experimentwise error rate, i.e. across multiple tests/models

I read the help in the confint_adjust() function in the api2lm package and the authors describe how to calculate adjusted CIs using the Bonferroni and Working-Hotellier methods. It all comes down to the multiplier that goes in front of the standard error for the point estimate. The confint_adjust() function only works on multiple comparisons within a single statistical model. But I wondered if the same procedure could be applied across tests/models. I used the following procedure to adjust the confidence interval for the cyl coefficient in the first model, based on the fact that I was conducting four further tests/models

Call number of comparisons/tests $k$ and control across all those tests at $\alpha = 0.05$.

a <- .05  # alpha level 
k <- 4 # number of comparisons

Now to get the CI we need the point estimate for the cyl coefficient, and the standard error of the estimate. We also need the residual degrees of freedom for the model.

est_test1 <- sum1$coefficients[2, "Estimate"] # estimate
ese_test1 <- sum1$coefficients[2, "Std. Error"] # se of estimate
residdf_test1 <- test1$df.residual # residual degreees of freedom for model

Now for the crucial Bonferroni multiplier: we adjust the two-tailed p-value passed into the quantile function for the t-distribution by multiplying the denominator in $p = 1 - \frac{\alpha}{2}$ by $k$, the number of comparisons, so $p = 1 - \frac{\alpha}{2k}$ with degrees of freedom taken from the residual degrees of freedom from the model.

m_test1_bonf <- qt(p = 1-a/(2*k),
                   df = residdf_test1) # multiplier for bonferroni method using residual degrees of freedom for the model

Now for the working-hotelling multiplier: we take the square root of the quantile function where $p = 1 - \alpha$, $df1$ = number of comparisons, $df2$ = residual degrees of freedom from the model.

m_test1_wh <- sqrt(k*qf(p = 1-a, 
                        df1 = k,
                        df2 = residdf_test1)) 

Now we construct the bonferroni-adjusted CI by taking the standard error of the estimate of the cyl coefficient, multiplying by our bonferroni multiplier and taking that value either side of the point estimate, like so

ci_test1_bonf <- est_test1 + c(-1,1)*m_test1_bonf*ese_test1 

Now calculate the Working-Hotellier-adjusted CI

ci_test1_wh <- est_test1 + c(-1,1)*m_test1_wh*ese_test1 

Mow compare these adjusted CIs to the unadjusted

data.frame(test = 1,
           unadj_lowCI = round(confint(test1)[2,1],2),
           unadj_highCI = round(confint(test1)[2,2],2),
           bonf_lowCI = round(ci_test1_bonf[1],2),
           bonf_highCI = round(ci_test1_bonf[2],2),
           wh_lowCI = round(ci_test1_wh[1],2),
           wh_highCI = round(ci_test1_wh[2],2))

# Output
#   test unadj_lowCI unadj_highCI bonf_lowCI bonf_highCI wh_lowCI wh_highCI
# 1    1       -3.53        -2.22      -3.73       -2.02    -3.93     -1.82

So our confidence intervals have widened with the correction procedure. My assumption, if the procedure is correct, is that I could go on and apply the same procedure to the remaining CIs for the disp, hp, and gear coefficients from the remaining three models.

Now (finally) my question

  1. Is this procedure kosher? Could I apply these corrections across models in the same way they are applied within models? My niggling doubt concerns the residual degrees of freedom for the different models, which would be different in different situations.

  2. If the procedure is kosher, is there a less conservative method than these two, something akin to the Benjamini-Hochberg procedure for example? One of the users that responded to this post suggested that it was relatively simple to calculate CIs using methods equivalent to the false discovery rate but did not elaborate on what those methods are.

  3. Once again, if the procedure is kosher, is there any software that would allow one to correct CIs across multiple models relatively easily?

Any help much appreciated. Will post a bounty on this in a few days.

p.s. I know there are lots of questions about adjustng CIs on CV but I have not been able to find one that answers mine. I have asked this question before here but was lazy and unsurprisingly got no response. I hope I haven't overcorrected.

$\endgroup$
7
  • 1
    $\begingroup$ Why do you think that adjustment of your confidence intervals is helpful? Are you using them for hypothesis testing, or for estimation? Or for characterisation? $\endgroup$ Jan 18 at 2:54
  • $\begingroup$ @Michael Lew, both hypothesis testing and estimation. Not sure what you mean by characterisation. I also do it because reviewers demand it. $\endgroup$
    – llewmills
    Jan 18 at 2:57
  • 1
    $\begingroup$ Have a look at this answer to a related question, and note that reading the chapter linked in a comment would be very helpful. stats.stackexchange.com/questions/630316/… $\endgroup$ Jan 18 at 3:07
  • 1
    $\begingroup$ Well, if you are intent on supplying 'corrections' then I suggest that you do not do it to the intervals, but, instead, supply a table of p-values that you can adjust however you like. I would also suggest that your inferences be unaffected by those 'corrections' unless your study is one of those rare ones where the 'corrections' are more helpful than is usually the case. $\endgroup$ Jan 18 at 3:36
  • 1
    $\begingroup$ Read the chapter, at least section 3.2: link.springer.com/chapter/10.1007/164_2019_286 $\endgroup$ Jan 18 at 4:21

2 Answers 2

8
$\begingroup$

You really can't adjust CIs in that many ways beyond a straight Bonferroni or Sidak adjustment. I think the others in p.adjust (excepting "none" of course) all deal with various ways of nesting tests, whereby the price of being less conservative is paid by not doing certain tests, depending on the outcomes of others. You can't do that with CIs because you've got to construct all of them. You can only do "simultaneous CI" control (see Section 5.4.3 in Oehlert's text. Also that chapter covers most of the p.adjust methods, and you will find that they are all Bonferroni-based adjustments for various nesting approaches.

$\endgroup$
1
  • $\begingroup$ Thanks @Russ Lenth. If there are a lot of tests in a paper controlling for the experimentwise error rate using Bonferroni would be unworkably conservative and lead to a lot of Type 2 error as Michael Lew mentioned. Thank you for the information. $\endgroup$
    – llewmills
    Jan 19 at 23:26
8
$\begingroup$

You did not give the motivation for fitting separate models. In the context of one model a solution that is better than arbitrary multiplicity adjustments is to use simultaneous confidence intervals. In the R rms package the contrast.rms function has an option to do this and provides simultaneous coverage over the whole set of requested contrasts. It also automatically accounts for linearly redundant contrasts.

$\endgroup$
2
  • 1
    $\begingroup$ the motivation for separate models is having separate outcomes and different analyses. I am familiar with using simultaneous CIs within a single model (is the default in R packages like glht) but is it possible to use simultaneous confidence intervals with separate models? $\endgroup$
    – llewmills
    Jan 19 at 23:21
  • 1
    $\begingroup$ It’s an excellent question that I hope someone can answer. You’d need to model the dependence structure of the different outcomes, otherwise you have to use the ad hoc multiplicity adjustments discussed earlier and they can be quite conservative. If modeling in a Bayesian way the easiest way to model dependence is through a shared random effect across outcomes. $\endgroup$ Jan 20 at 12:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.