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