一元配置分散分析を多要因に拡張した解析

　一元配置分散分析の項では1つのカテゴリカル変数を説明変数として組み込む解析として紹介した。本項目では、2つのカテゴリカル変数を説明変数とする解析として、二元配置分散分析Two-way ANOVAを紹介する。別にカテゴリカル変数がさらに増えてもよく、3つなら三元、4つなら四元……となる（ただし、目にすることは多くないと思う）。今後は、要因数に関わらず分散分析、と呼称する。さて、一元配置分散分析では要因内の水準間における平均値の比較として紹介した。もっと言えば、水準間の平均値が等しい、という帰無仮説を検討するための方法であった。統計モデリングの項で紹介した二元配置分散分析も同様に帰無仮説を有している（線形回帰も多重線形回帰も同様）。二元配置分散分析の統計モデルをおさらいすると、yを被説明変数、g1およびg2を二値のカテゴリカル変数（つまり、0 or 1の値をとる）として、交互作用を考慮する場合、

y ~ b0 + b1 × g1 + b2 × g2 + b3 × (g1×g2)

として、b0～b4を推定するものである。復習すると、ここの右辺を線形予測子と呼ぶのであった。ここで、帰無仮説は「b0 = 0、b1 = 0、……b4 = 0」となることである。これはいままでの検定の考え方に帰着することができる。b1 = 0とはg1に有意な効果がないことを意味し、g1に含まれる2水準がAとBとするなら「Aの平均 - Bの平均 = 0」すなわち「Aの平均 = Bの平均」となる。これは検定における平均値が等しい帰無仮説と同様である。

　g2を3水準含むカテゴリカル変数と考えたらどうだろうか？　g2の中の水準をX、Y、Zとして、水準Xを基準に、Yの効果の大きさ = 係数b2_y、Zの効果の大きさ = 係数b2_z、g1とYの交互作用 = 係数b3_y、g1とZの交互作用 = 係数b3_zとすれば、ダミー変数を考え方を用いて、

Xの場合：y ~ b0 + b1 × g1

Yの場合：y ~ b0 + b1 × g1 + b2_y + b3_y × g1

Zの場合：y ~ b0 + b1 × g1 + b2_z + b3_z × g1

として、b0、b1、b2_y、b2_z、b3_y、b3_zを推定する。ここで、g2側に着目したときに「b2_y = b2_z = 0」が帰無仮説となる。言い換えれば、「Xの平均 - Yの平均 = Xの平均 - Zの平均 = 0」すなわち「Xの平均 = Yの平均 = Zの平均」となり、3水準の一元配置分散分析の帰無仮説に帰着できる。

　2標本の平均値の検定や3水準の分散分析のときと、上記の状況で異なっている点は、交互作用の係数を推定している点である。上記の内容にならって、帰無仮説は同様に、交互作用項に有意な効果はない、すなわち「b3 = 0」であったり、「b3_y = b3_z = 0」である。

　そして、帰無仮説の大前提として、各水準の標本はすべて等分散の正規分布に従っていることが条件である。

検定統計量F

　一元配置分散分析と同様に二元配置分散分析の検定統計量は分散比のFである。計算は下記のように行う。基本的な考え方は一元配置分散分析と同じことがわかるだろう。全体のデータのばらつき（分散）を各説明変数（要因）で説明される部分と説明できない部分（残差）に分ける。異なるのは交互作用項に関する計算が増えていることくらいであろう。

交互作用の平均平方和の計算だけ特殊で、平均値からのずれの2乗和を計算するだけでなく、要因1と2の偏差平方和を引いている。これは、要因1×2の平均値からのずれの2乗和の中にはすでに要因1と2の偏差平方和が加算されているためである。

　そして、残差平均平方和に対する各要因や交互作用の平均平方和の比率を検定統計量Fとして計算する。各水準の標本が正規分布に従い、かつ、各水準の平均値が等しいときのF分布が知られているため、これを利用して検定を行うのである。

二元配置分散分析の検出力

　では、まずは検出力に関してシミュレートしてみよう。次のような設定で、解析に使うデータの生成を行う。

------------------------------------------------------

library(plotn)

library(multcomp)

cols <- col_genelator(number = 2, alpha = 1)

fils <- col_genelator(number = 2, alpha = 0.5)

nAX <- 30 #要因AXのサンプルサイズ

nAY <- 30 #要因AYのサンプルサイズ

nBX <- 30 #要因BXのサンプルサイズ

nBY <- 30 #要因BYのサンプルサイズ

meAX <- 70 #要因AXの平均

meAY <- 70 #要因AYの平均

meBX <- 50 #要因BXの平均

meBY <- 50 #要因BYの平均

s <- 15 #各要因の母標準偏差

mAX <- rnorm(nAX, mean = meAX, sd = s)#要因AかつXのデータの生成

mAY <- rnorm(nAY, mean = meAY, sd = s)#要因AかつYのデータの生成

mBX <- rnorm(nBX, mean = meBX, sd = s)#要因BかつXのデータの生成

mBY <- rnorm(nBY, mean = meBY, sd = s)#要因BかつYのデータの生成

d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図1の描画

------------------------------------------------------

設定から、要因A/Bは平均値に影響を与えているが、要因X/Yは影響を与えていないことを想定している。また、交互作用もなさそうであることがわかるだろう。このようなデータを10000回生成して、危険率5%の下での検出力を確認しよう。

------------------------------------------------------

#A > B

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 70

meAY <- 70

meBX <- 50

meBY <- 50

s <- 15

p1 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  mAX <- rnorm(nAX, mean = meAX, sd = s)

  mAY <- rnorm(nAY, mean = meAY, sd = s)

  mBX <- rnorm(nBX, mean = meBX, sd = s)

  mBY <- rnorm(nBY, mean = meBY, sd = s)

  d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

  fit1 <- aov(m ~ g1 * g2, data = d) #二元配置分散分析、交互作用を考慮する場合の表記

#g1 * g2 = g1 + g2 + g1:g2である。

#考慮しない場合は、m ~ g1 + g2とする。

  fit2 <- summary(fit1)

  p1 <- c(p1, fit2[[1]]$`Pr(>F)`[1])#g1のP値（要因A/B）

  p2 <- c(p2, fit2[[1]]$`Pr(>F)`[2])#g2のP値（要因X/Y）

  p3 <- c(p3, fit2[[1]]$`Pr(>F)`[3])#交互作用のP値（要因A/B）

}

histn(p1, xlab = "P value", ylab = "頻度") #図2の描画

histn(p2, xlab = "P value", ylab = "頻度") #図3の描画

histn(p3, xlab = "P value", ylab = "頻度") #図4の描画

sum(p1 < 0.05)#g1の5%以下のP値の数（要因A/B）

## [1] 10000

sum(p2 < 0.05)#g2の5%以下のP値の数（要因X/Y）

## [1] 476

sum(p3 < 0.05)#交互作用の5%以下のP値の数

## [1] 493

------------------------------------------------------

この時、要因A/Bに関しては100%の検出力であり、一方で要因X/Yおよび交互作用に関しては5%の危険率を維持できている。確かに、二元配置分散分析は要因の効果を分離して検出している。

　同様に要因X/Yだけに差があるとき……、としてもよいが、それは上記の議論と全く同じになるので、今度は要因A/Bも、要因X/Yも差があることを想定してみよう。ただし、交互作用はないものと考える。例えば、次のようにデータを生成する。

------------------------------------------------------

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 90

meAY <- 70

meBX <- 70

meBY <- 50

s <- 15

mAX <- rnorm(nAX, mean = meAX, sd = s)#要因AかつXのデータの生成

mAY <- rnorm(nAY, mean = meAY, sd = s)#要因AかつYのデータの生成

mBX <- rnorm(nBX, mean = meBX, sd = s)#要因BかつXのデータの生成

mBY <- rnorm(nBY, mean = meBY, sd = s)#要因BかつYのデータの生成

d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図5の描画

------------------------------------------------------

設定から、要因Bに対してAは平均が20増え、要因Yに対してXは平均が20増える。そして、要因BYに対して、要因AXは20 + 20 = 40程度、平均値が増加することを想定している。まさに、線形の関係にあることが強調されている。このようなデータを10000回生成して、危険率5%の下での検出力を確認しよう。

------------------------------------------------------

#A > B, X > Y

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 90

meAY <- 70

meBX <- 70

meBY <- 50

s <- 15

p1 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  mAX <- rnorm(nAX, mean = meAX, sd = s)

  mAY <- rnorm(nAY, mean = meAY, sd = s)

  mBX <- rnorm(nBX, mean = meBX, sd = s)

  mBY <- rnorm(nBY, mean = meBY, sd = s)

  d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

  fit1 <- aov(m ~ g1 * g2, data = d) #二元配置分散分析

  fit2 <- summary(fit1)

  p1 <- c(p1, fit2[[1]]$`Pr(>F)`[1])#g1のP値（要因A/B）

  p2 <- c(p2, fit2[[1]]$`Pr(>F)`[2])#g2のP値（要因X/Y）

  p3 <- c(p3, fit2[[1]]$`Pr(>F)`[3])#交互作用のP値（要因A/B）

}

histn(p1, xlab = "P value", ylab = "頻度") #図6の描画

histn(p2, xlab = "P value", ylab = "頻度") #図7の描画

histn(p3, xlab = "P value", ylab = "頻度") #図8の描画

sum(p1 < 0.05)#g1の5%以下のP値の数（要因A/B）

## [1] 10000

sum(p2 < 0.05)#g2の5%以下のP値の数（要因X/Y）

## [1] 10000

sum(p3 < 0.05)#交互作用の5%以下のP値の数

## [1] 504

------------------------------------------------------

この時、要因A/Bおよび要因X/Yのいずれも100%の検出力であり、一方で交互作用に関しては5%の危険率を維持できている。

　では今度は次のようにデータを生成する。何が起こっているか想像できるだろうか？

------------------------------------------------------

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 70

meAY <- 50

meBX <- 50

meBY <- 50

s <- 15

mAX <- rnorm(nAX, mean = meAX, sd = s)#要因AかつXのデータの生成

mAY <- rnorm(nAY, mean = meAY, sd = s)#要因AかつYのデータの生成

mBX <- rnorm(nBX, mean = meBX, sd = s)#要因BかつXのデータの生成

mBY <- rnorm(nBY, mean = meBY, sd = s)#要因BかつYのデータの生成

d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図9の描画

------------------------------------------------------

このとき、BX、AY、BYの時は平均値が変わらないのに、AXだけが平均値が異なっている。言い換えれば、要因A/Bや要因X/Yは単独では影響がないのに、要因AとXが組み合わさった時だけ、平均値が変わっている。こういう状態が、交互作用の存在を疑うべき状況だ。このようなデータを10000回生成して、危険率5%の下での検出力を確認しよう。

------------------------------------------------------

#A = B, X = Y, interaction

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 70

meAY <- 50

meBX <- 50

meBY <- 50

s <- 15

p1 <- NULL

p2 <- NULL

p3 <- NULL

pAY_AX <- NULL

pBX_AX <- NULL

pBY_AX <- NULL

pBX_AY <- NULL

pBY_AY <- NULL

pBY_BX <- NULL

for (i in 1:10000){

  mAX <- rnorm(nAX, mean = meAX, sd = s)

  mAY <- rnorm(nAY, mean = meAY, sd = s)

  mBX <- rnorm(nBX, mean = meBX, sd = s)

  mBY <- rnorm(nBY, mean = meBY, sd = s)

  d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

  d$g12 <- as.factor(paste0(d$g1, d$g2)) #交互作用の検討のためg1とg2を結合した文字列を作る

  fit1 <- aov(m ~ g1 * g2, data = d) #二元配置分散分析

  fit2 <- summary(fit1)

  fit3 <- aov(m ~ g12, data = d) #各要因組み合わせごとの一元配置分散分析

  fit4 <- summary(fit3)

  fit5 <- glht(fit3, linfct = mcp(g12 = "Tukey")) #各要因組み合わせごとのTukey法

  fit6 <- summary(fit5)

  p1 <- c(p1, fit2[[1]]$`Pr(>F)`[1])

  p2 <- c(p2, fit2[[1]]$`Pr(>F)`[2])

  p3 <- c(p3, fit2[[1]]$`Pr(>F)`[3])

  #以下Tukey-Kuramer法のP値

  pAY_AX <- c(pAY_AX, fit6$test$pvalues[1])

  pBX_AX <- c(pBX_AX, fit6$test$pvalues[2])

  pBY_AX <- c(pBY_AX, fit6$test$pvalues[3])

  pBX_AY <- c(pBX_AY, fit6$test$pvalues[4])

  pBY_AY <- c(pBY_AY, fit6$test$pvalues[5])

  pBY_BX <- c(pBY_BX, fit6$test$pvalues[6])

}

histn(p1, xlab = "P value", ylab = "頻度") #図10の描画

histn(p2, xlab = "P value", ylab = "頻度") #図11の描画

histn(p3, xlab = "P value", ylab = "頻度") #図12の描画

sum(p1 < 0.05)#g1の5%以下のP値の数（要因A/B）

## [1] 9538

sum(p2 < 0.05)#g2の5%以下のP値の数（要因X/Y）

## [1] 9523

sum(p3 < 0.05)#交互作用の5%以下のP値の数

## [1] 9535

------------------------------------------------------

すると、要因A/B、要因X/Yは単独では平均値が変わらなかったにもかかわらず、非常に高い検出力となった。これが、交互作用が存在するデータの解釈での注意点である。交互作用が有意な時、単純主効果も（見た目上）有意になりやすいのだ。交互作用を解析に組み込んでいなかった場合、本当は存在しない効果を見出してしまう可能性や深い考察が必要な効果を見逃してしまう可能性がある。

　ではこんな時、どうやって解析を続ければよいのだろう？　実はそのためのシミュレーションも併せて実行していることに気付いただろうか。まず、要因A/Bと要因X/Yの文字列を組み合わせた4水準（AX、AY、BX、BY）の説明変数を作成し、これを説明変数として一元配置分散分析を実行している。さらにその後、要因間の有意差を検出するためにTukey-Kramer法を合わせて実行している。要するに、2要因の解析から、各要因組み合わせごとの平均値差を検出するための解析に帰着しているのである。単純主効果だけを見ても全容がつかめないので、各要因組み合わせの平均値を解析する、というコンセプトだ。解析方法はこれだけに限らない。実験の設計によってはほかの解析方法もあり得るので、何を見たいのかをよく考えて解析を実行してみよう（後述）。では、以下のようにTukey-Kramer法による解析結果を図示してみよう。4水準の多重比較となるため、6つの比較の組み合わせがあることに注意する。

------------------------------------------------------

histn(pAY_AX, xlab = "P value", ylab = "頻度", breaks = "Sturges") #図13の描画

histn(pBX_AX, xlab = "P value", ylab = "頻度", breaks = "Sturges") #図14の描画

histn(pBY_AX, xlab = "P value", ylab = "頻度", breaks = "Sturges") #図15の描画

histn(pBX_AY, xlab = "P value", ylab = "頻度", breaks = "Sturges") #図16の描画

histn(pBY_AY, xlab = "P value", ylab = "頻度", breaks = "Sturges") #図17の描画

histn(pBY_BX, xlab = "P value", ylab = "頻度", breaks = "Sturges") #図18の描画

sum(pAY_AX < 0.05)

## [1] 9945

sum(pBX_AX < 0.05)

## [1] 9933

sum(pBY_AX < 0.05)

## [1] 9942

sum(pBX_AY < 0.05)

## [1] 118

sum(pBY_AY < 0.05)

## [1] 97

sum(pBY_BX < 0.05)

## [1] 109

------------------------------------------------------

解析の結果、データ生成時の仮定通り、要因組み合わせAXとそれ以外を比較した時、検出力は99%を超えている（図13～15）。その一方で、それ以外の比較では差がないことが正しいのであるが、これもちゃんと反映されて危険率は1%程度となっている。交互作用を考えるときは、今まで取り扱ってきた解析を総動員することで解釈を深めていくことができる。

二元配置分散分析の危険率

　では、次は次のようなデータを生成し、危険率を考えよう。

------------------------------------------------------

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 50

meAY <- 50

meBX <- 50

meBY <- 50

s <- 15

mAX <- rnorm(nAX, mean = meAX, sd = s)#要因AかつXのデータの生成

mAY <- rnorm(nAY, mean = meAY, sd = s)#要因AかつYのデータの生成

mBX <- rnorm(nBX, mean = meBX, sd = s)#要因BかつXのデータの生成

mBY <- rnorm(nBY, mean = meBY, sd = s)#要因BかつYのデータの生成

d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図19の描画

------------------------------------------------------

このようなデータを10000回生成して、危険率5%の下で分散分析をしよう。

------------------------------------------------------

#A = B, X = Y

nAX <- 30

nAY <- 30

nBX <- 30

nBY <- 30

meAX <- 50

meAY <- 50

meBX <- 50

meBY <- 50

s <- 15

p1 <- NULL

p2 <- NULL

p3 <- NULL

for (i in 1:10000){

  mAX <- rnorm(nAX, mean = meAX, sd = s)

  mAY <- rnorm(nAY, mean = meAY, sd = s)

  mBX <- rnorm(nBX, mean = meBX, sd = s)

  mBY <- rnorm(nBY, mean = meBY, sd = s)

  d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                  g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                  g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                         rep("X", times = nBX), rep("Y", times = nBY))

                  )#g1に要因A/B、g2に要因X/Yを格納してある。

  fit1 <- aov(m ~ g1 * g2, data = d) #二元配置分散分析

  fit2 <- summary(fit1)

  p1 <- c(p1, fit2[[1]]$`Pr(>F)`[1])#g1のP値（要因A/B）

  p2 <- c(p2, fit2[[1]]$`Pr(>F)`[2])#g2のP値（要因X/Y）

  p3 <- c(p3, fit2[[1]]$`Pr(>F)`[3])#交互作用のP値（要因A/B）

}

histn(p1, xlab = "P value", ylab = "頻度") #図20の描画

histn(p2, xlab = "P value", ylab = "頻度") #図21の描画

histn(p3, xlab = "P value", ylab = "頻度") #図22の描画

sum(p1 < 0.05)#g1の5%以下のP値の数（要因A/B）

## [1] 514

sum(p2 < 0.05)#g2の5%以下のP値の数（要因X/Y）

## [1] 517

sum(p3 < 0.05)#交互作用の5%以下のP値の数

## [1] 485

------------------------------------------------------

もう予想通りかと思うが、確かに要因A/B、要因X/Yおよび交互作用に関して、危険率を5%に抑えられている。二元配置分散分析は、このように複数の要因を説明変数として組み込み、平均値に影響を与える要因を分離するのに役立つわけである。

Rで行う二元配置分散分析

　では、次のようなデータに関して、二元配置分散分析を実行してみよう。要因は2つ存在し、まずは、何事もデータの図示からだ。

------------------------------------------------------

mAX <- c(15, 63, 86, 37, 55, 24, 10, 47, 69, 64, 

         58, 16, 51, 37, 64, 22, 39, 65, 57, 49, 

         63, 67, 67, 69, 52, 34, 26, 61, 67, 40)

mAY <- c(27, 82, 76, 42, 68, 40, 94, 48, 79, 46, 

         72, 39, 18, 58, 86, 56, 66, 56, 71, 56, 

         54, 64, 70, 70, 30, 30, 19, 58, 49, 55)

mBX <- c(49, 23, 19, 89, 39, 77, 12, 61, 32, 84,

         44, 43, 33, 48, 53, 36, 40, 69, 50, 58, 

         89, 52, 58, 37, 43, 20, 60, 70, 68, 82)

mBY <- c(22, 48, 57, 37, 42, 46, 73, 34, 35, 84, 

         10, 85, 105, 56, 26, 31, 57, 76, 47, 58, 

         45, 61, 35, 69, 89, 70, 53, 106, 27, 55)

nAX <- length(mAX)

nAY <- length(mAY)

nBX <- length(mBX)

nBY <- length(mBY)

d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                       rep("X", times = nBX), rep("Y", times = nBY))

                )#g1に要因A/B、g2に要因X/Yを格納してある。

head(d)

##    m g1 g2

## 1 15  A  X

## 2 63  A  X

## 3 86  A  X

## 4 37  A  X

## 5 55  A  X

## 6 24  A  X

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図23の描画

------------------------------------------------------

なんとなく、各要因の効果は弱そうだ。データのばらつきは、等しそうで、正規分布を仮定できなさそうなほどの偏りもない。つまり、仮定が満たされていそうなので、二元配置分散分析を実行しよう。

------------------------------------------------------

fit1 <- aov(m ~ g1 * g2, data = d) #二元配置分散分析

fit2 <- summary(fit1)

fit1

## Call:

##    aov(formula = m ~ g1 * g2, data = d)

## 

## Terms:

##                       g1       g2    g1:g2 Residuals

## Sum of Squares      4.80   780.30    90.13  51329.27

## Deg. of Freedom        1        1        1       116

## 

## Residual standard error: 21.03553

## Estimated effects may be unbalanced

fit2

##              Df Sum Sq Mean Sq F value Pr(>F)

## g1            1      5     4.8   0.011  0.917

## g2            1    780   780.3   1.763  0.187

## g1:g2         1     90    90.1   0.204  0.653

## Residuals   116  51329   442.5

------------------------------------------------------

上記のような解析結果が得られた。fit1とfit2は共通した要素が入っており、より見やすく調整されたものがfit2に入っている。fit2に格納した解析結果を確認すると、Pr(>F)と書かれている部分がいわゆるP値で、いずれも5%を上回っており、危険率5%のもとではいずれの要因および交互作用も有意な効果を持たないことがわかる。fit2の最左の列には要因が格納されていて、要因A/B = g1、要因X/Y = g2および要因間の交互作用 = g1:g2である。また、最後のResidualsが残差である。Dfは自由度、Sum Sqは偏差平方和、Mean Sqは平均平方和、F valueは分散比であるF値だ。今回は、各要因内の水準がともに2なので、各要因の自由度が1となり、偏差平方和と偏差平方和/自由度 = 平均平方和が一致している。

　では、各要素を手計算してみよう。各自由度は以下のように計算する。

------------------------------------------------------

dfT <- nAX + nBX + nAY + nBY - 1

dfAB <- 2 - 1

dfXY <- 2 - 1

dfIn <- dfAB * dfXY

dfR <- dfT - dfAB - dfXY - dfIn

dfT #全体の自由度

## [1] 119

dfAB #要因ABの自由度

## [1] 1

dfXY #要因XYの自由度

## [1] 1

dfIn #交互作用の自由度

## [1] 1

dfR #残差の自由度

## [1] 116

------------------------------------------------------

次に偏差平方和と平均平方和を計算しよう。各要因ごとの平均値を求めることから始める。

------------------------------------------------------

mT <- mean(c(mAX, mAY, mBX, mBY)) #全体の平均

m_mA <- mean(c(mAX, mAY)) #水準Aの平均

m_mB <- mean(c(mBX, mBY)) #水準Bの平均

m_mX <- mean(c(mAX, mBX)) #水準Xの平均

m_mY <- mean(c(mAY, mBY)) #水準Yの平均

m_mAX <- mean(mAX) #水準AかつXの平均

m_mAY <- mean(mAY) #水準AかつYの平均

m_mBX <- mean(mBX) #水準BかつXの平均

m_mBY <- mean(mBY) #水準BかつYの平均

ST <- sum((c(mAX, mAY, mBX, mBY) - mT)^2) #全体の平方和

SAB <- sum(((m_mA - mT)^2)*(nAX+nAY), ((m_mB - mT)^2)*(nBX+nBY)) #要因ABの平方和

SXY <- sum(((m_mX - mT)^2)*(nAX+nBX), ((m_mY - mT)^2)*(nAY+nBY)) #要因XYの平方和

SIn <- sum(((m_mAX - mT)^2)*nAX, ((m_mAY - mT)^2)*nAY, 

           ((m_mBX - mT)^2)*nBX, ((m_mBY - mT)^2)*nBY) - SAB - SXY #交互作用の平方和

SR <- ST - SAB - SXY - SIn

sT <- ST/dfT #全体の平均平方和

sAB <- SAB/dfAB #要因ABの平均平方和

sXY <- SXY/dfXY #要因XYの平均平方和

sIn <- SIn/dfIn #交互作用の平均平方和

sR <- SR/dfR  #残差の平均平方和

ST

## [1] 52204.5

SAB

## [1] 4.8

SXY

## [1] 780.3

SIn

## [1] 90.13333

SR

## [1] 51329.27

sT

## [1] 438.6933

sAB

## [1] 4.8

sXY

## [1] 780.3

sIn

## [1] 90.13333

sR

## [1] 442.4937

------------------------------------------------------

各平均平方和が求まったので、これらからF値を計算しよう。

------------------------------------------------------

FAB <- sAB/sR #要因ABに関するF値

FXY <- sXY/sR #要因XYに関するF値

FIn <- sIn/sR #交互作用に関するF値

FAB

## [1] 0.01084761

FXY

## [1] 1.763415

FIn

## [1] 0.2036941

------------------------------------------------------

最後、F値が求まったので、これをもとにP値を求める。

------------------------------------------------------

pf(df1 = dfAB, df2 = dfR, q = FAB, lower.tail = F) #要因ABに関するP値

## [1] 0.9172287

pf(df1 = dfXY, df2 = dfR, q = FXY, lower.tail = F) #要因XYに関するP値

## [1] 0.1868062

pf(df1 = dfIn, df2 = dfR, q = FIn, lower.tail = F) #交互作用に関するP値

## [1] 0.6525981

------------------------------------------------------

いずれも一致していることがわかるだろうか。無事、二元配置分散分析を実行できた。

　では、もう一つの例題を考えてみよう。次のようなデータだ。まずは図示しよう。

------------------------------------------------------

mAX <- c(81, 61, 69, 51, 27, 72, 70, 78, 29, 60,

         81, 66, 64, 71, 104, 50, 68, 84, 92, 71,

         46, 64, 72, 61, 38, 84, 54, 103, 73, 83)

mAY <- c(25, 26, 36, 20, 32, 55, 34, 3, 42, 34, 

         44, 53, 26, 28, 5, 10, 11, 9, 21, 47, 

         15, 15, 19, 18, 22, 11, -5, 38, 34, 35)

mBX <- c(33, 23, 26, 18, 0, 12, 2, 22, 20, 19, 

         41, 37, 27, 53, 42, 5, 17, 12, 8, 14,

         45, 24, 17, 32, 41, 36, -3, 22, 14, 32)

mBY <- c(21, 73, 58, 34, 40, 24, 65, 57, 60, 47,

         47, 41, 26, 81, 54, 31, 31, 65, 53, 65,

         76, 45, 71, 50, 88, 83, 64, 46, 66, 72)

nAX <- length(mAX)

nAY <- length(mAY)

nBX <- length(mBX)

nBY <- length(mBY)

d <- data.frame(m = c(mAX, mAY, mBX, mBY), 

                g1 = c(rep("A", times = nAX+nAY), rep("B", times = nBX+nBY)),

                g2 = c(rep("X", times = nAX), rep("Y", times = nAY),

                       rep("X", times = nBX), rep("Y", times = nBY))

                )

head(d)

##    m g1 g2

## 1 81  A  X

## 2 61  A  X

## 3 69  A  X

## 4 51  A  X

## 5 27  A  X

## 6 72  A  X

boxplotn(m ~ g1 + g2, data = d, ylab = "data", jitter.method = "swarm", 

         cex.dot = 1, pch.dot = c(16,16,17,17), col.dot = cols, 

         col.bor = cols, col.fill = fils) #図24の描画

------------------------------------------------------

何やら今度は交互作用が出ていそうな気配がしている。二元配置分散分析を行って確認してみよう。

------------------------------------------------------

fit1 <- aov(m ~ g1 * g2, data = d) #二元配置分散分析

fit2 <- summary(fit1)

fit2

##              Df Sum Sq Mean Sq F value Pr(>F)    

## g1            1   1802    1802   6.551 0.0118 *  

## g2            1    859     859   3.122 0.0799 .  

## g1:g2         1  40590   40590 147.568 <2e-16 ***

## Residuals   116  31907     275                   

## ---

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

------------------------------------------------------

確かに交互作用に有意な効果が認められている。こんな時は、単純主効果だけに目を向けるのではなく、要因の組み合わせに注意を払わなければならない。練習として、F値などの各要素を手計算してみよう。上記の方法にならえばすぐにできるはずだ。

　さて、こんな時の解析方針は個人的には2つあると思っている。一つ目はすでに紹介した、要因の組み合わせを説明変数として、一元配置分散分析に帰着する方法である。例えば、次のように実行する。

------------------------------------------------------

d$gs <- as.factor(paste0(d$g1, d$g2))

fit3 <- aov(m ~ gs, data = d) #要因組み合わせに対する一元配置分散分析

fit4 <- summary(fit3)

fit3

## Call:

##    aov(formula = m ~ gs, data = d)

## 

## Terms:

##                       gs Residuals

## Sum of Squares  43250.96  31907.17

## Deg. of Freedom        3       116

## 

## Residual standard error: 16.58499

## Estimated effects may be unbalanced

fit4

##              Df Sum Sq Mean Sq F value Pr(>F)    

## gs            3  43251   14417   52.41 <2e-16 ***

## Residuals   116  31907     275                   

## ---

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

fit5 <- glht(model = fit3, linfct = mcp(gs = "Tukey")) #要因組み合わせに対するTukey法

fit6 <- summary(fit5)

fit5

## 

##   General Linear Hypotheses

## 

## Multiple Comparisons of Means: Tukey Contrasts

## 

## 

## Linear Hypotheses:

##              Estimate

## AY - AX == 0   -42.13

## BX - AX == 0   -44.53

## BY - AX == 0   -13.10

## BX - AY == 0    -2.40

## BY - AY == 0    29.03

## BY - BX == 0    31.43

fit6

## 

##   Simultaneous Tests for General Linear Hypotheses

## 

## Multiple Comparisons of Means: Tukey Contrasts

## 

## 

## Fit: aov(formula = m ~ gs, data = d)

## 

## Linear Hypotheses:

##              Estimate Std. Error t value Pr(>|t|)    

## AY - AX == 0  -42.133      4.282  -9.839   <0.001 ***

## BX - AX == 0  -44.533      4.282 -10.400   <0.001 ***

## BY - AX == 0  -13.100      4.282  -3.059   0.0145 *  

## BX - AY == 0   -2.400      4.282  -0.560   0.9435    

## BY - AY == 0   29.033      4.282   6.780   <0.001 ***

## BY - BX == 0   31.433      4.282   7.340   <0.001 ***

## ---

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

## (Adjusted p values reported -- single-step method)

cld(fit5)

##  AX  AY  BX  BY 

## "c" "a" "a" "b"

------------------------------------------------------

以上のようなTukey-Kramer法を絡めた各要因組み合わせの多重比較により、要因組み合わせAXが最も平均値が大きく、次にBY、そしてAY = BXと続くことがわかる。これは図24と見比べても、則していると言えるだろう。このような解析が必要になるときは、すべての要因組み合わせを、いうなれば平等に見ているときだ。例えば、2つの薬剤の投与の有無にとって、被験者の状態がどのように変わるか、であれば、単独投与が効果があるのか、複数投与が良いのか、あるいはダメなのか、比較したくなるだろう。こんな時には、上記のように要因組み合わせを説明変数にするとよい。

　もう一つの方法は、一方の要因内でもう一方の要因の効果を検討する方法である。例えば、要因A内で要因X/Yの効果を検討し、要因B内で同様に要因X/Yの効果を検討する。例えば、次のようにする。

------------------------------------------------------

A <- d[d$g1 == "A",]

B <- d[d$g1 == "B",]

fit7 <- aov(m ~ g2, data = A) #要因A内で要因X/Yの効果を検討

fit8 <- summary(fit7)

fit7

## Call:

##    aov(formula = m ~ g2, data = A)

## 

## Terms:

##                       g2 Residuals

## Sum of Squares  26628.27  16436.73

## Deg. of Freedom        1        58

## 

## Residual standard error: 16.83425

## Estimated effects may be unbalanced

fit8

##             Df Sum Sq Mean Sq F value   Pr(>F)    

## g2           1  26628   26628   93.96 9.71e-14 ***

## Residuals   58  16437     283                     

## ---

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

fit9 <- aov(m ~ g2, data = B) #要因B内で要因X/Yの効果を検討

fit10 <- summary(fit9)

fit9

## Call:

##    aov(formula = m ~ g2, data = B)

## 

## Terms:

##                       g2 Residuals

## Sum of Squares  14820.82  15470.43

## Deg. of Freedom        1        58

## 

## Residual standard error: 16.33192

## Estimated effects may be unbalanced

fit10

##             Df Sum Sq Mean Sq F value   Pr(>F)    

## g2           1  14821   14821   55.56 5.05e-10 ***

## Residuals   58  15470     267                     

## ---

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

------------------------------------------------------

確かに、fit8やfit10の中身を見ると、要因A内で要因X/Y、要因B内で要因X/Yがともに有意な効果を持っていることがわかる。これも図24から見て取れるだろう。2標本の比較なんだからt検定を使うべきでは？と思われるかもしれないが、実は2標本の分散分析とt検定は全く同じP値を返す。途中仮定は違っても同じ仮定をもとに計算されるP値なので一致するのだ。このような解析が必要となるのは、例えば要因内の効果を重視しているときや比較の必要がない要因組み合わせが存在するときが該当する。例えば、植物の競争力の指標として植物の高さを測っているとしよう。乾いた条件での植物A種B種の高さの比較、湿った条件での植物A種B種の高さの比較、をしたいときがそうだ。このとき、興味があるのは植物A種とB種の競争なので、同じ条件にない植物（例えば、乾いた条件のA種vs湿った条件のB種）を比較する興味はない（なぜなら、競争することがないから）、ということだ。

　以上のように交互作用の解析には、これ1択というような正解はない。何を見たいかによって解析方法を考えていく必要がある。また、一元配置分散分析のときも話をしたが、必ずしも分散分析ののちにそのあとの解析をする必要はない。初めから、すべての条件を平等に見ているのであれば、そのように解析してもよいし、要因内のもう一つの要因の効果に興味があるなら、それを前提の解析をあらかじめしてもよい。一方で、要因を分解することでそれぞれの影響を定量化できる点や、交互作用があることを定量的に示す点に興味があるなら、分散分析する価値も十分にある。今回は、始めに紹介したように説明変数の係数を求めることなく解析を行った。これは、aov関数の出力の特徴で、内部的には計算されている（実際、Tukey-Kramer法をおこなうと係数が出てくる）。もし、係数を知りたい場合は、今後紹介するようにlm関数を用いればよい。