- R ver. 3.5.2
- RStudio ver. 1.1.463
- multcomp ver. 1.4.8
です.
そもそも多重比較とは何やねん, という方は, この記事をお読みください.
卒論などで, 分散分析の多重比較をしたとき, 有意差を表すアルファベットをつけることがあります. どういうことかというと, たとえば, 以下のサイトを参考にしてください.
- 多重比較のアルファベットの付け方 https://sites.google.com/site/hymd3a/statistics/numbering
- 多重比較の結果を表のみで記す、例のabc…の割り振り方 http://blog.kz-md.net/?p=566
こういった理屈は知っておかなければなりませんが, これを手計算で行うのは大変です.
R を使えば簡単に, 多重比較の有意差を表すアルファベットを得ることができます.
必要なのはパッケージ「multcomp」です.
インストールしていない場合は, インストールしておきましょう.
install.packages("multcomp")
library() で,「multcomp」を起動します.
library(multcomp)
たとえば, iris の萼片長 (Sepal.Length) の平均が, 品種ごとに差があるかどうかを確かめてみます. iris のデータ読み込んで,
data(iris)
品種と萼片長のデータのみを取り出します (別に取り出さなくてもできますが).
df_S.L <- subset(iris, select=c(Species,Sepal.Length))
df.S.L で, Tukey の多重比較の検定をしてみましょう(iris の萼片長が, Tukey の多重比較の検定に適しているかどうかは, 今回は不問にします). R のデフォルトのコマンドでもできます.
TukeyHSD(aov(Sepal.Length ~ Species , data = df_S.L))
結果は以下のように表示されます.
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Sepal.Length ~ Species, data = df_S.L)
$Species
diff lwr upr p adj
versicolor-setosa 0.930 0.6862273 1.1737727 0
virginica-setosa 1.582 1.3382273 1.8257727 0
virginica-versicolor 0.652 0.4082273 0.8957727 0
これで多重比較の結果を得ることができたわけですが, この結果から有意差を表すアルファベットをつけるのはけっこう大変です. ということで, multcomp パッケージの関数を使います.
ONWY <- aov(Sepal.Length ~ Species , data = df_S.L) MLTC <- glht(ONWY, linfct = mcp(Species = "Tukey")) cld(MLTC, level = 0.05, decreasing = TRUE)
この 1 行目の「data = 」に, iris の品種と萼片長のデータである「df_S.L」を入力して, あとは, 2 行目 3 行目はそのまま打ち込めば, 多重比較で有意差を表すアルファベットを得ることができます.
setosa versicolor virginica
"c" "b" "a"
しかし, これをそのまま csv ファイルなどへ書き出すことはできません. たとえば,
cld.S.L <- cld(MLTC, level = 0.05, decreasing = TRUE)
などとして,
write.csv(cld.S.L, file=“irisS.L_alphabet.csv")
などとしても, csv へ書き出すことはできません. これだとあまり便利ではないので, アルファベットだけを取り出したいと思います. それが
[["mcletters"]][["Letters”]]
です.
さきほど,「cld.S.L」に「cld(MLTC, level = 0.05, decreasing = TRUE)」を格納したので,
cld.S.L[["mcletters"]][["Letters”]]
とすれば,
setosa versicolor virginica
"c" "b" "a"
という結果を得ることができます.
これで得られた結果は, マトリックスへ変換できます.
matrix(cldS.L[["mcletters"]][["Letters"]])
これを例えば「cldS.L_mx」に格納して,
cldS.L_mx <- matrix(cldS.L[["mcletters"]][["Letters"]])
行と名前をつければ (「s.d」は「有意差」」の意味です),
rownames(cldS.L_mx) <- c("setlsa", "versicolor","virginica") colnames(cldS.L_mx) <- c("s.d”)
こういうのができます.
cldS.L_mx
s.d
setlsa "c"
versicolor "b"
virginica "a"
この状態だとデータフレームへ変換できるので,
as.data.frame(cldS.L_mx)
たとえば「cldS.L_df」に格納する, なんてすれば,
cldS.L_df <- as.data.frame(cldS.L_mx)
cldS.L_df を, csv へ書き出すことができます. これを応用すれば, 有意差のアルファベットが既に記入されている表を自動生成できたりできそうですね.