【統計】二元配置分散分析（Two-way ANOVA）

二元配置分散分析（Two-way ANOVA）について。

【目次】

帰無仮説、対立仮設
計算式等
計算例
プログラムコード
交互作用のイメージ
参考

一元配置分散分析は１要因に対する分析手法。
２要因の場合は、二元配置分散分析を用いる。
※本記事の二元配置分散分析は、対応なしの被験者間実験計画（ABSタイプ）。

帰無仮説、対立仮設

要因Aの主効果

帰無仮説 $H_{0}:$ 要因Aの水準間の主効果が同じ（平均が同じ）
対立仮設 $H_{1}:$ 要因Aの水準間のうち、少なくとも一つの水準で主効果に差がある（平均に差がある）

要因Bの主効果

帰無仮説 $H_{0}:$ 要因Bの水準間の主効果が同じ（平均が同じ）
対立仮設 $H_{1}:$ 要因Bの水準間のうち、少なくとも一つの水準で主効果に差がある（平均に差がある）

要因Aと要因Bの交互作用

帰無仮説 $H_{0}:$ 要因Aと要因Bに交互作用がある
対立仮設 $H_{1}:$ 要因Aと要因Bに交互作用がない

二元配置分散分析では、要因Aの水準間の平均の違い、要因Bの水準間の平均の違い、要因AとBの交互作用の有無を評価できる。

また、要因AとBに交互作用が無いと分かっている場合は交互作用を含めないモデルで分析可能。交互作用が無いとは、例えば「要因A1の時の要因B1とB2の差」と「要因A2の時の要因B1とB2の差」に違いがないと言うこと。

計算式等

※本計算は各水準のサンプルサイズが等しいことを前提とする。また、対応なしの被験者間実験計画（ABSタイプ）。

交互作用項を含めた分析
■ 平方和
$\quad S_{A} = \sum _{k=1} ^{a} n_{k} \cdot \left( \overline{y_{k}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{B} = \sum _{m=1} ^{b} n_{m} \cdot \left( \overline{y_{m}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{AB} = \sum _{k=1} ^{a} \sum _{m=1} ^{b} n_{km} \cdot \left( \overline{y_{km}} - \overline{y_{k}} - \overline{y_{m}} + \bar{\bar{y}} \right)^{2}$
$\quad S_{e} = \sum _{k=1} ^{a} \sum _{m=1} ^{b} \sum _{i=1} ^{n_{km}} \left( y_{km_{i}} - \overline{y_{km}} \right)^{2}$
$\quad S_{T} = \sum _{i=1} ^{n} \left( y_{i} - \bar{\bar{y}} \right)^{2}$

■ 分散分析表

変動因	平方和	自由度	平均平方	F値（F比）
水準間A	$S_{A}$	$f_{A}=a-1$	$V_{A}=\dfrac{S_{A}}{f_{A}}$	$F_{A}=\dfrac{V_{A}}{V_{e}}$
水準間B	$S_{B}$	$f_{B}=b-1$	$V_{B}=\dfrac{S_{B}}{f_{B}}$	$F_{B}=\dfrac{V_{B}}{V_{e}}$
交互作用AB	$S_{AB}$	$f_{AB}=a\cdot b-a-b+1$	$V_{AB}=\dfrac{S_{AB}}{f_{AB}}$	$F_{AB}=\dfrac{V_{AB}}{V_{e}}$
残差e	$S_{e}$	$f_{e}=n-a \cdot b$	$V_{e}=\dfrac{S_{e}}{f_{e}}$
全体T	$S_{T}$	$f_{T}=n-1$

※ k=要因Aの特定水準、m=要因Bの特定水準、i=特定の観測値、 $\overline{y_{k}}$ =要因Aの特定水準の平均値、 $\overline{y_{m}}$ =要因Bの特定水準の平均値、 $\bar{\bar{y}}$ =全体の平均値（総平均）、a=要因Aの水準数、b=要因Bの水準数

ここで平方和と自由度には次の関係がある。
$\quad S_{T} = S_{A} + S_{B} + S_{AB} + S_{e}$
$\quad f_{T} = f_{A} + f_{B} + f_{AB} + f_{e}$

交互作用項を含めない分析
■ 平方和　※ $S_{e}$ 以外同じ。交互作用項を残差項にプールしている。
$\quad S_{A} = \sum _{k=1} ^{a} n_{k} \cdot \left( \overline{y_{k}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{B} = \sum _{m=1} ^{b} n_{m} \cdot \left( \overline{y_{m}} - \bar{\bar{y}} \right)^{2}$
$\quad S_{e} = \sum _{k=1} ^{a} \sum _{m=1} ^{b} \sum _{i=1} ^{n_{km}} \left( y_{km_{i}} - \overline{y_{k}} - \overline{y_{m}} + \bar{\bar{y}} \right)^{2}$
$\quad S_{T} = \sum _{i=1} ^{n} \left( y_{i} - \bar{\bar{y}} \right)^{2}$

■ 分散分析表　※ $S_{e}$ 以外同じ。

変動因	平方和	自由度	平均平方	F値（F比）
水準間A	$S_{A}$	$f_{A}=a-1$	$V_{A}=\dfrac{S_{A}}{f_{A}}$	$F_{A}=\dfrac{V_{A}}{V_{e}}$
水準間B	$S_{B}$	$f_{B}=b-1$	$V_{B}=\dfrac{S_{B}}{f_{B}}$	$F_{B}=\dfrac{V_{B}}{V_{e}}$
残差e	$S_{e}$	$f_{e}=n-a-b+1$	$V_{e}=\dfrac{S_{e}}{f_{e}}$
全体T	$S_{T}$	$f_{T}=n-1$

ここで平方和と自由度には次の関係がある。
$\quad S_{T} = S_{A} + S_{B} + S_{e}$
$\quad f_{T} = f_{A} + f_{B} + f_{e}$

また、質的変数をダミー変数（０と１の数量変数）にして重回帰分析を行う（正規方程式を解く）ことで、各水準の条件付き期待値を算出できる（数量化理論１類）。
一元配置分散分析の記事を参照。

計算例

簡単なデータで計算方法を確認。
今回は薬剤（A薬、B薬、C薬の３水準）と前処置（有、無の２水準）の２要因において、評価指標Xの結果に差があるかを見る（架空データ）。交互作用項を含めた分析を行う。なお、評価指標Xは特定のものではないが、低いほど良いとしておく。

　A薬＋前処置無 {20,21,22,20,21}
　A薬＋前処置有 {18,16,17,16,17}
　B薬＋前処置無 {19,18,17,19,19}
　B薬＋前処置有 {19,18,18,17,17}
　C薬＋前処置無 {20,18,16,19,18}
　C薬＋前処置有 {15,14,17,16,18}

まず、全体・薬剤別・前処置別で平均を求め、総平均からのズレ（偏差平方和）を求める。

薬剤	前処置	N	平均	要因:薬剤によるズレ（偏差平方和）と要因:前処置によるズレ（偏差平方和）
全体	全体	30	18.0	---
A薬	全体	10	18.8	$\left(18.8-18.0\right)^{2}\cdot 10=6.4$
B薬	全体	10	18.1	$\left(18.1-18.0\right)^{2}\cdot 10=0.1$
C薬	全体	10	17.1	$\left(17.1-18.0\right)^{2}\cdot 10=8.1$
全体	無	15	19.1	$\left(19.1-18.0\right)^{2}\cdot 10≒19.267$
全体	有	15	16.9	$\left(16.9-18.0\right)^{2}\cdot 10≒19.267$

同じく、薬剤×前処置別で平均を求め、総平均・各要因の平均からのズレ（交互作用に関する変動）を求める。

薬剤	前処置	N	平均	総平均・各要因の平均からのズレ（交互作用に関する変動）
A薬	無	5	20.8	$\left(20.8-18.8-19.1+18.0\right)^{2}\cdot 5≒3.756$
A薬	有	5	16.8	$\left(16.8-18.8-16.9+18.0\right)^{2}\cdot 5≒3.756$
B薬	無	5	18.4	$\left(18.4-18.1-19.1+18.0\right)^{2}\cdot 5≒3.472$
B薬	有	5	17.8	$\left(17.8-18.1-16.9+18.0\right)^{2}\cdot 5≒3.472$
C薬	無	5	18.2	$\left(18.2-17.1-19.1+18.0\right)^{2}\cdot 5≒0.006$
C薬	有	5	16.0	$\left(16.0-17.1-16.9+18.0\right)^{2}\cdot 5≒0.006$

各表の最右列のズレを足し合わせて、各要因の水準間平方和と交互作用の平方和を計算する。
$\quad S_{A}=6.4+0.1+8.1=14.6$
$\quad S_{B}=19.267+19.267≒38.533$
$\quad S_{AB}=3.756+3.756+3.472+3.472+0.006+0.006=14.467$

また、今回は薬剤が３水準（A薬、B薬、C薬）、前処置が２水準（有、無）のため自由度は次のとおり。
$\quad f_{A}=3-1=2$
$\quad f_{B}=2-1=1$
$\quad f_{AB}=3\times2-3-2+1=2$

次に、個々の観測値について、各水準の平均からのズレ（群内の偏差平方）と総平均からのズレ（全体の偏差平方）を求める。
※ $S_{T} = S_{A} + S_{B} + S_{AB} + S_{e}$ なので $S_{T}$ か $S_{e}$ の１つ分かれば十分だが、計算の理解のため両方求める。

薬剤	前処置	評価指標	各水準の平均からのズレ（群内の偏差平方）	総平均からのズレ全体の偏差平方
A薬	無	20	$(20 - 20.8)^{2}=0.64$	$(20-18)^{2}=4$
A薬	無	21	$(21 - 20.8)^{2}=0.04$	$(21-18)^{2}=9$
A薬	無	22	$(22 - 20.8)^{2}=1.44$	$(22-18)^{2}=16$
A薬	無	20	$(20 - 20.8)^{2}=0.64$	$(20-18)^{2}=4$
A薬	無	21	$(21 - 20.8)^{2}=0.04$	$(21-18)^{2}=9$
A薬	有	18	$(18 - 16.8)^{2}=1.44$	$(18-18)^{2}=0$
A薬	有	16	$(16 - 16.8)^{2}=0.64$	$(16-18)^{2}=4$
A薬	有	17	$(17 - 16.8)^{2}=0.04$	$(17-18)^{2}=1$
A薬	有	16	$(16 - 16.8)^{2}=0.64$	$(16-18)^{2}=4$
A薬	有	17	$(17 - 16.8)^{2}=0.04$	$(17-18)^{2}=1$
B薬	無	19	$(19 - 18.4)^{2}=0.36$	$(19-18)^{2}=1$
B薬	無	18	$(18 - 18.4)^{2}=0.16$	$(18-18)^{2}=0$
B薬	無	17	$(17 - 18.4)^{2}=1.96$	$(17-18)^{2}=1$
B薬	無	19	$(19 - 18.4)^{2}=0.36$	$(19-18)^{2}=1$
B薬	無	19	$(19 - 18.4)^{2}=0.36$	$(19-18)^{2}=1$
B薬	有	19	$(19 - 17.8)^{2}=1.44$	$(19-18)^{2}=1$
B薬	有	18	$(18 - 17.8)^{2}=0.04$	$(18-18)^{2}=0$
B薬	有	18	$(18 - 17.8)^{2}=0.04$	$(18-18)^{2}=0$
B薬	有	17	$(17 - 17.8)^{2}=0.64$	$(17-18)^{2}=1$
B薬	有	17	$(17 - 17.8)^{2}=0.64$	$(17-18)^{2}=1$
C薬	無	20	$(20 - 18.2)^{2}=3.24$	$(20-18)^{2}=4$
C薬	無	18	$(18 - 18.2)^{2}=0.04$	$(18-18)^{2}=0$
C薬	無	16	$(16 - 18.2)^{2}=4.84$	$(16-18)^{2}=4$
C薬	無	19	$(19 - 18.2)^{2}=0.64$	$(19-18)^{2}=1$
C薬	無	18	$(18 - 18.2)^{2}=0.04$	$(18-18)^{2}=0$
C薬	有	15	$(15 - 16)^{2}=1$	$(15-18)^{2}=9$
C薬	有	14	$(14 - 16)^{2}=4$	$(14-18)^{2}=16$
C薬	有	17	$(17 - 16)^{2}=1$	$(17-18)^{2}=1$
C薬	有	16	$(16 - 16)^{2}=0$	$(16-18)^{2}=4$
C薬	有	18	$(18 - 16)^{2}=4$	$(18-18)^{2}=0$

このうち、各水準の平均からのズレ（群内の偏差平方）を足し合わせたものが残差平方和で、総平均からのズレ（全体の偏差平方）を足し合わせたものが全体平方和。
$\quad S_{e}≒0.64 + 0.04 + ... + 0 + 4≒30.4$
$\quad S_{T}≒4+9+...+4+0≒98$

また、全体のN数が30、薬剤が３水準（A薬、B薬、C薬）、前処置が２水準（有、無）のため自由度は次のとおり。
$\quad f_{e}=n-a\cdot b=30-3\cdot2=24$
$\quad f_{T}=n-1=30-1=29$

最後に、これら平方和・自由度を分散分析表にまとめて F値と p値を算出する。

変動因	平方和	自由度	平均平方	F値（F比）	p値
薬剤	$S_{A}=14.6$	$f_{A}=2$	$V_{A}=7.3$	$F_{A}≒5.7632$	$0.0090$
前処置	$S_{B}≒38.533$	$f_{B}=1$	$V_{B}≒38.533$	$F_{B}≒30.4211$	$0.0000$
交互作用	$S_{AB}≒14.467$	$f_{AB}=2$	$V_{AB}≒7.233$	$F_{AB}≒5.7105$	$0.0094$
残差e	$S_{e}=30.4$	$f_{e}=24$	$V_{e}≒1.267$
全体T	$S_{T}=98$	$f_{T}=29$

分散分析表のとおり、今回の架空データでは薬剤と前処置の交互作用が p<0.01 で有意であり、これら２要因の間には交互作用がある様子（何れかの組み合わせで相乗効果か相殺効果がある）。

なお、薬剤の主効果は p<0.01、前処置の主効果は p<0.001 でともに統計学的に有意だったが、交互作用があるため参考程度にとどめておく（独立していないため、薬剤と前処置を切り離して考えない方が良い）。

交互作用プロットを見ると、薬剤毎の前処置有無の平均の差が一定ではないため、交互作用が有意になったと考えられる。特にB薬は前処置に関する効果が薄いように見える。

f:id:cochineal19:20210220141315p:plain

※変数名：TRT01AN=薬剤、水準：1=A薬、2=B薬、3=C薬
　変数名：PRCFL=前処置有無、水準：0=無、1=有
　変数名：AVAL=評価指標X

二元配置分散分析の結果の解釈と進め方

今回の架空データでは交互作用ありだったが、二元配置分散分析の結果によって次のアクションに違いがあるため、備忘録がてらメモ。
※多重比較を事後検定と認識するのはどうだ、交互作用があったら単純主効果検定など色々な見解を見かけましたが、こちらのサイト様が分かりやすかったので、参考にさせていただきました。

１．交互作用が認められる
　⇒　薬剤と前処置が互いに影響を与えていると結論。
　　　必要であれば薬剤×前処置の６つの組み合わせで多重比較。
　　　各要因を切り離して考えない。

２．交互作用が認められない＆薬剤の主効果が有意
　⇒　薬剤に有意な変動があったと結論。
　　　薬剤が３水準なので、必要であれば多重比較。

３．交互作用が認められない＆前処置の主効果が有意

　⇒　前処置に有意な変動があったと結論。
　　　２水準なので多重比較する必要なし。

４．交互作用が認められない＆薬剤・前処置ともに主効果が有意

　⇒　薬剤・前処置にともに有意な変動があったと結論。
　　　加えて、薬剤×前処置に交互作用がなく互いに影響を与えていないと判断。
　　　必要であれば、薬剤が３水準なので、薬剤について多重比較。
　　　前処置は２水準なので完結。

Rでの実行：

> library(car) 
> ADS <- data.frame(
+   TRT01AN=c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2)
+   ,PRCFL=c(0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,1,1,1,1,1)
+   ,AVAL=c(20,21,22,20,21,18,16,17,16,17,19,18,17,19,19,19,18,18,17,17,20,18,16,19,18,15,14,17,16,18))
> ADS$前処置 <- factor(ifelse(ADS$PRCFL==0,"無","有"),levels=c("無","有"))
> ADS$投与群 <- factor(ifelse(ADS$TRT01AN==0,"A薬",ifelse(ADS$TRT01AN==1,"B薬","C薬")),levels=c("A薬","B薬","C薬"))
> fit1 <- lm(AVAL ~ 投与群 + 前処置 + 投与群:前処置, data=ADS)
> car::Anova(fit1, Type="II")
Anova Table (Type II tests)

Response: AVAL
              Sum Sq Df F value    Pr(>F)    
投与群        14.600  2  5.7632  0.009035 ** 
前処置        38.533  1 30.4211 1.134e-05 ***
投与群:前処置 14.467  2  5.7105  0.009363 ** 
Residuals     30.400 24                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SASでの実行：

data ADS;
 do TRT01AN = 1; do PRCFL = 0; do AVAL = 20,21,22,20,21; output; end; end; end;
 do TRT01AN = 1; do PRCFL = 1; do AVAL = 18,16,17,16,17; output; end; end; end;
 do TRT01AN = 2; do PRCFL = 0; do AVAL = 19,18,17,19,19; output; end; end; end;
 do TRT01AN = 2; do PRCFL = 1; do AVAL = 19,18,18,17,17; output; end; end; end;
 do TRT01AN = 3; do PRCFL = 0; do AVAL = 20,18,16,19,18; output; end; end; end;
 do TRT01AN = 3; do PRCFL = 1; do AVAL = 15,14,17,16,18; output; end; end; end;
run;
proc glm data=ADS;
 class TRT01AN PRCFL;
 model AVAL = TRT01AN PRCFL TRT01AN*PRCFL / ss1 ss2 ss3 e1 e2 e3 solution; 
run;

f:id:cochineal19:20210220140739p:plain

f:id:cochineal19:20210220140807p:plain

プログラムコード

Rのコード

library(car)
fit1 <- lm(AVAL ~ FACTOR1 + FACTOR2 + FACTOR1:FACTOR2, data=ADS) #-- FACTOR1:FACTOR2は交互作用項
car::Anova(fit1, Type="II") #-- Type2平方和

SASのコード

proc glm data=ADS1 outstat=OUTDS;
  class FACTOR1; /*class指定必須*/
  model AVAL = FACTOR1 FACTOR2 FACTOR1*FACTOR2 / ss1 ss2 ss3 e1 e2 e3 solution;
run;
/* option: ss1=Type1平方和, ss2=Type2平方和, ss3=Type3平方和 */
/* option: e#= Type#平方和のEstimable Functions（推定可能関数）を図示してくれる */
/* option: solution=回帰モデルの推定値を出力してくれる */

Pythonのコード