如何通过 R 中的交互进行 LSD 测试？

如何解决如何通过 R 中的交互进行 LSD 测试？

我有一个字段数据

sowing_date<- rep(c("Early","normal"),each=12)
herbicide<- rep (rep(c("No","Yes"),each=6),2)
nitrogen<- rep (rep(c("No",each=3),4)
Block<- rep(c("Block 1","Block 2","Block 3"),times=8)
Yield<- c(30,27,25,40,41,42,37,38,48,47,46,26,39,57,59,60)

DataA<- data.frame(sowing_date,herbicide,nitrogen,Block,Yield)

我进行了三向方差分析

anova3way <- aov (Yield ~ sowing_date + herbicide + nitrogen + 
            sowing_date:herbicide + sowing_date:nitrogen + 
            herbicide:nitrogen + sowing_date:herbicide:nitrogen + 
           factor(Block),data=DataA)
summary(anova3way)

3 个因素之间存在 3 种相互作用。所以，我想看看哪种组合显示出最大的收益。

我知道如何将平均差异与如下所示的单个因素进行比较，但在交互的情况下，我该怎么做？

library(agricolae)
LSD_Test<- LSD.test(anova3way,"sowing_date")
LSD_Test

例如，我想检查三向交互作用下的平均差异，以及每两个因素之间的交互作用。

例如，我想在 R 中得到这个 LSD 结果

你能告诉我我该怎么做吗？

非常感谢，

解决方法

需要一些手动工作的一种方法是将实验参数编码为 -1 和 1，以便正确分离 2 和 3 参数交互。
对值进行编码后，您可以从方差分析模型中提取剩余自由度和误差平方和，并将其传递给 LSD.test 函数。

见下面的例子：

sowing_date<- rep(c("Early","Normal"),each=12)
herbicide<- rep (rep(c("No","Yes"),each=6),2)
nitrogen<- rep (rep(c("No",each=3),4)
Block<- rep(c("Block 1","Block 2","Block 3"),times=8)
Yield<- c(30,27,25,40,41,42,37,38,48,47,46,26,39,57,59,60)
DataA<- data.frame(sowing_date,herbicide,nitrogen,Block,Yield)

anova3way <- aov (Yield ~ sowing_date * herbicide * nitrogen + 
                 factor(Block),data=DataA)
summary(anova3way)

#Encode the experiment's parameters as -1 and 1
DataA$codeSD <- ifelse(DataA$sowing_date == "Early",-1,1)
DataA$codeherb <- ifelse(DataA$herbicide == "No",1)
DataA$codeN2 <- ifelse(DataA$nitrogen == "No",1)

library(agricolae)

LSD_Test<- LSD.test(anova3way,c("sowing_date"))
LSD_Test

#Manually defining the treatment and specifying the
# degrees of freedom and Sum of the squares (Frin the resduals from the ANOVA)
print(LSD.test(y=DataA$Yield,trt=DataA$sowing_date,DFerror=14,MSerror=34.3))

#Example for a two parameter value
print(LSD.test(y=DataA$Yield,trt=(DataA$codeSD*DataA$codeherb),MSerror=34.3))

print(LSD.test(y=DataA$Yield,trt=(DataA$codeSD*DataA$codeherb*DataA$codeN2),MSerror=34.3))


#calaculate the means and std (as a check)
#DataA %>% group_by(sowing_date) %>% summarize(mean=mean(Yield),sd=sd(Yield))
#DataA %>% group_by(codeSD*codeherb*codeN2) %>% summarize(mean=mean(Yield),sd=sd(Yield))

您需要手动跟踪最终报告中与 -1 和 1 对应的运行/条件。

编辑： 所以我上面的回答显示了基于交互的整体效果。例如，除草剂和氮肥的相互作用如何影响产量。

根据您想要确定哪种组合提供最大收益的评论，您再次使用 LSD.test() 函数，但传递参数名称向量。

LSD_Test<- LSD.test(anova3way,c("sowing_date","herbicide","nitrogen"))
LSD_Test

从输出的组部分，您可以看到“正常”、“是”和“是”是最佳产量。 “组”列将标识唯一的结果集群。例如，最后 2 行提供了类似的产量。

...
$groups
                  Yield groups
Normal:Yes:Yes 58.66667      a
Early:Yes:Yes  47.00000      b
Normal:No:Yes  41.33333      c
Early:No:Yes   41.00000     cd
Normal:Yes:No  39.66667     cd
Early:Yes:No   38.33333      d
Early:No:No    27.33333      e
Normal:No:No   26.00000      e
...