使用multcomp::glht()生成K矩阵，进行交互式Tukey分析的2路方差分析。

Question

我想知道如何通过multcomp::glht()函数实现广义线性模型的同步推理。具体来说，我想对来自以下的2路方差分析进行完整的Tukey分析 此处. 他们对有交互作用的模型进行了Tukey的事后分析。

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

• 羊毛是一个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应该是怎样的，所以我可以手动计算出来。然而，这只是一个例子来熟悉这个包。我的真正的方差分析有一个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矩阵。有谁知道怎么做？

谢谢大家