用multcomp::glht()交互Tukey分析生成双向方差K矩阵

Question

我正在尝试通过multcomp::glht()函数来实现通用线性模型的同时推理。具体来说，我想对来自这里的双向方差分析进行一个完整的Tukey分析。他们对具有交互作用的模型进行图基的事后分析。

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

• 羊毛是一个2级因子(A，B)
• 张力是一个三级因子(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
• 甲:米-乙:h == 0

我知道对比矩阵K应该是怎样的，所以我可以手动计算出来。不过，这只是一个熟悉包的例子。我的实际变异系数有5级因子和另外10级因子，所以手动操作将是一种痛苦。谢谢

Answer 1

我已经通过电子邮件包找到了解决办法。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)

然而，在做了这么多的小片段之后，我仍然对用基或glht工具找到K矩阵感到好奇。有人知道怎么做吗？

谢谢