汎化性能が高いモデルが欲しいーpart2

　以前の説明を繰り返させてもらうが、我々は説明変数をモデルに組み込むことで、その説明変数を使って興味ある現象を予測しようとする。予測の言うからには、今手元にあるデータだけでなく、将来取れるデータでも同じような結果になっていないと意味がないだろう。予測精度の良いモデルとは、どのように構築すればよいのだろう。

　以前述べたことと同様に、説明変数を追加すれば追加するほど、誤差に相当する値、GLMでは残差逸脱度residual devianceが減少する。しかし、新たなデータを生成して、あらかじめ推定したモデルを使って逸脱度を計算すると、無駄な説明変数が増えるほど逸脱度が増加する。つまり、これは過学習である。ゆえに、予測性の良いモデルとは、必ずしも説明変数がたくさんあるモデルではない。ではどうやって変数を選択するかにあたって、本稿では、いうなれば罰則項付き逸脱度（あるいはその真意としては罰則項付き最大対数尤度）とも言える赤池情報量基準Akaike's Information Criterion, AICについて紹介する。

GLMにおける「モデルでは説明できない誤差」ー残差逸脱度residual deviance

　線形回帰において、モデルでは説明できない誤差は残差平方和で表されていた。そもそも、最小二乗法とはこの残差平方和を最小にするようにパラメータを選ぶ方法と紹介している。しかし、GLMでは残差平方和を用いない。そもそも残差平方和はデータのばらつきが正規分布に従うときに大きな意味を持つパラメータだからだ。GLMでは下記で定義する残差逸脱度residual deviance（あるいは単に逸脱度deviance）を「モデルでは説明できない誤差に相当する量」として使っている。

サンプルサイズがnのとき、可能なモデルで最大の最大対数尤度はどんな値かを考える。それは、n個のパラメータを使って作ったモデル、フルモデルの最大対数尤度logLfullである。yiはβiのパラメータから生成されたと考えて、このβiはある確率分布の平均に相当する値とする（例えば下図）。

あるモデルの最大対数尤度logLmodelを計算したとする。それとフルモデルとの最大対数尤度の差は可能なモデルの中で最大の最大対数尤度から興味あるモデルはどれくらい離れているか、の指標となる。これが残差逸脱度の基本的な考え方である。一般には、残差逸脱度は(あるモデルの最大対数尤度ーフルモデルの最大対数尤度)に-2を掛けたものである。(あるモデルの最大対数尤度) < (フルモデルの最大対数尤度)であるから、-2を掛ければ残差逸脱度は常に正の値になり、最大の最大対数尤度からの距離としてわかりやすい指標となる。-1ではなく、-2を掛けるのは、最大対数尤度の差（logの中にまとめれば尤度の比）に-2を掛けた値が、χ2乗分布に従い、この性質を使って検定できるためである（これを尤度比検定likelihood ratio testという）。残差逸脱度が小さいほど、あるモデルの最大対数尤度は、可能なモデルの中での最大の最大対数尤度に近い、つまり当てはまりが良いことになる。

　また、残差逸脱度の計算式を整理しなおして、残差逸脱度は下記の式で定義される、逸脱度残差deviance residualsの和として考えることもできる。逸脱度残差は、ある1つの観察値が、あるモデルで生成される確率とフルモデルで生成される確率の差に-2を掛けて定義される。正直名前の付け方がややこしすぎるので、もう少し考えてほしかったが仕方ない。

説明変数が多すぎるとモデルの予測性が下がるーpart2

　では、過学習が生じるさまを改めて確認しよう。以下のようにデータを生成する。

------------------------------------------------------

library(plotn)

b0 <- 1

n <- 30

y <- rpois(n = n, lambda = exp(b0))

histn(y)#図1の描画

------------------------------------------------------

これをポアソン回帰のGLMによって解析する。その残差逸脱度を確認しよう。

------------------------------------------------------

fit1 <- glm(y ~ 1, family = poisson(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.39444  -0.54159   0.07815   0.63126   1.61141  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)   1.0531     0.1078   9.767   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 27.049  on 29  degrees of freedom

## Residual deviance: 27.049  on 29  degrees of freedom

## AIC: 112.06

## 

## Number of Fisher Scoring iterations: 4

fit2$deviance

## [1] 27.04922

------------------------------------------------------

残差逸脱度は27.04程度である。今ここで、yとは無関係な説明変数x1~x3を生成し、それらを説明変数としてGLMに組み込む。

------------------------------------------------------

x1 <- runif(n = n, 0, 30)

x2 <- runif(n = n, 0, 30)

x3 <- runif(n = n, 0, 30)

fit3 <- glm(y ~ x1, family = poisson(link = "log")) #一般化線形モデル

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x1, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.1864  -0.7357   0.1894   0.3241   1.4817  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)   

## (Intercept)  0.72822    0.26120   2.788   0.0053 **

## x1           0.02142    0.01514   1.415   0.1572   

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 27.049  on 29  degrees of freedom

## Residual deviance: 25.048  on 28  degrees of freedom

## AIC: 112.06

## 

## Number of Fisher Scoring iterations: 4

fit5 <- glm(y ~ x1 + x2, family = poisson(link = "log")) #一般化線形モデル

fit6 <- summary(fit5)

fit6

## 

## Call:

## glm(formula = y ~ x1 + x2, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.2353  -0.7495   0.2173   0.3678   1.3968  

## 

## Coefficients:

##              Estimate Std. Error z value Pr(>|z|)   

## (Intercept)  0.770737   0.281258   2.740  0.00614 **

## x1           0.024511   0.017025   1.440  0.14996   

## x2          -0.006303   0.016053  -0.393  0.69460   

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 27.049  on 29  degrees of freedom

## Residual deviance: 24.893  on 27  degrees of freedom

## AIC: 113.91

## 

## Number of Fisher Scoring iterations: 5

fit7 <- glm(y ~ x1 + x2 + x3, family = poisson(link = "log")) #一般化線形モデル

fit8 <- summary(fit7)

fit8

## 

## Call:

## glm(formula = y ~ x1 + x2 + x3, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -2.01931  -0.63299   0.01787   0.55341   1.21707  

## 

## Coefficients:

##              Estimate Std. Error z value Pr(>|z|)   

## (Intercept)  1.003527   0.317794   3.158  0.00159 **

## x1           0.026963   0.017304   1.558  0.11918   

## x2          -0.009246   0.016260  -0.569  0.56958   

## x3          -0.017653   0.012072  -1.462  0.14366   

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 27.049  on 29  degrees of freedom

## Residual deviance: 22.728  on 26  degrees of freedom

## AIC: 113.74

## 

## Number of Fisher Scoring iterations: 4

fit4$deviance

## [1] 25.04792

fit6$deviance

## [1] 24.89287

fit8$deviance

## [1] 22.72821

------------------------------------------------------

ここで示したように、x1~x3はyとは独立に生成したにもかかわらず、説明変数に加えることで逸脱度は減少した。では、ここで構築したモデルを使って、新たに生成したyにおける残差逸脱度を計算してみよう。

------------------------------------------------------

#各モデルの係数を格納

b0_m1 <- fit2$coef[1,1]

b0_m2 <- fit4$coef[1,1]

b1_m2 <- fit4$coef[2,1]

b0_m3 <- fit6$coef[1,1]

b1_m3 <- fit6$coef[2,1]

b2_m3 <- fit6$coef[3,1]

b0_m4 <- fit8$coef[1,1]

b1_m4 <- fit8$coef[2,1]

b2_m4 <- fit8$coef[3,1]

b3_m4 <- fit8$coef[4,1]

#新規データの生成

y <- rpois(n = n, lambda = exp(b0))

x1 <- runif(n = n, 0, 30)

x2 <- runif(n = n, 0, 30)

x3 <- runif(n = n, 0, 30)

#逸脱度残差の計算

d1 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m1))))

d2 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m2 + b1_m2*x1))))

d3 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m3 + b1_m3*x1 + b2_m3*x2))))

d4 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m4 + b1_m4*x1 + b2_m4*x2 + b3_m4*x3))))

#残差逸脱度の計算

sum(d1)

## [1] 46.60078

sum(d2)

## [1] 48.08419

sum(d3)

## [1] 48.05897

sum(d4)

## [1] 54.17727

cr1 <- c(sum(d1), sum(d2), sum(d3), sum(d4))

------------------------------------------------------

確かに新規のデータに対して、モデルの説明変数が多いほど残差逸脱度は増加した。もう1回、データを生成して確認しよう。

------------------------------------------------------

y <- rpois(n = n, lambda = exp(b0))

x1 <- runif(n = n, 0, 30)

x2 <- runif(n = n, 0, 30)

x3 <- runif(n = n, 0, 30)

d1 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m1))))

d2 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m2 + b1_m2*x1))))

d3 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m3 + b1_m3*x1 + b2_m3*x2))))

d4 <- 2*(log(dpois(x = y, lambda = y)) - 

           log(dpois(x = y, lambda = exp(b0_m4 + b1_m4*x1 + b2_m4*x2 + b3_m4*x3))))

sum(d1)

## [1] 23.08797

sum(d2)

## [1] 27.7372

sum(d3)

## [1] 28.11644

sum(d4)

## [1] 28.18135

cr2 <- c(sum(d1), sum(d2), sum(d3), sum(d4))

------------------------------------------------------

やはり増加傾向にあることがわかる。では、初め、モデルを生成した時の残差逸脱度とともに図示しよう。

------------------------------------------------------

tr <- c(fit2$deviance,

        fit4$deviance,

        fit6$deviance,

        fit8$deviance)

barplotn(rbind(tr, cr1, cr2), beside = T, legend = T,

         ylab = "deviance", leg.lab = c("group1", "group2", "group3"), leg.sp = 4) #図2の図示

------------------------------------------------------

こんな感じで、過剰に複雑なモデルは必ずしもデータ一般に対して予測性は高いわけではないのである。ちなみに、この後のAICの話のつながりをよくするために断っておくと、ある共通のデータのモデル間で、どれが相対的に当てはまりが良いかを考えるだけなら、フルモデルの最大対数尤度は計算する必要はない。モデル間での最大対数尤度を比較すれば十分である（データが同じならフルモデルの最大対数尤度は共通なので、モデル間の残差逸脱度の差をとると消えてしまう項であるためである）。

赤池情報量基準AIC

　そこで汎化性能の高いモデルを選ぶための基準として使われる指標が、赤池情報量基準Akaike's Information Criterionである。AICは、下記のように、あるモデルの最尤推定パラメータ^β、最大対数尤度logL(^β)、推定したパラメータの個数kを使って定義される。

最大対数尤度が当てはまりの良さを示すなら、それに-2を掛けた値は、久保先生にならって「当てはまりの悪さ」を表す指標である。「対数尤度は小さい=当てはまりが悪い」ということは、「-2×対数尤度が大きい=当てはまりが悪い」ということである。

　対数尤度だけを考えると「当てはまりの悪さ」が小さい値のモデルを選ぶことを目指す。そこにAICは「当てはまりの悪さに、推定したパラメータの数×2だけ足す」。つまり、パラメータを推定すればするほど、当てはまりの悪さを悪化させる。各モデルでAICを計算し、最小のAICのモデルが最も予測性の良いモデルであるという判断基準を与える。この発想は、罰則項付き最小二乗法に近い。いうなればAICは罰則項付き最大対数尤度である。パラメータの追加によって、最大対数尤度がほとんど増加しないのであれば、そのパラメータは汎化性能を悪化させるだけ、と考えているのである。計算自体は非常に簡単で、罰則項付き最小二乗法とは比にならない。

AICは何の値を示している？

　確かに、当てはまりの良さ＝対数尤度とパラメータ数のバランスを考えるうえで、-2×対数尤度にパラメータ数を足す、というのはシンプルな発想である。だが、足し方にもいろいろあるはずだ。パラメータに掛ける値が2ではなく、1や4でも罰則としての定性的な効果は同じである。もちろん、掛ける数が変われば選ばれるモデルは変わる。AICでは、パラメータが1個増えるにつき、最大対数尤度が1以上増えないのであれば、モデルの改善とはみなさない、と考えている。もし、パラメータ数に掛ける数が1なら、最大対数尤度が0.5増加するだけでもモデルの改善と考えるし、4なら最大対数尤度が2増加しないとモデルの改善にはならない。AICはそれらの妥当な塩梅が、2であると主張している。これには意味があるはずだ。結論を言うと、AICは-2×(あるモデルの平均対数尤度の推定値)になっている。

平均対数尤度とAIC

　統計を行う、もっとも原初のモチベーションまで戻ってしまうが、そもそも私たちは、データの生成過程における真のモデルを知ることはできない。真のモデルを知ることができないということは、分布の実態も含め、当然、説明変数にかかるパラメータの真の値もわかりはしない。そこで、真のモデルから生成されたデータをもとにモデリングし、係数を最尤推定しているわけである。

実際、正しく解析を行うことができれば、推定値は真の値に限りなく近くなる。しかし、真の値にぴったり一致することは、やはりない。なぜかと言えば、標本が有限であるからである。こうして、真の値から推定値がずれてしまうと、当然であるが、対数尤度にも影響が出る。言い方を変えれば、推定値のずれによって対数尤度にバイアスが生じてしまう。これは母分散と標本分散の期待値が一致せず、それは母平均と標本平均のずれによって生じる、という統計で一番初めにならう問題に通ずるものである。

　対数尤度の話に戻ると、真のモデルでは、真の尤度関数が存在し、真のパラメータを頂点とする山形の形状を示すはずである。しかし、真のモデルがわからない以上、真の尤度関数を知る術はない。そこで、真の尤度関数に平行な形状の関数（つまり、真のパラメータを頂点とする関数）として、平均対数尤度（具体的計算は後述）を計算することを考える。もし、手元のデータに基づく最大対数尤度から平均対数尤度を計算できれば、平均対数尤度のより大きなモデルを選べば、それは真のモデルにより近いと考えることができる。今、真のモデルから生成された手元のデータをもとに最大対数尤度とパラメータ推定値を得る。パラメータはb0の1つのみとする。

------------------------------------------------------

b0 <- 1

n <- 30

y <- rpois(n = n, lambda = exp(b0))

fit1 <- glm(y ~ 1, family = poisson(link = "log")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

## 

## Call:

## glm(formula = y ~ 1, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -2.1448  -0.9665  -0.2024   0.4407   2.0264  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)   0.8329     0.1204   6.919 4.56e-12 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 29.931  on 29  degrees of freedom

## Residual deviance: 29.931  on 29  degrees of freedom

## AIC: 106.83

## 

## Number of Fisher Scoring iterations: 5

b0e <- coef(fit2)[1,1] #最尤推定値

histn(y) #図3の図示

------------------------------------------------------

このデータにおける尤度関数を図示しよう。対数尤度を計算する関数f1（1変数用）とf2（2変数用）を定義した。最尤推定値周りでの対数尤度関数の形状は下記の通り。併せて、今回は真のパラメータを知っているので、それも図示する。当たりまえだが、最尤推定値において尤度関数は最大を示す。

------------------------------------------------------

f1 <- function(y, b0t){

  sapply(b0t, function(tmp) sum(log(dpois(x = y, lambda = exp(tmp)))))

}

f2 <- function(x, y, b0t, b1t){

  mapply(function(tmp1, tmp2) {

    sum(log(dpois(x = y, lambda = exp(tmp1 + tmp2*x))))

    }, b0t, b1t)

}

db0 <- 0.3

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f1(y, seq(b0min, b0max, length = 2000))

plotn(seq(b0min, b0max, length = 2000), logL, xlab = "b0", ylab = "logL", type = "l") #図4の図示

overdraw(abline(v = b0e),

         axis(side = 1, at = b0e, labels = "^b0"), 

         abline(v = b0),

         axis(side = 1, at = b0, labels = "b0*"))

------------------------------------------------------

さて、ここで平均対数尤度は下記のようにシミュレーションに基づき計算できる。まず、（1）真のモデル（b0 = 1）に基づき、データを生成する。この時、サンプルサイズは上記で生成した時と同じにする。（2）新たに生成したデータと上記で推定した最尤推定値^b0 = 0.83を使って対数尤度を計算する。（3）これを複数回繰り返し、対数尤度の平均をとる。これが平均対数尤度である。下記では200回繰り返した。また、下記のスクリプトでは繰り返し計算をapply関数を使って並列的に処理する方法にしている（後述のシミュレーションが重たいため）。

------------------------------------------------------

r <- 200

y_new <- matrix(rpois(n = n*r, lambda = exp(b0)), ncol = n, byrow = T)

L_tmp <- apply(y_new, 1, f1, b0t = b0e)

L_m <- mean(L_tmp)#平均対数尤度

plotn(seq(b0min, b0max, length = 2000), logL, 

      xlab = "b0", ylab = "logL", type = "l",

      ylim = range(L_tmp))#図5の図示

overdraw(abline(v = b0e),

         axis(side = 1, at = b0e, labels = "^b0"), 

         abline(v = b0),

         axis(side = 1, at = b0, labels = "b0*"),

         points(rep(b0e, r), L_tmp, col = "#00000020", pch = 16),#各対数尤度

         points(b0e, L_m, pch = 21, bg = "white", cex = 1.2))#平均対数尤度

------------------------------------------------------

このように、あるモデルの最大対数尤度と平均対数尤度はずれることがわかる。後述でさらに解説するが、あるモデルの最大対数尤度は平均対数尤度よりも大きくなる傾向がある。それは、このモデルは手元の一つのデータを過学習しており、新たに生成されるデータを考慮していないからである。直感的に、真の値から最尤推定値がずれるほど平均対数尤度はより小さくなることが予測されるだろう。なぜならば、真の値から大きくずれた最尤推定値は手元のデータにだけうまくフィットするように決められた値であり、新たに生成されるデータの大多数にはフィットしない＝対数尤度は低くなるからである。

　さて、平均対数尤度は計算できたが、これは真のモデルを知っていたからであり、結局、手元のデータだけからは計算できない。平均対数尤度の性質を理解するために、もっと様々な状況の平均対数尤度を計算することを考える。つまり、最尤推定するところからシミュレートしてみよう。10000回繰り返し計算する。

------------------------------------------------------

b0 <- 1

n <- 100

r <- 200

b0es <- NULL

L_t <- NULL

L_m <- NULL

L_tmp <- NULL

for (i in 1:10000) {

  y <- rpois(n = n, lambda = exp(b0))#1つの標本を生成

  b0e <- log(mean(y))#説明変数が切片項だけの時は、最尤推定値はこのように計算できる

  b0es <- c(b0es, b0e)

  L_t <- c(L_t, f1(y, b0e))#1つの標本の最大対数尤度を計算

  y_new <- matrix(rpois(n = n*r, lambda = exp(b0)), ncol = n, byrow = T)#新たな標本をr個、生成

  L_tmp <- apply(y_new, 1, f1, b0t = b0e)#r個の対数尤度

  L_m <- c(L_m, mean(L_tmp))#平均対数尤度を計算

}

plotn(L_m ~ b0es)#図6の図示

#最尤推定値とその時の平均対数尤度を図示

------------------------------------------------------

すると、平均対数尤度は真のパラメータ付近を頂点とする形状になっていそうなことがわかる。下記のように二次関数で回帰してみる。

------------------------------------------------------

fit1 <- lm(L_m ~ b0es + I(b0es^2))

fit2 <- summary(fit1)

I <- fit2$coef[1,1]

L <- fit2$coef[2,1]

S <- fit2$coef[3,1]

-L/(2*S)#2次関数の軸の計算

## [1] 0.9982117

xx <- seq(min(b0es), max(b0es), length = 200)

yy <- I + L*xx + S*(xx)^2

plotn(L_m ~ b0es)#図7の図示

overdraw(points(xx, yy, type = "l", col = "red"))

------------------------------------------------------

2次回帰し、その軸を計算すると真のパラメータに近い値となった。このように、確かに平均対数尤度の関数は真のモデルの尤度関数と平行な関係にあることがわかるだろう。ここで、いくつかの最大対数尤度とその時の平均対数尤度の関係を抜き出してみる。

------------------------------------------------------

tmp <- data.frame(b0es, L_m, L_t)

plotn(tmp[1:10, 1], tmp[1:10, 2:3], mode = "m", legend = T, 

      pch = c(16,1), col.dot = "black", leg.lab = c("mean logL", "max logL"),

      leg.sp = 5)#図8の図示

for(i in 1:10){

  overdraw(

    arrows(tmp[i,1], tmp[i,3], tmp[i,1], tmp[i,2], angle = 15, length = 0.2)

  )

}

------------------------------------------------------

このように平均対数尤度に対して、あるモデルの最大対数尤度は大きい場合もあるが、小さい場合もある。だが、下記で計算するように、最大対数尤度は平均対数尤度と比較して1程度大きいことがわかる。この差は偶然ではなく、各人がシミュレートしても、b0の値を変えても、やはり最大対数尤度は平均対数尤度よりも1程度大きく、この差がひらくことも縮まることはない。

------------------------------------------------------

d_L <- L_t - L_m

mean(d_L)

## [1] 1.016453

var(d_L)

## [1] 44.03061

histn(d_L)#図9の図示

------------------------------------------------------

これがもし正しいならば、真のモデルがわからず平均対数尤度だってわからないにもかかわらず、今、手元のデータで最大対数尤度を求めたとき、最大対数尤度 - 1をすればそれは平均対数尤度の推定値＝真のモデルの尤度関数における位置情報がわかる指標となることを意味する。残る問題は、最大対数尤度からどんな値を引くかである。

　そこで次は説明変数を1つ追加して平均対数尤度を計算してみよう（推定するパラメータは2つ）。真のモデルではb0 = 1、b1 = 0.05である。計算は重たいが、PCに頑張ってもらおう。

------------------------------------------------------

b0 <- 1

b1 <- 0.05

n <- 100

r <- 200

L_t <- NULL

L_m <- NULL

L_tmp <- NULL

for (i in 1:10000) {

  x <- runif(n = n, 0, 30)#説明変数を作成

  y <- rpois(n = n, lambda = exp(b0 + b1*x))#説明変数に基づく標本の作製

  fit1 <- glm(y ~ x, family = poisson(link = "log"))

  fit2 <- summary(fit1)

  b0e <- fit2$coef[1,1]#切片の最尤推定値

  b1e <- fit2$coef[2,1]#x係数の最尤推定値

  L_t <- c(L_t, f2(x, y, b0e, b1e))#最大対数尤度

  x_new <- runif(n = n*r, 0, 30)#新たな説明変数を作成

  y_new <- rpois(n = n*r, lambda = exp(b0+b1*x_new))#新たな説明変数に基づく標本の作製

  x_new <- matrix(x_new, ncol = n, byrow = T)

  y_new <- matrix(y_new, ncol = n, byrow = T)

  d <- cbind(x_new, y_new)

  L_tmp <- apply(d, 1, FUN = function(v, b0t, b1t){

    f2(x = v[1:n], y = v[(n+1):(2*n)], b0t, b1t)

  }, b0t = b0e, b1t = b1e)#r個の対数尤度

  L_m <- c(L_m, mean(L_tmp))#平均対数尤度の計算

}

d_L <- L_t - L_m

mean(d_L)

## [1] 2.108326

var(d_L)

## [1] 53.42703

histn(d_L)#図10の図示

------------------------------------------------------

推定するパラメータが2つになった時、最大対数尤度と平均対数尤度の差の平均は2程度となった。そう、実は推定するパラメータの数をkとすれば、最大対数尤度logL(^β)と平均対数尤度E(logL(^β))は下記の式の関係がある。同時にここで、AICの定義が登場する。つまり、AICは平均対数尤度の推定値を-2倍したものである。

以上の議論から、AICが最小のモデルを選ぶという方法は、平均対数尤度が最大のモデルを選ぶ＝真のモデルから生成されたあらゆるデータに対して平均的に当てはまりのよいモデルを選択することに他ならないのである。

ネストされたモデルのAIC比較

　とはいえ、一つ不安になる要因として、いくら何でも最大対数尤度と平均対数尤度の差のばらつきが大きすぎないか、という点である。実際、上記の例では、パラメータが1つの場合は分散が44、2つの場合は53もある。安心して使うためには、もう少し分散が小さい状況を想定できるほうが嬉しいだろう。

　通常、AICはネストnestされたモデル間の比較で用いる。あるモデルAが、あるモデルBに含まれている状況を、AがBにネストされている、という。具体的には、切片パラメータをβ0とし、以降、説明変数xiに対応するパラメータをβiとすると、モデルAはβ0、モデルBはβ0とβ1を推定するモデルであれば、モデルAはモデルBにネストされている。モデルCはβ0～β2を推定するモデルなら、AとBはともにCにネストされている。モデルDではβ0とβ2を推定するモデルだとすれば、AはDにネストされるし、AとDはともにCにネストされているが、BとDはたがいにネストの関係ではない。

上記のモデルAとモデルBをAICで比較する状況を考える。AICの性質から、モデルAにおける最大対数尤度 - 平均対数尤度の期待値は1となり、モデルBにおける最大対数尤度 - 平均対数尤度の期待値は2となる。今、モデルAおよびBの対数尤度の差の期待値のばらつきが大きすぎることが懸念されている。

　ここで、さらに「モデルAの対数尤度の差」と「モデルBの対数尤度の差」の差をとることを考える。もし、β1の推定によって対数尤度がほとんど改善しない、つまりモデルAの対数尤度logL(y|^β0)≒モデルBの対数尤度logL(y|^β0, ^β1)のとき、「対数尤度の差」の差は1になることが期待される。これは(モデルBのAIC) - (モデルAのAIC)の1/2量でもある。このような「対数尤度の差」の差はばらつきが小さいことが知られる。つまり、ネスト関係にあるモデルのAICの差のばらつきは小さいため、あまり問題にならないだろう、と考えられている。確認してみよう。

　今、切片パラメータのモデルと切片および説明変数xの係数パラメータを推定する2つのモデルを考える。同じデータに関して、それぞれでモデル推定し、対数尤度の差を計算する。さらに、モデル間で「対数尤度の差」の差を計算する。上記のモデルAの対数尤度logL(y|^β0)≒モデルBの対数尤度logL(y|^β0, ^β1)という状況は、つまり真のモデルではβ1 = 0であることを想定していることに注意する。

------------------------------------------------------

b0 <- 1

b1 <- 0

n <- 100

r <- 200

L_t_n <- NULL

L_m_n <- NULL

L_t_x <- NULL

L_m_x <- NULL

L_tmp_n <- NULL

L_tmp_x <- NULL

for (i in 1:10000) {

  x <- runif(n = n, 0, 30)#説明変数を作成

  y <- rpois(n = n, lambda = exp(b0 + b1*x))#説明変数に基づく標本の作製

  b0e_n <- log(mean(y))#1パラメータモデルの切片パラメータの最尤推定値

  fit1 <- glm(y ~ x, family = poisson(link = "log"))

  fit2 <- summary(fit1)

  b0e_x <- fit2$coef[1,1]#2パラメータモデルの切片パラメータの最尤推定値

  b1e_x <- fit2$coef[2,1]#2パラメータモデルの係数パラメータの最尤推定値

  L_t_n <- c(L_t_n, f1(y, b0e_n))#1パラメータモデルの最大対数尤度

  L_t_x <- c(L_t_x, f2(x, y, b0e_x, b1e_x))#2パラメータモデルの最大対数尤度

  x_new <- runif(n = n*r, 0, 30)#新たな説明変数の生成

  y_new <- rpois(n = n*r, lambda = exp(b0+b1*x_new))#新たな標本の生成

  x_new <- matrix(x_new, ncol = n, byrow = T)

  y_new <- matrix(y_new, ncol = n, byrow = T)

  d <- cbind(x_new, y_new)

  L_tmp_n <- apply(y_new, 1, f1, b0t = b0e_n)#1パラメータの最尤推定値をもとにした対数尤度

  L_tmp_x <- apply(d, 1, FUN = function(v, b0t, b1t){

    f2(x = v[1:n], y = v[(n+1):(2*n)], b0t, b1t)

    }, b0t = b0e_x, b1t = b1e_x)#2パラメータの最尤推定値をもとにした対数尤度

  L_m_n <- c(L_m_n, mean(L_tmp_n))#1パラメータの平均対数尤度

  L_m_x <- c(L_m_x, mean(L_tmp_x))#2パラメータの平均対数尤度

}

d_L_n <- L_t_n - L_m_n #1パラメータの対数尤度の差

d_L_x <- L_t_x - L_m_x #2パラメータの対数尤度の差

d_L <- d_L_x - d_L_n #モデル間の「対数尤度の差」の差

mean(d_L) #「対数尤度の差」の差の平均

## [1] 1.012558

var(d_L) #「対数尤度の差」の差の分散

## [1] 2.031392

histn(d_L, xlab = "logL") #図11の作図

------------------------------------------------------

すると期待通り、「対数尤度の差」の差の期待値は1になり、そのばらつきも対数尤度の差と比較しても極めて小さくなった。これは、結局、モデルAの対数尤度の差のばらつきとモデルBの対数尤度の差のばらつきが、ちょうどいい具合にキャンセルしたことを意味する。それゆえ、ネストされたモデル間でのAICの比較に限って言えば、AICの差（モデル間の距離）はどんな状況でもそれなりの精度で予測できるということである。

Rで行うAICに基づくモデル選択

　では、RでAICによるモデル選択を実行してみよう。とは言えどもglm関数では自動でAICを計算してくれているので、いくつか検討したいモデルを作成した後は、AICを比較するだけで終わる。もし、説明変数がたくさんあって、比較するべきモデルが大量にある場合は、MASSパッケージのstepAIC関数によって自動でAICが最も小さなモデルの候補（あくまで候補）を選んでくれる（後述）。

　以下のようなデータにおいて、説明変数x1とx2をモデルに組み込むかを考えてみよう。まずはデータを図示する。

------------------------------------------------------

x1 <- c(19.7, -2.9, 5.9, -4.7, 12.7, -1.8, 19.5, 

        -4.1, -3.3, 12, 5.8, 2.1, 2.4, 12.4, 0.8, 

        -9.1, -3.3, 18.1, 8.7, 2.8, -9.9, 2.7,

        1.5, -2.8, 19.6, 4.6, 4.4, -6, 19.6, -2.9)

x2 <- c(18.3, 24.9, 11.8, 11.9, 9.5, 12.9, 29.4,

        29.3, 4.3, 3.3, 15.5, 8.2, 10, 19.4, 12.1, 

        13, 26.7, 23.4, 19.5, 28.8, 1.9, 24.2, 2.8, 

        5.4, 20.3, 19.7, 11.9, 26.2, 27.1, 10.6)

y <- c(7, 3, 4, 3, 7, 1, 6, 6, 1, 4, 3, 7,

       5, 2, 4, 1, 1, 6, 3, 2, 4, 4, 5, 5,

       6, 2, 3, 1, 11, 2)

plotn(y ~ x1) #図12の作図

plotn(y ~ x2) #図13の作図

------------------------------------------------------

データの様子を見るに、yとx1は関係が強そうだが、yとx2は関係が弱くてx2を説明変数に組み込む意義は薄いように思える。では、切片だけのモデル、切片とx1のモデル、切片とx2のモデル、切片、x1、x2のモデルでAICを比較してみよう。

------------------------------------------------------

fit1 <- glm(y ~ 1, family = poisson(link = "log"))

fit2 <- glm(y ~ x1, family = poisson(link = "log"))

fit3 <- glm(y ~ x2, family = poisson(link = "log"))

fit4 <- glm(y ~ x1 + x2, family = poisson(link = "log"))

summary(fit1)

## 

## Call:

## glm(formula = y ~ 1, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -1.78255  -1.09280   0.01671   0.83586   2.89355  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  1.37793    0.09167   15.03   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 40.324  on 29  degrees of freedom

## Residual deviance: 40.324  on 29  degrees of freedom

## AIC: 134.81

## 

## Number of Fisher Scoring iterations: 5

summary(fit2)

## 

## Call:

## glm(formula = y ~ x1, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.5811  -0.8538  -0.2288   0.6593   1.7163  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept) 1.157213   0.117114   9.881  < 2e-16 ***

## x1          0.038637   0.009954   3.882 0.000104 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 40.324  on 29  degrees of freedom

## Residual deviance: 25.516  on 28  degrees of freedom

## AIC: 122

## 

## Number of Fisher Scoring iterations: 4

summary(fit3)

## 

## Call:

## glm(formula = y ~ x2, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -1.8783  -1.1020  -0.0041   0.7571   2.7488  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept) 1.289164   0.200102   6.443 1.17e-10 ***

## x2          0.005455   0.010800   0.505    0.614    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 40.324  on 29  degrees of freedom

## Residual deviance: 40.068  on 28  degrees of freedom

## AIC: 136.56

## 

## Number of Fisher Scoring iterations: 5

summary(fit4)

## 

## Call:

## glm(formula = y ~ x1 + x2, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##      Min        1Q    Median        3Q       Max  

## -1.56354  -0.77943  -0.09019   0.50594   2.01213  

## 

## Coefficients:

##             Estimate Std. Error z value Pr(>|z|)    

## (Intercept)  1.29441    0.20374   6.353 2.11e-10 ***

## x1           0.04191    0.01087   3.856 0.000115 ***

## x2          -0.00962    0.01197  -0.804 0.421673    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## (Dispersion parameter for poisson family taken to be 1)

## 

##     Null deviance: 40.324  on 29  degrees of freedom

## Residual deviance: 24.868  on 27  degrees of freedom

## AIC: 123.36

## 

## Number of Fisher Scoring iterations: 5

------------------------------------------------------

上記のようにGLMを行い、それを格納したオブジェクトをsummary関数に通せば、AICを確認できる。以下のように、AIC関数に複数のオブジェクトを入力すれば、AICだけを同時に確認できる。

------------------------------------------------------

AIC(fit1, fit2, fit3, fit4)

##      df      AIC

## fit1  1 134.8114

## fit2  2 122.0041

## fit3  2 136.5563

## fit4  3 123.3556

------------------------------------------------------

この結果に基づけば、AIC最小のモデルはfit2、つまり切片とx1のモデルである。ゆえに、モデルとしての変数はx1だけを採用すると判断する。

　モデルが4つだけならこのようにしても簡単だが、もっと多いときは大変だ。そこでMASSパッケージのstepAIC関数を使えば、自動で最小AICのモデルの候補を選んでくれる。stepAICのアルゴリズムは、最も多くの説明変数を含むモデルから一つずつ説明変数の除去（forward方向）、あるいは切片モデルから一つずつ説明変数の追加（backward方向）、あるいはその両方（both方向）を行いながらAICを計算して最小のモデルをさがす。あくまで最小のAICモデルの候補である理由は、すべての説明変数の組み合わせについて検討しているわけではないからである。

　stepAICに最も多くの説明変数の含むモデルを入力し、direction引数でAICの計算の方向（forward, backward, both）を決めたら、実行してみよう。

------------------------------------------------------

library(MASS)

stepAIC(fit4, direction="both")

## Start:  AIC=123.36

## y ~ x1 + x2

## 

##        Df Deviance    AIC

## - x2    1   25.516 122.00

## <none>      24.868 123.36

## - x1    1   40.068 136.56

## 

## Step:  AIC=122

## y ~ x1

## 

##        Df Deviance    AIC

## <none>      25.516 122.00

## + x2    1   24.868 123.36

## - x1    1   40.324 134.81

## 

## Call:  glm(formula = y ~ x1, family = poisson(link = "log"))

## 

## Coefficients:

## (Intercept)           x1  

##     1.15721      0.03864  

## 

## Degrees of Freedom: 29 Total (i.e. Null);  28 Residual

## Null Deviance:       40.32 

## Residual Deviance: 25.52     AIC: 122

------------------------------------------------------

最終的な結果として、やはりx1を説明変数に含むモデルが選択された。以上でAICの比較により、汎化性能が高いと思われるモデルを構築できた。さて、重要な注意点を述べる。罰則項付き最小二乗法の時と同様、あくまで未知のデータに対してもある程度、予測精度を担保したモデルを構築する手段であり、真のモデルを推定するものではない。得られたモデルが真のモデルであると信じることは、とんでもない間違いを産む可能性があるので、十分注意しよう。ただし、もしAICを比較するモデルの中に正しいモデルが含まれている場合、データが無限にとれれば正しいモデルを選択することは可能である。

追記：カルバック・ライブラー情報量Kullback–Leibler divergence