使用 multcomp::glht() 通过交互 Tukey 分析为 2 路方差分析生成 K 矩阵

Question

我试图弄清楚如何通过 multcomp::glht() 函数实现广义线性模型的同时推理。具体来说，我想对来自here 的 2-way ANOVA 示例进行完整的 Tukey 分析。他们对具有交互作用的模型进行 Tukey 事后分析。

mod <- lm(breaks ~ wool * tension, data = warpbreaks)

• wool 是 2 级因子 (A, B)
• 张力是一个 3 级因子（L、M、H）

在小插图的示例中，他们检查了每个羊毛级别内张力级别的平均值差异。但是我有兴趣学习如何寻找每个可能的级别组合之间的差异。 按照示例，使用如下代码：

tmp <- expand.grid(tension = unique(warpbreaks$tension),
                    wool = unique(warpbreaks$wool))
X <- model.matrix(~ wool * tension, data = tmp)
glht(mod, linfct = X)
Tukey <- contrMat(table(warpbreaks$tension), "Tukey")
K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(warpbreaks$wool)[1], rownames(K1), sep = ":")
K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$wool)[2], rownames(K2), sep = ":")
K <- rbind(K1, K2)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))
summary(glht(mod, linfct = K %*% X))

他们得到类似下面的东西：

Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = breaks ~ wool * tension, data = warpbreaks)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)   
A:M - L == 0 -20.5556     5.1573  -3.986  0.00131 **
A:H - L == 0 -20.0000     5.1573  -3.878  0.00187 **
A:H - M == 0   0.5556     5.1573   0.108  0.99996   
B:M - L == 0   0.5556     5.1573   0.108  0.99996   
B:H - L == 0  -9.4444     5.1573  -1.831  0.30795   
B:H - M == 0 -10.0000     5.1573  -1.939  0.25535   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

那么，我怎样才能得到正确的对比度矩阵，以便进行如下比较：

• A:M - B:M == 0
• A:M - B:H == 0

我知道对比度矩阵 K 应该如何，所以我可以手动计算出来。然而，这只是一个熟悉这个包的例子。我真正的 ANOVA 有一个 5 级因子和另一个 10 级因子，所以手动操作会很痛苦。 谢谢

Answer 1

我通过 emmeans 包找到了解决方案。函数 lsm 似乎可以满足我的要求。特别是：

library(emmeans)        # for lsm
model_glht <- glht(mod1, lsm(pairwise ~ wool : tension), by = NULL)
summary(model_glht)

我终于得到了想要的比较：

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = breaks ~ wool * tension, data = warpbreaks)

Linear Hypotheses:
               Estimate Std. Error t value Pr(>|t|)    
A,L - B,L == 0  16.3333     6.4767   2.522   0.1285    
A,L - A,M == 0  20.5556     6.3052   3.260   0.0216 *  
A,L - B,M == 0  15.7778     6.4135   2.460   0.1465    
A,L - A,H == 0  20.0000     6.5399   3.058   0.0369 *  
A,L - B,H == 0  25.7778     5.8918   4.375   <0.001 ***
B,L - A,M == 0   4.2222     4.1239   1.024   0.9002    
B,L - B,M == 0  -0.5556     4.2877  -0.130   1.0000    
B,L - A,H == 0   3.6667     4.4746   0.819   0.9591    
B,L - B,H == 0   9.4444     3.4589   2.730   0.0809 .  
A,M - B,M == 0  -4.7778     4.0239  -1.187   0.8294    
A,M - A,H == 0  -0.5556     4.2225  -0.132   1.0000    
A,M - B,H == 0   5.2222     3.1261   1.671   0.5384    
B,M - A,H == 0   4.2222     4.3826   0.963   0.9210    
B,M - B,H == 0  10.0000     3.3391   2.995   0.0431 *  
A,H - B,H == 0   5.7778     3.5759   1.616   0.5738    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

然而，在对小插图进行了这么多工作之后，我仍然对通过 base 或 glht 工具找到 K 矩阵感到好奇。有人知道怎么做吗？

谢谢