emmeans包多重比较误差校正的应用时机与设置方法

我最初在Cross Validated发布了该问题，但由于内容纯为软件语法相关，认为更适合在Stack Overflow提问，本问题是相关帖子的后续。我构建了多项logistic回归，分析受访者使用合法（licit）/非法（illicit）医用大麻治疗不同病症（疼痛pain、睡眠问题sleep、心理健康/物质使用问题mhsu、其他病症allOther）的对数几率差异。

演示数据 #

df <- tibble(mcType = factor(rep(c("licit", "illicit"),
                                 times = c(534,1207))),
             cond = factor(c(rep(c("pain","mhsu","allOther","sleep"), 
                                 times = c(280,141,82,31)),
                             rep(c("pain","mhsu","allOther","sleep"), 
                                 times = c(491,360,208,148))),
                           levels = c("pain","sleep","mhsu","allOther")))

不同大麻类型对应的各病症报告比例 #

mcType  cond         n   tot  perc
<fct>   <fct>    <int> <int> <dbl>
1 illicit pain       491  1207 40.7 
2 illicit sleep      148  1207 12.3 
3 illicit mhsu       360  1207 29.8 
4 illicit allOther   208  1207 17.2 
5 licit   pain       280   534 52.4 
6 licit   sleep       31   534  5.81
7 licit   mhsu       141   534 26.4 
8 licit   allOther    82   534 15.4 

多项logistic回归拟合代码及输出 #

我使用nnet包的multinom()函数拟合了模型，代码及输出如下：

library(nnet)
summary(mm <- multinom(cond ~ mcType,
                       data = df))

输出结果：

Coefficients:
  (Intercept) mcTypelicit
sleep     -1.1992431  -1.0014884
mhsu      -0.3103369  -0.3756443
allOther  -0.8589398  -0.3691759

Std. Errors:
  (Intercept) mcTypelicit
sleep     0.09377333   0.2112368
mhsu      0.06938587   0.1244098
allOther  0.08273132   0.1503720

Residual Deviance: 4327.814 
AIC: 4339.814 

emmeans简单效应检验过程 #

参考相关使用说明，emmeans默认会应用误差校正，可通过adjust =参数指定校正方法。我先在构建emmeans对象时指定adjust="bonferroni"，代码及输出如下：

# 检验不同病症下大麻类型的效应，先构建emmeans对象
library(emmeans)
(em_mcTypeByCond <- emmeans(object = mm,
                            specs = ~mcType|cond,
                            adjust = "bonferroni"))

输出结果：

cond = pain:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.4068 0.01414  6   0.3648   0.4488
 licit   0.5243 0.02161  6   0.4602   0.5885

cond = sleep:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.1226 0.00944  6   0.0946   0.1506
 licit   0.0581 0.01012  6   0.0280   0.0881

cond = mhsu:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.2983 0.01317  6   0.2592   0.3374
 licit   0.2641 0.01908  6   0.2074   0.3207

cond = allOther:
 mcType    prob      SE df lower.CL upper.CL
 illicit 0.1723 0.01087  6   0.1401   0.2046
 licit   0.1535 0.01560  6   0.1072   0.1999

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates

后续我尝试使用其他误差校正方法（如"BH"、"fdr"、"westfall"、"holm"等）时无法生效，怀疑是校正应用时机错误，即在检验前就指定了校正参数。因此我尝试在pairs()函数中指定adjust参数，检验不同大麻类型在各病症下的概率差异，代码及输出如下：

(mcTypeByCond_test <- pairs(em_mcTypeByCond,
                            adjust = "bonferroni"))

输出结果：

cond = pain:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit  -0.1175 0.0258  6 -4.551  0.0039 

cond = sleep:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit   0.0646 0.0138  6  4.665  0.0034 

cond = mhsu:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit   0.0342 0.0232  6  1.476  0.1905 

cond = allOther:
 contrast        estimate     SE df t.ratio p.value
 illicit - licit   0.0188 0.0190  6  0.987  0.3616 

核心疑问 #

当前输出中没有提示使用了何种误差校正方法，且我需要对全部4组两两比较的误差进行统一控制。请问p值校正的参数需要在哪个阶段、如何设置才能生效？

内容的提问来源于stack exchange，提问作者llewmills