1

在我的工作中，我经常不得不使用 Anova 和 Tukey 检验进行不同的治疗比较，以确定一个因素实验中的多个治疗中的哪一个在统计上彼此不同。

我附上的代码产生了两个单独的数字：一个带有值的处理分布（示例 graph1），另一个带有 Tukey 测试结果，显示哪一对处理彼此之间存在显着差异（示例 graph2）。

我过去所做的是查看 Tukey 结果并手动编辑第一个图表，其中包含表示统计等效组组的字母（例如 graph3）。我一直在研究不同的 r 库，以寻找自动生成类似于图 3 的东西的方法，该图总结了这些分组，但还没有找到方法。有没有人有什么建议？

PS- 如果下面的图形例程有点麻烦，我很抱歉，但它本质上是我开发的一组更全面的函数的一个片段，用于测试数据分布、有条件地应用相关测试并生成输出表格和图形。

我为制作前两个图表而编写的代码如下。t?usp=分享

Group=c("G1","G1","G1","G1","G2","G2","G2","G2","G3","G3","G3","G3")
Vals=c(runif(4),runif(4)+0.5,runif(4)+0.1)
data=data.frame(Group)
data=cbind(data, Vals)  
anova_results=aov(Vals~Group,data=data)
anova_results2=anova(anova_results)[1, ]
anova_significance=anova_results2[1,"Pr(>F)"]
significant=anova_significance[1]<=0.05
if (significant==1) {
  Tukey_results=TukeyHSD(anova_results,"Group")
  Tukey_results=Tukey_results$Group
}  
plot(data$Group, data$Vals) 
if (significant==1) {
  plot(TukeyHSD(anova_results,"Group"), las=1)   
}
4

1 回答 1

1

Roman Lustrik 对上述评论的建议是正确的。我最终找到了两种基于两个相关库的替代方法。运行问题中发布的代码后，创建分组处理图运行：

#simple looking solution
library(multcomp)
tuk <- glht(anova_results, linfct = mcp(Group = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))
plot(tuk.cld)
par(opar)

#more fancy looking solution using the multcompView library with a lot of ways to 
#alter the plot appearance as necessary
library(multcompView)
multcompBoxplot(Vals~Group,data=data)

# Now, the solution below is my favorite solution as the text direction of the groups 
#work very well if you have many treatments that you are comparing
opar <- par()  
par(oma = c(6, 0, 0, 0)) #extra space for extra large treatment names
xzx <-multcompBoxplot(Vals~Group,data=data,sortFn=median,  decreasing=FALSE,    
                      horizontal=FALSE,
                      plotList=list(
                        boxplot=list(fig=c(0,  1,  0,  1),  las=3,
                                     cex.axis=1.5),                      
                        multcompLetters=list(
                          fig=c(0.87,  0.97,  0.115,  0.923), #0.1108, 0.9432 Top of     
#page 18 manual for very convoluted explanation (c(y bottom, y top,x L, x R))
                          type='Letters') ) )
par(opar)
于 2013-02-10T20:03:00.543 回答