FDR-consistent confidence intervals for emmeans contrasts

The problem

Adjusting p-values for the false discovery rate (FDR) with the Benjamini-Hochberg procedure is common when many comparisons are made. In ‘emmeans’, this is done with adjust = "fdr". However, FDR adjustments are defined for tests, not for confidence intervals, so confint() replaces them with a Bonferroni adjustment. The intervals are then more conservative than the p-values, and the two can lead to different conclusions.

Consider the chickwts data (base R): chick weights after six weeks on six feed supplements, giving 15 pairwise comparisons.

library(fcrci)
library(emmeans)
#> Welcome to emmeans.
#> Caution: You lose important information if you filter this package's results.
#> See '? untidy'

fit <- lm(weight ~ feed, data = chickwts)
prs <- pairs(emmeans(fit, ~ feed), adjust = "fdr")

The FDR-adjusted p-values:

test(prs)
#>  contrast              estimate   SE df t.ratio p.value
#>  casein - horsebean      163.38 23.5 65   6.957 <0.0001
#>  casein - linseed        104.83 22.4 65   4.682 <0.0001
#>  casein - meatmeal        46.67 22.9 65   2.039  0.0570
#>  casein - soybean         77.15 21.6 65   3.576  0.0012
#>  casein - sunflower       -5.33 22.4 65  -0.238  0.8125
#>  horsebean - linseed     -58.55 23.5 65  -2.493  0.0228
#>  horsebean - meatmeal   -116.71 24.0 65  -4.870 <0.0001
#>  horsebean - soybean     -86.23 22.7 65  -3.797  0.0007
#>  horsebean - sunflower  -168.72 23.5 65  -7.184 <0.0001
#>  linseed - meatmeal      -58.16 22.9 65  -2.540  0.0225
#>  linseed - soybean       -27.68 21.6 65  -1.283  0.2187
#>  linseed - sunflower    -110.17 22.4 65  -4.920 <0.0001
#>  meatmeal - soybean       30.48 22.1 65   1.379  0.1991
#>  meatmeal - sunflower    -52.01 22.9 65  -2.271  0.0360
#>  soybean - sunflower     -82.49 21.6 65  -3.823  0.0007
#> 
#> P value adjustment: fdr method for 15 tests

And the confidence intervals of the same contrasts:

confint(prs)
#>  contrast              estimate   SE df lower.CL upper.CL
#>  casein - horsebean      163.38 23.5 65     91.8    235.0
#>  casein - linseed        104.83 22.4 65     36.6    173.1
#>  casein - meatmeal        46.67 22.9 65    -23.1    116.5
#>  casein - soybean         77.15 21.6 65     11.4    142.9
#>  casein - sunflower       -5.33 22.4 65    -73.6     62.9
#>  horsebean - linseed     -58.55 23.5 65   -130.1     13.0
#>  horsebean - meatmeal   -116.71 24.0 65   -189.7    -43.7
#>  horsebean - soybean     -86.23 22.7 65   -155.4    -17.0
#>  horsebean - sunflower  -168.72 23.5 65   -240.3    -97.1
#>  linseed - meatmeal      -58.16 22.9 65   -127.9     11.6
#>  linseed - soybean       -27.68 21.6 65    -93.4     38.1
#>  linseed - sunflower    -110.17 22.4 65   -178.4    -41.9
#>  meatmeal - soybean       30.48 22.1 65    -36.9     97.8
#>  meatmeal - sunflower    -52.01 22.9 65   -121.8     17.8
#>  soybean - sunflower     -82.49 21.6 65   -148.2    -16.7
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: bonferroni method for 15 estimates

Note the footer: the intervals were adjusted with the Bonferroni method, the only indication that the requested FDR adjustment was not applied to them. For example, horsebean - linseed has an FDR-adjusted p-value of about 0.023, yet its interval includes zero. A reader comparing the two tables would not know whether that difference is significant.

FCR-adjusted intervals

Benjamini and Yekutieli (2005) defined the false coverage-statement rate (FCR), the interval counterpart of the FDR: the expected proportion of intervals that fail to cover their parameter among the intervals reported for the selected (significant) parameters. They showed that if \(R\) of \(m\) hypotheses are rejected by the Benjamini-Hochberg procedure at level \(q\), constructing marginal intervals at confidence level

\[1 - \frac{R\,q}{m}\]

controls the FCR at level \(q\) when the test statistics are independent. Since the procedure rejects exactly the hypotheses with \(p_i \le Rq/m\), these intervals exclude the null value exactly for the rejected hypotheses.

confint_fcr() applies this adjustment to an ‘emmeans’ object, reading the FDR method from the adjust setting used to create it:

res <- confint_fcr(prs)
res
#>  contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
#>  casein - horsebean      163.38 23.5 65   113.28   213.49   6.957 <0.0001
#>  casein - linseed        104.83 22.4 65    57.06   152.61   4.682 <0.0001
#>  casein - meatmeal        46.67 22.9 65    -2.17    95.52   2.039  0.0570
#>  casein - soybean         77.15 21.6 65    31.12   123.19   3.576  0.0012
#>  casein - sunflower       -5.33 22.4 65   -53.11    42.44  -0.238  0.8125
#>  horsebean - linseed     -58.55 23.5 65  -108.66    -8.44  -2.493  0.0228
#>  horsebean - meatmeal   -116.71 24.0 65  -167.84   -65.58  -4.870 <0.0001
#>  horsebean - soybean     -86.23 22.7 65  -134.68   -37.78  -3.797  0.0007
#>  horsebean - sunflower  -168.72 23.5 65  -218.82  -118.61  -7.184 <0.0001
#>  linseed - meatmeal      -58.16 22.9 65  -107.01    -9.31  -2.540  0.0225
#>  linseed - soybean       -27.68 21.6 65   -73.71    18.36  -1.283  0.2187
#>  linseed - sunflower    -110.17 22.4 65  -157.94   -62.39  -4.920 <0.0001
#>  meatmeal - soybean       30.48 22.1 65   -16.67    77.63   1.379  0.1991
#>  meatmeal - sunflower    -52.01 22.9 65  -100.85    -3.16  -2.271  0.0360
#>  soybean - sunflower     -82.49 21.6 65  -128.52   -36.45  -3.823  0.0007
#> 
#> Confidence level used: 0.9633 
#> Conf-level adjustment: FCR (Benjamini & Yekutieli, 2005), q = 0.05 
#> P value adjustment: fdr (Benjamini-Hochberg) method for 15 tests

Now horsebean - linseed has an interval that excludes zero, matching its p-value, and casein - meatmeal (FDR-adjusted p = 0.057) has an interval that includes zero. The footer reports the confidence level used, \(1 - 11 \times 0.05 / 15 \approx 0.963\). The details are available with fcr_info():

fcr_info(res)
#>   family  m  R     level
#> 1    all 15 11 0.9633333

Because the intervals are computed by ‘emmeans’ itself (with adjust = "none" and the FCR level), degrees of freedom, back-transformations and one-sided intervals behave exactly as in ‘emmeans’. The common shortcut emmeans(fit, pairwise ~ feed, adjust = "fdr") also works: confint_fcr() then uses its contrasts component.

Ratios and families defined by by variables

The pigs data from ‘emmeans’ records the concentration of leucine in pigs fed three protein sources at four protein percentages. On the log scale, differences become ratios after back-transformation. Here each percentage is a separate family of three comparisons, as in ‘emmeans’:

fit2 <- lm(log(conc) ~ source * factor(percent), data = pigs)
prs2 <- pairs(emmeans(fit2, ~ source | percent), adjust = "fdr")
res2 <- confint_fcr(prs2, type = "response")
res2
#> percent = 9:
#>  contrast    ratio     SE df lower.CL upper.CL null t.ratio p.value
#>  fish / soy  0.746 0.0792 17    0.583    0.954    1  -2.763  0.0200
#>  fish / skim 0.730 0.0774 17    0.571    0.933    1  -2.971  0.0200
#>  soy / skim  0.978 0.0929 17    0.785    1.219    1  -0.233  0.8188
#> 
#> percent = 12:
#>  contrast    ratio     SE df lower.CL upper.CL null t.ratio p.value
#>  fish / soy  0.779 0.0740 17    0.625    0.971    1  -2.629  0.0264
#>  fish / skim 0.716 0.0679 17    0.574    0.891    1  -3.526  0.0078
#>  soy / skim  0.918 0.0872 17    0.737    1.144    1  -0.898  0.3818
#> 
#> percent = 15:
#>  contrast    ratio     SE df lower.CL upper.CL null t.ratio p.value
#>  fish / soy  0.791 0.0839 17    0.632    0.989    1  -2.212  0.0409
#>  fish / skim 0.625 0.0727 17    0.489    0.799    1  -4.043  0.0025
#>  soy / skim  0.790 0.0839 17    0.632    0.989    1  -2.216  0.0409
#> 
#> percent = 18:
#>  contrast    ratio     SE df lower.CL upper.CL null t.ratio p.value
#>  fish / soy  0.753 0.1010 17    0.527    1.076    1  -2.112  0.0594
#>  fish / skim 0.540 0.0725 17    0.378    0.772    1  -4.587  0.0008
#>  soy / skim  0.717 0.1180 17    0.464    1.110    1  -2.020  0.0594
#> 
#> Confidence level used per family: 
#>   percent = 9: 0.9667 
#>   percent = 12: 0.9667 
#>   percent = 15: 0.95 
#>   percent = 18: 0.9833 
#> Conf-level adjustment: FCR (Benjamini & Yekutieli, 2005), q = 0.05 
#> P value adjustment: fdr (Benjamini-Hochberg) method for 3 tests 
#> Intervals are back-transformed from the log scale 
#> Tests are performed on the log scale

Each family receives its own confidence level, depending on how many of its hypotheses were rejected. An interval for a ratio excludes 1 exactly when the FDR-adjusted p-value is at most 0.05. When a family has no rejections, the level is computed with \(R = 1\) (the Bonferroni level), so that no interval excludes the null value. Benjamini and Yekutieli construct no intervals in that case; this fallback is a convention of the package that keeps the table complete.

fcr_info(res2)
#>         family m R     level
#> 1  percent = 9 3 2 0.9666667
#> 2 percent = 12 3 2 0.9666667
#> 3 percent = 15 3 3 0.9500000
#> 4 percent = 18 3 1 0.9833333

To treat all 12 comparisons as a single family, use by = NULL, as you would in summary():

fcr_info(confint_fcr(prs2, type = "response", by = NULL))
#>   family  m R level
#> 1    all 12 6 0.975

Dependent contrasts: the Benjamini-Yekutieli version

The FCR guarantee of Benjamini and Yekutieli (2005) for the level above assumes independent test statistics. Contrasts estimated from the same model are usually correlated, and t statistics that share a variance estimate are dependent even when the estimates are not. If a guarantee under any dependence structure is required, use the Benjamini-Yekutieli procedure, either through adjust = "BY" in ‘emmeans’ or with method = "BY":

fcr_info(confint_fcr(prs, method = "BY"))
#>   family  m R     level
#> 1    all 15 8 0.9919636

The level now includes the factor \(\sum_{i=1}^{m} 1/i\), which makes both the p-values and the intervals more conservative. The intervals and p-values still agree with each other.

What the intervals mean

  • For the significant contrasts, the intervals control the FCR (on average, at most a proportion \(q\) of them fail to cover the true value) when the test statistics are independent ("BH") or under any dependence ("BY"). For correlated contrasts with "BH", this is not a formal guarantee.
  • Intervals are also shown for the non-significant contrasts, so that the table is complete and consistent with the p-values. The FCR criterion makes no coverage statement about these intervals.
  • A single, unadjusted 95% interval should not be read from this table; use confint(prs, adjust = "none") for that.

Using the level elsewhere

fcr_level() works with any vector of p-values, so the same adjustment can be applied to results from other packages, as long as the intervals and p-values are based on the same test statistic:

p <- c(0.0004, 0.003, 0.012, 0.04, 0.21, 0.67)
fcr_level(p, q = 0.05)
#> [1] 0.975

References

Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B, 57(1), 289-300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

Benjamini, Y. and Yekutieli, D. (2005). False discovery rate-adjusted multiple confidence intervals for selected parameters. Journal of the American Statistical Association, 100(469), 71-81. https://doi.org/10.1198/016214504000001907

Jung, K., Friede, T. and Beissbarth, T. (2011). Reporting FDR analogous confidence intervals for the log fold change of differentially expressed genes. BMC Bioinformatics, 12, 288. https://doi.org/10.1186/1471-2105-12-288