多重比較をRで実施する方法

はじめに

Rは細々続けています。

2群以上実験や複数水準のデータを取得したとき、多重比較をしたくなることがあると思います。

総当たりの多重比較は実験デザインとしてあまりよくないとはわかっていつつも、そういう場面があることもしばしばです。ということで、今回は多重比較をRでやる方法についてです。 参考は以下のページです。(ANOVAの部分をすっ飛ばしていますが)

ANOVA in R: The Ultimate Guide - Datanovia

最初に使用するパッケージを読み込んでおきます。

library(tidyverse)
library(ggpubr)
library(rstatix)
library(gt)
library(datarium)
library(emmeans)

前提

以下の解析はtidyverseのパッケージを用いています。 実際の解析で自分のデータを使う際にはtidy dataに変換してから実行が必要です。 tidy dataについてはたくさん解説があるので、そちらを参照 (例えばこちら) 私たちのR - 17  整然データ構造

一元配置デザインの場合

要するに3群以上の平均値について比較したい場合です。教科書的には (?) よく、ANOVAしてpost hoc みたいな言い方をしますが、別にANOVAは必須ではないと思われます。

使用するサンプルデータはplant growthで、3条件 (control, treatment1, treatment2) での植物の成長を見ています。 組み込まれているので、

data("PlantGrowth")

で読み込めます。

まずはグラフをチェックしてみます。(ggpubrだとこういうとき楽)

ggboxplot(PlantGrowth, x = "group", y = "weight", color = "group")

箱ヒゲ図

なんとなく群間に差がありそうな気がします。

Tukey-HSD

よく使うTukey-HSDはいかで実行可能です。

PlantGrowth %>% tukey_hsd(weight ~ group)

ちなみにgtパッケージを使うと以下のような整った表形式で表示可能です。

PlantGrowth %>% tukey_hsd(weight ~ group) %>% gt()

tukey

Tukey-HSDは論文でもよくみますが、独立、正規性、等分散性の3条件が必要なので、適応可能な範囲は比較的狭いと思われます。

Bonferroniとその変法たち

Bonferroniもよく見かける手法です。これも解説はいくらでもネットにあるので飛ばしますが、検定の回数で有意水準を割るというのが簡単な解説です。(逆にp-valueに検定の回数をかけて表示することも) そこから派生したHolmやShafferの方法は割る数の方を少しずつ減らしたりしている方法になります。 計算方法からわかるように、Bonferroniやその変法はかなり適応範囲が広いです。

pairwise_t_test(data=PlantGrowth,
                weight~group, pool.sd = TRUE,
                p.adjust.method = "holm") %>% gt()

ここではHolmをやってみて、gt()で表形式で出してみました。

holm

pairwise_t_testの引数は最初はデータ、2つ目はfomula形式になっています。 fomula形式は 独立変数 ~ 従属変数 という形式での記載になっております。この記載が変数間の関係性が明示的で私は好みです。 3つ目はpool.sdで、(間違っているかもしれませんが) 検定統計量を計算する際の分散としてプールされた分散を使用しているという意味だと思っております。 最後のp.adjust.methodで調整方法を選べて、ここではholmにしていますが、他にも"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"が選択可能です。(あまり使わないので、bonferroniとholm以外わかりませんが)

ちなみに違いを見るためにBonferroniも出してみます。

bonf

最もp-valueが小さい比較のp-valueは同じですが、他で違いがあります。holmはp-valueを小さい順に並べて、1つずつp-valueにかける値を減らします。そのため最初に評価される、最も小さいp-valueは同一になります、

二元配置デザインの場合

二元配置の場合、多重比較をやりたくなる場面としては交互作用があって、全群比較をしたくなった場合だと思います。 ここでは上記と異なるデータで、jobsatisfactionというデータを用います。

データの見た目

データはid列とgender列、education列とscore列からなります。 まずは簡単に図示してみます。

ggboxplot(
  jobsatisfaction, x = "gender", y = "score",
  color = "education_level"
)

jobsatisfaction

今回はANOVAの結果も載せておきます。

jobsatisfaction %>% anova_test(score ~ gender * education_level) %>% gt()

2 way ANOVA

interactionが有意になっています。 この場合、全群での比較をしたいわけですが、1元配置のときのpairwise_t_testはそのままだと使えません。 一度データをグループに分けてから使用する必要があります。

jobsatisfaction %>% 
  tukey_hsd(score~education_level + gender,
                  pool.sd = TRUE,
                  p.adjust.method = "holm")

multiple comp

三元以上の場合

おそらく同様にできると思います。 ただ三元配置以上は交互作用の解釈がかなり難しくなるので、可能であれば実験デザインの段階で避けるべきでしょうね(少なくとも私がやっている生物系の実験の場合には)

感想

Rでの多重比較について書いてみました。 以前はSPSSを使っていたのですが、SPSSではGUIだったのでRとはだいぶん感覚が違います。Rの方が手法を勉強する気にはなるので良いと思っております。