【統計の勉強】多群比較の手順④〜検定手順をif文に起こす

【追記：2020年5月29日】コントロールの有無のところを今後の流れも含めて検討し直しています。

こんにちは。実験データを自分で解析できるようになりたい！と思って統計プログラミング言語であるRを勉強している、生命科学系研究者のえいこです。

最近は多群比較について勉強しています。

今回は、多群比較の勉強の第四回目。前回書いたフローチャートを元にif文を書いてみようと思います。

多群比較の手順の勉強をして、私が自己判断で書いたフローチャートはこちら。

これを元にif文を作っていきます。

前回の記事で、データの正規性を確かめているのでその後からコードを書いていきます。

Shapiro-wilk検定のp値によってif分岐をする

Shapiro-wilk検定の帰無仮説は「正規性に従う」なので、有意水準を0.05とした場合、p値が0.05を下回ると帰無仮説は棄却されて「正規分布ではない」という判断をします。

Shapiro-wilk検定のp値の取り出して、if文を書いていくのは2群検定でもやっているので、これを参考にして書いていきます。

3群比較なのでpの値が3つ出てきますが、同じです。

1つでもpの値が0.05よりも小さい場合は「ノンパラメトリック検定」に行くようにif文を作ります。

if (shapiro[["control"]][["p.value"]]<0.05
   |shapiro[["A"]][["p.value"]]<0.05
   |shapiro[["B"]][["p.value"]]<0.05){
 #non parametric
 print("non parametric")
 print("Kraskal wallis")
 } else {
 #parametric
   print("parametric")
   print("Levene test") #分散の検定
 }

という感じで、最初はif文の中に検定方法を直接入れるのではなく、どんな検定をしていくのかを”print()”関数を使って出力して行くようにしています。

これで、第一層目のif分岐が完成しました！

だけど、少しだけ第二層目の準備も書いてあるのです。

ノンパラメトリックの場合は、「Kruskal wallis 検定」で中央値に差があるかどうか検定していきます。

パラメトリック検定の場合は、「Levene 検定」を使って等分散性の確認をしていきます。

Levene検定の結果からさらにif分岐をする

ノンパラメトリックの「Kruskal wallis 検定」はいずれにしてもやることになるので置いておきます。

パラメトリック検定の場合、「Levene 検定」で等分散性を確認できなかった場合は「Welchのt検定をBonferroni補正する」という手順に進んでいきます。

では、Rで「Levene 検定」ってどうすれば良いのでしょう？

Rによるルビーン検定

上のサイトを参考にして、「Levene 検定」をしていきます。ここでは、2種類のパッケージを使って「Levene 検定」をしていますが、ほかのサイトなどを見ると”car”パッケージを使っている人の方が多そうなので、”car”を使うことにしました。

まずは、パッケージをインストールして”car”を呼び出します。

install.packages("car")
library(car)

で、「Levene 検定」するための関数は、“leveneTest()”（center=meanを忘れずに！）

levene<-leveneTest(data$value~data$group, center=mean)
levene

として実行すると...

leveneTestの結果はdata.frameの形で保存されていて、pの値は1行目の3列目にあることがわかります。このpの値を取り出すのは、こんなコードを書いておけば良いのです。

dpval<-levene[1,3]

“dpval”という変数の中に、leveneTestのpの値が入るようになりました。

“dpval”の値が有意水準よりも、大きいか小さいかでif分岐していけばOKです。

で、書いたコードがこちら

if(dpval<0.05){
     print("Dunnett")
     print("Tukey-Kramer")
     }
   } else {
     print("Bonferroni")
   }

最後に、比較対象群の有無で検定を分けていきたいと思います。

比較対照群の有無でif文を分岐する

ノンパラメトリック検定の場合、比較対照があればSteel、比較対照群がなければSteel-Dwassを採用しようと思っています。

パラメトリック検定の場合、比較対照があればDunnett、比較対照群がなければTukey-Kramerを採用することにします。

【2020年5月29日　追記】
コントロールの有無に関しては、こちらの記事でコントロール有無をif文分岐で検索しているので、こちらの値を使うことにします。

controlの有無を手動で宣言しなくても良いコードに変更しています！

ただ、コントロール有無で検定方法を分岐するのは必要。

まずは、比較対照群（control）があるかどうかをまず宣言しておきます。

#with control or not?#
control <- 1 #in case of w/o control, write "0"

比較対照がある場合は、1、ない場合は0と書くようにしておきます。

で、どうやってif分岐すると...

パラメトリック検定の場合は、

  if(contro==1){
   print("steel")
 } else if (control==0){
   print("Steel-Dwass")
 }

ノンパラメトリック検定の場合は、

 if(control==1){
   print("Dunnett")
   } else if (control==0) {
   print("Tukey-Kramer")
   }

このコードを該当する場所に書いておきます。

完成したコードで、作ったデータを使ってif文を回してみます。

すると...

と出力されました。

次回は、それぞれの検定に関数を調べて今回作ったif文に入れていこうと思います。

それでは、また！