農場の区画ごとに異なる肥料を施して収量を比較したい。工場の生産ラインで複数の条件を変えて品質への影響を調べたい。このような場面で、限られた資源のもとで最も多くの情報を引き出すための方法論が実験計画法（Design of Experiments, DoE）です。
ロナルド・フィッシャー（R.A. Fisher）が1920年代にロサムステッド農業試験場で体系化した実験計画法は、現代の統計学と因果推論の基盤となっています。ランダム化比較試験（RCT）はその最も基本的な形であり、本記事ではRCTを超えたより洗練されたデザインを解説します。
実験計画法は以下の分野で活用されます。
- 農学: 肥料、品種、灌漑方法の効果比較
- 製造業: 品質管理、プロセス最適化（タグチメソッド）
- 医学: 臨床試験のデザイン（クロスオーバー試験、要因計画）
- テクノロジー: A/Bテストの拡張（多変量テスト）
本記事の内容
- フィッシャーの3原則（反復、ランダム化、局所管理）
- 完全無作為化デザイン（CRD）
- 乱塊法（RCBD）
- ラテン方格法
- 要因計画法
- Pythonによる各デザインのシミュレーションと分散分析
前提知識
- ランダム化比較試験（RCT） — RCTの理論と設計
- 平均処置効果（ATE） — ATEの推定
フィッシャーの3原則
フィッシャーが体系化した実験計画法の核心は、以下の3つの原則です。これらは「良い実験」を設計するための普遍的なガイドラインであり、分野を問わず適用されます。
反復（Replication）
同じ処置を複数の実験単位に適用することです。反復がなければ、処置効果と偶然の変動を区別できません。
例えば、肥料Aと肥料Bの効果を比較するとき、各肥料をそれぞれ1つの区画にしか施さなければ、収量の差が「肥料の効果」なのか「区画固有の土壌条件の違い」なのか判断できません。各肥料を5つの区画に施せば、処置ごとの平均収量を計算し、区画間のばらつきも推定できます。
反復の数 $n$ を増やすと、平均値の標準誤差は $\sigma / \sqrt{n}$ に比例して減少するため、推定の精度が向上し、検定の検出力が高まります。必要な反復数は、事前に検出力分析（power analysis）で決定することが推奨されます。
ここで「反復」と「擬似反復」（pseudoreplication）の違いに注意しましょう。真の反復は、処置を独立した実験単位に適用することです。もし1つの区画に肥料を施して5箇所で測定しても、それは5回の反復ではなく1回の反復（5回の測定）です。擬似反復を真の反復と取り違えると、自由度が過大に見積もられ、偽陽性率が跳ね上がります。Hurlbert（1984）はこの問題を系統的に指摘し、生態学の実験における擬似反復の蔓延を批判しました。
ランダム化（Randomization）
処置の割り当てをランダムに行うことです。RCTで詳しく解説したように、ランダム化は未観測の交絡変数の影響を期待値的にゼロにし、因果推定のバイアスを除去します。
ランダム化が重要な理由を農業試験で考えましょう。もし研究者が「良い区画に新しい肥料、悪い区画に古い肥料」を割り当てれば、新しい肥料の効果が過大評価されます。ランダム化は、このような系統的な割り当てバイアスを防ぎます。
フィッシャーの深い洞察は、ランダム化が推測統計の基盤を提供するという点にあります。ランダム化により、帰無仮説のもとでの検定統計量の分布が決定され、p値の計算が正当化されます。これは後述する並べ替え検定の考え方に直結しています。
局所管理（Local Control / Blocking）
実験単位間の変動を制御するために、類似した実験単位をまとめてブロック（block）とする方法です。ブロック内ではランダム化を行いますが、ブロック間の変動は統計的に除去されます。
ブロッキングの直感は「比較は同じ条件のもとで行う」ということです。農業試験で農場に日当たりの良い区画と悪い区画がある場合、日当たりの同じ区画をまとめてブロックとし、各ブロック内ですべての処置を施します。こうすることで、日当たりの差（ブロック間変動）は処置効果の推定に影響しなくなります。
数学的には、ブロッキングは全体の変動を「ブロック間変動」と「ブロック内変動」に分解し、ブロック間変動を分析から除去することに対応します。これにより、誤差項の分散が小さくなり、処置効果の推定精度と検定の検出力が向上します。
ブロッキングの効果を定量的に評価するために、相対効率（relative efficiency, RE）という指標がよく使われます。
$$ RE = \frac{MS_{\text{error(CRD)}}}{MS_{\text{error(RCBD)}}} $$
$RE > 1$ であれば、ブロッキングにより誤差が減少し、効率が向上したことを意味します。
フィッシャーのランダム化推論
ランダム化検定の基本的考え方
フィッシャーが提唱したランダム化の原理は、推測統計の基盤そのものを提供します。これは後世の漸近理論に基づくアプローチとは異なる、有限標本で厳密に成り立つ推論方法です。
ランダム化検定（permutation test / randomization test）の基本的な考え方は以下の通りです。
-
シャープ帰無仮説: 帰無仮説を「すべての実験単位に対して処置効果がゼロ」（$H_0: Y_i(1) = Y_i(0)$ for all $i$）と定式化します。これはフィッシャーのシャープ帰無仮説（Fisher’s sharp null hypothesis）と呼ばれます。
-
参照分布の構成: シャープ帰無仮説のもとでは、各個体のアウトカムは処置の割り当てに関係なく一定です。したがって、処置の割り当てを「ランダムにやり直す」と、帰無仮説のもとでの検定統計量の分布が得られます。
-
p値の計算: 実際に観測された検定統計量が、参照分布の中でどの程度極端かを計算し、p値とします。
具体的に定式化しましょう。$n$ 個の実験単位を $k$ 群にランダムに割り当てるCRDを考えます。帰無仮説のもとで、処置の割り当て $\bm{W} = (W_1, \ldots, W_n)$ のすべての並べ替えに対する検定統計量の値を計算します。全部で $\binom{n}{n_1, n_2, \ldots, n_k}$ 通りの並べ替えがあります。
$$ \begin{equation} p\text{-value} = \frac{\#\{\bm{W}^* : T(\bm{W}^*) \geq T(\bm{W}_{\text{obs}})\}}{\text{Total permutations}} \end{equation} $$
この推論方法はフィッシャーのランダム化推論と呼ばれ、正規性やその他の分布仮定を一切必要としません。ランダム化という実験操作そのものが推論を正当化するのです。
ANOVAとの関係
ランダム化検定はANOVAのF検定と密接に関係しています。正規性の仮定が成り立つとき、ランダム化検定とF検定は漸近的に同じ結果を与えます。しかし、正規性が疑わしい場合や標本が小さい場合には、ランダム化検定の方が分布仮定に依存しない分、頑健です。
実務的には、ランダム化検定は並べ替えの数が膨大になりうるため、すべてを列挙する代わりにモンテカルロ近似（ランダムに多数の並べ替えを抽出）が用いられます。
フィッシャーのランダム化推論は、実験計画法の理論的基盤であると同時に、現代のノンパラメトリック検定やブートストラップ法にもつながる重要な概念です。
この理論的基盤を踏まえて、具体的な実験デザインを順に見ていきましょう。
完全無作為化デザイン（CRD）
CRD（Completely Randomized Design）は最も単純な実験デザインです。$n$ 個の実験単位を $k$ 個の処置群にランダムに割り当てます。ブロッキングは行わず、すべての実験単位を均質と見なします。
CRDは実験単位間の異質性が小さい場合（例: 実験室で管理された条件、同一ロットの材料）に適しています。逆に、実験単位間に大きな変動がある場合は、後述するRCBDの方が効率的です。
統計モデル
$$ \begin{equation} Y_{ij} = \mu + \tau_i + \epsilon_{ij}, \quad i = 1, \ldots, k, \quad j = 1, \ldots, n_i \end{equation} $$
ここで $\mu$ は全体平均、$\tau_i$ は第 $i$ 処置の効果、$\epsilon_{ij} \sim N(0, \sigma^2)$ は誤差項です。パラメータの識別のため、$\sum_{i=1}^k \tau_i = 0$ または $\tau_1 = 0$（基準カテゴリ）の制約を置きます。
一元配置分散分析（One-way ANOVA）
処置間の差を検定するために、一元配置ANOVAを用います。ANOVAの基本的なアイデアは、データの全変動を「処置間変動」と「処置内変動（誤差）」に分解することです。
帰無仮説 $H_0: \tau_1 = \tau_2 = \cdots = \tau_k$（すべての処置効果が等しい）のもとで、$F$ 統計量を構成します。
まず、平方和の分解を行います。
$$ \underbrace{\sum_{i=1}^k \sum_{j=1}^{n_i} (Y_{ij} – \bar{Y}_{\cdot\cdot})^2}_{SS_{\text{total}}} = \underbrace{\sum_{i=1}^k n_i (\bar{Y}_{i\cdot} – \bar{Y}_{\cdot\cdot})^2}_{SS_{\text{between}}} + \underbrace{\sum_{i=1}^k \sum_{j=1}^{n_i} (Y_{ij} – \bar{Y}_{i\cdot})^2}_{SS_{\text{within}}} $$
ここで $\bar{Y}_{i\cdot}$ は第 $i$ 処置群の平均、$\bar{Y}_{\cdot\cdot}$ は全体平均です。
$F$ 統計量は平均平方の比として定義されます。
$$ F = \frac{SS_{\text{between}} / (k – 1)}{SS_{\text{within}} / (n – k)} = \frac{MS_{\text{between}}}{MS_{\text{within}}} $$
$H_0$ のもとで $F \sim F(k-1, n-k)$ に従います。$F$ 値が大きいほど、処置間の差が誤差に対して大きいことを意味し、$H_0$ が棄却されやすくなります。
ANOVAの仮定と診断
ANOVAのF検定が正確であるためには、以下の仮定が必要です。
- 独立性: 各観測値が独立であること。実験のランダム化により保証される
- 正規性: 各処置群内で誤差が正規分布に従うこと。中心極限定理により、サンプルサイズが十分大きければ多少の非正規性は問題にならない
- 等分散性（homoscedasticity）: すべての処置群で分散が等しいこと。分散が不等な場合、F検定のサイズ（第1種の過誤率）が名目水準から乖離する
等分散性の診断にはバートレット検定（Bartlett’s test）やレベネ検定（Levene’s test）が用いられます。レベネ検定は正規性の逸脱に対してバートレット検定よりもロバストです。等分散性が疑わしい場合は、ウェルチのANOVA（Welch’s ANOVA）を使うことが推奨されます。
ANOVAの結果は通常、ANOVA表として整理されます。
| 要因 | 平方和 | 自由度 | 平均平方 | F値 |
|---|---|---|---|---|
| 処置間 | $SS_B$ | $k-1$ | $MS_B$ | $F = MS_B / MS_W$ |
| 処置内 | $SS_W$ | $n-k$ | $MS_W$ | |
| 合計 | $SS_T$ | $n-1$ |
乱塊法（RCBD）
動機
CRDでは実験単位間の変動がすべて誤差項に含まれます。もし実験単位間に既知の変動要因（例: 農場の区画の日当たり、患者の年齢層、工場の生産ライン）がある場合、これをブロック変数として制御することで推定精度を向上させることができます。
RCBDの直感は「リンゴとリンゴを比較する」ことです。異なるブロック（条件の異なるグループ）間で処置を比較するのではなく、同じブロック内で処置を比較します。ブロック間の差は除去されるため、処置効果の推定がより精密になります。
統計モデル
$$ \begin{equation} Y_{ij} = \mu + \tau_i + \beta_j + \epsilon_{ij}, \quad i = 1, \ldots, k, \quad j = 1, \ldots, b \end{equation} $$
ここで $\beta_j$ は第 $j$ ブロックの効果、$b$ はブロック数です。各ブロック内ですべての処置が1回ずつ施されます（したがって合計の観測数は $n = kb$）。
このモデルは加法的であることに注意してください。すなわち、処置効果 $\tau_i$ はブロックによらず一定であると仮定しています。もし処置効果がブロックによって異なる（交互作用が存在する）場合、この加法的モデルは不適切です。
平方和の分解と分散分析
RCBDでは全変動を3つの成分に分解します。
$$ SS_{\text{total}} = SS_{\text{treatment}} + SS_{\text{block}} + SS_{\text{error}} $$
各成分は以下の通りです。
$$ SS_{\text{treatment}} = b \sum_{i=1}^k (\bar{Y}_{i\cdot} – \bar{Y}_{\cdot\cdot})^2, \quad df = k – 1 $$
$$ SS_{\text{block}} = k \sum_{j=1}^b (\bar{Y}_{\cdot j} – \bar{Y}_{\cdot\cdot})^2, \quad df = b – 1 $$
$$ SS_{\text{error}} = SS_{\text{total}} – SS_{\text{treatment}} – SS_{\text{block}}, \quad df = (k-1)(b-1) $$
処置効果の検定は次の $F$ 統計量で行います。
$$ F_{\text{treatment}} = \frac{SS_{\text{treatment}} / (k – 1)}{SS_{\text{error}} / ((k-1)(b-1))} = \frac{MS_{\text{treatment}}}{MS_{\text{error}}} $$
ブロック効果を分離することで、$SS_{\text{error}}$ がCRDの $SS_{\text{within}}$ よりも小さくなります。CRDでは $SS_{\text{within}} = SS_{\text{block}} + SS_{\text{error}}$ であるため、ブロック効果が大きいほど誤差の減少が大きく、検定の検出力が劇的に向上します。
RCBDの仮定
RCBDが適切に機能するためには以下の仮定が必要です。
- 加法性: 処置効果とブロック効果が加法的（交互作用がない）
- 独立性: 誤差 $\epsilon_{ij}$ が独立
- 等分散性: 誤差の分散がすべてのセルで等しい
- 正規性: 誤差が正規分布に従う
加法性の仮定が疑わしい場合は、タキーの加法性検定（Tukey’s test for non-additivity）で確認できます。
ラテン方格法
アイデア
RCBDでは1つのブロック変数を制御しましたが、2つの変動要因を同時に制御したい場合はどうでしょうか。例えば、農業試験で農場の「位置（東西）」と「土壌の種類（南北）」の両方が収量に影響するかもしれません。
ラテン方格法（Latin square design）は、2つのブロック変数（行と列）を同時に制御するデザインです。$k \times k$ の格子において、各行と各列にすべての処置がちょうど1回ずつ現れるように配置します。この配置は数独パズルの構造に似ています。
統計モデル
$$ \begin{equation} Y_{ijk} = \mu + \tau_i + \rho_j + \gamma_k + \epsilon_{ijk} \end{equation} $$
ここで $\tau_i$ は第 $i$ 処置の効果、$\rho_j$ は第 $j$ 行の効果、$\gamma_k$ は第 $k$ 列の効果です。
例（4×4ラテン方格）
| 列1 | 列2 | 列3 | 列4 | |
|---|---|---|---|---|
| 行1 | A | B | C | D |
| 行2 | B | C | D | A |
| 行3 | C | D | A | B |
| 行4 | D | A | B | C |
各処置（A, B, C, D）は各行に1回ずつ、各列に1回ずつ現れています。
平方和の分解
$$ SS_{\text{total}} = SS_{\text{treatment}} + SS_{\text{row}} + SS_{\text{col}} + SS_{\text{error}} $$
誤差の自由度は $(k-1)(k-2)$ です。RCBDの $(k-1)(b-1)$ に比べて自由度が減少する点に注意が必要です。自由度が小さくなるため、処置数が少ない場合（$k \leq 3$）には検出力が低くなる可能性があります。
ラテン方格法の利点と制限
利点: 2つのブロック変数を同時に制御できるため、CRDやRCBDよりもさらに誤差を減らせる可能性があります。
制限: 処置数 $k$、行数、列数がすべて等しくなければなりません。また、自由度が小さくなるため、処置数が少ない場合には推奨されません。さらに、処置と行・列の間に交互作用がないという仮定が必要です。
グレコ・ラテン方格
3つのブロック変数を同時に制御したい場合は、グレコ・ラテン方格（Graeco-Latin square）を使います。2つのラテン方格を重ね合わせて、各セルに2つの処置文字（ギリシャ文字とラテン文字）が入るようにします。4つ目のブロック変数まで制御が可能ですが、条件が厳しくなるため、実用上は3つのブロック変数までが現実的です。
不完備ブロック計画
不完備ブロック計画の必要性
RCBDでは各ブロック内ですべての処置が1回ずつ施される必要がありますが、処置数 $k$ が大きい場合、1つのブロック内にすべての処置を収めることが物理的に不可能なことがあります。
例えば、ワインの官能評価で10種類のワインを比較したいが、1人のパネリストが1回の試飲セッションで評価できるのは3種類までである場合、各ブロック（パネリスト×セッション）にはすべての処置を含められません。このような場合に使われるのが不完備ブロック計画（incomplete block design）です。
釣り合い型不完備ブロック計画（BIBD）
最も整った不完備ブロック計画は釣り合い型不完備ブロック計画（Balanced Incomplete Block Design, BIBD）です。BIBDは以下の条件を満たします。
- $k$ 個の処置、$b$ 個のブロック、各ブロックに $r’$ 個の処置（$r’ < k$）
- 各処置はちょうど $r$ 個のブロックに現れる
- 任意の2つの処置が同じブロックに現れる回数 $\lambda$ が一定
これらのパラメータは次の関係式を満たします。
$$ \begin{equation} br’ = kr, \quad \lambda(k – 1) = r(r’ – 1) \end{equation} $$
BIBDの利点は、すべての処置対が同じ回数だけ直接比較されるため、処置効果の推定値の分散がすべての対で等しくなることです。ただし、整数制約のためにBIBDが存在しないパラメータの組み合わせもあります。
BIBDの分散分析ではブロック内平方和を用いた調整済み処置効果を計算する必要があり、CRDやRCBDよりも複雑になります。
反復不完備ブロック計画
BIBDの条件を満たせない場合でも、各処置対の比較回数をなるべく均等にした部分的に釣り合い型のデザインが使われます。最適デザイン理論に基づいてA最適性やD最適性を基準に計画を構築することも一般的です。
不完備ブロック計画は、実験の物理的制約を考慮しながら統計的効率を最大化するための工夫であり、大規模な実験（農業試験、薬物スクリーニング、官能評価）で広く利用されています。
不完備ブロック計画という発展的な話題を見たところで、複数の因子を同時に扱う要因計画法に移りましょう。
要因計画法
ここまでの実験デザイン（CRD、RCBD、ラテン方格）は、1つの因子の効果を調べるものでした。しかし実際の実験では、複数の因子が同時にアウトカムに影響することが多いです。例えば、化学反応の収率は温度と触媒の種類の両方に依存するかもしれません。
要因計画法（factorial design）は、複数の因子の効果と交互作用を同時に効率的に評価するデザインです。
$2^k$ 要因計画
最も基本的な要因計画は、$k$ 個の因子がそれぞれ2水準（高/低、あり/なし）を持つ $2^k$ 要因計画です。$2^k$ 通りの処置の組み合わせがあります。
例えば、因子A（温度: 高/低）と因子B（圧力: 高/低）の $2^2 = 4$ 通りの組み合わせを調べます。
| 条件 | A（温度） | B（圧力） |
|---|---|---|
| (1) | 低 | 低 |
| a | 高 | 低 |
| b | 低 | 高 |
| ab | 高 | 高 |
$2^k$ 要因計画の重要な利点は、1因子ずつ変える実験（One-Factor-At-a-Time, OFAT）と比べて効率が良いことです。OFATでは因子Aの効果を調べた後に因子Bの効果を調べますが、要因計画ではすべてのデータがすべての因子の効果推定に使われます。
主効果と交互作用
2因子の要因計画モデルは以下のように書けます。
$$ \begin{equation} Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk} \end{equation} $$
$\alpha_i$ は因子Aの主効果（main effect）、$\beta_j$ は因子Bの主効果、$(\alpha\beta)_{ij}$ はAとBの交互作用効果（interaction effect）です。
主効果は、他の因子の水準を平均化したときの、ある因子の各水準間の差です。例えば、因子Aの主効果は「温度を高にしたとき、圧力の水準を平均した収率の変化」です。
交互作用効果は、一方の因子の効果が他方の因子の水準によって異なる場合に生じます。例えば、「低圧では温度を上げると収率が上がるが、高圧では温度を上げると収率が下がる」という場合、温度と圧力に交互作用があります。
交互作用が存在する場合、因子Aの効果は因子Bの水準によって異なるため、主効果だけを見ても不十分です。交互作用の有無は実験結果の解釈に大きく影響するため、要因計画の最も重要な特徴は交互作用を検出できることにあります。
一部実施要因計画
因子数 $k$ が大きくなると、$2^k$ の全組み合わせを実施するのはコストが高くなります。例えば $k = 7$ では $2^7 = 128$ 条件が必要です。
一部実施要因計画（fractional factorial design）は、全組み合わせの一部（通常 $2^{k-p}$ 通り、$p \geq 1$）のみを実施するデザインです。これにより実験回数を大幅に削減できますが、一部の効果が交絡（confounding / aliasing）し、区別できなくなるというトレードオフがあります。
例えば、$2^{7-4}$ の一部実施計画では $2^3 = 8$ 条件で7因子の主効果を推定できますが、主効果と2因子交互作用が交絡する場合があります。解像度（resolution）という概念で交絡の程度を評価し、一般に解像度IIIは主効果が2因子交互作用と交絡し、解像度Vでは主効果と2因子交互作用がともに推定可能です。
効果の推定とスパースの原理
一部実施要因計画が実用的に機能する背景には、効果のスパース性原理（effect sparsity principle / Pareto principle）があります。これは「多くの因子のうち、実際にアウトカムに大きな影響を与えるのはごく少数の因子である」という経験則です。
この原理に基づけば、効果の階層性原理（effect hierarchy principle）も成り立ちます。すなわち、低次の効果（主効果）は高次の効果（高次交互作用）よりも重要である可能性が高いです。したがって、高次交互作用が無視できるという仮定のもとで、一部実施要因計画は効率的に主効果と低次交互作用を推定できます。
プラケット・バーマン計画
$k$ 個の因子をわずか $k+1$ 回の実験で主効果のスクリーニングを行いたい場合に有用なのがプラケット・バーマン計画（Plackett-Burman design）です。$N = 4m$（$m$ は正の整数）回の実験で $N-1$ 個の因子の主効果を推定できます。例えば、$N = 12$ 回の実験で11因子の主効果をスクリーニングできます。ただし解像度はIIIであるため、主効果と2因子交互作用が交絡する可能性があります。
Pythonシミュレーション
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# --- 完全無作為化デザイン（CRD）vs 乱塊法（RCBD）---
k = 4 # 処置数
b = 5 # ブロック数
n_per = 5 # CRDの各処置のサンプル数
# 真の処置効果
tau = np.array([0, 2, 4, 1]) # 処置1-4の効果
block_effect = np.array([0, 3, 6, 1, 8]) # ブロック効果（大きな変動）
# --- CRDのデータ生成 ---
Y_crd = []
T_crd = []
for i in range(k):
for j in range(n_per):
block = np.random.choice(b)
y = 10 + tau[i] + block_effect[block] + np.random.normal(0, 1)
Y_crd.append(y)
T_crd.append(i)
Y_crd = np.array(Y_crd)
T_crd = np.array(T_crd)
# CRDのANOVA
groups_crd = [Y_crd[T_crd == i] for i in range(k)]
F_crd, p_crd = stats.f_oneway(*groups_crd)
# --- RCBDのデータ生成 ---
Y_rcbd = np.zeros((k, b))
for i in range(k):
for j in range(b):
Y_rcbd[i, j] = 10 + tau[i] + block_effect[j] + np.random.normal(0, 1)
# RCBDのANOVA（手動計算）
grand_mean = Y_rcbd.mean()
SS_treat = b * np.sum((Y_rcbd.mean(axis=1) - grand_mean)**2)
SS_block = k * np.sum((Y_rcbd.mean(axis=0) - grand_mean)**2)
SS_total = np.sum((Y_rcbd - grand_mean)**2)
SS_error = SS_total - SS_treat - SS_block
df_treat = k - 1
df_block = b - 1
df_error = (k - 1) * (b - 1)
MS_treat = SS_treat / df_treat
MS_error = SS_error / df_error
F_rcbd = MS_treat / MS_error
p_rcbd = 1 - stats.f.cdf(F_rcbd, df_treat, df_error)
print(f"=== CRD ===")
print(f"F統計量: {F_crd:.4f}, p値: {p_crd:.4f}")
print(f"\n=== RCBD ===")
print(f"F統計量: {F_rcbd:.4f}, p値: {p_rcbd:.6f}")
print(f"MS_treat: {MS_treat:.4f}")
print(f"MS_error: {MS_error:.4f}")
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) CRD
ax = axes[0]
bp = ax.boxplot(groups_crd, labels=[f"Trt {i+1}" for i in range(k)],
patch_artist=True)
colors = ["lightblue", "lightsalmon", "lightgreen", "lightyellow"]
for patch, color in zip(bp["boxes"], colors):
patch.set_facecolor(color)
ax.set_ylabel("Outcome Y", fontsize=12)
ax.set_title(f"CRD: F={F_crd:.2f}, p={p_crd:.4f}", fontsize=11)
ax.grid(True, alpha=0.3)
# (b) RCBD
ax = axes[1]
x = np.arange(k)
w = 0.15
for j in range(b):
ax.bar(x + j * w, Y_rcbd[:, j], w, label=f"Block {j+1}", alpha=0.7)
ax.set_xticks(x + w * (b - 1) / 2)
ax.set_xticklabels([f"Trt {i+1}" for i in range(k)], fontsize=10)
ax.set_ylabel("Outcome Y", fontsize=12)
ax.set_title(f"RCBD: F={F_rcbd:.2f}, p={p_rcbd:.6f}", fontsize=11)
ax.legend(fontsize=7, ncol=3)
ax.grid(True, alpha=0.3, axis="y")
# (c) 検出力の比較
ax = axes[2]
n_sims = 1000
crd_significant = 0
rcbd_significant = 0
for _ in range(n_sims):
# CRD
Y_c = []
T_c = []
for i in range(k):
for j in range(n_per):
blk = np.random.choice(b)
y = 10 + tau[i] + block_effect[blk] + np.random.normal(0, 1)
Y_c.append(y)
T_c.append(i)
Y_c = np.array(Y_c)
T_c = np.array(T_c)
_, p = stats.f_oneway(*[Y_c[T_c==i] for i in range(k)])
if p < 0.05:
crd_significant += 1
# RCBD
Y_r = np.zeros((k, b))
for i in range(k):
for j in range(b):
Y_r[i, j] = 10 + tau[i] + block_effect[j] + np.random.normal(0, 1)
gm = Y_r.mean()
sst = b * np.sum((Y_r.mean(axis=1) - gm)**2)
ssb = k * np.sum((Y_r.mean(axis=0) - gm)**2)
sse = np.sum((Y_r - gm)**2) - sst - ssb
f = (sst / (k-1)) / (sse / ((k-1)*(b-1)))
p = 1 - stats.f.cdf(f, k-1, (k-1)*(b-1))
if p < 0.05:
rcbd_significant += 1
power_crd = crd_significant / n_sims
power_rcbd = rcbd_significant / n_sims
bars = ax.bar(["CRD", "RCBD"], [power_crd, power_rcbd],
color=["steelblue", "salmon"], alpha=0.7, edgecolor="black")
for bar, val in zip(bars, [power_crd, power_rcbd]):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.02,
f"{val:.3f}", ha="center", fontsize=12)
ax.set_ylabel("Power", fontsize=12)
ax.set_title(f"Power comparison ({n_sims} simulations)", fontsize=11)
ax.set_ylim(0, 1.1)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("experimental_design.png", dpi=150, bbox_inches="tight")
plt.show()
シミュレーション結果から以下のことが確認できます。
-
CRDではブロック間の変動が誤差に含まれるため、箱ひげ図の幅が広く、F統計量とp値が比較的弱い結果になっています
-
RCBDではブロック効果を除去できるため、処置間の差がより明確に検出されます。F統計量が大幅に改善し、p値が小さくなっています
-
検出力の比較（右図）: RCBDの検出力がCRDを大きく上回っています。ブロック効果が大きい場合、RCBDの効率向上は劇的です
分割区計画と交差計画
分割区計画（Split-Plot Design）
要因計画法では各因子の水準を自由に組み合わせられることを前提としていましたが、実際には一部の因子を変えるのが困難な場合があります。例えば、製造ラインの温度設定を変えるには大きなコストと時間がかかるため、温度を変えるたびにすべての他の因子の組み合わせを試すのが合理的です。
分割区計画（split-plot design）は、変更が困難な因子（whole-plot factor）と容易な因子（sub-plot factor）を区別して効率的に実験を行うデザインです。
分割区計画では、まず大きな実験単位（whole plot）にwhole-plot因子を割り当て、次にその中の小さな実験単位（sub-plot）にsub-plot因子を割り当てます。この結果、2つのレベルの誤差項が存在し、分散分析ではwhole-plot誤差とsub-plot誤差を別々に推定する必要があります。
whole-plot因子の検定にはwhole-plot誤差を使い、sub-plot因子と交互作用の検定にはsub-plot誤差を使います。一般にwhole-plot誤差の方が大きい（反復が少ない）ため、whole-plot因子の検出力はsub-plot因子よりも低くなります。
交差（クロスオーバー）計画
交差計画（crossover design）は、各被験者が複数の処置を順番に受けるデザインです。時期（period）と処置（treatment）の両方の効果を推定できます。
最も単純な2処置2時期の交差計画（AB/BA交差計画）では、被験者を2群に分け、第1群は「処置A → 処置B」の順に、第2群は「処置B → 処置A」の順に処置を受けます。各被験者が自分自身の対照となるため、被験者間変動を除去でき、統計的効率が高くなります。
交差計画における主な注意点は持ち越し効果（carryover effect）です。第1期の処置の効果が第2期に持ち越される場合、第2期の結果が汚染されます。持ち越し効果を最小化するために、期間間に洗い出し期間（washout period）を設けるのが一般的です。
持ち越し効果がないとき、交差計画の統計モデルは以下の通りです。
$$ \begin{equation} Y_{ijk} = \mu + \pi_i + \tau_{d(i,j)} + \lambda_j + \epsilon_{ijk} \end{equation} $$
ここで $\pi_i$ は被験者 $i$ の効果、$\tau_{d(i,j)}$ は時期 $j$ に被験者 $i$ が受けた処置の効果、$\lambda_j$ は時期 $j$ の効果です。
交差計画は臨床試験（慢性疾患の対症療法の比較など）で特に有用ですが、疾患の状態が時間とともに変化する場合や、処置の効果が持続的な場合には適用が困難です。
さまざまなデザインの応用を見たところで、実験結果の解析における多重比較の問題に進みましょう。
多重比較法
ANOVAで $H_0$（すべての処置効果が等しい）が棄却された場合、「どの処置とどの処置の間に差があるのか」を特定する必要があります。この問題を多重比較（multiple comparisons）と呼びます。
多重性の問題
$k$ 個の処置間で全ペアを比較すると $\binom{k}{2} = k(k-1)/2$ 通りの検定が必要です。各検定を有意水準 $\alpha = 0.05$ で行うと、少なくとも1つの偽陽性が生じる確率（族単位誤り率, FWER）は $1 – (1 – \alpha)^{k(k-1)/2}$ まで膨らみます。$k = 5$ では10通りの比較で、FWERは約40%にもなります。
ボンフェローニ補正
最も単純な方法は、各比較の有意水準を $\alpha / m$（$m$: 比較の数）に引き下げることです。保守的ですが、あらゆる状況で適用でき、FWERを $\alpha$ 以下に制御します。
テューキーのHSD検定
Tukey’s HSD（Honestly Significant Difference）は、全ペア比較に特化した検定法です。すべてのペアの平均値の差を、スチューデント化された範囲統計量（studentized range statistic）$q$ の分布と比較します。ボンフェローニ補正よりも検出力が高いという利点があります。
$$ q = \frac{\bar{Y}_{i\cdot} – \bar{Y}_{j\cdot}}{\sqrt{MS_E / n_i}} $$
ダネット検定
Dunnett’s testは、対照群との比較に特化した方法です。$k – 1$ 個の処置群をそれぞれ対照群と比較する場合に、全ペア比較よりも検出力が高くなります。新薬の効果をプラセボと比較する臨床試験などで使われます。
応答曲面法
応答曲面法の概要
要因計画法がカテゴリカルな因子や2水準の因子の効果を調べるのに適しているのに対し、応答曲面法（Response Surface Methodology, RSM）は連続量の因子の最適な組み合わせを見つけることを目的とします。RSMはGeorge E. P. BoxとK. B. Wilsonによって1951年に提案されました。
RSMの基本的な流れは以下の通りです。
- スクリーニング段階: 一部実施要因計画やプラケット・バーマン計画を使って、重要な因子を絞り込む
- 一次モデルの当てはめ: $2^k$ 要因計画を実施し、一次（線形）モデルを当てはめる。最急上昇法（steepest ascent）で応答の改善方向を探索する
- 二次モデルの当てはめ: 最適点の近傍で中心複合計画（CCD）やBox-Behnken計画を実施し、二次（曲面）モデルを当てはめる
- 最適点の同定: 二次モデルの停留点を計算し、応答を最大化（または最小化）する因子の組み合わせを求める
二次モデル
RSMの中心は、2因子の場合、以下の二次モデルです。
$$ \begin{equation} Y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{11} x_1^2 + \beta_{22} x_2^2 + \beta_{12} x_1 x_2 + \epsilon \end{equation} $$
このモデルの停留点は $\nabla Y = \bm{0}$ を解くことで求まります。停留点が最大値、最小値、鞍点のいずれであるかは、$\beta_{11}$, $\beta_{22}$, $\beta_{12}$ から構成されるヘッセ行列の固有値の符号で判定されます。
中心複合計画（CCD）
二次モデルを推定するための標準的なデザインが中心複合計画（Central Composite Design, CCD）です。CCDは3つの要素で構成されます。
- 要因点: $2^k$ 要因計画の各点（立方体の頂点）
- 軸点（axial points / star points）: 各因子軸上の $\pm \alpha$ の位置にある $2k$ 個の点
- 中心点: 原点に置かれた $n_c$ 個の反復
$\alpha$ の値を適切に選ぶことで、回転可能性（rotatability）——原点からの距離が同じ点での予測分散が等しくなる性質——を達成できます。$\alpha = \sqrt[4]{2^k}$ が回転可能な設計を与えます。
RSMは製造プロセスの最適化、食品の配合設計、化学反応条件の探索など、実用的な最適化問題で広く使われています。
検出力分析とサンプルサイズの決定
実験を始める前に「どのくらいのサンプルサイズが必要か」を決めることは、実験計画の重要なステップです。
検出力の定義
検出力（power）$1 – \beta$ は、対立仮説が真であるときに帰無仮説を正しく棄却する確率です。一般に、検出力 $\geq 0.8$ が目標とされます。
検出力に影響する要素
- 効果量（effect size）: 処置間の差の大きさ。効果量が大きいほど検出しやすい
- 有意水準 $\alpha$: 有意水準を小さくすると検出力は下がる
- サンプルサイズ $n$: サンプルサイズが大きいほど検出力は上がる
- 誤差の分散 $\sigma^2$: 分散が小さいほど検出力は上がる（ブロッキングが有効な理由）
CRDの一元配置ANOVAの場合、非心度パラメータは次のように計算されます。
$$ \lambda = \frac{n \sum_{i=1}^k \tau_i^2}{\sigma^2} $$
この非心度パラメータと $F$ 分布の非心分布を用いて検出力を計算し、目標の検出力を達成するサンプルサイズ $n$ を求めます。
まとめ
- フィッシャーの3原則（反復・ランダム化・局所管理）は実験計画の基盤
- CRDは最も単純だが、ブロック内変動を制御できない
- RCBDはブロック変数で既知の変動を制御し、推定精度と検出力を向上させる
- ラテン方格法は2つのブロック変数を同時に制御する
- 要因計画法は複数因子の主効果と交互作用を効率的に推定する
次のステップとして、以下の記事も参考にしてください。
- ネイマン・ピアソンの基本補題 — 仮説検定の最適性理論
- ワルド検定 — パラメータの検定理論