Zenn
🌱

【R】multcompで多重比較結果をアルファベット表示させる

2022/09/09に公開
R
統計学
多重比較
tech

はじめに

目的

以前はExcelを使ったT検定を行なっていましたが、多群間サンプルでT検定を繰り返し使うことがよろしくないと叱責いただいたので、統計学の勉強をし直すことになりました。
T検定を繰り返してはいけない!←これ、めちゃくちゃ分かりやすかったです。

…

さて、原理を理解したとはいえ、この計算はやりたくないのでRを使います。
今回は、multcompというpackageの、glht関数を使うと簡単に多群間サンプルの分散分析(よく論文のfigureで見るやつ）ができるらしいので、やってみました。
とりあえず、データの読み込みから、検定結果を表示するところまでをできるようになってみようと思います。

実験条件

今回のfigでは、植物が特定のストレスを受けた際に、"強い" or "弱い"を、定量的に、統計検定をもって示したいです。(エコタイプA〜Dの4群を用意)

本実験では、ストレス処理後に、植物の持つ緑色色素(クロロフィル)が特異的に吸光する波長を検出することで、その色素含量を定量しました。

スクリーンショット 2022-09-07 16.18.41.png

要は、植物が緑色であるほど生存していることを示していて、その"ストレス耐性"に、統計学的に有意な差があるのか、ということです。

以下、今回の検定で使っていくデータになるので、ご活用ください。

test_abc.csv
ecotype,chlorophyll
A,0.257875877
A,0.2109656862
A,0.2064517446
B,0.8902901736
B,0.6678915832
B,1.146381712
C,0.2005284342
C,0.2509381304
C,0.0866050339
D,0.7237973545
D,0.5866545809
D,0.6178009144

要求package

  • multcomp
  • car

動作環境

  • MacBook Air(M1, 2020), macOS Monterey
  • RStudio, "Ghost Orchid" Release for macOS Mozilla/5.0 (Macintosh; Intel Mac OS X 12_0_1)

1. ワーキングディレクトリの設定とパッケージのロード

僕は忘れがちなので、この2つは必ず初めにやります。

ワーキングディレクトリの設定.R
setwd("~/Desktop/R/figure")
#検定したいデータのcsvファイルが入っているディレクトリをワーキングディレクトリに設定

自分の場合は、Desktopにおいてある、Rフォルダの下のfigureフォルダに、test_abc.csvを置いていたので、そこにsetwdで仕事場を指定。続いて、パッケージのロードのためにlibrary()を実行します。
　未インストールの方はinstall.package()であらかじめインストールしておきましょう。（インストールは初回のみ）

パッケージのロード.R
library(multcomp) #多重比較につかう
library(car) #多重比較前の等分散性の検定につかう

2. csvファイルの読み込みと加工

次に実験データの入っているcsvファイルをRStudioに読み込ませます。Excelから読み込む方法もあるらしいけど、僕は未修です。

csvファイルの読み込み.R
d <- read.csv("test_abc.csv", header = TRUE)

上の操作では、read.csv()を用いて読み込んだcsvファイルを、dという文字に代入する操作をしています。（dはなんの文字を指定しても良いです。僕の場合は、短い方が楽なのでdataの頭文字を指定しています。）

header = TRUEはデータの１行目がサンプル名とか、数値データじゃない時に指定してあげます。１行目から数値データの時にはheader = FALSEを指定します。

dを入力、実行してあげると、Consoleに表が出力されます。(実行結果↓)

　続いて、ちょっとだけデータの加工をします。とはいえ、目に見える加工はしません。この後に使う関数がecotypeを認識できるデータ型に変更するために、データの持つ属性を編集します。

データ型の加工.R
d$ecotype <- factor(d$ecotype)
attach(d)

3. ついに検定していきます！

まず統計学によって多重比較をする際の流れについてですが、

このような流れで行います。
（図表は以下から引用させていただきました。）
統計検定を理解せずに使っている人のために III - 池田郁男

初めに、パラメトリック検定/ノンパラメトリック検定のどちらをすべきなのか、多群間での分散について検討できるLevene's test（F検定の３群以上バージョンみたいな検定試験）をやってみます。

等分散性の検定: Levene's test

まず、得られたサンプルのデータの母分散が等分散するか、不当分散するかについて検討します。

僕たち研究者は、目の前にあるサンプルの分散や平均を調べたいのではなく、集団全体の分散や平均について考えたいのだ（by えらい人）

等分散性の検定.R
leveneTest(chlorophyll, ecotype, center = mean)

実行結果です↓

levene検定の帰無仮説H0は「各群の分散は等しい」です。検定結果Prのcolumnを見てみると、
p値 = 0.1551 > 0.05(有意水準5%)
となり、帰無仮説H0は棄却されません。
つまり、ecotype A~Dは等分散性を仮定して、一元配置分散分析(ANOVA)で平均の差があるのか、検定をしていきます。

one way ANOVA（一元配置分散分析）

次に、one way ANOVA(analysis of variance)で、多群のデータの平均に差があるのかを検定します。

一元配置分散分析.R
anova(lm(chlorophyll ~ ecotype))

実行結果です↓

今回の検定における仮説は、

  • 帰無仮説H0は「4データ群の平均は等しい」
  • 対立仮説H1は「少なくとも1つのペアの平均に差がある」

です。検定結果Prのcolumnを見てみると、

p値 = 0.0004176 < 0.05(有意水準5%)

となり、帰無仮説H0は棄却され、対立仮説Hnを受容します。よって、平均の少なくとも1つのペアに有意差があることが分かります。しかし、それがどのデータ群間のものなのかは分かりません。

4. 多重比較

前項のANOVAで、データ群のどこかのペアの平均には差があることが示されたので、多重比較をしていきます。多重比較にもいろいろ方法はあるみたいですが、今回は全ての群間の平均値について比較をするTukeyの方法で検定します。

Tukey’s T-test

早速、コードを入力していきます。

TukeyのT検定.R
TukeyHSD(aov(chlorophyll ~ ecotype))

実行しましょう↓

検定結果のp値が、A~Dの各データ群総当たりで算出されています。

TukeyのT検定の結果をアルファベットとして出力する

これで終わりではありません。よく論文のfigureで目にする、棒グラフの上についているアルファベットで検定結果を出力できるようにしていきます。

アルファベット.R
res <- aov(chlorophyll ~ ecotype, data = d) #ANOVAをresという函に格納
tuk <- glht(res, linfct = mcp(ecotype = "Tukey")) #多重比較の結果をtukという函へ格納
cld(tuk, decreasing = F) #アルファベットとして検定結果を出力

実行しましょう↓

無事、出力できました〜

結果の解釈

（わかると思うけど念の為…）
このグラフをみると、エコタイプAとC、BとDがそれぞれ同程度のストレス耐性を保持していて、A,C・B,D間のストレス耐性の度合いには有意に差があることがわかりました〜
（もっとabcdにふれるデータセットを用意すればよかったと少しだけ後悔）

おわりに

今回は、Rで多重比較の結果をアルファベット出力しました。次回は棒グラフと同時に出力するscriptについて記事を書こうかな。(めちゃくちゃ書くの大変だった)

最近僕らのラボ内で 「正しい検定を」 の波が来ています。この記事が誰かの助けになれたら嬉しいです。

References

今回のR script

test_abc.R
setwd("~/Desktop/R/figure")
library(multcomp)
library(psych)

d <-read.csv("test_abc.csv", header = TRUE)
d$ecotype <- factor(d$ecotype)
attach(d)

leveneTest(chlorophyll, ecotype, center = mean) 
anova(lm(chlorophyll ~ ecotype))
TukeyHSD(aov(chlorophyll ~ ecotype))

res <- aov(chlorophyll ~ ecotype, data = d)
tuk <- glht(res, linfct = mcp(ecotype = "Tukey"))
cld(tuk, decreasing = F)

Discussion