chapter04 GLMのモデル選択 -AICとモデルの予測の良さ

4.1 データは一つ, モデルはたくさん

• 第3章で扱ったモデルに対して, 4つのモデルを設定
  • x モデル
  • f モデル
  • x + f モデル
  • 一定モデル(切片\(\beta1\)のみ)
• 前章では最大対数尤度でモデルの良さを測定したが, 果たしてそれは正しいのか, というのがテーマ

4.2 統計モデルのあてはまりの悪さ : 逸脱度

• 逸脱度(deviance)(あてはまりの悪さ)は以下の式で表される
  • \(D = -2logL\) (ただし\(logL\)は対数尤度)

#まずは各モデルを作成
setwd("~/Desktop/Statistics_Study/StatisticalModeling/kubobook_2012/03poisson")
library(readr)
d<- read_csv("data3a.csv")

## 
## ─ Column specification ────────────────────────────
## cols(
##   y = col_double(),
##   x = col_double(),
##   f = col_character()
## )

#xモデル
fit <- glm(y ~ x, data = d, family = poisson(link = "log"))
#fモデル
fit.f <- glm(y ~ f, data = d, family = poisson)
#x + f モデル
fit.all <- glm(y ~ x + f, data = d, family = poisson)
#一定モデル(nullモデル)
fit.null <-glm(y ~ 1, data = d, family = poisson)

#各モデルの要約
library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

stargazer(fit, fit.f, fit.all, fit.null, type = "text", style = "all", ci = TRUE,
          star.cutoffs = NA, omit.table.layout = 'n',
          align = TRUE)

## 
## ===========================================================================================
##                                                 Dependent variable:                        
##                         -------------------------------------------------------------------
##                                                          y                                 
##                               (1)              (2)              (3)              (4)       
## -------------------------------------------------------------------------------------------
## x                            0.076                             0.080                       
##                          (0.006, 0.145)                    (0.007, 0.153)                  
##                            t = 2.125                         t = 2.162                     
##                            p = 0.034                         p = 0.031                     
## fT                                            0.013            -0.032                      
##                                          (-0.127, 0.153)  (-0.178, 0.114)                  
##                                             t = 0.179        t = -0.430                    
##                                             p = 0.859        p = 0.668                     
## Constant                     1.292            2.052            1.263            2.058      
##                          (0.579, 2.005)   (1.952, 2.151)   (0.539, 1.988)   (1.988, 2.128) 
##                            t = 3.552        t = 40.463       t = 3.417        t = 57.586   
##                            p = 0.0004       p = 0.000        p = 0.001        p = 0.000    
## -------------------------------------------------------------------------------------------
## Observations                  100              100              100              100       
## Log Likelihood              -235.386         -237.627         -235.294         -237.643    
## Akaike Inf. Crit.           474.773          479.255          476.587          477.286     
## Residual Deviance       84.993 (df = 98) 89.475 (df = 98) 84.808 (df = 97) 89.507 (df = 99)
## Null Deviance (df = 99)      89.507           89.507           89.507           89.507     
## ===========================================================================================

• 結果より, xモデルの最大対数尤度は\(-235.4\)なので, 逸脱度:\(D = 470.8\)となる
  • しかしどこにも\(470.8\)はない…
• ここで残差逸脱度(Residual Deviance)は以下の式で定義される
  • \(D - (ポアソン分布モデルで可能な最小逸脱度)\)
  • \(最小逸脱度 = フルモデルの逸脱度\) →今回であればデータ数100個をあてはめたモデル
  • 残差逸脱度は以下のように求まる

#最小逸脱度
MD <-sum(log(dpois(d$y, lambda = d$y)))
#xモデルの逸脱度
x_D <-(-2*(-235.386))
#xモデルの残差逸脱度
x_RD <- x_D -(-2 * MD)
#結果のxモデルのResidual Devianceと一致
x_RD

## [1] 84.99249

• つまり残差逸脱度は385.8を基準とする「あてはまりの悪さの相対値」
• 逸脱度が大きいほど「あてはまりの悪いモデル」→ここではnullモデルが該当
• 結果を見るとパラメータが増えるほど残差逸脱度は小さくなっているのがわかる

4.3 モデル選択基準AIC

• ここではモデル選択基準の1つであるAICを紹介
  \(AIC = -2(logL - k)= D + 2k(Dは逸脱度, kはパラメータ)\)
  • このAICが一番小さいモデルが良いモデル
  • 結果の表からはxモデルが一番AICが小さいとわかる

4.4 AICを説明するためのまた別の問題

• 種子数\(y\)と体サイズ\(x\)において2つのモデルが考えられる
  1. \(logl = \beta1\)(一定モデル \(k=1\))
  2. \(logl = \beta1 + \beta2x\)(xモデル \(k=2\))
• この2つのモデルはネストしている

4.5 なぜAICでモデル選択してよいのか?

4.5.1 統計モデルの予測の良さ : 平均対数尤度

• ※ここでは真のモデルがわかっている
• 最大対数尤度 \(=\) たまたま得られた観測データへのあてはまり
• 予測の良さを評価するには同じ方法で別のデータを取ってきた時の精確さを見る
• シミュレーションを繰り返して(ここでは200回)得られた対数尤度の平均 \(=\) 平均対数尤度

4.5.2 最大対数尤度のバイアス補正

• 本章での一定モデルの最大対数尤度は\(logL = -120.6\), 一方平均対数尤度は\(E(logL) = -122.9\)
  • つまりたまたま得られた一定モデルはあてはまりの良さを過大評価しているということになる(0に近いほうがいいから)
    
  • 毎回\(logL>E(logL)\)というわけではない
  • 上記のシミュレーションを1セットとして12回繰り返す
  • 最大対数尤度と平均対数尤度の差の分布を\(b = logL - E(logL)\)とする→バイアス
    • ここではバイアス\(b\)の平均が1.01になった
    • あてはまりの良さである最大対数尤度のほうが平均的に1くらい大きい
• ???とはいえ現実のデータ分析では真のモデルはわからないから, \(E(logL)\)はどうなるか不明
• 次のように式変形をする\(E(logL) = logL - b\)
  • \(logL\)と平均的な\(b\)がわかれば推定できそう→バイアス補正
  • パラメータをk個もつモデルの平均対数尤度の推定量は\(logL-k\)(シミュレーションの結果とほぼ一致!)
• これに\(-2\)をかけるとAICになる
  AIC = -2*(logL - 1)

4.5.3 ネストしているGLM間のAIC比較

• バイアスである\(b\)にバラつきはないのか?
  • \((xモデルのlogL) - (一定モデルのlogL) \fallingdotseq 0.5\)
    • パラメータが1増えると最大対数尤度が0.5大きくなる→あてはまりが良くなる
  • \((xモデルのE(logL)) - (一定モデルのE(logL)) \fallingdotseq -0.5\)
    • パラメータが1増えると平均対数尤度は0.5小さくなる→あてはまりが悪くなる -したがって, \(logL\)の平均バイアスは1増加
  • また, パラメータが1増えることでバイアスは平均1増加(ネストしているモデルだとバラつきは少ない)

4.6 この章のまとめ

• AICは良い予測をするモデルである(\(\neq\)あてはまりのよいモデル)
• AICはあくまでも自分で作ったモデルの中で比較して良いモデルを選んでいるだけ