R:当交互作用不显著时，如何显示2向方差分析的均值分离字母？

Question

我正在比较一个肥料实验，其中我有一个响应变量(生长速率)和两个独立变量(基因型和速率)。我在R中使用lsmeans、car和multcompView软件包进行了双向方差分析。我的交互效应不显著，但我的主要效应变量基因型和比率显著。我希望查看表示均值分离的差异的字母。我该怎么做呢？

以下是我尝试过的方法：

数据集

demo <- data.frame(genotype=c(1,
                              1,
                              1,
                              1,
                              2,
                              2,
                              2,
                              2,
                              3,
                              3,
                              3,
                              3,
                              2,
                              2,
                              1,
                              3,
                              3,
                              3,
                              1,
                              3,
                              2,
                              2,
                              1,
                              1,
                              1,
                              2,
                              2,
                              1,
                              3,
                              3,
                              3,
                              3,
                              2,
                              1,
                              2,
                              1),
                   rate = c(0,
                            1,
                            2,
                            3,
                            0,
                            1,
                            2,
                            3,
                            0,
                            1,
                            2,
                            3,
                            1,
                            3,
                            0,
                            0,
                            2,
                            1,
                            2,
                            3,
                            0,
                            2,
                            3,
                            1,
                            1,
                            0,
                            2,
                            2,
                            3,
                            0,
                            2,
                            1,
                            3,
                            0,
                            1,
                            3),
                   response = c(0.636008,
                                0.520401,
                                0.511821,
                                0.200255,
                                0.535433,
                                0.272521,
                                0.192943,
                                0.238736,
                                0.568422,
                                0.497654,
                                0.433995,
                                0.316064,
                                0.336663,
                                0.187805,
                                0.696257,
                                0.517592,
                                0.244133,
                                0.469655,
                                0.277319,
                                0.188577,
                                0.534307,
                                0.204349,
                                0.263632,
                                0.651846,
                                0.46279,
                                0.499629,
                                0.150601,
                                0.168777,
                                0.227221,
                                0.518858,
                                0.35837,
                                0.516537,
                                0.221283,
                                0.753765,
                                0.446882,
                                0.301356))

demo$genotype <- as.factor(demo$genotype)
demo$rate <- as.factor(demo$rate)


#Load packages
library(car)
library(multcompView)
library(lsmeans)

我之前检查了交互术语，它们并不重要。下面是主要的效果模型：

mod1 <- lm(response ~ genotype + rate, data = demo)
Anova(mod1, type = "III")

#Results
Response: response
             Sum Sq Df  F value    Pr(>F)    
(Intercept) 2.50289  1 389.5155 < 2.2e-16 ***
genotype    0.11256  2   8.7589  0.001009 ** 
rate        0.70042  3  36.3345 4.106e-10 ***
Residuals   0.19277 30                       
---

看一下基因型均值之间的差异：

genotype <- lsmeans(mod1, pairwise ~ genotype)
genotype

#Results:
Results are averaged over the levels of: rate 
Confidence level used: 0.95 

$contrasts
 contrast estimate     SE df t.ratio p.value
 1 - 2      0.1353 0.0327 30  4.133  0.0008 
 1 - 3      0.0489 0.0327 30  1.495  0.3075 
 2 - 3     -0.0863 0.0327 30 -2.638  0.0340 

Results are averaged over the levels of: rate 
P value adjustment: tukey method for comparing a family of 3 estimates

rate的过程相同：

rate <- lsmeans(mod1, pairwise ~ rate)
rate

#Results:
Results are averaged over the levels of: genotype 
Confidence level used: 0.95 

$contrasts
 contrast estimate     SE df t.ratio p.value
 0 - 1      0.1206 0.0378 30 3.191   0.0165 
 0 - 2      0.3020 0.0378 30 7.992   <.0001 
 0 - 3      0.3461 0.0378 30 9.160   <.0001 
 1 - 2      0.1814 0.0378 30 4.801   0.0002 
 1 - 3      0.2256 0.0378 30 5.969   <.0001 
 2 - 3      0.0442 0.0378 30 1.168   0.6509 

Results are averaged over the levels of: genotype 
P value adjustment: tukey method for comparing a family of 4 estimates

现在，我尝试通过基因型输出来获取字母，但得到了一条错误消息：

cld(genotype, 
    alpha = 0.05, 
    Letters = letters)


#Error message thrown here:
Error in cld(genotype, alpha = 0.05, Lettes = letters) : 
  could not find function "cld"

有没有更好的方法来获得分隔符，或者我只是犯了一个简单的错误？

Answer 1

解决方案是使用agricolae包：

install.packages("agricolae")
library(agricolae)

mod1 <- lm(response ~ genotype + rate, data = demo)
Anova(mod1, type = "III")

genotype <- HSD.test(mod1, "genotype")
genotype

#Part of the output:
response groups
1 0.4536856      a
3 0.4047565      a
2 0.3184293      b

rate <- HSD.test(mod1, "rate") 
rate

#Part of the output:
response groups
0 0.5844746      a
1 0.4638832      b
2 0.2824787      c
3 0.2383254      c