モデル選択の実験 (AIC vs AICc / R の AICcmodavg パッケージ)

  • ????????????????????

前回の「モデル選択の実験 (AIC vs LOOCV)」の続きです。

小標本の場合は、AIC じゃなくて AICc を使うといいよとのことなので、今回は、前回同様の方法で AIC と AICc を比較してみた。

真のモデルやその他のモデルの設定などは前回と全く同様。図の見方も前回と同様です。LOOCV の結果も並べてみたかったが計算量の関係で断念。

いろんな意味で手抜き気味です。あしからず。

AICc

「正規ノイズの線形モデル」のケースでは以下で定義される情報量基準。

\begin{align}
AICc = AIC + \frac{2k(k+1)}{n-k-1}
\end{align}
k : パラメータ数、n : データ数。

GLM のケースではこの定式化は使えないらしい[1]。 そんなこんなで GLM の場合は R の AICcmodavg パッケージをつかう。

glm.fit <- glm(...)
AICc(glm.fit)

とすることで AICc の値を良きに計らって計算してくれる。

実験

真のモデルからデータを生成し、推定、モデル選択を行う、という作業を 10000 回繰り返した。

次のグラフはデータ数 n=20 のとき。確かに、AICc のほうが Model 2 の選択率が高そうだ。真のモデルは 図の中で Model 2 となっているやつ。

n=40 とするとこんなかんじ。

n=60, 80, 100 の様子はそんなに変わらない。 n が大きくなっても正解率が 1 に近づいていく様子はない(一致性をもたない、らしい)

ここまでの結果のまとめ

横軸にデータ数、縦軸に正解率をとってグラフを描いてみた。

  • 常に AICc のほうが好成績。
  • \(n\to\infty\) で明らかに AIC=AICc なので、グラフのもっともっと右端では2つの線は一致する予定。
  • AIC はデータ数を増やしても正解率が改善していく様子はなし
  • AICc の正解率が途中から低下してくる。何故?どのように考えればいい?

実験に使った R スクリプト

#-------------------------------------------------
# AIC vs. AICc
#
# [なにをするの?]
#  -- モデル選択手法である AIC と AICc を比較してみる
#  -- 使うモデルは単純なポアソン回帰 (R の glm 関数を使う)
#  -- 
 
library(mvtnorm)
library(AICcmodavg)
 
generate.data <- function(n){
  X <- rmvnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),nc=2))
  mu <- exp(0.5*X[,1]+1)
  Y <- rpois(n,mu)
  d <- as.data.frame(cbind(X,Y))
  colnames(d) <- c("x1","x2","y")
  d
}
 
#------------------
# estimating GLM
#------------------
 
estimate.models <- function(dat,wgt=rep(1,dim(dat)[1])){
  models <- list()
 
  # Model 0
  models[[1]] <- glm(y~1,data=dat,family=poisson,weights=wgt)
 
  # Model 1
  models[[2]] <- glm(y~x1,data=dat,family=poisson,weights=wgt)
 
  # Model 2
  models[[3]] <- glm(y~x2,data=dat,family=poisson,weights=wgt)
 
  # Model 3
  models[[4]] <- glm(y~x1+x2,data=dat,family=poisson,weights=wgt)
 
  models
}
 
#--------------------------
# Model selection by AIC
#--------------------------
model.select.aic <- function(d){
  models <- estimate.models(d)
  which.min(sapply(models,AIC))  # Which model has minimum AIC?
}
 
#---------------------------
# Model selection by AICc
#---------------------------
model.select.aicc <- function(d){
  models <- estimate.models(d)
  # AICcmodavg package provides AICc!!
  which.min(sapply(models,AICc))  # Which model has minimum AICc?
}

#----------------------------
# Experiment / AIC vs. AICc
#----------------------------
library(tabplot)

ns <- 20*(1:5)
correct <- matrix(0,length(ns),2)
for( j in 1:length(ns) ){
  n <- ns[j]
  n.experiments <- 10000
  results <- matrix(0,n.experiments,5)
  colnames(results) <- c("AIC.selected","AICc.selected",
                         "p.val.Intersept","p.val.x1","p.val.x2")
  for(i in 1:n.experiments){
    d <- generate.data(n)
    p.values <- summary(estimate.models(d)[[4]])$coef[,4]
    results[i,] <- c(model.select.aic(d),model.select.aicc(d),p.values)
  }
  
  results.df <- transform(data.frame(results),
                          AIC.selected=factor(AIC.selected,levels=1:4),
                          AICc.selected=factor(AICc.selected,levels=1:4))
  levels(results.df[,1]) <- paste("Model",levels(results.df[,1]))
  levels(results.df[,2]) <- paste("Model",levels(results.df[,2]))

  correct[j,] <- colSums(results.df[,1:2]=="Model 2")
  
  tblobj <- tableplot(results.df,sortCol=5,nBins=100,scales="lin")
  plot(tblobj,title=paste("Data size =",n),showTitle=T,fontsize=12)
  dev2bitmap(paste("AicVsAICc_n_",n,".jpg",sep=""),
             width=10,height=10,gaa=4,taa=4)
}

correct.df <- data.frame(N=c(ns,ns),
                         CorrectRatio=c(correct[,1],correct[,2])/10000,
                         Method=rep(c("AIC","AICc"),each=length(ns)))
ggplot(aes(x=N,y=CorrectRatio),data=correct.df,colour=Method) +
  geom_line(aes(colour=Method),size=2) +
  geom_point(aes(colour=Method),size=5) +
  xlab("\nNumber of data used for estimation") +
  ylab("Correct ratio\n") +
  opts(legend.position=c(0.8,0.2),
       legend.background=theme_rect(fill="white"))
ggsave("AICvsAICc.png",height=4,width=4)
はてなブックマーク - モデル選択の実験 (AIC vs AICc / R の AICcmodavg パッケージ)
Pocket

  1. [1] でも実は間違えてこの定式化でやってみたけどそれなりに良い結果 (AIC よりもよい結果) が得られた。
  • 竹澤

    AICc は正規分布を仮定して導出されたものなので、ポアソン分布に対して優れた
    結果を与えるかどうかははっきりしません。この点に対して、久保さんのように、
    正規分布の場合以外はAICcを使うべきではない、と考えている人もいます。
    一方、正規分布以外のときも経験的にAICcがうまく機能することが多いことが
    分かっているのだから、AICcを使うべきと考えている人(We believe there are compelling reasons to use AICc and QAICc as effective methods for general use. http://warnercnr.colostate.edu/~anderson/PDF_files/AIC Myths and Misunderstandings.pdf)もいます。