ランダム化比較試験

View slides in full screen

多重比較の問題をシミュレーションで理解する

p 値の性質を理解する

ここでは多重比較の問題をシミュレーションで示す。我々が統計的仮説検定に使用する代表的な指標として p 値が挙げられるが、 p 値とは（ざっくりいうと）「帰無仮説が正しいと仮定した場合、今回得られた検定統計量と同じか、それ以上に帰無仮説から離れている確率」を意味する。 p 値が小さいということは「帰無仮説が正しいというには、今回の結果は異常すぎる。だから、帰無仮説が間違っていると考えた方がより自然だよね」ということになる。その異常さの指標が p 値であり、社会科学では通常、0.05（5%）を基準とする。この基準は「有意水準」（または、\(\alpha\)）と呼ばれる¹。\(\alpha\) = 0.05の場合、 p 値が0.05を下回ったら我々は「仮説は支持された」と表現する。しかし、 p 値の意味を考えると、 p 値が0でない限り、実際に因果効果がないにも関わらず、因果効果があるかのような結果がたまたま得られる可能性は存在する。

# 結果を再現するためにシードを固定する
set.seed(19861009)
# 平均 = 171.4、標準偏差 = 5.8の正規分布から500個の数字を無作為に抽出
group1 <- rnorm(500, mean = 171.4, sd = 5.8)
# group1と全く同じgroup2を作成
group2 <- group1

group1とgroup2は全く同じものだから、平均値も同じく約171.6である。

mean(group1)

[1] 171.5796

mean(group2)

[1] 171.5796

2群の平均値の差が0であることを帰無仮説とする検定をしてみると、 p = 1の結果が得られる。ざっくりいうと「帰無仮説が正しいとした場合、今回のような結果はめっちゃあり得る」ことをだ。

t.test(group1, group2)

    Welch Two Sample t-test

data:  group1 and group2
t = 0, df = 998, p-value = 1
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.7398454  0.7398454
sample estimates:
mean of x mean of y 
 171.5796  171.5796 

p 値だけ抽出したい場合、t.test()$p.valueを使えば良い。

t.test(group1, group2)$p.value

[1] 1

この N = 1,000を母集団とする。通常、我々は母集団を直接観察することはできず、母集団の一部である標本から母集団を推論することが多い。これからは各群で100個ずつ無作為に抽出したものを標本（ n = 200）とし、この標本を使って2群の差の検定を行う。

ノート参考）sample()関数

　第1引数のベクトルからsize個の行を無作為に抽出する関数である。sample(my_vector, size = 100)はmy_vectorから無作為に100個の要素を抽出するコードである。replace = TRUEを追加すると復元抽出を行う（既定値はFALSEであり、非復元抽出を行う）。

group1_sub_1 <- sample(group1, size = 100, replace = TRUE)
group2_sub_1 <- sample(group2, size = 100, replace = TRUE)
t.test(group1_sub_1, group2_sub_1)

    Welch Two Sample t-test

data:  group1_sub_1 and group2_sub_1
t = 0.16184, df = 190.31, p-value = 0.8716
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.613623  1.902077
sample estimates:
mean of x mean of y 
 171.4832  171.3389 

今回の p 値は約0.872である。p > 0.05であるため、母集団において2群間に差があるとは言えない。これをもう一度やってみよう。

group1_sub_2 <- sample(group1, size = 100, replace = TRUE)
group2_sub_2 <- sample(group2, size = 100, replace = TRUE)
t.test(group1_sub_2, group2_sub_2)

    Welch Two Sample t-test

data:  group1_sub_2 and group2_sub_2
t = 0.18275, df = 194.05, p-value = 0.8552
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.386072  1.669172
sample estimates:
mean of x mean of y 
 172.0572  171.9156 

今回は p \(\simeq\) 0.855であり、今回も p > 0.05であるため、母集団において2群間に差があるとは言えない。同じ結論であるが、 p 値には変動が生じた。もし、これを無数に繰り返したら、何回かは p \(\leq\) 0.05の結果が得られるかも知れない。

以上の作業を1,000回繰り返してみる。ここではfor()関数を利用したループ構文を使うが、for()関数の使い方については『私たちのR』の「Rプログラミングの基礎」を参照すること。

trials     <- 1000            # 試行回数
result_vec <- rep(NA, trials) # 結果を格納する長さ1000の空ベクトルを用意する

# iを1ずつ増やしながら「{}」内の作業を繰り返す（iの初期値は1）
# iがtrialsに達するまで繰り返す（ここだとtrials = 1000）
for (i in 1:trials) {
  # group1から無作為に100個の要素を反復抽出すし、group1_subに格納
  group1_sub <- sample(group1, size = 100, replace = TRUE)
  # group2から無作為に100個の要素を反復抽出すし、group2_subに格納
  group2_sub <- sample(group2, size = 100, replace = TRUE)
  
  # 結果ベクトルのi番目に平均値の差の検定のp値を格納
  result_vec[i] <- t.test(group1_sub, group2_sub)$p.value
}

結果が格納されているresult_vecの中身を確認してみよう。計1,000個の数字が入っているので、すべて出力せず、head()関数を使って最初の50ケースだけを出力してみよう²。

head(result_vec, n = 50)

 [1] 0.190327968 0.222601686 0.166202227 0.706861635 0.507185738 0.904233707 0.308284022 0.587728946
 [9] 0.584220793 0.771474434 0.583631579 0.351888684 0.821629231 0.250568478 0.474959883 0.268788816
[17] 0.543049151 0.748209635 0.575553911 0.838333618 0.610160479 0.460881573 0.086647600 0.380820889
[25] 0.461473898 0.589726513 0.114475163 0.743054055 0.565888352 0.543396088 0.757256967 0.007551659
[33] 0.363392991 0.519455208 0.308931568 0.967461390 0.253462409 0.863548096 0.451520105 0.889560868
[41] 0.475251439 0.567767986 0.549038458 0.861901764 0.517104060 0.124166328 0.973664109 0.782367208
[49] 0.121350999 0.207622189

これらの数字が0.05以下がどうか一つ一つ確認するのは面倒だから、関係演算子を使ってみよう。result_vec <= 0.05を入力すると、result_vecの要素が0.05以下かどうかを一つずつ判定し、TRUEかFALSEを返す。ここでは紙幅の関係上（?）、最初の50個のみ出力するが、result_vec <= 0.05のみでも良い。

head(result_vec <= 0.05, n = 50)

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[17] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
[33] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[49] FALSE FALSE

続いて、table()関数を使用し、TRUEとFALSEの数を度数分布表で出力してみよう。

ノート参考）table()関数

　第1引数のベクトルの度数分布表を返す関数である。連続変数（間隔尺度、または比率尺度で測定された変数）に対してはほぼ使わないが、離散変数（いわゆる名目変数、順序変数）の分布を知るには便利な関数だ。

table(result_vec <= 0.05)

FALSE  TRUE 
  957    43 

計1,000回の試行中、43回の試行において p \(\leq\) 0.05の結果が得られた。%で換算すると4.3%である。母集団においてはgroup1とgroup2の間に差がないにも関わらず、標本を用いた推論では約4.3%の割合で「group1とgroup2の間に差がある」と判定されたのである。これが1,000回でなく、100万回、1億回なら5%に近づく。つまり、\(\alpha\) = 0.05を採用した場合、帰無仮説が正しくても、約5%の割合で帰無仮説が棄却され、対立仮説（\(\simeq\) 自分が推したい仮説）が支持される結果になる。

多重比較（多重検定）とは何か

p 値は直感的な指標ではないが、その意味を理解した上で使ったら大変便利な指標ではある。しかし、以上の説明は2群間の比較（ここではgroup1とgroup2）で有効な話である。たとえば、「啓発メッセージは投票率を上げる」という仮説を検証したいとしよう。そこで啓発メッセージ4種を用意し、5つのグループに無作為に割り当てられた被験者グループに異なる啓発メッセージを呈示したとする（啓発メッセージを呈示していないグループも含めて5群である）。この場合、統制群が1つ、処置群が4つとなる。

啓発メッセージの効果を測定するためには何と何を比較すべきか。それは「処置群の投票率 - 統制群の投票率」である。しかし、問題は、今回のケースだとその処置群が4つもあるということだ。つまり、我々は「処置群1の投票率 - 統制群の投票率」、「処置群2の投票率 - 統制群の投票率」、「処置群3の投票率 - 統制群の投票率」、「処置群4の投票率 - 統制群の投票率」を計算する。つまり、比較の回数は4回である。そして、一つでも差分が統計的有意に0でないという結果が得られたら「啓発メッセージは投票率を上げる!」と結論づけることになる。3つのp値が0.05より大きくても1つでも p \(\leq\) 0.05で仮説を支持してしまう。

これはちょっと卑怯ではないだろうか。実際そうだ。そして、これには「多重比較の問題」（multiple comparisons problme）という名前が付いている。つづいて、多重比較の問題がどれほど深刻な問題かをシミュレーションで確認してみよう。

多重比較の問題を理解する

まず、N = 1,000の架空の母集団（pop_df）を作成する。1,000名をランダムに5つのグループ（1つの統制群を含む）に分け、異なる刺激を与えた後、応答変数（y）を測定したデータである。母集団における処置効果（処置群1〜4のy - 統制群のy）はほぼ0のデータである。

ノート参考）tibble()関数

　表形式オブジェクト（\(\simeq\) data.frame型）を作成する関数である。R内臓のdata.frame()関数とほぼ同じ関数である。この関数を使用するためには事前に{tibble}パッケージを読み込んでおく必要があるが、{tidyverse}を読み込むと自動的に読み込まれる。

pacman::p_load(tidyverse)
set.seed(19861008)
pop_df <- tibble(treat = rep(c("Control", "Treatment1", "Treatment2",
                               "Treatment3", "Treatment4"), each = 200),
                 y     = rnorm(1000))

head(pop_df)

# A tibble: 6 × 2
  treat         y
  <chr>     <dbl>
1 Control -0.0863
2 Control  0.396 
3 Control -1.33  
4 Control  0.574 
5 Control  0.152 
6 Control  0.555 

線形回帰モデル（lm()）を使用し、母集団における処置効果（真の因果効果; パラメーター）を確認してみよう。

pop_fit <- lm(y ~ treat, data = pop_df)
summary(pop_fit)

Call:
lm(formula = y ~ treat, data = pop_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2781 -0.6693 -0.0179  0.6565  4.3429 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)     -0.06227    0.07209  -0.864    0.388
treatTreatment1  0.03238    0.10195   0.318    0.751
treatTreatment2  0.02159    0.10195   0.212    0.832
treatTreatment3  0.09748    0.10195   0.956    0.339
treatTreatment4  0.04322    0.10195   0.424    0.672

Residual standard error: 1.019 on 995 degrees of freedom
Multiple R-squared:  0.001025,  Adjusted R-squared:  -0.002991 
F-statistic: 0.2553 on 4 and 995 DF,  p-value: 0.9065

処置群1の処置効果の点推定値は-0.021、 p 値は1であり、統計的有意な処置効果は確認できない。処置群2、3、4も同様であり、処置効果は確認できない。

しかし、通常、母集団のサイズは非常に大きく、直截観察することはできない。だから我々は母集団の一部である標本を抽出し、標本から母集団を推論する。ここでは母集団からn = 500の標本を抽出し、その標本を使用して処置効果を推定してみよう。

ノート参考）slice_sample()関数

　表形式オブジェクトからn個の行を無作為に抽出する関数である。slice_sample(pop_df, n = 100)はpop_dfから無作為に100行を抽出するコードである。replace = TRUEを追加すると復元抽出を行う（既定値はFALSEであり、非復元抽出を行う）。この関数を使用するためには事前に{dplyr}パッケージを読み込んでおく必要があるが、{tidyverse}を読み込むと自動的に読み込まれる。

sample_fit_1 <- lm(y ~ treat, data = slice_sample(pop_df, n = 500, replace = TRUE))
summary(sample_fit_1)

Call:
lm(formula = y ~ treat, data = slice_sample(pop_df, n = 500, 
    replace = TRUE))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9410 -0.7202 -0.0139  0.6824  3.5299 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)     -0.04174    0.10691  -0.390    0.696
treatTreatment1 -0.13840    0.14972  -0.924    0.356
treatTreatment2 -0.14629    0.15729  -0.930    0.353
treatTreatment3  0.04986    0.14801   0.337    0.736
treatTreatment4 -0.02827    0.14902  -0.190    0.850

Residual standard error: 1.064 on 495 degrees of freedom
Multiple R-squared:  0.005267,  Adjusted R-squared:  -0.002771 
F-statistic: 0.6553 on 4 and 495 DF,  p-value: 0.6234

真の処置効果とズレ（誤差; エラー）はあるが、今回も全ての処置効果の推定値の p 値が0.05以上であり、統計的に有意な処置効果は確認できない。推定値の p 値のみを抽出するためにはsummary(lmオブジェクト名)$coefficients[-1, 4]と入力する。

summary(sample_fit_1)$coefficients[-1, 4]

treatTreatment1 treatTreatment2 treatTreatment3 treatTreatment4 
      0.3557134       0.3528137       0.7363812       0.8496285 

これらの p 値が0.05未満かを確認する。

summary(sample_fit_1)$coefficients[-1, 4] <= 0.05

treatTreatment1 treatTreatment2 treatTreatment3 treatTreatment4 
          FALSE           FALSE           FALSE           FALSE 

一つでも統計的に処置効果が認められたら、つまり、一つでもTRUEが確認されたら「仮説は支持された」と解釈する場合、仮説の支持有無はany()関数で判定できる。

ノート参考）any()関数

　TRUEとFALSEのみで構成されたベクトルに対し、一つでもTRUEがあったらTRUEを返す関数である。ベクトルのすべての要素がFALSEの場合のみ、FALSEを返す。

any(summary(sample_fit_1)$coef[-1, 4] <= 0.05)

[1] FALSE

今回はすべての p 値が0.05以上だという結果が得られた。つまり、処置群1〜4すべて統計的有意な処置効果は確認できず、仮説は支持されなかった（棄却されたわけではない）。

この作業を1,000回繰り返してみよう。ここではTreatment 1から4までの各処置効果の p 値も別途格納しておこう。

trials     <- 1000 # 試行回数
result_vec <- rep(NA, 1000) # 結果を格納する空ベクトル
t1_vec <- rep(NA, 1000) # 処置群1の処置効果のp値を格納する空ベクトル
t2_vec <- rep(NA, 1000) # 処置群2の処置効果のp値を格納する空ベクトル
t3_vec <- rep(NA, 1000) # 処置群3の処置効果のp値を格納する空ベクトル
t4_vec <- rep(NA, 1000) # 処置群4の処置効果のp値を格納する空ベクトル

for (i in 1:trials) {
  temp_data <- slice_sample(pop_df, n = 500, replace = TRUE)
  temp_fit  <- lm(y ~ treat, data = temp_data)
  t1_vec[i] <- summary(temp_fit)$coefficients[2, 4] # 処置群1の処置効果のp値
  t2_vec[i] <- summary(temp_fit)$coefficients[3, 4] # 処置群2の処置効果のp値
  t3_vec[i] <- summary(temp_fit)$coefficients[4, 4] # 処置群3の処置効果のp値
  t4_vec[i] <- summary(temp_fit)$coefficients[5, 4] # 処置群4の処置効果のp値
  result_vec[i] <- any(summary(temp_fit)$coefficients[-1, 4] <= 0.05)
}

まず、それぞれの処置に統計的有意な処置効果が確認できた割合を確認してみよう。

table(t1_vec <= 0.05)

FALSE  TRUE 
  954    46 

計1000回の施行中、46回の試行に統計的有意な処置効果が確認された。処置群1の真の処置効果は0であるが、約5%の割合で統計的有意な処置効果が得られてしまう。これは処置群2、3、4でも確認できる（実際に確認してみよう）。

しかし、我々は一つでも統計的有意な処置効果が得られたら「仮説は支持された」と判定する³。その結果はresult_vecに格納されている。

table(result_vec)

result_vec
FALSE  TRUE 
  798   202 

計1000回の試行中、202回の試行において一つ以上の統計的有意な処置効果が確認され、仮説が支持される結果が得られた。真の処置効果は0であるが、このように処置効果が統計的有意に0でない状況が約20%の割合で出現することを意味する。この20%ということは我々が統計的有意性検定に使用した基準（\(\alpha\)）に比較回数の積である（後述するが、厳密には積ではない）。今回は\(\alpha\) = 0.05、比較回数は4回⁴だから、0.05 \(\times\) 4 = 0.2である。厳密には20%ではないが（後述）、重要なのは比較の回数と間違った判定を下ろす可能性は比例することだ。

重要参考）厳密には20%じゃない!

　今回の例でいると、仮説が支持される結果が得られる割合の期待値は20%でなく、厳密には約18.5%である。計算式は1 - (1 - \(\alpha\))^比較回数であり、今回は1 - (1 - 0.05)⁴ \(\simeq\) 0.185になる。比較回数が2回なら約9.75%、3回なら約14.3%、5回なら22.6%であり、本文に登場した積を使った計算はあくまでも簡易的な計算方法である。

ボンフェローニ補正をしてみよう

4回の比較を繰り返す場合、1回のみの比較に比べ、「仮説は支持された（帰無仮説は棄却された）」と判定される確率が約4倍高くなることが分かった。これを補正するためにはどうすれば良いだろうか。様々な手法が提唱されているが、ここでは最も簡単な方法であるボンフェローニ補正（Bonferroni correction）を施してみよう。ボンフェローニ補正は2つの考え方がある。比較回数を m とした場合、

• \(\alpha\)を\(\frac{1}{\mathsfit{m}}\)にし、 p \(\leq \alpha\)の場合、「5%水準で統計的有意である」と判定する。
• p 値を m 倍し、 p \(\leq \alpha\)の場合、「5%水準で統計的有意である」と判定する。

2つは同じことであるが、ここでは両方やってみよう。まず、\(\alpha\)を\(\frac{1}{\mathsfit{m}}\)にする方法である。今回の例だと、\(\alpha\) = 0.0125を採用し、 p 値が0.0125以下であれば「5%水準で統計的に有意である」と解釈する。処置群1〜4の処置効果の推定値の p 値はt1_vec、t2_vec、t3_vec、t4_vecに格納されている。一つでも p 値が0.0125以下であれば仮説は支持されたと解釈されるため、ペアごとに最も小さい p 値のみを抽出してみよう。

重要参考）pmin()関数

　pmin()は2つ以上の同じ長さのベクトルを受け付け、ペアごとに最小値を返す関数である。例えば、xの中身が1、3、5で、yの中身が2、2、4の場合、pmax(x, y)のようにコードを書く。まず、xの1番目の要素（=1）とyの1番目の要素（=2）を比較し、最小値である1を返す。続いて、xの2番目の要素（=3）とyの2番目の要素（=2）を比較し、最小値である1を返す。最後に、xの3番目の要素（=5）とyの3番目の要素（=4）を比較し、最小値である4を返す。

minimum_p <- pmin(t1_vec, t2_vec, t3_vec, t4_vec)

head(minimum_p, n = 10)

 [1] 0.195778366 0.578161172 0.234322363 0.089455091 0.284272075 0.058462114 0.363639984 0.544207581
 [9] 0.001115049 0.037441882

各試行における最も小さい p 値が格納されていることが分かる。続いて、table()関数で p 値が0.0125以下の試行数を出力してみる。

table(minimum_p <= 0.0125)

FALSE  TRUE 
  923    77 

1,000回の試行中、202回の試行において、1つ以上の処置効果が統計的有意と判定され、仮説は支持されたと解釈される。

続いて、 p 値を m 倍にするやりかたでやってみよう。先ほどのコードと同じだが、t1_vec、…、t4_vecを4倍にしておく必要がある。

t1_vec_bc <- t1_vec * 4
t2_vec_bc <- t2_vec * 4
t3_vec_bc <- t3_vec * 4
t4_vec_bc <- t4_vec * 4

minimum_p <- pmin(t1_vec_bc, t2_vec_bc, t3_vec_bc, t4_vec_bc)
head(minimum_p, n = 10)

 [1] 0.783113462 2.312644689 0.937289452 0.357820366 1.137088299 0.233848454 1.454559937 2.176830326
 [9] 0.004460195 0.149767528

p 値を m 倍したのであれば、判定の基準は\(\alpha\)のままで良い。

table(minimum_p <= 0.05)

FALSE  TRUE 
  923    77 

同じ結果が得られた。今回は約7.7%になったが、試行回数を増やせば増やすほど5%へ近づくと期待される。

多重比較の問題を補正する方法はボンフェローニ法以外にもホルム法（Holm）、BH法（Benjamini-Hochberg）、ダネット検定（Dunnett’s test）などがあり、実は今回のケースだと、複数に比較にわたって一つの統制群が共有されているため、ダネット検定が最も適切と考えられる。ただし、ボンフェローニ補正は他の補正法よりも保守的、つまり統計的有意な結果が得られにくく、分析する側からみれば最も都合の悪い手法だから、ボンフェローニ補正を使えば文句は言われないだろう。

重要参考）厳密には5%じゃない!

　今回の例でいると、ボンフェローニ補正後、仮説が支持される結果が得られる割合の期待値は5%でなく、厳密には約4.9%である。計算式は1 - (1 - \(\frac{\alpha}{\textsf{比較回数}}\))^比較回数であり、今回は1 - (1 - 0.0125)⁴ \(\simeq\) 0.049、つまり約4.9%になる。

ノート参考）ダネット検定のやりかた

　ダネット検定をRで実装するには{DescTools}か{multcomp}パッケージを使用する。{DescTools}のDunnettTest()関数はlm()関数と使い方が似ているため、ここでは{DescTools}を使用した例を紹介する。lm()関数とは違って、control引数を使って統制群のラベルを指定することができる。

pacman::p_load(DescTools)
sample_500 <- slice_sample(pop_df, n = 500, replace = TRUE)
summary(lm(y ~ treat, data = sample_500))

Call:
lm(formula = y ~ treat, data = sample_500)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9273 -0.7047 -0.0179  0.6375  4.4385 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
(Intercept)     -0.26559    0.09658  -2.750  0.00618 **
treatTreatment1  0.14014    0.14620   0.959  0.33826   
treatTreatment2  0.12720    0.14376   0.885  0.37667   
treatTreatment3  0.30967    0.14941   2.073  0.03872 * 
treatTreatment4  0.18186    0.14941   1.217  0.22410   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.075 on 495 degrees of freedom
Multiple R-squared:  0.008941,  Adjusted R-squared:  0.0009327 
F-statistic: 1.116 on 4 and 495 DF,  p-value: 0.3479

DunnettTest(y ~ treat, data = sample_500, control = "Control")

  Dunnett's test for comparing several treatments with a control :  
    95% family-wise confidence level

$Control
                        diff      lwr.ci    upr.ci   pval    
Treatment1-Control 0.1401398 -0.22002134 0.5003009 0.7580    
Treatment2-Control 0.1272025 -0.22693768 0.4813427 0.8050    
Treatment3-Control 0.3096694 -0.05838828 0.6777271 0.1289    
Treatment4-Control 0.1818623 -0.18619538 0.5499200 0.5763    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

　通常の線形回帰モデルから得られたp値よりも大きいp値が得られたことが分かる。

脚注

1. 有意確率とは別の概念である。そもそも有意確率という表現は避けるべきである。

2. n = 30を省略すると最初の6個まで出力される

3. どういう時に「仮説は支持された」と判定するかをあらかじめ明確にしておく必要がある。

4. Control vs. Treatment 1 / Control vs. Treatment 2 / Control vs. Treatment 3 / Control vs. Treatment 4