渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

モデル選択とAICとcross validationの関係を大雑把に実験してみる

Stanの開発者でもある統計学界の重鎮、Andrew Gelmanがこんなブログ記事をupしていました。

ちなみに@さんがこの件についてこんなコメントをされてました。

WAICとcross validationの関係については渡辺澄夫先生の本にも当然のように記述があって、そこで密林でポチって読んでみたのですが、いかんせん僕の数学力が弱すぎてまだ全然理解できてない有様ですorz

ベイズ統計の理論と方法

ベイズ統計の理論と方法

というザマでこの件について自分では何も勉強出来てないんですが(汗)、そもそも論としてモデル選択においてAICとcross validationってどんな関係性にあるんだっけ?というのをふと思ったので、大雑把に実験してみようかと思います。


データセットの準備


既に挙動がはっきりしているものということで、テニス四大大会データセットを用います。

'men.txt'を学習データに、'women.txt'をholdoutデータに使います。

# Training
> dm <- read.table('men.txt',sep='\t',header=T)
> dm <- dm[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]
> dm$Result<-as.factor(dm$Result)

# Holdout
> dw <- read.table('women.txt',sep='\t',header=T)
> dw <- dw[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]

何も考えずにロジスティック回帰モデルを推定して、AIC / LOO ACC / holdout ACCを出してみる


とりあえず素のロジスティック回帰をやった上で、AIC / LOO ACC / holdout ACCを出してみましょう。

> dm.glm <- glm(Result~., dm, family=binomial)
> summary(dm.glm)

Call:
glm(formula = Result ~ ., family = binomial, data = dm)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.91243  -0.14289  -0.00641   0.08835   3.04139  

Coefficients: (2 not defined because of singularities)
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.378486   5.064233  -0.272 0.785468    
FSP.1        0.116125   0.065947   1.761 0.078258 .  
FSW.1        0.198896   0.043094   4.615 3.92e-06 ***
SSP.1              NA         NA      NA       NA    
SSW.1        0.290083   0.069708   4.161 3.16e-05 ***
ACE.1       -0.040074   0.048374  -0.828 0.407430    
DBF.1       -0.136684   0.082271  -1.661 0.096636 .  
WNR.1        0.074055   0.029644   2.498 0.012485 *  
UFE.1       -0.007523   0.030277  -0.248 0.803779    
BPC.1        0.393523   0.102120   3.854 0.000116 ***
BPW.1        0.313412   0.069018   4.541 5.60e-06 ***
NPA.1        0.012960   0.042315   0.306 0.759398    
NPW.1       -0.014493   0.036409  -0.398 0.690578    
FSP.2       -0.065368   0.066044  -0.990 0.322288    
FSW.2       -0.263955   0.047725  -5.531 3.19e-08 ***
SSP.2              NA         NA      NA       NA    
SSW.2       -0.247359   0.063257  -3.910 9.22e-05 ***
ACE.2        0.026067   0.038823   0.671 0.501937    
DBF.2        0.117291   0.084828   1.383 0.166759    
WNR.2       -0.032345   0.028113  -1.151 0.249925    
UFE.2       -0.024302   0.028906  -0.841 0.400490    
BPC.2       -0.329055   0.111769  -2.944 0.003239 ** 
BPW.2       -0.351452   0.087206  -4.030 5.57e-05 ***
NPA.2       -0.036414   0.059266  -0.614 0.538935    
NPW.2        0.047514   0.042770   1.111 0.266603    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 680.42  on 490  degrees of freedom
Residual deviance: 150.87  on 468  degrees of freedom
AIC: 196.87

Number of Fisher Scoring iterations: 8

> res <- rep(0, nrow(dm))
> for (i in 1:nrow(dm)){
+     train <- dm[-i,]
+     test <- dm[i,]
+     model <- glm(Result~., train, family=binomial)
+     pred <- round(predict(model, newdata=test, type='response'),0)
+     res[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0)
+ }
> sum(res)/nrow(dm)
[1] 0.9205703
> sum(diag(table(dw$Result, round(predict(dm.glm, newdata=dw[,-1], type='response'),0))))/nrow(dw)
[1] 0.9269912

AIC = 196.87, LOO ACC = 0.921, holdout ACC = 0.927という結果になりました。


ステップワイズ法でAIC最小のモデルを選んで、AIC / LOO ACC / holdout ACCを出してみる


普通に{MASS}パッケージのstepAICを使うだけです。交互作用は想定していないのでご注意ください。

> library(MASS)
> stepAIC(dm.glm)

# 中略 #

        Df Deviance    AIC
<none>       153.86 185.86
- DBF.2  1   155.88 185.88
- FSP.2  1   156.11 186.11
- NPW.2  1   156.16 186.16
- DBF.1  1   156.25 186.25
- FSP.1  1   158.54 188.54
- UFE.2  1   161.87 191.87
- WNR.1  1   165.91 195.91
- SSW.2  1   176.82 206.82
- SSW.1  1   178.59 208.59
- BPW.2  1   180.71 210.71
- FSW.1  1   181.31 211.31
- BPC.2  1   182.61 212.61
- BPC.1  1   187.76 217.76
- BPW.1  1   189.80 219.80
- FSW.2  1   200.18 230.18

Call:  glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + 
    BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + 
    BPW.2 + NPW.2, family = binomial, data = dm)

Coefficients:
(Intercept)        FSP.1        FSW.1        SSW.1        DBF.1        WNR.1        BPC.1        BPW.1  
   -0.67573      0.12962      0.18332      0.29546     -0.11832      0.05049      0.38807      0.32564  
      FSP.2        FSW.2        SSW.2        DBF.2        UFE.2        BPC.2        BPW.2        NPW.2  
   -0.09027     -0.25356     -0.26049      0.11584     -0.04035     -0.39068     -0.31849      0.02078  

Degrees of Freedom: 490 Total (i.e. Null);  475 Residual
Null Deviance:	    680.4 
Residual Deviance: 153.9 	AIC: 185.9

> dm.glm.aic <- glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + 
+                       BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + 
+                       BPW.2 + NPW.2, family = binomial, data = dm)
> summary(dm.glm.aic)

Call:
glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + 
    BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + 
    BPW.2 + NPW.2, family = binomial, data = dm)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8678  -0.1372  -0.0069   0.1008   3.0421  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.67573    5.02335  -0.135  0.89299    
FSP.1        0.12962    0.06182   2.097  0.03602 *  
FSW.1        0.18332    0.03890   4.713 2.44e-06 ***
SSW.1        0.29546    0.06685   4.420 9.87e-06 ***
DBF.1       -0.11832    0.07914  -1.495  0.13490    
WNR.1        0.05049    0.01536   3.287  0.00101 ** 
BPC.1        0.38807    0.08934   4.344 1.40e-05 ***
BPW.1        0.32564    0.06178   5.271 1.36e-07 ***
FSP.2       -0.09027    0.06117  -1.476  0.14003    
FSW.2       -0.25356    0.04424  -5.732 9.94e-09 ***
SSW.2       -0.26049    0.06132  -4.248 2.16e-05 ***
DBF.2        0.11584    0.08149   1.422  0.15516    
UFE.2       -0.04035    0.01475  -2.735  0.00624 ** 
BPC.2       -0.39068    0.09954  -3.925 8.67e-05 ***
BPW.2       -0.31849    0.06928  -4.598 4.28e-06 ***
NPW.2        0.02078    0.01396   1.488  0.13665    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 680.42  on 490  degrees of freedom
Residual deviance: 153.86  on 475  degrees of freedom
AIC: 185.86

Number of Fisher Scoring iterations: 8

> res.aic <- rep(0, nrow(dm))
> for (i in 1:nrow(dm)){
+     train <- dm[-i,]
+     test <- dm[i,]
+     model <- glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + 
+             BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + 
+             BPW.2 + NPW.2, family = binomial, data = train)
+     pred <- round(predict(model, newdata=test, type='response'),0)
+     res.aic[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0)
+ }
> sum(res.aic)/nrow(dm)
[1] 0.9287169
> sum(diag(table(dw$Result, round(predict(dm.glm.aic, newdata=dw[,-1], type='response'),0))))/nrow(dw)
[1] 0.9314159

AIC = 185.86, LOO ACC = 0.929, holdout ACC = 0.931という結果になりました。


L1正則化モデル


{glmnet}でalpha = 1としてL1正則化ロジスティック回帰でやってみます。ちなみにLOOは2つ考え方があるかと思うのですが、一応両論併記になるように「あらかじめ決め打ちでL1正則化モデルを固定して回す」のと「都度L1正則化モデルを10-folds CVでパラメータ決定して回す」のとを使い分けています。ただしholdout ACCは10-fold CVでパラメータ決定したものにしてあります。

> library(glmnet)
> dm.cv.glmnet <- cv.glmnet(as.matrix(dm[,-1]), as.integer(dm[,1])-1, family="binomial", alpha=1)
> plot(dm.cv.glmnet)
> coef(dm.cv.glmnet, s=dm.cv.glmnet$lambda.min)
25 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  3.533402e-01
FSP.1        3.805604e-02
FSW.1        1.179697e-01
SSP.1       -3.275595e-05
SSW.1        1.475791e-01
ACE.1        .           
DBF.1       -8.934231e-02
WNR.1        3.628403e-02
UFE.1       -7.839983e-03
BPC.1        3.758665e-01
BPW.1        2.064167e-01
NPA.1        .           
NPW.1        .           
FSP.2       -2.924528e-02
FSW.2       -1.568441e-01
SSP.2        .           
SSW.2       -1.324209e-01
ACE.2        1.233763e-02
DBF.2        4.032510e-02
WNR.2       -2.071361e-02
UFE.2       -6.114823e-06
BPC.2       -3.648171e-01
BPW.2       -1.985184e-01
NPA.2        .           
NPW.2        1.340329e-02

# 一旦glmを使って直す
> dm.glm.l1 <- glm(formula = Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + DBF.1 + WNR.1 + 
+         UFE.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + ACE.2 + DBF.2 + WNR.2 + UFE.2 + BPC.2 
+            BPW.2 + NPW.2, family = binomial, data = dm)
> summary(dm.glm.l1)

Call:
glm(formula = Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + DBF.1 + 
    WNR.1 + UFE.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + ACE.2 + 
    DBF.2 + WNR.2 + UFE.2 + BPC.2 + BPW.2 + NPW.2, family = binomial, 
    data = dm)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8984  -0.1487  -0.0076   0.1013   3.1776  

Coefficients: (1 not defined because of singularities)
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.430594   5.007621  -0.286 0.775121    
FSP.1        0.131931   0.063677   2.072 0.038276 *  
FSW.1        0.183218   0.039567   4.631 3.65e-06 ***
SSP.1              NA         NA      NA       NA    
SSW.1        0.290247   0.068486   4.238 2.25e-05 ***
DBF.1       -0.129578   0.082095  -1.578 0.114479    
WNR.1        0.067105   0.026901   2.495 0.012612 *  
UFE.1       -0.008339   0.028484  -0.293 0.769697    
BPC.1        0.375908   0.088707   4.238 2.26e-05 ***
BPW.1        0.313854   0.063923   4.910 9.11e-07 ***
FSP.2       -0.080219   0.064170  -1.250 0.211257    
FSW.2       -0.257784   0.046546  -5.538 3.05e-08 ***
SSW.2       -0.250462   0.062454  -4.010 6.06e-05 ***
ACE.2        0.031084   0.038078   0.816 0.414312    
DBF.2        0.098090   0.082402   1.190 0.233895    
WNR.2       -0.032804   0.027685  -1.185 0.236053    
UFE.2       -0.017172   0.027139  -0.633 0.526898    
BPC.2       -0.350231   0.100279  -3.493 0.000478 ***
BPW.2       -0.318826   0.074215  -4.296 1.74e-05 ***
NPW.2        0.020829   0.014036   1.484 0.137833    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 680.42  on 490  degrees of freedom
Residual deviance: 152.04  on 472  degrees of freedom
AIC: 190.04

Number of Fisher Scoring iterations: 8

> res.l1 <- rep(0, nrow(dm))
> for (i in 1:nrow(dm)){
+     train <- dm[-i,]
+     test <- dm[i,]
+     model <- glm(formula = Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + DBF.1 + WNR.1 + 
+             UFE.1 + BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + ACE.2 + DBF.2 + WNR.2 + UFE.2 + BPC.2 + 
+             BPW.2 + NPW.2, family = binomial, data = train)
+     pred <- round(predict(model, newdata=test, type='response'),0)
+     res.l1[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0)
+ }
> sum(res.l1)/nrow(dm)
[1] 0.9205703
> sum(diag(table(dw$Result, round(predict(dm.glm.l1, newdata=dw[,-1], type='response'),0))))/nrow(dw)
[1] 0.9269912

# {glmnet}縛り
> res.l1 <- rep(0, nrow(dm))
> for (i in 1:nrow(dm)){
+     train <- dm[-i,]
+     test <- dm[i,]
+     model <- cv.glmnet(as.matrix(train[,-1]),as.integer(train[,1])-1,family="binomial",alpha=1)
+     pred <- round(predict(model, newx=as.matrix(test[,-1]), type='response', s=model$lambda.min),0)
+     res.l1[i] <- ifelse(pred==as.integer(test[1])-1, 1, 0)
+ }
> sum(res.l1)/nrow(dm)
[1] 0.9205703
> sum(diag(table(dw$Result, round(predict(dm.cv.glmnet, newx=as.matrix(dw[,-1]), s=dm.cv.glmnet$lambda.min, type='response'),0))))/nrow(dw)
[1] 0.9336283

AIC = 190.04, LOO ACC = 0.921, holdout ACC = 0.934という結果になりました。

3つのモデルのパフォーマンスを比べてみる


単にプロットしてみると以下のような結果が得られます。

f:id:TJO:20161027153521p:plain


全体としてみればもちろんAICが小さければ小さいほどACCは2つとも改善しているのが分かるかと思いますが、LOOとholdoutとで微妙に結果が割れています。これについては@先生が本来のWAICとLOOとの比較についてコメントされているのと話が似ているのかなと。

つまり、LOOというのは結局のところ学習データの内部に完結させているものなので、たまたま学習データ自体に変なばらつきがあって毎回抽出するごとに推定されたモデルと真のデータ生成過程とが大きく乖離しているようなケースではLOOは悪くなるけれども、holdoutのような大きめのデータセットに対しては「均される」のでパフォーマンスが良くなるみたいなことはあるのかもしれないな、とはちょっと思ってます。


とはいえ、本当は間を埋めるような微妙なAICの値のモデルを幾つか選んできてサンプルサイズを増やした方が良かったのでしょうが、とりあえずザッとやってみましたという程度なので今回はこの辺にしときます。。。


追記


何と@先生から言及いただいてしまいました。以下そのまま貼ります。