具有三明治估计的回归的成对比较（在 R 中）

机器算法验证 回归 数据可视化 p 值 箱形图 三明治
2022-04-10 19:08:14

简短的问题

我在 R 中运行回归并制作了响应变量的箱线图，并按其中一个预测变量分组。在这个箱线图上，我想添加一些关于统计模型的信息。你会建议我提供什么信息（以及如何显示它（这不是编程问题））？

展开的问题

我有几个预测变量：两个分类的、非序数的预测变量和一个连续的预测变量（下面用 R 编码）

set.seed(81)
pred1 = rep(c('Car', 'Bike', 'Train', 'Airplane'), 6)
pred2 = rep(c('High', 'Low', 'Middle'), 8)
pred3 = rnorm(24)
resp = c(rnorm(12, sd = 1), rnorm(12, sd = 5))

resp是响应变量。我用三明治估计进行了回归：

require(sandwich)
require(lmtest)


m = aov(resp ~ pred1 + pred2)
coeftest(m, sandwich)

t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.49642    0.73911 -0.6716  0.51034  
pred1Bike    1.55917    1.16568  1.3376  0.19769  
pred1Car     1.23873    1.24080  0.9983  0.33135  
pred1Train   2.50882    0.91468  2.7428  0.01338 *
pred2Low     0.11613    1.00540  0.1155  0.90932  
pred2Middle  0.51476    0.90924  0.5661  0.57829  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我为以下组绘制了箱线图pred1：

require(ggplot2)
ggplot(data.frame(pred1, resp), aes(x=pred1, y=resp)) + geom_boxplot()

在此处输入图像描述

在此图上，我想添加一些字母以指示统计上相似的组（p.value < 0.05），如此处所述。像这样的东西：

ggplot(data.frame(pred1, resp), aes(x=pred1, y=resp)) + geom_boxplot() + annotate('text', x=1:4, y=6, label=c('a','b','a','b'), size = 8, color='red')

在此处输入图像描述

我的问题是：

我怎样才能找到这些 p.values 与我的稳健回归进行成对比较？我可以做以下m简单的 aov 模型：

TukeyHSD(m)

但以下不起作用：

TukeyHSD(coeftest(m, sandwich))

我可能会误解这些成对比较是什么，以及我目前的结果是什么意思！如果您有这种感觉，请告诉我！我的问题的目的是让我了解在箱线图上显示统计模型结果的最佳方式是什么。

注意：变量pred2和pred3用于提取我不希望影响的方差的某些部分pred1（如pred1，pred2并且pred3在我的情况下是相关的）。因此，我想最好不要运行简单的成对 t 检验来获得我想在每个箱线图顶部添加的 p.values。

1个回答

multcomp一个解决方案实际上是在关于包的书第 4.6 节中作为示例给出的：

Bretz, F.、Hothorn, T. 和 Westfall, PH (2011)。使用 R 进行多重比较。佛罗里达州博卡拉顿：CRC 出版社。

只需稍微调整您的代码（所有内容都需要在一个 data.frame 中，而不是四处浮动）：

require(multcomp)
require(sandwich)

set.seed(81)
pred3 = rnorm(24)
df <- data.frame(pred1 = rep(c('Car', 'Bike', 'Train', 'Airplane'), 6), pred2 = rep(c('High', 'Low', 'Middle'), 8)
, resp = c(rnorm(12, sd = 1), rnorm(12, sd = 5)))

m <- aov(resp ~ pred1 + pred2, df)

tukey <- glht(m, linfct = mcp(pred1 = "Tukey") , vcov = sandwich)

summary(tukey, test = adjusted())

##          Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = resp ~ pred1 + pred2, data = df)
## 
## Linear Hypotheses:
##                       Estimate Std. Error t value Pr(>|t|)  
## Bike - Airplane == 0     1.559      1.166    1.34    0.547  
## Car - Airplane == 0      1.239      1.241    1.00    0.748  
## Train - Airplane == 0    2.509      0.915    2.74    0.058 .
## Car - Bike == 0         -0.320      1.422   -0.23    0.996  
## Train - Bike == 0        0.950      1.149    0.83    0.838  
## Train - Car == 0         1.270      1.225    1.04    0.726  
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## (Adjusted p values reported -- single-step method)

请注意，glht默认情况下使用单步法作为 alpha 误差累积的调整。如果你想要别的东西，你需要把它喂给adjusted()

例如，使用 Bonferroni-Holm 校正（我倾向于使用它，因为我不明白单步实际上做了什么）：

summary(tukey, test = adjusted("holm"))

如果你不想要我不推荐的 alpha 错误更正，这也是可能的：

summary(tukey , test = adjusted("none"))
其它你可能感兴趣的问题
上一篇如何估计反双曲正弦变换的 theta？ 下一篇指数分布的拟合优度检验