如何从具有交互的GLM获得Tukey紧凑字母显示

Question

我用一个广义线性模型分析了一组数据，这个模型在三向交互（factorA，factorB，factorC）和第四个连续因子（factorD）中有三个分类因素，它们只是简单地添加到模型中。我试图从模型中获取一组Tukey字母组（即紧凑字母显示），但尚未找到成功包含交互的方法。我对包括factorD只是三个在交互中不感兴趣。

我已经得到了Tukey调整的成对比较：

lsmeans(my.glm, factorA*factorB*factorC)

但我无法弄清楚如何制作一个紧凑的字母显示器。它可以使用multcomp包来完成，但我只能找到使用该包进行主效应的方法，而不是交互。

那么我尝试了agricolae包，因为这篇文章（https://stats.stackexchange.com/questions/31547/how-to-obtain-the-results-of-a-tukey-hsd-post-hoc-test-in-a-table-showing-groupe）讨论了应该有效。但是，按照该答案中的说明导致HSD.test的非功能性响应。具体来说，我可以让主要效果测试工作正常，例如HSD.test(my.glm,"factorA")但我无法让互动发挥作用。我试过这个：

intxns<-with(my.data, interaction(factorA,factorB,factorC))
HSD.test(my.glm,"intxns",group=TRUE)

但是得到一个错误，表明HSD.test函数没有将“intxns”识别为有效对象，它看起来像这样（同样，我检查了intxns对象，它看起来很好，行数与残差数相匹配）我的glm）：

Name: inxtns
 factorA factorB factorC factorD

如果我只是将废话放入HSD.test函数调用中的factor字段，我会得到同样的错误。我检查了inxtns对象，它看起来很好，行数与残留数量匹配.agricolae笔记实际上并不涵盖HSD.test中交互的使用，但我认为它可以工作。

有谁知道如何让HSD.test与交互一起工作？或者，您是否还有其他功能可以为具有交互作用的glm生成紧凑的字母显示？

我已经在这个问题上工作了很多天并且找不到解决办法，希望我没有错过任何明显的东西。

谢谢！

Answer 1

我不知道你是如何指定你的glm模型的，但是对于HSD.test，它希望将特定的治疗名称与glm公式中指定的相同名称以及数据框相匹配。这就是为什么你的主要效果factorA会起作用，而不是三方互动。对于交互的多个比较测试，我发现最简单的方法是单独生成交互并将它们作为附加列添加到数据框中。然后可以使用编码交互的新变量来指定glm模型。

例如，

set.seed(42)
glm.dat <- data.frame(y = rnorm(1000), factorA = sample(letters[1:2], 
   size = 1000, replace = TRUE),
 factorB = sample(letters[1:2], size = 1000, replace = TRUE),
 factorC = sample(letters[1:2], size = 1000, replace = TRUE))

# Generate interactions explicitly and add them to the data.frame
glm.dat$factorAB <- with(glm.dat, interaction(factorA, factorB))
glm.dat$factorAC <- with(glm.dat, interaction(factorA, factorC))
glm.dat$factorBC <- with(glm.dat, interaction(factorB, factorC))
glm.dat$factorABC <- with(glm.dat, interaction(factorA, factorB, factorC))

# General linear model 
 glm.mod <- glm(y ~ factorA + factorB + factorC +  factorAB + factorAC + 
   factorBC + factorABC, family = 'gaussian', data = glm.dat) 

# Multiple comparison test

library(agricolae)
comp <- HSD.test(glm.mod, trt = "factorABC", group = TRUE)

给

comp$groups giving

    trt        means M
 1 a.a.a  0.070052189 a
 2 a.b.b  0.035684571 a
 3 b.a.a  0.020517535 a
 4 b.b.b -0.008153257 a
 5 a.b.a -0.036136140 a
 6 a.a.b -0.078891136 a
 7 b.a.b -0.080845419 a
 8 b.b.a -0.115808772 a