ある統計モデルのパラメータが「特定の値であるかどうか」を判定したいとき、私たちはどのような方法を使えるでしょうか。たとえば、新しい治療法の効果を調べるとき、「治療効果パラメータ $\theta$ はゼロか（＝効果なし）、それともゼロでないか」を検定したい場面が頻繁に登場します。
このような仮説検定に対する古典的なアプローチとして、尤度比検定（Likelihood Ratio Test）、ワルド検定（Wald Test）、そしてスコア検定（Score Test）の三つが知られています。これらは「三大漸近検定」と呼ばれ、大標本のもとで同等の性能を持つことが知られていますが、それぞれの計算手順には大きな違いがあります。
スコア検定は、C. R. ラオ（C. R. Rao）によって1948年に提案されたことからラオ検定（Rao Test）とも呼ばれます。この検定の最大の特徴は、帰無仮説のもとでのみパラメータを推定すればよいという計算上の利点です。対立仮説のもとでの最尤推定（MLE）を必要としないため、モデルが複雑で対立仮説のもとでのMLEが困難な場合に特に有用です。
スコア検定を理解すると、以下のような幅広い応用が開けます。
- 一般化線形モデル（GLM）: ロジスティック回帰やポアソン回帰で、変数の追加が有意かどうかを帰無仮説のもとの推定のみで判定できます
- 生存時間分析: Cox比例ハザードモデルにおいて、共変量の効果をスコア検定で評価する手法は臨床試験で広く使われています
- 遺伝統計学: ゲノムワイド関連研究（GWAS）では、数百万の遺伝子座を検定する必要があり、スコア検定の計算効率が不可欠です
- 計量経済学: ラグランジュ乗数検定（LM検定）として知られ、制約付きモデルの残差から検定統計量を構成します
本記事の内容
- スコア関数とフィッシャー情報量の直感的理解
- スコア検定統計量の数学的な導出（省略なし）
- 尤度比検定・ワルド検定との幾何学的な関係
- 多次元パラメータへの拡張
- Pythonによる実装と数値シミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
スコア関数の直感的な意味
対数尤度の「勾配」としてのスコア
スコア検定を理解するために、まずスコア関数（score function）の直感的な意味を押さえましょう。
観測データ $x_1, x_2, \dots, x_n$ が確率分布 $f(x; \theta)$ に従うとき、対数尤度関数は次のように定義されます。
$$ \ell(\theta) = \sum_{i=1}^{n} \ln f(x_i; \theta) $$
対数尤度関数は「パラメータ $\theta$ がこの値だったら、手元のデータがどれだけ『もっともらしい』か」を表す指標です。スコア関数はこの対数尤度関数の $\theta$ に関する微分（勾配）です。
$$ \begin{equation} S(\theta) = \frac{\partial \ell(\theta)}{\partial \theta} \end{equation} $$
スコア関数の値は、「パラメータを少しだけ変えたとき、対数尤度がどれだけ変わるか」を教えてくれます。山登りに例えると、対数尤度関数が山の地形で、スコア関数はその地点での「傾き」です。
スコア関数がゼロになる点 — 最尤推定量
山の頂上では傾きがゼロになるのと同じように、対数尤度の最大値を与える点（最尤推定量 $\hat{\theta}_{\text{MLE}}$）ではスコア関数がゼロになります。
$$ S(\hat{\theta}_{\text{MLE}}) = 0 $$
これは最尤推定量を求めるための基本方程式（スコア方程式、尤度方程式）です。逆に言えば、スコア関数がゼロから大きくずれているパラメータ値は、データに対して「もっともらしくない」ことを意味します。
スコア検定の基本アイデア
ここがスコア検定の核心です。帰無仮説 $H_0: \theta = \theta_0$ のもとでスコア関数 $S(\theta_0)$ を計算してみましょう。
- もし $\theta_0$ が真の値に近ければ、$S(\theta_0)$ はゼロに近い値を取るはずです。なぜなら、最尤推定量 $\hat{\theta}_{\text{MLE}}$ は真の値の近くにあり、$S(\hat{\theta}_{\text{MLE}}) = 0$ だからです
- もし $\theta_0$ が真の値から離れていれば、$S(\theta_0)$ はゼロから大きくずれます。対数尤度の山の中腹にいるようなもので、傾きが急になります
つまり、$S(\theta_0)$ の値が「大きい」か「小さい」かを見ることで、$\theta_0$ が妥当かどうかを判定できるのです。これがスコア検定のアイデアです。
ただし、$S(\theta_0)$ の「大きさ」を評価するには、その分散（ばらつき）を考慮する必要があります。そこで登場するのがフィッシャー情報量です。
フィッシャー情報量とスコアの分散
スコア関数の期待値と分散
スコア関数には、2つの基本的な性質があります。
性質1: 期待値がゼロ
真のパラメータ $\theta_0$ のもとで、スコア関数の期待値はゼロです。
$$ \begin{equation} E_{\theta_0}[S(\theta_0)] = 0 \end{equation} $$
この性質の証明は以下の通りです。$\int f(x; \theta)\,dx = 1$ の両辺を $\theta$ で微分すると、
$$ \int \frac{\partial f(x; \theta)}{\partial \theta}\,dx = 0 $$
となります。ここで $\frac{\partial f}{\partial \theta} = f \cdot \frac{\partial \ln f}{\partial \theta}$ を使うと、
$$ \int \frac{\partial \ln f(x; \theta)}{\partial \theta} f(x; \theta)\,dx = E\left[\frac{\partial \ln f}{\partial \theta}\right] = E[S(\theta)] = 0 $$
が得られます。つまり、「正しいパラメータ値でのスコア関数は平均的にゼロ」ということです。
性質2: 分散がフィッシャー情報量
スコア関数の分散は、フィッシャー情報量（Fisher information）$I(\theta)$ に等しくなります。
$$ \begin{equation} \text{Var}_{\theta_0}[S(\theta_0)] = E_{\theta_0}[S(\theta_0)^2] = I(\theta_0) \end{equation} $$
期待値がゼロなので、分散は二次モーメントに一致します。フィッシャー情報量は、正則条件のもとで次のようにも表されます。
$$ \begin{equation} I(\theta) = -E\left[\frac{\partial^2 \ell(\theta)}{\partial \theta^2}\right] \end{equation} $$
直感的には、フィッシャー情報量は「対数尤度関数の曲がり具合（曲率）の期待値」です。対数尤度の山が鋭く尖っていれば情報量が大きく、なだらかであれば情報量が小さくなります。鋭い山は「パラメータの推定精度が高い」ことを意味し、なだらかな山は「推定が不確実」であることを意味します。
n 個のデータに対する情報量
$n$ 個の独立同一分布（i.i.d.）の観測値に対しては、スコア関数とフィッシャー情報量はそれぞれ次のようになります。
$$ S_n(\theta) = \sum_{i=1}^{n} \frac{\partial \ln f(x_i; \theta)}{\partial \theta}, \quad I_n(\theta) = n \cdot I_1(\theta) $$
ここで $I_1(\theta)$ は1個のデータに対するフィッシャー情報量です。情報量はデータ数 $n$ に比例して増加するため、データが多いほどパラメータの推定精度が高くなります。
これらの準備ができたので、いよいよスコア検定統計量を厳密に導出しましょう。
スコア検定統計量の導出
検定問題の設定
検定問題を次のように設定します。
$$ H_0: \theta = \theta_0 \quad \text{vs.} \quad H_1: \theta \neq \theta_0 $$
データ $x_1, \dots, x_n$ は $f(x; \theta)$ から i.i.d. で生成されるとします。
スコア関数の漸近正規性
中心極限定理により、十分大きな $n$ のもとでは、スコア関数 $S_n(\theta_0)$ は漸近的に正規分布に従います。帰無仮説 $H_0: \theta = \theta_0$ が正しいとき、
$$ \begin{equation} S_n(\theta_0) \xrightarrow{d} N\left(0,\, I_n(\theta_0)\right) \end{equation} $$
これは、$S_n(\theta_0)$ が $n$ 個のスコア $\frac{\partial \ln f(x_i; \theta_0)}{\partial \theta}$ の和であり、各項が独立で期待値ゼロ、分散 $I_1(\theta_0)$ を持つことから、中心極限定理が適用できるためです。
標準化されたスコア検定統計量
漸近正規性から、スコア関数を標準化した量を考えます。
$$ \frac{S_n(\theta_0)}{\sqrt{I_n(\theta_0)}} \xrightarrow{d} N(0, 1) $$
この標準化された量の2乗をとると、スコア検定統計量（Score test statistic）が得られます。
$$ \begin{equation} T_S = \frac{S_n(\theta_0)^2}{I_n(\theta_0)} \end{equation} $$
帰無仮説のもとで、$T_S$ は漸近的に自由度1の $\chi^2$ 分布に従います。
$$ \begin{equation} T_S \xrightarrow{d} \chi^2(1) \quad \text{under } H_0 \end{equation} $$
これがスコア検定統計量です。$T_S$ の値が大きい（$\chi^2(1)$ 分布の上側確率が有意水準以下になる）とき、帰無仮説を棄却します。
情報量の推定
実際の計算では、フィッシャー情報量 $I_n(\theta_0)$ を解析的に求められない場合があります。そのときは以下のいずれかの推定量を使います。
方法1: 期待フィッシャー情報量
$$ I_n(\theta_0) = -E_{\theta_0}\left[\frac{\partial^2 \ell(\theta)}{\partial \theta^2}\bigg|_{\theta=\theta_0}\right] $$
方法2: 観測フィッシャー情報量
$$ \hat{I}_n(\theta_0) = -\frac{\partial^2 \ell(\theta)}{\partial \theta^2}\bigg|_{\theta=\theta_0} $$
方法3: スコア関数の外積推定（OPG推定量）
$$ \tilde{I}_n(\theta_0) = \sum_{i=1}^{n} \left[\frac{\partial \ln f(x_i; \theta_0)}{\partial \theta}\right]^2 $$
これら3つの推定量は、正則条件のもとで大標本では同じ値に収束しますが、有限標本では異なる値を取ることがあります。観測情報量は計算が容易で、期待情報量は解析的な式が利用できる場合に便利です。
ここまでの議論は1次元パラメータの場合でした。次に、多次元パラメータの場合へ拡張しましょう。
多次元パラメータへの拡張
ベクトルパラメータの場合
実際の統計モデルでは、パラメータは $\bm{\theta} = (\theta_1, \theta_2, \dots, \theta_p)^T$ のように多次元であることが一般的です。この場合、スコア関数はベクトル、フィッシャー情報量は行列になります。
スコアベクトル:
$$ \begin{equation} \bm{S}_n(\bm{\theta}) = \frac{\partial \ell(\bm{\theta})}{\partial \bm{\theta}} = \left(\frac{\partial \ell}{\partial \theta_1}, \frac{\partial \ell}{\partial \theta_2}, \dots, \frac{\partial \ell}{\partial \theta_p}\right)^T \end{equation} $$
フィッシャー情報行列:
$$ \begin{equation} \bm{I}_n(\bm{\theta}) = -E\left[\frac{\partial^2 \ell(\bm{\theta})}{\partial \bm{\theta} \partial \bm{\theta}^T}\right] \end{equation} $$
$\bm{I}_n$ は $p \times p$ の正定値対称行列です。
多次元スコア検定統計量
帰無仮説 $H_0: \bm{\theta} = \bm{\theta}_0$ に対するスコア検定統計量は、スコアベクトルの二次形式として定義されます。
$$ \begin{equation} T_S = \bm{S}_n(\bm{\theta}_0)^T \bm{I}_n(\bm{\theta}_0)^{-1} \bm{S}_n(\bm{\theta}_0) \end{equation} $$
帰無仮説のもとで、$T_S$ は漸近的に自由度 $p$ の $\chi^2$ 分布に従います。
$$ T_S \xrightarrow{d} \chi^2(p) \quad \text{under } H_0 $$
1次元の場合の $S^2/I$ を、ベクトルと行列の言葉で自然に拡張した形になっています。情報行列の逆行列 $\bm{I}_n^{-1}$ を掛けることで、スコアベクトルの各成分の分散とそれらの間の相関を適切に考慮しています。
複合仮説の場合 — パラメータの一部を検定する
最も実用的な状況は、パラメータの一部だけを検定する場合です。パラメータを $\bm{\theta} = (\bm{\psi}^T, \bm{\lambda}^T)^T$ と分割し、$\bm{\psi}$ が興味のあるパラメータ（$q$ 次元）、$\bm{\lambda}$ が局外パラメータ（$p – q$ 次元）であるとします。
検定問題は次のようになります。
$$ H_0: \bm{\psi} = \bm{\psi}_0 \quad \text{vs.} \quad H_1: \bm{\psi} \neq \bm{\psi}_0 $$
ここで $\bm{\lambda}$ は帰無仮説のもとで自由（未知）です。
この場合、帰無仮説のもとでの最尤推定は $\bm{\psi} = \bm{\psi}_0$ に固定したうえで $\bm{\lambda}$ のみを推定します。この制約付き最尤推定量を $\hat{\bm{\lambda}}_0$ とし、$\tilde{\bm{\theta}}_0 = (\bm{\psi}_0^T, \hat{\bm{\lambda}}_0^T)^T$ とします。
スコア検定統計量は、$\bm{\psi}$ に対するスコア（$\bm{\lambda}$ に対するスコアは $\hat{\bm{\lambda}}_0$ で消去されている）を用いて構成されます。情報行列を以下のようにブロック分割します。
$$ \bm{I}_n = \begin{pmatrix} \bm{I}_{\psi\psi} & \bm{I}_{\psi\lambda} \\ \bm{I}_{\lambda\psi} & \bm{I}_{\lambda\lambda} \end{pmatrix} $$
$\bm{\psi}$ の周辺的なスコア検定統計量は、$\bm{\psi}$ に関する有効スコア（efficient score）を用いて次のように表されます。
$$ \begin{equation} T_S = \bm{S}_{\psi}^T \left(\bm{I}_{\psi\psi} – \bm{I}_{\psi\lambda}\bm{I}_{\lambda\lambda}^{-1}\bm{I}_{\lambda\psi}\right)^{-1} \bm{S}_{\psi} \end{equation} $$
ここで $\bm{S}_{\psi}$ は $\tilde{\bm{\theta}}_0$ で評価した $\bm{\psi}$ に関するスコアベクトルです。括弧内のシューア補行列 $\bm{I}_{\psi\psi} – \bm{I}_{\psi\lambda}\bm{I}_{\lambda\lambda}^{-1}\bm{I}_{\lambda\psi}$ は、$\bm{\lambda}$ の影響を除去した後の $\bm{\psi}$ に関する有効情報量です。
この検定統計量は帰無仮説のもとで漸近的に $\chi^2(q)$ に従います。
重要なのは、この計算全体が帰無仮説のもとでの推定のみで完結するということです。対立仮説のもとで $\bm{\psi}$ と $\bm{\lambda}$ を同時に推定する必要がありません。
次に、スコア検定を尤度比検定・ワルド検定と比較し、三者の幾何学的な関係を明らかにしましょう。
三大漸近検定の比較 — 幾何学的な理解
三つの検定統計量
三大漸近検定の検定統計量を並べて比較します。いずれも帰無仮説 $H_0: \theta = \theta_0$（1次元の場合）に対する検定です。
尤度比検定（LR検定）:
$$ \begin{equation} T_{LR} = 2\left[\ell(\hat{\theta}) – \ell(\theta_0)\right] \end{equation} $$
対数尤度の高さの差を見ます。帰無仮説のもとでのモデルと制約なしのモデルの間で、対数尤度がどれだけ改善するかを評価します。
ワルド検定:
$$ \begin{equation} T_W = I_n(\hat{\theta})(\hat{\theta} – \theta_0)^2 \end{equation} $$
最尤推定量 $\hat{\theta}$ と帰無仮説の値 $\theta_0$ の水平距離を見ます。推定値が帰無仮説の値からどれだけ離れているかを、推定精度（情報量）で標準化して評価します。
スコア検定:
$$ \begin{equation} T_S = \frac{S_n(\theta_0)^2}{I_n(\theta_0)} \end{equation} $$
帰無仮説の点 $\theta_0$ での対数尤度の傾き（勾配） を見ます。勾配がゼロに近ければ帰無仮説の値が尤度の頂上に近く、勾配が急であれば帰無仮説の値は頂上から離れていると判断します。
対数尤度曲線上の幾何学的解釈
3つの検定統計量は、対数尤度曲線 $\ell(\theta)$ 上の異なる幾何学的量を測定しています。
- LR検定: 曲線上の2点 $(\theta_0, \ell(\theta_0))$ と $(\hat{\theta}, \ell(\hat{\theta}))$ の垂直距離（高さの差）
- ワルド検定: 2点 $\theta_0$ と $\hat{\theta}$ の水平距離
- スコア検定: $\theta_0$ における曲線の接線の傾き
もし対数尤度関数が完全な二次関数（放物線）であれば、これら3つの量は完全に同値になります。正規分布の平均の検定がまさにこの例です。一般的な場合には、対数尤度関数は厳密な二次関数ではないため3つの検定は異なる値を返しますが、大標本のもとでは対数尤度関数が二次近似で十分よく近似されるため、漸近的に同値になります。
三者の漸近的同値性
帰無仮説のもとで、3つの検定統計量はすべて漸近的に同じ $\chi^2$ 分布に従います。
$$ T_{LR},\; T_W,\; T_S \xrightarrow{d} \chi^2(1) \quad \text{under } H_0 $$
さらに、対立仮説のもとでも、これら3つの検定は同じ漸近検出力を持ちます。つまり、十分に大きな標本サイズでは、どの検定を使っても実質的に同じ結論が得られます。
有限標本での違い
大標本では同等とはいえ、有限標本では3つの検定統計量の値が異なることがあります。一般的な経験則として、次の不等式が成り立つことが多いです。
$$ T_W \geq T_{LR} \geq T_S $$
これは、ワルド検定が最も「自由」な（帰無仮説を棄却しやすい）検定であり、スコア検定が最も「保守的」な検定であることを意味します。ただし、この順序は必ずしも常に成り立つわけではなく、パラメータ化の方法にも依存します。
各検定の計算コストの比較
3つの検定の実用上の最大の違いは、計算に必要なものです。
| 検定 | 必要な計算 | 計算コスト |
|---|---|---|
| 尤度比検定 | $\hat{\theta}$（制約なしMLE）と $\theta_0$ での対数尤度 | 制約なしMLEの計算が必要 |
| ワルド検定 | $\hat{\theta}$（制約なしMLE）と情報行列 | 制約なしMLEの計算が必要 |
| スコア検定 | $\theta_0$ でのスコアと情報量 | 帰無仮説のもとでの計算のみ |
スコア検定は、対立仮説のもとでの最尤推定を必要としないため、次のような場面で計算上の大きな優位性を持ちます。
- 制約なしMLEの計算が困難な場合: 混合モデルや非線形モデルでは、MLEの数値計算が収束しにくいことがあります
- 多数の検定を同時に行う場合: ゲノムワイド関連研究では数百万の遺伝子座を検定するため、各検定の計算コストを最小限にする必要があります
- 帰無モデルの推定が一度で済む場合: 複数の共変量を一つずつ追加する変数選択の場面では、帰無モデル（基本モデル）の推定を1回行うだけで、すべての候補変数に対するスコア検定が計算できます
スコア検定の理論的基盤を理解したところで、次に具体的な例を通じて計算手順を確認しましょう。
具体例: 正規分布の平均の検定
問題設定
最もシンプルな例として、正規分布 $N(\mu, \sigma^2)$ の平均 $\mu$ に対するスコア検定を導出します。分散 $\sigma^2$ は既知とします。
$$ H_0: \mu = \mu_0 \quad \text{vs.} \quad H_1: \mu \neq \mu_0 $$
データ $x_1, \dots, x_n \overset{\text{i.i.d.}}{\sim} N(\mu, \sigma^2)$ が与えられています。
スコア関数の計算
対数尤度関数は次の通りです。
$$ \ell(\mu) = -\frac{n}{2}\ln(2\pi\sigma^2) – \frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i – \mu)^2 $$
$\mu$ で微分してスコア関数を求めます。
$$ S_n(\mu) = \frac{\partial \ell}{\partial \mu} = \frac{1}{\sigma^2}\sum_{i=1}^{n}(x_i – \mu) = \frac{n(\bar{x} – \mu)}{\sigma^2} $$
ここで $\bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i$ は標本平均です。
フィッシャー情報量の計算
対数尤度の二次微分は次のようになります。
$$ \frac{\partial^2 \ell}{\partial \mu^2} = -\frac{n}{\sigma^2} $$
これはデータに依存しない定数なので、期待値を取っても同じ値です。したがって、フィッシャー情報量は次のようになります。
$$ I_n(\mu) = -E\left[\frac{\partial^2 \ell}{\partial \mu^2}\right] = \frac{n}{\sigma^2} $$
スコア検定統計量
帰無仮説のもとでスコア関数を評価します。
$$ S_n(\mu_0) = \frac{n(\bar{x} – \mu_0)}{\sigma^2} $$
スコア検定統計量は次のようになります。
$$ T_S = \frac{S_n(\mu_0)^2}{I_n(\mu_0)} = \frac{\left[\frac{n(\bar{x} – \mu_0)}{\sigma^2}\right]^2}{\frac{n}{\sigma^2}} $$
分子を展開して整理します。分子は $\frac{n^2(\bar{x} – \mu_0)^2}{\sigma^4}$、分母は $\frac{n}{\sigma^2}$ なので、
$$ T_S = \frac{n^2(\bar{x} – \mu_0)^2}{\sigma^4} \cdot \frac{\sigma^2}{n} = \frac{n(\bar{x} – \mu_0)^2}{\sigma^2} $$
つまり、
$$ \begin{equation} T_S = \frac{n(\bar{x} – \mu_0)^2}{\sigma^2} = \left(\frac{\bar{x} – \mu_0}{\sigma/\sqrt{n}}\right)^2 = z^2 \end{equation} $$
これは、標準的な $z$ 検定統計量の2乗に他なりません。正規分布の平均の検定という最もシンプルな場合では、スコア検定は $z$ 検定と完全に一致します。
この例では、対数尤度関数が $\mu$ の完全な二次関数であるため、スコア検定、尤度比検定、ワルド検定の3つすべてが同じ検定統計量を与えます。
正規分布の例は教育的でしたが、やや単純すぎました。次に、スコア検定が真に力を発揮するポアソン分布の例を見てみましょう。
具体例: ポアソン分布のスコア検定
問題設定
ポアソン分布 $\text{Poisson}(\lambda)$ の母数 $\lambda$ に対する検定を考えます。
$$ H_0: \lambda = \lambda_0 \quad \text{vs.} \quad H_1: \lambda \neq \lambda_0 $$
データ $x_1, \dots, x_n \overset{\text{i.i.d.}}{\sim} \text{Poisson}(\lambda)$ が与えられています。
導出
ポアソン分布の確率質量関数は $P(X = x) = \frac{\lambda^x e^{-\lambda}}{x!}$ であり、対数尤度は次のようになります。
$$ \ell(\lambda) = \sum_{i=1}^{n}\left[x_i \ln \lambda – \lambda – \ln(x_i!)\right] = n\bar{x}\ln\lambda – n\lambda – \sum_{i=1}^{n}\ln(x_i!) $$
スコア関数を求めるために $\lambda$ で微分します。
$$ S_n(\lambda) = \frac{n\bar{x}}{\lambda} – n = \frac{n(\bar{x} – \lambda)}{\lambda} $$
フィッシャー情報量は、対数尤度の二次微分の負の期待値です。
$$ \frac{\partial^2 \ell}{\partial \lambda^2} = -\frac{n\bar{x}}{\lambda^2} $$
帰無仮説のもとで $E[\bar{x}] = \lambda_0$ なので、
$$ I_n(\lambda_0) = -E\left[\frac{\partial^2 \ell}{\partial \lambda^2}\bigg|_{\lambda=\lambda_0}\right] = \frac{n}{\lambda_0} $$
スコア検定統計量は次のようになります。
$$ T_S = \frac{S_n(\lambda_0)^2}{I_n(\lambda_0)} = \frac{\left[\frac{n(\bar{x} – \lambda_0)}{\lambda_0}\right]^2}{\frac{n}{\lambda_0}} = \frac{n(\bar{x} – \lambda_0)^2}{\lambda_0} $$
この結果は、ポアソン分布の場合の $\chi^2$ 型のスコア検定統計量になっています。$\bar{x}$ が $\lambda_0$ から離れるほど $T_S$ が大きくなり、帰無仮説を棄却する方向に作用します。
ポアソン分布の場合、対数尤度関数は $\lambda$ の二次関数ではないため、スコア検定、尤度比検定、ワルド検定は異なる値を返します。しかし、$n$ が大きければ3つの値は近づきます。
理論的な理解が深まったところで、Pythonを使って三つの検定統計量を実際に計算し、比較してみましょう。
Pythonによる実装と数値シミュレーション
三大漸近検定の比較
ポアソン分布の例について、スコア検定・尤度比検定・ワルド検定の3つの検定統計量を計算し、それらの分布を帰無仮説のもとでシミュレーションします。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
np.random.seed(42)
# --- パラメータ設定 ---
lambda_0 = 5.0 # 帰無仮説のもとでのλ
n = 20 # サンプルサイズ
n_sim = 10000 # シミュレーション回数
alpha_level = 0.05 # 有意水準
# --- シミュレーション ---
T_score = np.zeros(n_sim)
T_lr = np.zeros(n_sim)
T_wald = np.zeros(n_sim)
for i in range(n_sim):
# 帰無仮説のもとでデータを生成
x = np.random.poisson(lambda_0, size=n)
x_bar = np.mean(x)
# スコア検定統計量
T_score[i] = n * (x_bar - lambda_0)**2 / lambda_0
# 尤度比検定統計量
if x_bar > 0:
T_lr[i] = 2 * n * (x_bar * np.log(x_bar / lambda_0) - (x_bar - lambda_0))
else:
T_lr[i] = 2 * n * lambda_0
# ワルド検定統計量
if x_bar > 0:
T_wald[i] = n * (x_bar - lambda_0)**2 / x_bar
else:
T_wald[i] = n * lambda_0
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# χ²(1)分布の理論曲線
chi2_x = np.linspace(0, 15, 200)
chi2_pdf = stats.chi2.pdf(chi2_x, df=1)
titles = ["Score test", "Likelihood ratio test", "Wald test"]
data_list = [T_score, T_lr, T_wald]
for ax, T, title in zip(axes, data_list, titles):
ax.hist(T, bins=60, density=True, alpha=0.6, color="steelblue",
edgecolor="white", label="Simulated")
ax.plot(chi2_x, chi2_pdf, "r-", linewidth=2, label=r"$\chi^2(1)$")
# 棄却率の計算
reject_rate = np.mean(T > stats.chi2.ppf(1 - alpha_level, df=1))
ax.set_title(f"{title}\nReject rate = {reject_rate:.3f} "
f"(nominal = {alpha_level})", fontsize=11)
ax.set_xlabel("Test statistic", fontsize=11)
ax.set_ylabel("Density", fontsize=11)
ax.set_xlim(0, 15)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("score_test_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
上のヒストグラムから、以下のことが読み取れます。
-
3つの検定統計量はいずれも $\chi^2(1)$ 分布によく近似されている。シミュレーションによるヒストグラム（青）と理論的な $\chi^2(1)$ 分布の密度関数（赤線）がよく一致しています。これは、$n = 20$ の標本サイズでも漸近近似がかなり良いことを示しています
-
棄却率は名目水準 $0.05$ にほぼ一致している。3つの検定とも第1種の過誤率が $5\%$ に近い値を示しています。もし棄却率が名目水準から大きくずれていれば、漸近近似が不十分であることを意味しますが、ポアソン分布の場合は $n = 20$ で十分です
-
3つの分布は似ているが完全には一致しない。有限標本では検定統計量の値が異なるため、特に分布の裾の部分で微妙な違いがあります。これは先ほどの理論的な議論と整合しています
検出力の比較
次に、対立仮説のもとでの検出力（帰無仮説を正しく棄却する確率）を3つの検定で比較します。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
np.random.seed(42)
lambda_0 = 5.0
n = 20
n_sim = 10000
alpha_level = 0.05
chi2_crit = stats.chi2.ppf(1 - alpha_level, df=1)
# 真のλの範囲（対立仮説）
lambda_true_list = np.linspace(3.0, 8.0, 30)
power_score = []
power_lr = []
power_wald = []
for lambda_true in lambda_true_list:
reject_s = 0
reject_lr = 0
reject_w = 0
for _ in range(n_sim):
x = np.random.poisson(lambda_true, size=n)
x_bar = np.mean(x)
# スコア検定
Ts = n * (x_bar - lambda_0)**2 / lambda_0
if Ts > chi2_crit:
reject_s += 1
# 尤度比検定
if x_bar > 0:
Tlr = 2 * n * (x_bar * np.log(x_bar / lambda_0) - (x_bar - lambda_0))
else:
Tlr = 2 * n * lambda_0
if Tlr > chi2_crit:
reject_lr += 1
# ワルド検定
if x_bar > 0:
Tw = n * (x_bar - lambda_0)**2 / x_bar
else:
Tw = n * lambda_0
if Tw > chi2_crit:
reject_w += 1
power_score.append(reject_s / n_sim)
power_lr.append(reject_lr / n_sim)
power_wald.append(reject_w / n_sim)
fig, ax = plt.subplots(figsize=(9, 6))
ax.plot(lambda_true_list, power_score, "b-o", markersize=4, linewidth=1.5,
label="Score test")
ax.plot(lambda_true_list, power_lr, "r-s", markersize=4, linewidth=1.5,
label="LR test")
ax.plot(lambda_true_list, power_wald, "g-^", markersize=4, linewidth=1.5,
label="Wald test")
ax.axhline(alpha_level, color="gray", linestyle="--", linewidth=0.8,
label=f"Significance level = {alpha_level}")
ax.axvline(lambda_0, color="gray", linestyle=":", linewidth=0.8)
ax.set_xlabel(r"True $\lambda$", fontsize=12)
ax.set_ylabel("Power (rejection rate)", fontsize=12)
ax.set_title(f"Power comparison (n={n}, " + r"$\lambda_0$" + f"={lambda_0})",
fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)
plt.tight_layout()
plt.savefig("score_test_power.png", dpi=150, bbox_inches="tight")
plt.show()
検出力の比較グラフから、以下の重要な知見が得られます。
-
帰無仮説が正しい $\lambda = 5.0$ では、3つの検定とも棄却率が $0.05$ に近い。灰色の点線（$\lambda = 5.0$）で3つの曲線が水平の破線（$\alpha = 0.05$）に交わっていることが確認できます
-
対立仮説のもとでは、3つの検定の検出力はほぼ同等。特に、帰無仮説の値 $\lambda_0 = 5.0$ から十分離れた領域では、3つの曲線がほぼ重なっています。これは漸近的同値性の数値的な確認です
-
$\lambda_0$ の近傍では微妙な差がある。帰無仮説の値の近く（検出力が低い領域）で3つの検定にわずかな差が見られます。尤度比検定が最もバランスの取れた性能を持ち、ワルド検定はやや積極的（棄却しやすい）、スコア検定はやや保守的（棄却しにくい）な傾向があります
ロジスティック回帰でのスコア検定
最後に、実用的な例として、ロジスティック回帰における共変量の有意性をスコア検定で判定します。
import numpy as np
from scipy import stats
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
np.random.seed(42)
# --- データ生成 ---
n = 200
x1 = np.random.randn(n) # 第1共変量（有意）
x2 = np.random.randn(n) # 第2共変量（検定対象）
beta_true = np.array([0.5, 1.0, 0.3]) # [切片, β1, β2]
# ロジスティック回帰の真のモデル
eta = beta_true[0] + beta_true[1] * x1 + beta_true[2] * x2
prob = 1 / (1 + np.exp(-eta))
y = np.random.binomial(1, prob)
# --- 帰無モデル（x2を含まない）の最尤推定 ---
# H0: β2 = 0 のもとでβ0, β1を推定
def neg_loglik_null(params):
b0, b1 = params
eta_null = b0 + b1 * x1
p = 1 / (1 + np.exp(-eta_null))
p = np.clip(p, 1e-10, 1 - 1e-10)
return -np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))
from scipy.optimize import minimize
res_null = minimize(neg_loglik_null, x0=[0.0, 0.0], method="BFGS")
b0_hat, b1_hat = res_null.x
# --- 帰無仮説のもとでのスコア ---
eta_null = b0_hat + b1_hat * x1
p_null = 1 / (1 + np.exp(-eta_null))
# β2に関するスコア（x2の係数に対応）
score_beta2 = np.sum((y - p_null) * x2)
# β2に関する有効フィッシャー情報量
# I_22 - I_21 @ I_11^{-1} @ I_12
W = p_null * (1 - p_null) # 重み
X_null = np.column_stack([np.ones(n), x1]) # 帰無モデルの計画行列
# 情報行列のブロック
I_11 = X_null.T @ np.diag(W) @ X_null # 2x2
I_12 = X_null.T @ (W * x2) # 2x1
I_22 = np.sum(W * x2**2) # scalar
# シューア補行列
I_eff = I_22 - I_12 @ np.linalg.solve(I_11, I_12)
# スコア検定統計量
T_score = score_beta2**2 / I_eff
p_value = 1 - stats.chi2.cdf(T_score, df=1)
print("=" * 50)
print("ロジスティック回帰におけるスコア検定")
print("=" * 50)
print(f"帰無仮説: β2 = 0 (x2の効果なし)")
print(f"帰無モデル推定: β0 = {b0_hat:.4f}, β1 = {b1_hat:.4f}")
print(f"スコア S(β2) = {score_beta2:.4f}")
print(f"有効情報量 I_eff = {I_eff:.4f}")
print(f"スコア検定統計量 T_S = {T_score:.4f}")
print(f"p値 = {p_value:.6f}")
print(f"有意水準5%で: {'棄却' if p_value < 0.05 else '棄却しない'}")
print(f"(真のβ2 = {beta_true[2]})")
この実装結果から、以下のことが確認できます。
-
スコア検定は帰無モデルの推定のみで実行できている。フルモデル（$x_1$ と $x_2$ の両方を含むモデル）を一度も推定していません。帰無モデル（$x_2$ を含まない）の推定結果を使って、$x_2$ の係数 $\beta_2$ に関するスコアと情報量を計算しています
-
真の $\beta_2 = 0.3$ に対して、スコア検定は帰無仮説（$\beta_2 = 0$）を棄却する方向に作用する。$n = 200$ の標本サイズで、$\beta_2 = 0.3$ の効果を検出できるかどうかは、効果の大きさとノイズの比率に依存します
-
有効情報量はシューア補行列として計算されている。$\beta_0$ と $\beta_1$ の影響を除去した後の $\beta_2$ に関する情報量は、情報行列のブロック構造を用いて正確に計算されています
スコア検定の発展的な話題
ラグランジュ乗数検定との同値性
計量経済学では、スコア検定はラグランジュ乗数検定（Lagrange Multiplier test, LM検定）と呼ばれています。この名前は、帰無仮説が「等式制約付き最適化」の制約として表現できることに由来します。
帰無仮説 $\theta = \theta_0$ のもとでの対数尤度の最大化は、制約 $\theta = \theta_0$ のもとでの $\ell(\theta)$ の最大化と同じです。ラグランジュ乗数法を使うと、制約付き最適化の停留条件は次のようになります。
$$ \frac{\partial \ell}{\partial \theta}\bigg|_{\theta=\theta_0} = \lambda \cdot \frac{\partial g}{\partial \theta} $$
ここで $g(\theta) = \theta – \theta_0 = 0$ が制約であり、$\lambda$ がラグランジュ乗数です。$\frac{\partial g}{\partial \theta} = 1$ なので、ラグランジュ乗数は $\lambda = S(\theta_0)$、すなわちスコア関数の値そのものです。
したがって、「ラグランジュ乗数が大きい ⟺ 制約が対数尤度を大きく損なっている ⟺ 帰無仮説は正しくない」という論理でスコア検定が構成されます。スコア検定とLM検定は全く同じ検定統計量を持ち、名前だけが異なります。
ロバストスコア検定
標準的なスコア検定は、データが仮定された確率分布に正しく従うことを前提としています。しかし、実際のデータでは分布の仮定が完全に正しいとは限りません。このようなモデルの誤特定（model misspecification）に対して頑健な変種がロバストスコア検定です。
ロバストスコア検定では、フィッシャー情報量の代わりに、スコア関数のサンドイッチ型の分散推定量を用います。
$$ \begin{equation} T_{RS} = \bm{S}_n^T \left(\bm{V}_S\right)^{-1} \bm{S}_n \end{equation} $$
ここで $\bm{V}_S$ はスコアベクトルのロバスト分散推定量です。この検定は、分布の仮定が誤っていても漸近的に正しいサイズを維持します。
スコア検定と情報幾何学
より深い理論的な観点として、スコア検定は情報幾何学（information geometry）の枠組みで自然に理解できます。パラメータ空間上にフィッシャー計量を導入すると、スコア検定統計量は帰無仮説の点 $\theta_0$ における対数尤度のリーマン勾配の大きさの2乗に対応します。
この幾何学的視点は、三大漸近検定を統一的に理解する上で非常に示唆的です。パラメータ空間の曲がり（曲率）がゼロの場合（指数型分布族のアフィンパラメータ化）、三つの検定は完全に一致します。曲率がゼロでない場合に初めて三者に差が生じるのです。
まとめ
本記事では、スコア検定（ラオ検定）の理論を直感的な理解から厳密な導出まで解説しました。
- スコア関数は対数尤度の勾配であり、最尤推定量ではゼロになります。帰無仮説の値でスコア関数がゼロから大きくずれていれば、帰無仮説は正しくないと判断できます
- スコア検定統計量 $T_S = S(\theta_0)^2 / I(\theta_0)$ は、帰無仮説のもとで漸近的に $\chi^2$ 分布に従います。この検定統計量はスコア関数の分散（フィッシャー情報量）で標準化することで構成されます
- 三大漸近検定（スコア検定・尤度比検定・ワルド検定）は対数尤度の異なる幾何学的量（傾き・高さの差・水平距離）を測定しますが、大標本のもとでは漸近的に同値です
- スコア検定の最大の利点は、帰無仮説のもとでの推定のみで計算が完結することです。対立仮説のもとでの最尤推定が困難な場合に特に有用です
- 多次元パラメータへの拡張では、スコアベクトルと情報行列の二次形式として検定統計量が構成されます。局外パラメータの影響はシューア補行列を通じて除去されます
- Pythonでのシミュレーションにより、3つの検定統計量の $\chi^2$ 分布への適合と検出力の漸近的同値性を数値的に確認しました
スコア検定は、統計的推測の理論体系における重要な柱であるとともに、ゲノミクスや大規模データ分析など現代の応用でも不可欠なツールです。
次のステップとして、以下の記事も参考にしてください。