【汎用性抜群】尤度比検定をわかりやすく解説します

尤度比検定(likelihood ratio test)

尤度比検定とは、汎用性の高い統計モデルの検定です。

専門的な用語抜きに説明すると、尤度比検定とは二つのモデルのうち、観測データをよりよく説明するのはどちらだろうか、という問いに答えるための統計的な手法です。

ここで一方のモデルの尤もらしさをもう一方のモデルの尤もらしさで割ってあげることで、統計量を作ります。

統計的仮説検定に関して興味がある方は青の統計学の以下のコンテンツをご覧ください。

尤度とは

尤度(likelihood)は、統計学においてデータがある確率モデルに従って生成されたと仮定した(確率モデルのパラメーターが与えられた)ときの、そのデータが観測される確率を表します。

例えば、平均値 ${\mu}$、標準偏差 ${\sigma}$が与えられた正規分布に従ってデータが生成されたと仮定した場合、尤度は、次の式で計算されます

$$L(\mu,\sigma)= \frac{1}{\sigma \sqrt{2\pi}^2}exp(- \frac{(x-\mu)^2}{2\sigma^2})$$

ここで、${x}$は観測されたデータ、${n}$は観測されたデータの数です。

もっと具体に言えば、「コインを${100}$回投げて、${50}$回表になる時の確率」を尤度${L(p)}$、コインを一回投げて表になる確率のことをパラメーター${p}$と当てはめています。

最尤法とはこの尤度が最大になるためのパラメータを推定量として求めようという手法です。これは、観測された事象は、「起きやすい確率が最大のものが起きている」という考えから派生しています。

もっと尤度についておさらいしたい方は、青の統計学のこちらの記事をご覧ください。

【尤度って？】尤度関数と最尤推定量の解説と例題

繰り返しになりますが、尤度比検定（Likelihood ratio test, LRT）は統計学において、2つの統計モデルの適合度を比較するために使用される手法です。

モデルの設定

この検定は、ネストされたモデル（すなわち、一方のモデルが他方のモデルに含まれるパラメータの制約版である場合）の適合度を比較する際によく使われます。

では基本的なアプローチを説明します。

まず、完全モデル（制約なしモデル）と縮小モデル（制約ありモデル）を設定します。

完全モデルはパラメータがより多く、一方の縮小モデルはそのパラメータの一部に制約（例えば、あるパラメータが0であるなど）が課されます。

なので、以下のように帰無仮説を設定します

• 帰無仮説H₀: 縮小モデルが適切である（縮小モデルで十分）
• 対立仮説H₁: 完全モデルが必要である（縮小モデルでは不十分）

完全モデル

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

・\(\mathbf{y}\)は\(n×1\)の目的変数ベクトルです。

・\(\mathbf{X}\)は、\(n×(p+1)\)の設計行列です。ちなみに、説明変数の値を含んでします。1列目はすべて1でバイアス項\(\beta_0\)のためのものです、

・\(\beta\)はパラメータベクトルで、\((p＋1)×1\)です。

・\(\epsilon\)は\(n×1\)の誤差項ベクトルです。

縮小モデル

もし \(\beta_2​=0\) という制約を設ける場合、このパラメータをモデルから除外します。

行列 \(\mathbf{X}\) から該当する列 \(x_2\) を削除し、パラメータベクトル\(\beta\)からも\(\beta_2\) を除きます。この場合のモデルは次のように表されます

$$\mathbf{y}=\mathbf{X}_{-2}\beta_{-2}+\epsilon$$

・\(\mathbf{X}_{-2}\)は、\(\mathbf{X}\)から\(\beta_2\)に対応する列を除いた\(n×p\)の設計行列です。

・\(\beta_{-2}\)は\(\beta_2\)を除いた\(p×1\)のパラメータベクトルです。

応用例 回帰分析での変数選択を例に取ると

縮小モデル：y = β₀ + β₁x₁
完全モデル：y = β₀ + β₁x₁ + β₂x₂ + β₃x₃
この場合、df = 2（β₂とβ₃の2つのパラメータ差）となります。

尤度比検定統計量の作成

尤度比検定では、これら2つのモデルの尤度を比較します。

尤度比 \(\Lambda\) は次のように定義されます

$$\Lambda=\frac{L(\beta_{-2}|\mathbf{X}_2,\mathbf{y})}{L(\beta|\mathbf{X},\mathbf{y})}$$

分子が縮小モデルの尤度で、分母が完全モデルの尤度ですね。

さて、尤度比から尤度比検定統計量を作りましょう。

$${\lambda=-2log\Lambda=-2log\frac{L(\beta_{-2}|\mathbf{X}_2,\mathbf{y})}{L(\beta|\mathbf{X},\mathbf{y})}}$$

尤度比検定量\(\lambda\)は、大標本（大きなサンプルサイズ）の下で、自由度が完全モデルと縮小モデル間のパラメータ数の差に等しいカイ二乗分布に従います。

λ = -2(対数尤度[縮小モデル] - 対数尤度[完全モデル])

この性質は大標本の理論の一部ウィルクスの定理に基づいています。

$${\lambda\sim \chi^2_{df}}$$

ここで、\(df\)は自由度で、完全モデルと縮小モデルの間のパラメータ数の差です。

たとえば、完全モデルが5個のパラメータを持ち、縮小モデルが3個のパラメータを持つ場合、自由度は

${df=5−3=2}$

この性質により、${\lambda}$ を ${\chi^2}$分布の臨界値と比較することで、帰無仮説の棄却の可否を判断できます。

さて、次は有意水準${\alpha}$(検定で誤って帰無仮説を棄却する確率)を設定します。

一般的には \alpha = 0.05が使われます。

補足｜カイ二乗分布

さてカイ二乗分布とは、確率変数\(Z_n\)が標準正規分布に従うとき、各確率変数の平方和が従う確率分布のことです。

尤度比検定の漸近理論で使います。

尤度比検定の汎用性の高さは、サンプル数が十分大きい時には、尤度比検定統計量の対数に2をかけたものがカイ２乗分布に従う性質にあります。

$$Z_{1},…,Z_{n}\sim N(0,1)$$

$$Z_{1}^2+…+Z_{n}^2 \sim\chi^2(d)$$

また、カイ二乗分布は以下のような確率密度関数に従います。

$$f(x;n)=\frac{1}{2^{\frac{n}{2}}}x^{\frac{n}{2}-1}exp(-\frac{x}{2})$$

ここで、 \(x\) は自由度 \(n\)のカイ二乗分布に従う乱数、\(n\)は自由度です。

適合度検定における\(\chi^2\)統計量を見てみましょう。

$${\chi^2=\sum\frac{(O_i-E_i)^2}{E_i}}$$

統計量全体がカイ二乗分布に従うのは、独立した正規分布の二乗和として表されるからです。

詳しくは、青の統計学のこちらの記事をご覧ください。

補足｜カイ二乗分布とガンマ分布

補足ですが、カイ二乗分布はガンマ分布の特別な形としても理解することができます。

ガンマ分布は形状パラメータ\(\frac{k}{2}\)と尺度パラメータ\(2\)のガンマ分布がカイ二乗分布に一致します。

この関係から、カイ二乗分布の確率密度関数は以下のように表されます。

$$f(x;k)=\frac{1}{2^{\frac{k}{2}Γ(\frac{k}{2})}}x^{\frac{k}{2}-1}e^{-\frac{x}{2}}$$

ここで \(x≥0\) であり、\(Γ\)はガンマ関数です。

尤度比の使い道、KL divergenceについて

尤度比の使い方として、KL divergenceをご紹介します。

二つの確率分布\(q(\textbf{x})\)と\(p(\textbf{x})\)がどれくらい異なるかを測るための非対称な距離尺度です。

これは情報理論の中で一般的に使われ、特に確率分布や統計モデル間の相違度を測定するのに便利です。

$$KL[q(\textbf{x})|p(\textbf{x})]=\int q(\textbf{x})log\frac{q(\textbf{x})}{p(\textbf{x})}d\textbf{x}$$

上のように尤度比の対数の期待値をとっていることがわかります。

KLダイバージェンスがゼロ、つまり最小値になるのはPとQが同じ分布である場合のみです。

ただ、二つの確率分布を並列に扱っているわけではなく、一つの確率分布(ここでは\(p(x)\))を真の確率分布と仮定した場合に、別の確率分布を(ここでは\(q(x)\))用いて情報を符号化するときの“情報損失”を測定します。

KLダイバージェンスは尤度比の対数を一種の”期待値”として捉えることができます。

つまり、Pが真の分布であるとき、それに基づいて生成される各データ点において、Pに基づく尤度の対数とQに基づく尤度の対数との差（これが尤度比の対数に相当します）の”Pによる期待値”がKLダイバージェンスとなります。

そのため、KLダイバージェンスは統計的推論における尤度比の概念を、全確率分布に対する尺度として一般化したものと考えることができます。

詳しくは以下のコンテンツをご覧ください。

【AICで使う】KL divergence(カルバック-ライブラー情報量)をわかりやすく解説|python

CODE｜一般化線形モデル（GLM）における尤度比検定

ちょっとした検定を行います。

ポアソン回帰モデルを用いて、ある説明変数が応答変数に対して有意な影響を持つかどうかを尤度比検定で評価します。

ポアソン回帰モデルは、一般化線形モデルの一つです。

詳しくはこちらのコンテンツをご覧ください。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import chi2

# データ生成
np.random.seed(0)
data = np.random.poisson(lam=30, size=1000)
covariate = np.random.normal(size=1000)
data_frame = pd.DataFrame({'Response': data, 'Covariate': covariate})

# 完全モデル（Covariateを含む）
full_model = smf.glm('Response ~ Covariate', data=data_frame, family=sm.families.Poisson()).fit()

# 縮小モデル（Covariateを除外）
reduced_model = smf.glm('Response ~ 1', data=data_frame, family=sm.families.Poisson()).fit()

# 尤度比検定
lr_stat = reduced_model.deviance - full_model.deviance
lr_df = full_model.df_model - reduced_model.df_model
lr_pvalue = chi2.sf(lr_stat, lr_df)  # カイ二乗分布からのp値を計算

print(f"LR statistic: {lr_stat}")
print(f"p-value: {lr_pvalue}")
print(f"Degrees of freedom: {lr_df}")

LR statistic: 1.9353660252947975

p-value: 0.16417256842907785

Degrees of freedom: 1

次元を一つ落としただけなので、自由度は1です。

p値の部分を補足をすると、カイ二乗分布の生存関数（1 – CDF）を計算しています。

これは、検定統計量がカイ二乗分布で観察された値よりも大きくなる確率を返し、これによってp値が得られます。

有意水準とp値について復習したい方はこちらをご覧ください。

FOLLOW ME ！

@blue_statistics