2

我已经使用广义线性模型分析了一组数据，该模型在三向交互中具有三个分类因子（因子A、因子B、因子C）和简单地添加到模型中的第四个连续因子（因子D）。我正在尝试从模型中获取一组 Tukey 字母组（即紧凑型字母显示），但还没有找到成功包含交互的方法。我对包含因子D 不感兴趣，只是交互中的三个。

我已经得到了 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 函数调用的因子字段中，我会得到同样的错误。我检查了 inxtns 对象，它看起来不错，并且行数与残差数相匹配。agricolae注释实际上并未涵盖 HSD.test 中交互的使用，但我认为它可以工作。

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

我已经为此工作了好几天，但一直未能找到解决方案，希望我没有遗漏一些明显的东西。

谢谢！

4

1 回答 1

2

我不知道您是如何指定 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
于 2013-08-01T02:11:22.097 回答