【汎用性抜群】尤度比検定をわかりやすく解説します
尤度比検定とは、汎用性の高い統計モデルの検定です。
その汎用性の高さは、サンプル数が十分大きい時には、尤度比検定統計量の対数に2をかけたものがカイ２乗分布に従う性質にあります。
python中心に解説したコンテンツは以下になります。
【python】尤度比検定で統計モデルの比較をしよう|統計的仮説検定
尤度比検定(likelihood ratio test)
まず、専門的な用語抜きに説明すると、尤度比検定とは二つのモデルのうち、観測データをよりよく説明するのはどちらだろうか、という問いに答えるための統計的な手法です。
ここで一方のモデルの尤もらしさをもう一方のモデルの尤もらしさで割ってあげることで、統計量を作ります。
ここの「尤もらしさ」が「最大尤度」になります。
尤度についておさらいしたい方はこちらをご覧ください。
「大学卒業経験と就業年数は、個人の年収とどのような関係があるのか」
以下のようなデータを扱います。被験者は100人いるとします。大卒50人と高卒50人を無作為に選びました。
x：個人iの就業年数。整数型(int)
y：個人iの直近の年収。単位は円。整数型(int)
c：個人iが大卒資格があるかどうか。Yなら大卒、Nなら高卒です。因子型(factor)
データを見てみましょう。
data <- read.csv("income.csv")
data
data$x
data$y
data$c
前回のコンテンツで説明した通り、データ名$カラム名で、列のデータだけ出力できます。
前回の記事：【GLM】一般化線形モデルを解説(統計モデル編)
このデータを使って、大卒経験と就業年数が、個人iの年収にどのように関係してくるのかを分析します。
ちなみに、高卒50人を「コントロール群」と呼びます。
尤度比検定の準備
今回は、以下のようなモデルを比較したいと思います。今回は、ある個人\(i\)においての年収が\(y_i\)である確率\(p(yi | λ_i)\)は、ポアソン分布に従うとします。
今回は、一般化線形モデル(GLM)を扱うので事前に勉強したい方はは以下のコンテンツを見ておくことをお勧めします。
モデル１：定数項と就業年数\(x_i\)で、yを表したモデル。パラメーター数は2。
$$logλ_{i} =β_{1}+β_{2}x_{i}$$
このモデルでは、大卒経験は無視しています。
ちなみに、\(λ_i\)は個人iの平均年収です。
$$λ_{i} = exp(β_{1}+β_{2}x_{i})$$
まず、帰無仮説として「大卒資格の有無は個人の年収に関係がない」を設定します。
数式に表すと、
$$H0:β_{3}=0$$
今回棄却したい(否定したい)仮説は、モデル1です。
モデル２：定数項と就業年数xiと大卒経験有無cで、yを表したモデル。パラメーター数は3。
$$logλ_{i}=β_{1}+β_{2}x_{i}+β_{3}c$$
このモデルは、大卒資格の有無を考慮しています。
帰無仮説であるモデル1に対して、モデル２を対立仮説と言います。
cの係数である\(β_3\)は0ではないので、大卒資格は+であれ-であれ年収と連動があるということです。
*ちなみに、ある係数を0にするとモデル同士が一致するモデルをネストしているといいます。
尤度比検定はネストしているモデル同士を比較検討する統計モデル検定の一つです。
まず、モデル1(帰無仮説)とモデル２(対立仮説)の最大対数尤度を求めてみましょう。
CODE
ひとまずグラフを書いてみましょう。
as.factorで因子型に変換しています。
plot(data$x, data$y, pch = c(21,19)[as.factor(data$c)])
legend("topleft", legend = c("N", "Y"),pch = c(21,19))
ここでは、因子である大卒資格データもグラフに含めています。
大卒資格のない白丸のいくつかは、かなり年収が低い印象を受けます。
【モデル1(帰無仮説)の推定】
fit <- glm(y ~ x, data = data, family = poisson(link = log))
fit
logLik(fit)
glm()でGLMによる回帰分析を行います。今回はx+cとせずに、xのみです。
interceptは定数項のことです。推定結果は以下のようになります。
$$logλ_{i}=6.75917+0.01111x_{i}$$
> logLik(fit)
'log Lik.' -2445.33 (df=2)最大対数尤度は-2445.33になりました。
【モデル2(対立仮説)の推定】
fit.college <- glm(y ~ x+c, data = data, family = poisson(link = log))
fit.college
logLik(fit.college)
モデル1のコードと比べると、x+cになっています。
出力結果は以下の通りです。
cの係数は0.05966になりました。式にすると以下の通りになります。
正の数なので、「大卒資格は年収にプラスの影響がある」ということになります。
$$logλ_{i}=6.70690+0.01395x_{i}+0.05966c_{i}$$
> logLik(fit.college)
'log Lik.' -2406.12 (df=3)対立仮説についても、最大対数尤度を見ます。-2406.12でした。
尤度比検定統計量について
一般には尤度比検定統計量の帰無仮説のもとでの確率分布は分かりません。
そこで尤度比検定は、帰無仮説の最大対数尤度から対立仮説の最大対数尤度の差に２をかけたものを検定統計量として使います。
$$D=2L=2(logL_{1}^*-logL_{2}^*)$$
また、最大尤度の比の対数に2をかけたものも同値です。
対数なので、割り算は引き算になります。
これは、【良いモデルとは】AIC(赤池情報量基準)について。で登場した「逸脱度」です。
$$2log\frac{対立仮説の最大尤度}{帰無仮説の最大尤度}$$
そして、この検定統計量は自由度1のカイ二乗分布に従います。
自由度は帰無仮説の制約数によります。
つまり、最尤推定するパラメータの数ですね。
$$H0:β_{3}=0$$
今回は、β3の値が制約されているだけなので1です。β3=β4=0なら自由度2です。
-2*(logLik(fit)-logLik(fit.college))
'log Lik.' 78.4204実際に計算してみると、検定統計量は78.4204でした。
実際にカイ二乗分布表を見てみましょう。
有意になっておりますね。
モデルに入れる説明変数の取捨選択に役に立ちました。
→因果があるまでは言えません。
有意やp値については、以下のコンテンツに議論があります。
【仮説検定】p値をゼロから解説(第一種の過誤,第二種の過誤,検出力)
| 自由度m α | 0.99 | 0.975 | 0.95 | 0.05 | 0.025 | 0.01 |
| 1 | 0.00016 | 0.00098 | 0.0039 | 3.84 | 5.02 | 6.63 |
| 2 | 0.020 | 0.051 | 0.10 | 5.99 | 7.38 | 9.21 |
| 3 | 0.11 | 0.22 | 0.35 | 7.81 | 9.35 | 11.34 |
| 4 | 0.30 | 0.48 | 0.71 | 9.49 | 11.14 | 13.28 |
| 5 | 0.55 | 0.83 | 1.15 | 11.07 | 12.83 | 15.09 |
カイ二乗分布
カイ二乗分布とは、確率変数Z_nが標準正規分布に従うとき、各確率変数の平方和が従う確率分布のことです。
尤度比検定の漸近理論で使います。
$$Z_{1},…,Z_{n}~N(0,1)$$
$$Z_{1}^2+…+Z_{n}^2~χ^2(d)$$
統計的仮説検定に関して興味がある方は以下のコンテンツをご覧ください。
【仮説検定】p値をゼロから解説(第一種の過誤,第二種の過誤,検出力)
