何时以及如何在 emmeans 包中应用纠错

When and how to apply error correction in emmeans package

我最初是在 cross--validated 上发布的,但我认为它可能更适合 SO,因为它纯粹是关于软件语法的。

这是 this post 的后续问题。我 运行 一项多项逻辑回归检验了受访者的对数比值差异,表明他们治疗了 运行ge 不同的医疗条件(painsleep、mental-health/substance 将 (mhsu) 和所有其他条件 (allOther)) 与 licitillicit 医用大麻一起使用。

这是玩具数据

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 

为了查看受访者根据他们报告的大麻类型表明每种情况的相对比例是否存在差异,我使用 运行 多项式逻辑回归 multinom() nnet包。下面输出,

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


# output
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 

I 运行 简单效果测试,使用 emmeans 包。在 this blog post 中,作者建议 emmeans 包默认应用纠错,但您可以通过 adjust = 参数控制它。

# testing effect of mc type at each level of condition. first create emmeans object
library(emmeans)
(em_mcTypeByCond <- emmeans(object = mm,
                            specs = ~mcType|cond,
                            adjust = "bonferroni"))

# output  
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() 函数中应用调整参数(测试两种​​大麻之间每种情况的概率差异)

(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 

但是如您所见,它没有提供任何消息告诉我应用了哪种类型的纠错(我假设 none,并尝试了几种不同的方法)。我还想控制所有四个成对比较的误差。

所以我需要知道我需要如何以及在什么阶段提出指定调整 p 值的参数。

非常感谢任何帮助

P 值调整应用于每个 by 组,并且每个组中只有一个比较 - 因此没有多重性。不做调整时不显示调整的注释。

要对所有结果应用调整,您需要在显示结果时从考虑中删除 by 变量:

summary(pairs(...), by = NULL, adjust = "bonf")