GLM関数群のAIC算出法まとめ

East_scrofaさんのところでまとめられているように、GLMおよびGLMMの関数群にはAICを求める対数尤度の基準が２通りあるようで、違う求め方で求めたAICは単純に比べられない。そこで、よくつかわれるGLM, GLMMの関数群のAICの求め方を補足としてまとめる。

従来のAICの求め方をしている関数群。
glm(), glm.nb(MASSパッケージ), negbin(aodパッケージ), betabin(aodパッケージ)

「residual devianceを計算する基準となっているモデルの対数尤度を0として出したもの」
lmer(lme4パッケージ), glmmML(glmmMLパッケージ)

下にnegbin, betabin, glm.nbの検証プログラムを載せておく。

#negbin
n1<-100
a1<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-factor(seq(6))
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-negbin(y1~1+x1,random=~ID1,data=tdata1)
a1[i1]<- AIC(res1)@istats$AIC
rd1[i1]<-deviance(res1)
}
a1-rd1 #AICとresidual devianceの差

#betabin
n1<-100
a1<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-factor(seq(6))
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
y2<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-betabin(cbind(y1,y2)~1+x1,random=~ID1,data=tdata1)
a1[i1]<-AIC(res1)@istats$AIC
rd1[i1]<-deviance(res1)
}

a1-rd1 #AICとresidual devianceの差

#glm.nb
n1<-100
a1<-numeric(n1)
rd1<-numeric(n1)
x1<-seq(6)
ID1<-factor(seq(6))
for (i1 in 1:n1){
y1<-floor(runif(6,min=1,max=20))
tdata1<-data.frame(x1,ID1,y1,y2)
res1<-glm.nb(y1~1+x1,data=tdata1)
a1[i1]<- res1$aic
rd1[i1]<-res1$deviance
}
a1-rd1 #AICとresidual devianceの差