1

私の仕事では、分散分析とテューキー検定を使用してさまざまな治療の比較を行い、1つの因子実験の複数の治療のどれが互いに統計的に異なるかを判断する必要があります。

私が添付したコードは、2つの別々の図を生成します。1つは値の処理分布（グラフ1の例）で、もう1つは処理のどのペアが互いに有意に異なるかを示すテューキー検定結果（グラフ2の例）です。

私が過去に行ったことは、テューキーの結果を見て、統計的に同等のグループのグループを示す文字で最初のグラフを手動で編集することです（グラフ3の例）。そのようなグループ化を要約したグラフ3に似たものを自動的に生成する方法について、さまざまなrライブラリを調べてきましたが、まだ方法が見つかりませんでした。誰か提案はありますか？

PS-以下のグラフルーチンが少し面倒な場合は申し訳ありませんが、これは基本的に、データ分散をテストし、関連するテストを条件付きで適用し、出力テーブルと図を作成するために開発した、はるかに包括的な関数セットのフラグメントです。

最初の2つのグラフを作成するために作成したコードは次のとおりです。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

上記のコメントに関するRomanLustrikの提案は的を射ていました。私は2つの関連するライブラリに基づいてそれを行うための2つの代替方法を見つけることになりました。質問に投稿されたコードを実行した後、グループ化された治療プロットを作成するには、次のコマンドを実行します。

#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 に答える