chapter05 GLMの尤度比検定と検定の非対称性
Ei Taguchi
本章では尤度比検定(likelyhood ratio test)について説明
- ネストしているモデル間の比較に使用
5.1 統計学的な検定のわくぐみ
- 検定とモデル選択の共通点
- モデル設計→パラメータの最尤推定計算
- 統計モデルの検定(違うところ) = Neyman-Pearson検定
- 帰無仮説の危険率を評価(有意水準を設定)
- 帰無仮説が棄却できるか判断
- AICによるモデル選択(違うところ)
- AICの評価
- 予測の良いモデルを選ぶ
5.2 尤度比検定の例題 : 逸脱度の差を調べる
- 第3章の種子数のデータを使用
#各自データを読み込んでください
head(d)## # A tibble: 6 x 3
## y x f
## <dbl> <dbl> <chr>
## 1 6 8.31 C
## 2 6 9.44 C
## 3 6 9.5 C
## 4 12 9.07 C
## 5 10 10.2 C
## 6 4 8.32 C- 一定モデル : 種子数の平均が体サイズxに依存しないモデル(傾き \(= 0\))
fit.null <-glm(y ~ 1, data = d, family = poisson)- xモデル : 種子数の平均が体サイズxに依存するモデル(傾き \(\neq 0\))
fit.x <- glm(y ~ x, data = d, family = poisson(link = "log"))stargazer(fit.x, fit.null, type = "text", style = "all", ci = TRUE,
star.cutoffs = NA, omit.table.layout = 'n',
align = TRUE)##
## =========================================================
## Dependent variable:
## ---------------------------------
## y
## (1) (2)
## ---------------------------------------------------------
## x 0.076
## (0.006, 0.145)
## t = 2.125
## p = 0.034
## Constant 1.292 2.058
## (0.579, 2.005) (1.988, 2.128)
## t = 3.552 t = 57.586
## p = 0.0004 p = 0.000
## ---------------------------------------------------------
## Observations 100 100
## Log Likelihood -235.386 -237.643
## Akaike Inf. Crit. 474.773 477.286
## Residual Deviance 84.993 (df = 98) 89.507 (df = 99)
## Null Deviance (df = 99) 89.507 89.507
## =========================================================2つのモデルを比較すると, 逸脱度(Residual Deviance)の差は\(4.5\)くらい
- しかし, パラメータ数が多い方が逸脱度は小さくなる
尤度比 : \[\begin{align} \frac{L_1}{L_2} = \frac{一定モデルの最大尤度 : exp(-237.6)}{xモデルの最大尤度 : exp(-235.4)}\\ \end{align}\]
- このままでは使わず\(log\)をとり\(-2\)をかける
逸脱度の差 : \[\begin{align} Δ_{1, 2} & = -2 × log\frac{L_1}{L_2}\\ &= -2 × (logL_1 - logL_2)\\ \end{align}\]
\(D_1 = -2logL_1\), \(D_2 = -2logL_2\)とすると, \[\begin{align} ΔD_{1, 2} = D_1 - D_2 \\ \end{align}\] となる
今回, \(ΔD_{1, 2} = 4.5\) となった
- この値がモデルのあてはまりの良さが改善されたと判断していい値かどうかを調べる
5.3 2種類の過誤と統計学的な検定の非対称性
本例題での帰無仮説と対立仮説
- 帰無仮説 : 一定モデル(パラメータ数 : \(k=1\), \(\beta_2 = 0\))
- 対立仮説 : xモデル(パラメータ数 : \(k=2\), \(\beta_2 ≠ 0\))
第一種の過誤 : 本当は帰無仮説が正しいのに棄却されてしまう(対立仮説が正しいと判断)
第二種の過誤 : 本当は対立仮説が正しいのに帰無仮説を採択してしまう
Neyman-Pearson検定では第一種の過誤の検討にだけ専念
検定の流れ
- 帰無仮説が正しいと仮定
- 一定モデルの\(\hat{\beta_1} = 2.06\)を真のモデルとほぼ同じと考える
- 真のモデルからデータを何度も生成し, \(\beta_2 = 0\)と\(\beta_2 ≠ 0\)のモデルをあてはめ, たくさんの\(ΔD_{1, 2}\)から\(ΔD_{1, 2}\)の分布を知る
- これによって, \(ΔD_{1, 2} \geqq 4.5\)となる確率Pを評価
\(ΔD_{1, 2} = 4.5\)がありえない値と評価されれば帰無仮説が棄却される
5.4 帰無仮説を棄却するための有意水準
- \(ΔD_{1, 2} \geqq 4.5\)となる確率P = P値(p value)
- P値が十分小さければ帰無仮説を棄却できる
- 有意水準\(\alpha\)を設定(基本は0.05)
5.4.1 方法(1) 汎用性のあるパラメトリックブートストラップ法
- パラメトリックブートストラップ法
- シュミレーションによってP値を計算
#一定モデルとxモデルの逸脱度の差
fit.null$deviance - fit.x$deviance## [1] 4.513941- 一定モデルによる推定が\(\hat{\beta_1} = 2.06\)だったので平均種子数は\(\exp(2.058) = 7.83\)となる
- 真のモデルから生成されるデータは, 平均7.83の100個のポアソン乱数
exp(2.058)## [1] 7.830294#標本平均は7.83, データdに追加する
d$y.rnd <- rpois(100, lambda = mean(d$y))#一定モデルとxモデルをこの新データに当てはめてみる
fit.null.rnd <- glm(y.rnd ~ 1, data = d, family = poisson)
fit.x.rnd <- glm(y.rnd ~ x, data = d, family = poisson)
#逸脱度の差
fit.null.rnd$deviance - fit.x.rnd$deviance## [1] 0.4153217- このステップを1000回ほど繰り返すと逸脱度の差\(ΔD_{1, 2}\)の分布ができる
get.dd <- function(d)
{
#dのデータの数
n.sample <- nrow(d)
#yの平均
y.mean <- mean(d$y)
#n.sampleだけ乱数発生
d$y.rnd <- rpois(n.sample, lambda = y.mean)
fit.null.rnd <- glm(y.rnd ~ 1, data = d, family = poisson)
fit.x.rnd <- glm(y.rnd ~ x, data = d, family = poisson)
#逸脱度の差
fit.null.rnd$deviance - fit.x.rnd$deviance
}
pb <- function(d, n.bootstrap)
{
#get.ddをn.bootstrap回実施(回数を指定)
replicate(n.bootstrap, get.dd(d))
}#1000回シミュレーション
dd12 <- pb(d, n.bootstrap = 1000)
summary(dd12)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.1109 0.4806 1.0611 1.3920 15.6523#ヒストグラムにする
dd <- as.data.frame(dd12)
ggplot(data = dd) +
aes(x = dd12) +
geom_histogram(alpha = 0.5, bins = 100, color = "blue", fill = "blue") +
geom_vline(aes(xintercept = 4.5), color = "#990000", linetype = "dashed")
#p値(4.5以上となる数÷1000)
sum(dd12 >= 4.5)/1000## [1] 0.04#95パーセンタイル値
quantile(dd$dd12, 0.95)## 95%
## 3.84572この値以下はよくある差(つまりこれ以上だとレア)
尤度比検定の結論としては\(p\leqq0.05\)となるので, 帰無仮説は棄却される
- よって今回はxモデルを採択
5.4.2 方法(2) \(X^2\)分布を使った近似計算法
- \(ΔD_{1, 2}\)は近似的に自由度1の\(X^2\)分布に従う
- ただしサンプルサイズが十分に大きいことが前提
anova(fit.null, fit.x, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 99 89.507
## 2 98 84.993 1 4.5139 0.03362 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#ちなみに自由度1のカイ二乗確率点
qchisq(0.95, 1)## [1] 3.841459#自由度1のカイ2乗分布
curve( dchisq( x, df = 1 ), xlim = c( 0, 10 ), lwd = 3 )
5.5 「帰無仮説を棄却できない」は「差がない」ではない
- もし帰無仮説が棄却できなかった場合どうするのか
- 「対立仮説が正しい」ではない
- 「どっちも正しいとも正しくないともいえない」 \(=\) 判断を保留にする
- 第二種の過誤の確率\(P_2\)については検出力 : \(1-P_2\) で定義される
- 帰無仮説が棄却されべき時に正しく帰無仮説が棄却される確率
- 基本的には検出力は0.8に設定されることが多い
- https://bellcurve.jp/statistics/course/12767.html
5.6 検定とモデル選択, そして推定された統計モデルの解釈
- モデル選択や検定は文脈によって解釈も異なる
5.7 この章のまとめと参考文献
- 尤度比検定はどのような統計モデルでも使える
- 分布に関しては平均ならt分布, 分散ならF分布などを使う
- とにかくなんでも\(P<0.05\)とすればいいわけではない