二元配置分散分析の理論と実践：数理モデルからMATLABによる実装まで

1. はじめに

分散分析（Analysis of Variance, ANOVA）は、統計学者であるロナルド・エイルマー・フィッシャー（Ronald Aylmer Fisher）によって開拓された、現代統計学の重要な概念の一つです。フィッシャーがロザムステッド農業試験場において、1840年代から蓄積された膨大な農業データを分析する中で、複数の要因（例えば、異なる肥料や作物の品種）がもたらす影響を同時に、かつ系統的に分離・評価する必要に迫られたことにより、この手法は開発されました。この歴史的背景は、分散分析が単なる抽象的な数学的理論ではなく、現実世界の複雑な研究課題から生まれた実践的なツールであることを示しています。フィッシャーの研究は、統計的分析と適切な実験計画法の間に不可分な関係があることを明らかにし、その後の科学的研究のあり方に深遠な影響を与えました 。

本記事で詳述する二元配置分散分析（Two-Way ANOVA）は、このフィッシャーの思想を体現する統計手法であり、一つの連続的な従属変数（応答変数）に対して、二つの独立したカテゴリカル変数（要因）が及ぼす影響を評価するために用いられます。この分析手法の際立った特徴は、三つの異なる帰無仮説を同時に検定できる点にあります。それはすなわち、第一の要因がもたらす「主効果」、第二の要因がもたらす「主効果」、そして二つの要因が互いに影響を及ぼし合う「交互作用効果」です 。これにより、研究者は要因の独立した影響だけでなく、それらが組み合わさったときに生じる相乗的な効果をも明らかにすることが可能となります。

本記事の構成は以下の通りです。まず第2章では、二元配置分散分析の理論的基礎として、主効果と交互作用の基本概念、その数理モデル、そして分析の妥当性を担保するための前提条件について解説します。第3章では、仮説検定の枠組み、分散分析の核心である平方和の分解、そして分散分析表の構築と解釈、さらには効果量や多重比較といった統計的推論の具体的な手順を詳しく解説します。第4部では、これらの理論を実践に移すため、汎用的な数値計算ソフトウェアであるMATLABを用いて、「植物の成長」という具体的なデータセットを用いた総合的な分析事例をステップ・バイ・ステップで示します。本記事を通じて、読者が二元配置分散分析の理論的背景を深く理解し、それを自らの研究において正確に適用、解釈、そして報告するための知識を習得することを目指します。

2. 二元配置分散分析の理論的基礎

2.1. 主効果と交互作用

要因が2つ以上存在し、それらが結果に影響を与えていると考えられる時に二元分散分析を利用します。例えば、肥料の種類と、日射量が植物の草丈に影響を与えていると考えられるときに、本当にそれら要因が草丈に影響を与えているかを知るために利用します。

この二元配置分散分析を理解する上で、その核心となる以下の3つの効果を正確に把握することが不可欠です。

• 主効果（Main Effect）: 1つ目の要因が与える影響
• 主効果（Main Effect）: 2つ目の要因が与える影響
• 交互作用（Interaction Effect）: 2つの要因が組み合わさることで生じる相乗効果や相殺効果

以下の図のように交互作用を考慮するモデルでは、1つのデータは「全体平均」、「要因Aの主効果」、「要因Bの主効果」、「交互作用」、そして「誤差（説明できない個々のばらつき）」の各要素で構成されていると考えます。

2.1.1. 主効果

主効果とは、もう一方の要因の影響を平均化（あるいは無視）したときに、一つの要因が従属変数に与える独立した影響を指します 。例えば、植物の成長に関する実験で、「肥料の種類」と「日照量」という二つの要因を考える場合、「肥料の種類」の主効果とは、日照量の各水準（例えば「高」「中」「低」）における成長量の平均値をとり、肥料の種類ごと（例えば「肥料A」「肥料B」「対照群」）の成長量を比較することに相当します。これにより、「全体として、肥料の種類は植物の成長に影響を与えるか」という問いに答えることができます。同様に、「日照量」の主効果は、肥料の種類の違いを平均化し、日照量の水準間で成長量に差があるかを評価するものです。

2.1.2. 交互作用

交互作用は二元配置分散分析において最も重要かつ洞察に富む概念であり、一方の要因が従属変数に与える影響のパターンが、もう一方の要因の水準によって変化する状況を指します 。言い換えれば、二つの要因の効果が単純な足し算（加法的）で説明できず、それらが組み合わさることで特有の効果が生じる場合、交互作用が存在すると言えます 。先の例で言えば、「肥料Aは日照量が多い条件下で最も効果的だが、日照量が少ない条件下では肥料Bの方が効果的である」といった状況がこれに該当します。もし交互作用が存在しない（すなわち効果が加法的な）場合、肥料Aの効果は、日照量が多いか少ないかに関わらず一定であるはずです。

以下の図は、各観測値に含まれる変動要素を、交互作用、ならびに各主効果に段階的に分解していくプロセスを示したものです。このプロセスにより、観測データ間の差異が、交互作用および主効果の各要素によって構成されていることが視覚的に示唆されます。

2.1.3. 交互作用の優先的解釈

分析結果を解釈する際の鉄則は、常に交互作用効果から検討を開始することです。もし交互作用が統計的に有意であると判断された場合、主効果を単独で解釈することは極めて誤解を招きやすいです。なぜなら、有意な交互作用は、要因Aの効果が要因Bの水準に依存していることを意味するため、「要因Aの全体的な効果」という主効果の概念そのものが、状況を単純化しすぎており、実態を正確に反映していないからです。この場合、解釈の中心は「要因Bのどの水準において、要因Aの各水準間にどのような差があるか」という、より具体的な比較（単純主効果の検定やセルごとの平均値の比較）に移るべきです。逆に、交互作用が有意でない場合は、二つの要因が独立して（加法的に）従属変数に影響していると解釈でき、それぞれの主効果を独立に評価することが正当化されます。

2.1.4. 交互作用の視覚的理解

交互作用の有無は、交互作用プロット（Interaction Plot）を作成することで直感的に把握することができます 。このプロットは、横軸に一方の要因の水準を、縦軸に従属変数の平均値をとり、もう一方の要因の水準ごとに線で結んで描画されます。もし、プロットされた線が互いにほぼ平行であれば、それは交互作用が存在しない（あるいは非常に小さい）ことを示唆します。この場合、一方の要因の水準が変化したときの影響（線の傾き）が、もう一方の要因の水準によらず一定であることを意味するからです 。対照的に、線が平行でなく、交差したり、傾きが大きく異なったりする場合には、交互作用の存在が強く示唆されます 。

2.2. 数理モデル

二元配置分散分析の統計的挙動は、その根底にある数理モデルによって厳密に記述されます。このモデルを理解することは、分析の前提条件や結果の解釈を深く把握するために不可欠です。

2.2.1. 分散分析と一般線形モデル

まず重要な視点として、分散分析は一般線形モデル（General Linear Model, GLM）の特殊なケースとして定式化できるという点が挙げられます 。一般線形モデルは、通常、行列を用いて次のように表現されます。

$$
\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}
$$

ここで、$\mathbf{y}$は観測された従属変数のベクトル、$\mathbf{X}$は計画行列（Design Matrix）、$\boldsymbol{\beta}$は推定されるパラメータ（係数）のベクトル、そして$\boldsymbol{\epsilon}$は誤差項のベクトルです。分散分析の文脈では、要因（カテゴリカル変数）の各水準が、計画行列$\mathbf{X}$内で「ダミー変数」または「指示変数」と呼ばれる0と1（あるいは-1, 0, 1）からなる列ベクトルとして表現されます 。この統一的な視点により、分散分析は回帰分析と理論的に結びつき、より柔軟で拡張性の高い枠組みの中で理解することが可能となります 。

具体的な数値例による説明

この関係を理解するために、要因が1つである場合の具体的な例を考えてみましょう。ある学習方法の効果を測定するために、3つの異なるグループ（A, B, C）にそれぞれ2人ずつ被験者を割り当て、テストを実施したとします。結果は以下のようになりました。

• グループA (統制群/対照群): 5点, 7点 - 実験的介入を行わない基準となるグループ
• グループB (手法B): 8点, 10点 - 学習手法Bを適用したグループ
• グループC (手法C): 4点, 6点 - 学習手法Cを適用したグループ

このデータを一般線形モデル$y=X\beta+\epsilon$ に当てはめてみます。

1. 従属変数のベクトル y

まず、観測された全てのテスト得点をベクトル y として表現します。

\mathbf{y} =
\begin{pmatrix}
5 \\
7 \\
8 \\
10 \\
4 \\
6
\end{pmatrix}

1. パラメータのベクトル β と計画行列 X

次に、パラメータ β と計画行列 X を定義します。ここでは、グループAを「基準（リファレンス）」とするダミーコーディングを用います。

パラメータ β は、以下の3つの要素で構成されます。

• β0: 切片。基準であるグループAの平均値。
• β1: グループBの平均値とグループAの平均値の差。
• β2: グループCの平均値とグループAの平均値の差。

\boldsymbol{\beta} =
\begin{pmatrix}
\beta_0 \\
\beta_1 \\
\beta_2
\end{pmatrix}

計画行列 X は、各観測値がどのパラメータに対応するかを「0」と「1」で示します。

• 1列目: 全て1であり、切片 β0（全体の基準）に対応します。
• 2列目: 観測値がグループBのものであれば1、そうでなければ0となるダミー変数です。これは β1 に対応します。
• 3列目: 観測値がグループCのものであれば1、そうでなければ0となるダミー変数です。これは β2 に対応します。

したがって、計画行列 X は以下のようになります。

\mathbf{X} =\begin{pmatrix}1 & 0 & 0 \\1 & 0 & 0 \\1 & 1 & 0 \\1 & 1 & 0 \\1 & 0 & 1 \\1 & 0 & 1\end{pmatrix}

1. モデルの結合

これらを y=Xβ+ϵ の式にまとめると、以下のようになります。

\begin{pmatrix}5 \\7 \\8 \\10 \\4 \\6\end{pmatrix}=\begin{pmatrix}1 & 0 & 0 \\1 & 0 & 0 \\1 & 1 & 0 \\1 & 1 & 0 \\1 & 0 & 1 \\1 & 0 & 1\end{pmatrix}\begin{pmatrix}\beta_0 \\\beta_1 \\\beta_2\end{pmatrix}+\begin{pmatrix}\epsilon_1 \\\epsilon_2 \\\epsilon_3 \\\epsilon_4 \\\epsilon_5 \\\epsilon_6\end{pmatrix}

この行列計算を展開すると、各グループの期待値（誤差を除いた予測値）がどのようにパラメータ $\beta$で表現されるかが分かります。

• グループAの観測値: $1 \cdot \beta_0 + 0 \cdot \beta_1 + 0 \cdot \beta_2 = \beta_0$
• グループBの観測値: $1 \cdot \beta_0 + 1 \cdot \beta_1 + 0 \cdot \beta_2 = \beta_0 + \beta_1$
• グループCの観測値: $1 \cdot \beta_0 + 0 \cdot \beta_1 + 1 \cdot \beta_2 = \beta_0 + \beta_2$

この計算から、パラメータは以下を意味していることがわかります。

• $\beta_0$はグループAの平均点そのものです。（データ平均: $(5+7)/2=6$）
• $\beta_0+\beta_1$ はグループBの平均点です。（データ平均: $(8+10)/2=9$）
• $\beta_0+\beta_2$ はグループCの平均点です。（データ平均: $(4+6)/2=$5）

最小二乗法などを用いてこのデータからパラメータを推定すると、$\hat{\beta}_0=6$, $\hat{\beta}_1=3$, $\hat{\beta}_2=-1$ となります。

分散分析における「グループ間に差があるか」という問いは、一般線形モデルの文脈では「$\beta_1$ または $\beta_2$ が0と有意に異なるか（$\beta_1$=0 かつ $\beta_2$=0 という帰無仮説が棄却されるか）」という問いに変換されます。このように、計画行列 $X$の作り方を工夫することで、分散分析は回帰分析と同じ枠組みで扱うことができるのです。

2.2.2. 加法モデル（主効果モデル）

もし、二つの要因間に交互作用が存在しないと理論的に仮定できる場合、あるいはデータから交互作用が有意でないと判断された場合には、交互作用項$(\alpha\beta)_{ij}$をモデルから除外した、より単純な加法モデル（Additive Model）を先ほどと書き表し方を変えることで以下のように表すことができます 。

$$
Y_{ijk} = \mu + \alpha_i + \beta_j + \epsilon_{ijk}
$$

• $Y_{ijk}$: 要因Aの第$i$水準（$i = 1, 2, \dots, a$）、要因Bの第$j$水準（$j = 1, 2, \dots, b$）における、$k$番目の観測値（$k = 1, 2, \dots, n$）。$a$は要因Aの水準数、$b$は要因Bの水準数、$n$は各セル（水準の組み合わせ）内の繰り返し数（サンプルサイズ）です。
• $\mu$: 全体平均（Grand Mean）。すべての観測値の母平均を表します。
• $\alpha_i$: 要因Aの第$i$水準における主効果。全体平均$\mu$からのズレを表します。通常、$\sum_{i=1}^{a} \alpha_i = 0$という制約が課されます 。
• $\beta_j$: 要因Bの第$j$水準における主効果。全体平均$\mu$からのズレを表します。同様に、$\sum_{j=1}^{b} \beta_j = 0$という制約が課されます 。
• $\epsilon_{ijk}$: 実験誤差。個々の観測値に影響を与えるランダムな変動を表します。この誤差項は、互いに独立であり、平均が0、分散が$\sigma^2$の正規分布$N(0, \sigma^2)$に従うと仮定されます。これが分散分析の重要な前提条件となります 。

このモデルは、要因Aの効果（$\alpha_i$）が要因Bのどの水準においても一定であり、同様に要因Bの効果（$\beta_j$）が要因Aのどの水準においても一定であることを仮定しています。つまり、二つの要因の効果は単純に足し合わされるだけで、組み合わせによる特別な効果は生じません。このモデルの母平均は$\mu_{ij} = \mu + \alpha_i + \beta_j$となり、交互作用プロットでは線が完全に平行になる状況に対応します 。

2.2.3. 交互作用モデル（完全モデル）

二つの要因（要因A、要因B）とその交互作用を含む、完全な二元配置分散分析のモデルは以下の数式で表されます 。

Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk}

• $(\alpha\beta)_{ij}: 要因Aの第$i$水準と要因Bの第$j$水準の組み合わせによって生じる交互作用効果。主効果だけでは説明できない、特定の組み合わせに固有の効果を表します。これも同様に、各行、各列について和が0になる以下のような制約が課されます。

\sum_{i=1}^{a} (\alpha\beta)_{ij} = 0 \quad \text{for all } j, \quad \text{and} \quad \sum_{j=1}^{b} (\alpha\beta)_{ij} = 0 \quad \text{for all } i

このモデルは、各セル（要因AとBの特定の水準の組み合わせ）が、それぞれ固有の母平均$\mu_{ij} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij}$を持つことを許容しており、最も柔軟なモデルです。

2.3. 前提条件と頑健性

二元配置分散分析から得られるF検定の結果が統計的に妥当であるためには、その数理モデルが依拠するいくつかの重要な前提条件が満たされている必要があります。これらの条件は、モデルの誤差項$\epsilon_{ijk}$の性質に関するものです 。

1. 観測の独立性 (Independence of Observations)

   各観測値は、他のどの観測値からも統計的に独立している必要があります。これは、ある被験者の測定値が他の被験者の測定値に影響を与えない、あるいはある測定時点の結果が次の測定時点の結果に影響しないことを意味します。この前提は、統計的検定によって確認するものではなく、主に実験計画の段階で保証されるべき最も重要な条件です。例えば、被験者を各グループに無作為に割り当てること（ランダム化）は、この独立性を確保するための基本的な手続きです。この前提が破られると（例えば、同じ被験者から繰り返し測定を行う場合）、第一種の過誤（帰無仮説が真であるにもかかわらず棄却してしまう誤り）の確率が著しく増大する可能性があります。

2. 分散の均一性 (Homogeneity of Variances / Homoscedasticity)

   比較するすべてのグループ（すなわち、要因Aと要因Bのすべての水準の組み合わせからなる各セル）において、誤差の分散が等しい必要があります。つまり、$\text{Var}(\epsilon_{ijk}) = \sigma^2$がすべての$i, j$に対して一定でなければなりません。下図の左側のように分散が等しいと考えられる場合にはこの前提は満たされますが、一方で左図のように分散の均一性が等しいと考えられない場合にはこの前提は満たされません。この前提は、通常、レーベンの検定（Levene's test）によって検定されます 。この前提が満たされない場合、特に各セルのサンプルサイズが不均一（アンバランス）であると、F検定の結果が歪められ、信頼性が低下します。
   

3. 残差の正規性 (Normality of Residuals):

   各グループにおける誤差項（あるいは、モデルをデータにあてはめた後の残差）が正規分布に従う必要があります。これは、各グループのデータそのものが正規分布に従うこととほぼ同義ですが、より正確にはモデルから予測される値と実測値の差である残差の分布を評価します。

   この前提は、シャピロ-ウィルク検定（Shapiro-Wilk test）のような統計的検定や、下図のようなQ-Qプロット（Quantile-Quantile plot）と呼ばれるグラフ的手法を用いて評価されます。

   

3. 仮説検定と統計的推論

二元配置分散分析の目的は、データに基づいて統計的な推論を行い、要因の効果に関する結論を導き出すことです。このプロセスは、仮説の設定、データの変動を分解する平方和の計算、そしてその結果をまとめた分散分析表の解釈という一連の手順を経て行われます。

2.1. 仮説の設定

二元配置分散分析では、前述の通り、三組の帰無仮説（$H_0$）と対立仮説（$H_1$）を同時に検定します。

1. 交互作用効果に関する仮説

   • 帰無仮説 $H_0$: 要因Aと要因Bの間に交互作用は存在しません。

   
   H_0: (\alpha\beta)_{ij} = 0 \quad \text{for all } i, j
   
   

   • 対立仮説 $H_1$: 要因Aと要因Bの間に交互作用が存在します。

     
     H_1: \text{At least one } (\alpha\beta)_{ij} \neq 0
     
     

2. 要因Aの主効果に関する仮説

   • 帰無仮説 $H_0$: 要因Aの主効果は存在しません。すなわち、要因Aのすべての水準における母平均は等しいです。

   H_0: \mu_{1.} = \mu_{2.} = \dots = \mu_{a.} \quad (\text{or equivalently, } \alpha_i = 0 \text{ for all } i)
   
   

   ここで、$\mu_{i.}$は要因Aの第$i$水準における母平均を表します。

   • 対立仮説 $H_1$: 要因Aの主効果が存在します。すなわち、要因Aの少なくとも一つの水準で母平均が他と異なります。

   H_1: \text{At least one } \mu_{i.} \text{ is different} \quad (\text{or equivalently, at least one } \alpha_i \neq 0)
   
   

3. 要因Bの主効果に関する仮説

   • 帰無仮説 $H_0$: 要因Bの主効果は存在しません。すなわち、要因Bのすべての水準における母平均は等しいです。

   ここで、$\mu_{.j}$は要因Bの第$j$水準における母平均を表します。

   
   H_0: \mu_{.1} = \mu_{.2} = \dots = \mu_{.b} \quad (\text{or equivalently, } \beta_j = 0 \text{ for all } j)
   
   

   • 対立仮説 $H_1$: 要因Bの主効果が存在します。すなわち、要因Bの少なくとも一つの水準で母平均が他と異なります。

   
   H_1: \text{At least one } \mu_{.j} \text{ is different} \quad (\text{or equivalently, at least one } \beta_j \neq 0)
   
   

これらの仮説を検定するために、分散分析ではデータの総変動を各要因に起因する変動へと分解していきます。

2.2. 平方和の分解

分散分析の数学的な核心は、データの総変動（Total Sum of Squares, SST）を、各要因の主効果、交互作用、そして誤差に起因する変動へと分割することにあります 。この「平方和の分解」により、各変動要因が全体のばらつきに対してどれだけ寄与しているかを定量的に評価できます。

まず、以下の記法を導入します。

• $Y_{ijk}$: 要因Aの第$i$水準、要因Bの第$j$水準における$k$番目の観測値
• $\bar{Y}_{ij.}$: セル$(i,j)$の平均値
• $\bar{Y}_{i..}$: 要因Aの第$i$水準の平均値（行平均）
• $\bar{Y}_{.j.}$: 要因Bの第$j$水準の平均値（列平均）
• $\bar{Y}_{...}$: 全体平均（Grand Mean）
• $a$: 要因Aの水準数
• $b$: 要因Bの水準数
• $n$: 各セル内の繰り返し数（バランスデザインの場合）

データの総変動は、各観測値と全体平均との差の平方和である総平方和（SST）で表されます。

SST = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{ijk} - \bar{Y}_{...})^2

この総平方和は、代数的な操作により、以下の4つの要素に直交分解されます（詳細な証明は付録A.1を参照）。

SST = SSA + SSB + SSAB + SSE

各構成要素は以下の通りです。

• 要因Aの平方和（SSA）: 要因Aの主効果による変動。各行平均と全体平均の差の平方和で計算されます。

SSA = bn \sum_{i=1}^{a} (\bar{Y}{i..} - \bar{Y}{...})^2

• 要因Bの平方和（SSB）: 要因Bの主効果による変動。各列平均と全体平均の差の平方和で計算されます。

SSB = an \sum_{j=1}^{b} (\bar{Y}{.j.} - \bar{Y}{...})^2

• 交互作用の平方和（SSAB）: 要因AとBの交互作用による変動。セル平均から主効果の影響を取り除いた残りの変動を表します。

SSAB = n \sum_{i=1}^{a} \sum_{j=1}^{b} (\bar{Y}{ij.} - \bar{Y}{i..} - \bar{Y}{.j.} + \bar{Y}{...})^2

• 誤差平方和（SSE）: いずれの要因によっても説明されない、グループ内のランダムな変動（残差変動）。各観測値と所属するセル平均の差の平方和で計算されます。

SSE = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{ijk} - \bar{Y}_{ij.})^2

この分解により、例えばSSAがSSTに対して大きな割合を占めていれば、要因Aがデータのばらつきに大きく寄与していることが示唆されます。

2.3. 分散分析表の構築と解釈

平方和の計算結果は、標準的な分散分析表（ANOVA Table）にまとめられます。この表は、各変動要因の評価を系統的に行い、仮説検定の結果を明瞭に提示するためのものです。

表：標準的な二元配置分散分析表の構造

変動要因	平方和 (SS)	自由度 (df)	平均平方 (MS)	F値 (F-statistic)	p値 (p-value)
要因A	$SSA$	$df_A = a-1$	$MSA = \frac{SSA}{df_A}$	$F_A = \frac{MSA}{MSE}$	$P(F > F_A)$
要因B	$SSB$	$df_B = b-1$	$MSB = \frac{SSB}{df_B}$	$F_B = \frac{MSB}{MSE}$	$P(F > F_B)$
交互作用	$SSAB$	$df_{AB} = (a-1)(b-1)$	$MSAB = \frac{SSAB}{df_{AB}}$	$F_{AB} = \frac{MSAB}{MSE}$	$P(F > F_{AB})$
誤差	$SSE$	$df_E = ab(n-1)$	$MSE = \frac{SSE}{df_E}$
合計	$SST$	$df_T = abn-1$

分散分析表の各列の意味は以下の通りです。

• 自由度（Degrees of Freedom, df）: 各平方和を計算するために用いられる、独立な情報の数。自由度は、水準数やサンプルサイズから決定されます。
• 平均平方（Mean Square, MS）: 平方和を対応する自由度で割った値（$MS = SS/df$）。これは、各変動要因に起因する分散の不偏推定量と解釈されます。特に、$MSE$は、要因の効果がない場合のデータのばらつき（誤差分散$\sigma^2$）の推定量となります。
• F値（F-statistic）: 仮説検定のための検定統計量です。各要因（主効果または交互作用）の平均平方を誤差の平均平方で割った比率（$F = MS_{effect} / MSE$）です 。このF値は、「要因による変動が、偶然による変動（誤差）と比べてどれだけ大きいか」を示す指標となります。帰無仮説が真である場合、このF値はF分布に従います。
• p値（p-value）: 計算されたF値、あるいはそれよりも極端な値が、帰無仮説の下で偶然生じる確率です。このp値が、あらかじめ設定した有意水準（通常は$\alpha = 0.05$）よりも小さい場合、帰無仮説を棄却し、対立仮説を採択します。すなわち、その要因の効果は「統計的に有意である」と結論付けられます。

2.4. 効果量と多重比較

分散分析表から得られるp値は、効果の有無を判断するための二元的な情報（有意か、非有意か）しか提供しません。しかし、研究においては、その効果が「どの程度の大きさ」を持つのか、すなわち実践的な重要性を評価することが不可欠です。このために効果量と多重比較が用いられます。

効果量

効果量（Effect Size）は、要因が従属変数の変動をどの程度説明しているかを示す標準化された指標です 。p値がサンプルサイズの影響を大きく受けるのに対し、効果量はサンプルサイズに依存しない効果の大きさそのものを表すため、異なる研究間での結果の比較を可能にします。

分散分析において最も一般的に用いられる効果量の一つがパーシャル・イータ二乗（Partial Eta-Squared, $\eta_p^2$）です。これは、他の要因や交互作用の影響を除いた後で、特定の要因が説明する変動の割合を示します。計算式は以下の通りです。

$$
\eta_p^2 = \frac{SS_{\text{effect}}}{SS_{\text{effect}} + SS_{\text{error}}}
$$

ここで、$SS_{\text{effect}}$は関心のある要因の平方和（$SSA$, $SSB$, または$SSAB$）、$SS_{\text{error}}$は誤差平方和（$SSE$）です。$\eta_p^2$の値は0から1の範囲をとり、例えば$\eta_p^2 = 0.14$であれば、その要因が従属変数の変動の14%を説明していると解釈できます。慣例的な解釈の目安として、0.01が「小さい効果」、0.06が「中程度の効果」、0.14が「大きい効果」とされることがありますが、これらの基準は分野や文脈によって異なるため、絶対的なものではありません。

多重比較

分散分析の結果、ある要因の主効果が有意であると結論付けられた場合、その要因に3つ以上の水準が存在するならば、次に「どの水準とどの水準の間に具体的な差があるのか」を特定する必要があります。F検定はあくまで「少なくとも一つの水準の平均が他と異なる」ことしか示してくれないからです 。この目的のために行われるのが多重比較（Multiple Comparisons）または事後検定（Post-Hoc Tests）です。

多重比較には多くの手法が存在しますが、すべての水準のペアを比較する場合に広く用いられるのがテューキーのHSD（Honestly Significant Difference）検定です。この検定は、複数のペアワイズ比較を同時に行っても、全体の第一種の過誤の確率（ファミリーワイズエラー率）を特定の水準（例えば0.05）に制御する特徴を持ちます。これにより、検定を繰り返すことで偽陽性が増加する問題を回避できます。有意な交互作用がある場合は、主効果の多重比較ではなく、全セルの平均値間での多重比較を行うことがより適切です。

4. MATLABによる具体例：植物の成長データ分析

この章では、具体的な研究シナリオに基づき、二元配置分散分析の理論をMATLABを用いて実践する包括的な分析事例を示します。データ準備から分析、結果の解釈、そして結論の報告まで、一連のプロセスを解説します。

4.1. 研究シナリオとデータセット

研究目的
ある植物学者が、2種類の肥料（肥料A、肥料B）と日照条件が植物の成長（草丈）に与える影響を評価したいと考えています。特に、肥料の効果が日照条件によって変化するか（交互作用）を検証することが目的です。

• 要因A: 肥料の種類 (Fertilizer)
  • 3水準: 'Control'（対照群：肥料なし）、'Fertilizer_A'、'Fertilizer_B'
• 要因B: 日照量 (Sunlight)
  • 4水準: 'None'（日照なし）、'Low'（低い）、'Medium'（中程度）、'High'（高い）

これは3✖️4の要因計画であり、各組み合わせ（セル）ごとに5つの植物を栽培し、合計60個の観測値が得られたバランスの取れたデザインです。

データセット
分析に使用するデータは、MATLAB上で生成し、anovan関数が要求するロングフォーマットのテーブル形式に準備します。

clear; clc; close all;

%% 1. データローディングと準備
% --- データ生成 (再現性のための模擬データ) ---
rng(123); % 乱数シードを固定
Height = [15.5+randn(5,1)*0.5; 20.1+randn(5,1)*0.8; 22.4+randn(5,1)*0.7; 23.1+randn(5,1)*0.9;... % Control
          16.0+randn(5,1)*0.6; 25.5+randn(5,1)*1.0; 29.8+randn(5,1)*0.8; 35.2+randn(5,1)*1.1;... % Fertilizer_A
          15.8+randn(5,1)*0.7; 22.3+randn(5,1)*0.9; 32.5+randn(5,1)*1.2; 32.8+randn(5,1)*1.0];   % Fertilizer_B

FertilizerLevels = {'Control', 'Fertilizer_A', 'Fertilizer_B'};
SunlightLevels = {'None', 'Low', 'Medium', 'High'};

Fertilizer = repelem(categorical(FertilizerLevels'), 20, 1);
Sunlight = repmat(repelem(categorical(SunlightLevels'), 5, 1), 3, 1);

% データテーブルの作成
plant_data = table(Height, Fertilizer, Sunlight);
disp('--- データの一部を表示: ---');
disp(head(plant_data));     

4.2. 前提条件の検定

分析の妥当性を確認するため、分散の均一性と残差の正規性を検定します。

分散の均一性（レーベンの検定）vartestn関数を使用して、12個の全セル（3*4）間で分散が等しいかを確認します。

%% 2. 前提条件の検定
% --- 分散の均一性検定 (レーベンの検定) ---
disp('--- 分散の均一性検定 (レーベンの検定) ---');

% 2つの要因を文字列として結合し、単一のグループ化変数を作成します。
% 例: "Control-None", "Fertilizer_A-Low" など
% この方法は、grp2idxの複雑な入力要件を回避する、より直接的な解決策です。
combined_groups = string(plant_data.Fertilizer) + "-" + string(plant_data.Sunlight);
[p_levene, ~] = vartestn(plant_data.Height, combined_groups,...
    'TestType', 'LeveneAbsolute', 'Display', 'off');

fprintf('レーベンの検定のp値: %.4f\n', p_levene);
if p_levene >= 0.05
    disp('p >= 0.05 のため、分散の均一性の仮定は満たされていると判断します。');
else
    disp('p < 0.05 のため、分散の均一性の仮定は満たされていません。');
end

レーベンの検定のp値: 0.1850であり、p >= 0.05 のため、分散の均一性の仮定は満たされていると判断します。

残差の正規性
まずanovanを一度実行して残差を抽出し、その残差に対して正規性を評価します。視覚的な確認のためにQ-Qプロットを、統計的な検定のためにリリフォース検定を用います。

% --- 残差の正規性検定 ---
[~, ~, stats_temp] = anovan(plant_data.Height, {plant_data.Fertilizer, plant_data.Sunlight},...
    'model', 'interaction', 'display', 'off');
residuals = stats_temp.resid;

disp('--- 残差の正規性検定 ---');
% Q-Qプロットによる視覚的確認
figure;
qqplot(residuals);
title('残差のQ-Qプロット');
grid on;

% リリフォース検定による統計的確認
[h_lillie, p_lillie] = lillietest(residuals);
fprintf('リリフォース検定のp値: %.4f\n', p_lillie);
if h_lillie == 0
    disp('p >= 0.05 のため、残差の正規性の仮定は棄却されません。');
else
    disp('p < 0.05 のため、残差は正規分布に従わない可能性があります。');
end

リリフォース検定のp値: 0.1612であり、p >= 0.05 のため、残差の正規性の仮定は棄却されません。

4.3. 二元配置分散分析の実行と結果の解釈

前提条件が満たされていることを確認した上で、anovan関数を用いて主効果と交互作用を含むモデルを分析します。

disp('--- 二元配置分散分析の実行 ---');
[~, tbl_anova, stats_anova] = anovan(plant_data.Height, {plant_data.Fertilizer, plant_data.Sunlight},...
    'model', 'interaction',...
    'varnames', {'Fertilizer', 'Sunlight'},...
    'display', 'on'); % 'on'に設定して分散分析表を表示

分散分析表の解釈
出力された分散分析表から、各効果のp値を確認します。

• 交互作用 (Fertilizer*Sunlight): p値が有意水準（例: 0.05）を下回る場合、肥料の効果が日照量の水準によって異なる、あるいはその逆の関係があることを強く示唆します。
• 主効果 (Fertilizer, Sunlight): 交互作用が有意な場合、主効果を単独で解釈するには注意が必要です。解釈は交互作用の文脈で行うべきです。

効果量の計算
パーシャル・イータ二乗（eta_p2）を計算し、各効果が従属変数の変動をどの程度説明するかを評価します。

disp('--- 効果量の計算 (パーシャル・イータ二乗) ---');
SS_Fertilizer = tbl_anova{2,2};
SS_Sunlight = tbl_anova{3,2};
SS_Interaction = tbl_anova{4,2};
SS_Error = tbl_anova{5,2};

eta_p_sq_Fertilizer = SS_Fertilizer / (SS_Fertilizer + SS_Error);
eta_p_sq_Sunlight = SS_Sunlight / (SS_Sunlight + SS_Error);
eta_p_sq_Interaction = SS_Interaction / (SS_Interaction + SS_Error);

fprintf('Fertilizerのパーシャル・イータ二乗: %.3f\n', eta_p_sq_Fertilizer);
fprintf('Sunlightのパーシャル・イータ二乗: %.3f\n', eta_p_sq_Sunlight);
fprintf('Fertilizer*Sunlightのパーシャル・イータ二乗: %.3f\n', eta_p_sq_Interaction);

Fertilizerのパーシャル・イータ二乗: 0.955
Sunlightのパーシャル・イータ二乗: 0.988
Fertilizer*Sunlightのパーシャル・イータ二乗: 0.926

4.4. 結果の可視化と事後解析

交互作用プロット
有意であった交互作用をinteractionplotで視覚的に確認します。プロットされた線の非平行性（特に交差）は、交互作用の存在を示す視覚的な証拠となります。

disp('--- 交互作用プロットの作成 ---');
figure;
interactionplot(plant_data.Height, {plant_data.Fertilizer, plant_data.Sunlight},...
    'varnames', {'Fertilizer', 'Sunlight'});
title('FertilizerとSunlightの交互作用プロット');
ylabel('平均草丈 (cm)');

交互作用の事後解析（多重比較）
交互作用が有意であったため、multcompare関数を用いて、どのセルの平均値間に差があるのかを具体的に特定します。ここでは全12セル間の平均値をテューキーのHSD法で比較します。

disp('--- 多重比較検定 (テューキーのHSD法) ---');
multcompare(stats_anova, 'Dimension', [1 2], 'CType', 'hsd');

この検定により、「日照量が'High'の条件下では、'Fertilizer_A'は'Control'よりも有意に成長を促進した」といった、より詳細な知見を得ることができます。

5. 引用文献

• MathWorks. (n.d.-a). N-Way ANOVA. MATLAB Documentation. Retrieved from https://jp.mathworks.com/help/stats/n-way-anova.html
• MathWorks. (n.d.-b). anovan. MATLAB Documentation. Retrieved from https://www.mathworks.com/help/stats/anovan.html
• MathWorks. (n.d.-c). Two-Way ANOVA. MATLAB Documentation. Retrieved from https://www.mathworks.com/help/stats/two-way-anova.html
• Scheirer, J., Ray, W. S., & Hare, N. (1976). The Analysis of Ranked Data Derived from Completely Randomized Factorial Designs. Biometrics, 32(2), 429–434.
  https://doi.org/10.2307/2529511
• Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3/4), 591–611.
  https://doi.org/10.1093/biomet/52.3-4.591

6. 付録

A.1. 平方和の分解に関する代数的証明

二元配置分散分析における総平方和（SST）の分解、すなわち $SST = SSA + SSB + SSAB + SSE$が成立することの証明は、各観測値の全体平均からの偏差 $(Y_{ijk} - \bar{Y}_{...})$を段階的に分解し、その平方和を展開した際に生じる交差積項（cross-product terms）がすべて0になることを示すことで行われます [91, 92]。

まず、個々の偏差を以下のように分解します。

Y_{ijk} - \bar{Y}{...} = (\bar{Y}_{i..} - \bar{Y}_{...}) + (\bar{Y}_{.j.} - \bar{Y}_{...}) + (\bar{Y}_{ij.} - \bar{Y}_{i..} - \bar{Y}_{.j.} + \bar{Y}_{...}) + (Y_{ijk} - \bar{Y}_{ij.})

この式は、右辺の項を整理すると左辺と等しくなる恒等式です。右辺の各項はそれぞれ、要因Aの主効果、要因Bの主効果、交互作用効果、そして誤差に対応しています。

総平方和 $SST = \sum_{i}\sum_{j}\sum_{k}(Y_{ijk} - \bar{Y}{...})^2$ は、この分解式の右辺を二乗し、すべての観測値にわたって合計をとることで得られます。

SST = \sum{i,j,k} [(\bar{Y}{i..} - \bar{Y}{...}) + (\bar{Y}{.j.} - \bar{Y}{...}) + (\bar{Y}{ij.} - \bar{Y}{i..} - \bar{Y}{.j.} + \bar{Y}{...}) + (Y_{ijk} - \bar{Y}{ij.})]^2

この式を展開すると、各項の二乗和（これが$SSA, SSB, SSAB, SSE$に他なりません）と、異なる項同士の積（交差積）の和が現れます。例えば、要因Aの主効果の項と誤差項の交差積は以下のようになります。

2 \sum{i,j,k} (\bar{Y}{i..} - \bar{Y}{...})(Y_{ijk} - \bar{Y}{ij.})

この交差積項が0になることを示します。まず、$k$に関する和をとります。

2 \sum_{i,j} (\bar{Y}{i..} - \bar{Y}{...}) \sum_{k} (Y_{ijk} - \bar{Y}{ij.})

ここで、$\sum{k} (Y_{ijk} - \bar{Y}_{ij.})$ は、セル$(i,j)$内の観測値とそのセル平均との偏差の和であるため、定義により0になります。したがって、この交差積項全体も0になります。

同様の代数的計算により、他のすべての交差積項（主効果Aと主効果B、主効果と交互作用、交互作用と誤差など）の和も0になることが証明できます 。その結果、総平方和は各効果の平方和の単純な和として表現されます。

\sum_{i,j,k}(Y_{ijk} - \bar{Y}{...})^2 = \sum{i,j,k}(\bar{Y}{i..} - \bar{Y}{...})^2 + \sum_{i,j,k}(\bar{Y}{.j.} - \bar{Y}{...})^2 + \sum_{i,j,k}(\bar{Y}{ij.} - \bar{Y}{i..} - \bar{Y}{.j.} + \bar{Y}{...})^2 + \sum_{i,j,k}(Y_{ijk} - \bar{Y}_{ij.})^2

SST = SSA + SSB + SSAB + SSE

これは、データの総変動が、互いに直交する（相関のない）変動成分に完全に分解されることを数学的に示しています。

A.2. 本稿で使用したMATLABコード全文

第4部の具体例で使用したMATLABコードの完全版を以下に示します。

%% 第4章 MATLABによる総合的な分析事例：植物の成長データ分析
% このスクリプトは、二元配置分散分析の全プロセスを網羅しています。

clear; clc; close all;

%% 1. データローディングと準備
% --- データ生成 (再現性のための模擬データ) ---
rng(123); % 乱数シードを固定
Height = [15.5+randn(5,1)*0.5; 20.1+randn(5,1)*0.8; 22.4+randn(5,1)*0.7; 23.1+randn(5,1)*0.9;... % Control
          16.0+randn(5,1)*0.6; 25.5+randn(5,1)*1.0; 29.8+randn(5,1)*0.8; 35.2+randn(5,1)*1.1;... % Fertilizer_A
          15.8+randn(5,1)*0.7; 22.3+randn(5,1)*0.9; 32.5+randn(5,1)*1.2; 32.8+randn(5,1)*1.0];   % Fertilizer_B

FertilizerLevels = {'Control', 'Fertilizer_A', 'Fertilizer_B'};
SunlightLevels = {'None', 'Low', 'Medium', 'High'};

Fertilizer = repelem(categorical(FertilizerLevels'), 20, 1);
Sunlight = repmat(repelem(categorical(SunlightLevels'), 5, 1), 3, 1);

% データテーブルの作成
plant_data = table(Height, Fertilizer, Sunlight);
disp('--- データの一部を表示: ---');
disp(head(plant_data));

%% 2. 前提条件の検定
% --- 分散の均一性検定 (レーベンの検定) ---
disp('--- 分散の均一性検定 (レーベンの検定) ---');

% 2つの要因を文字列として結合し、単一のグループ化変数を作成します。
% 例: "Control-None", "Fertilizer_A-Low" など
% この方法は、grp2idxの複雑な入力要件を回避する、より直接的な解決策です。
combined_groups = string(plant_data.Fertilizer) + "-" + string(plant_data.Sunlight);
[p_levene, ~] = vartestn(plant_data.Height, combined_groups,...
    'TestType', 'LeveneAbsolute', 'Display', 'off');

fprintf('レーベンの検定のp値: %.4f\n', p_levene);
if p_levene >= 0.05
    disp('p >= 0.05 のため、分散の均一性の仮定は満たされていると判断します。');
else
    disp('p < 0.05 のため、分散の均一性の仮定は満たされていません。');
end

% --- 残差の正規性検定 ---
% 一時的にANOVAを実行して残差を取得
[~, ~, stats_temp] = anovan(plant_data.Height, {plant_data.Fertilizer, plant_data.Sunlight},...
    'model', 'interaction', 'display', 'off');
residuals = stats_temp.resid;

disp('--- 残差の正規性検定 ---');
% Q-Qプロットによる視覚的確認
figure('Name', '残差のQ-Qプロット');
qqplot(residuals);
title('残差のQ-Qプロット');
grid on;

% リリフォース検定による統計的確認
[h_lillie, p_lillie] = lillietest(residuals);
fprintf('リリフォース検定のp値: %.4f\n', p_lillie);
if h_lillie == 0
    disp('p >= 0.05 のため、残差の正規性の仮定は棄却されません。');
else
    disp('p < 0.05 のため、残差は正規分布に従わない可能性があります。');
end

%% 3. 二元配置分散分析の実行
disp('--- 二元配置分散分析の実行 ---');
[~, tbl_anova, stats_anova] = anovan(plant_data.Height, {plant_data.Fertilizer, plant_data.Sunlight},...
    'model', 'interaction',...
    'varnames', {'Fertilizer', 'Sunlight'},...
    'display', 'on'); % 'on'に設定して分散分析表を表示

%% 4. 効果量の計算 (パーシャル・イータ二乗)
disp('--- 効果量の計算 (パーシャル・イータ二乗) ---');
SS_Fertilizer = tbl_anova{2,2};
SS_Sunlight = tbl_anova{3,2};
SS_Interaction = tbl_anova{4,2};
SS_Error = tbl_anova{5,2};

eta_p_sq_Fertilizer = SS_Fertilizer / (SS_Fertilizer + SS_Error);
eta_p_sq_Sunlight = SS_Sunlight / (SS_Sunlight + SS_Error);
eta_p_sq_Interaction = SS_Interaction / (SS_Interaction + SS_Error);

fprintf('Fertilizerのパーシャル・イータ二乗: %.3f\n', eta_p_sq_Fertilizer);
fprintf('Sunlightのパーシャル・イータ二乗: %.3f\n', eta_p_sq_Sunlight);
fprintf('Fertilizer*Sunlightのパーシャル・イータ二乗: %.3f\n', eta_p_sq_Interaction);

%% 5. 結果の可視化と事後解析
% --- 交互作用プロットの作成 ---
disp('--- 交互作用プロットの作成 ---');
figure('Name', '交互作用プロット');
interactionplot(plant_data.Height, {plant_data.Fertilizer, plant_data.Sunlight},...
    'varnames', {'Fertilizer', 'Sunlight'});
title('FertilizerとSunlightの交互作用プロット');
ylabel('平均草丈 (cm)');

% --- 多重比較検定 (テューキーのHSD法) ---
disp('--- 多重比較検定 (テューキーのHSD法) ---');
% 'Dimension', [1 2] で全セルの組み合わせを比較
figure('Name', '多重比較');
multcompare(stats_anova, 'Dimension', [1 2], 'CType', 'hsd');