A way of thinking

筆者個人の思考過程です。意見には個人差があります。

Nakagawa, S., Schielzeth, H., 2013. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133-142.

以下は,作業用メモです。AICの値ではモデルの相対比較はできても,そのモデルがどれくらいデータを説明できていそうかはわからないので,AICで議論しているとモデルの当てはまりの良さを示してよ,的な査読コメントや意見を言われることがあるし,その気持ちはよく分かります。ということで,これに沿って*1,GLMM用にR2を計算してみました。しかし,同じパラメータ数のモデルでもAICとR2(GLMM)ではモデルの”良さ”のランキングが変わってしまいました。これでなぜなんでしょう?結構小さいAICの差で議論しているせいもあるのかもしれんが,結局逸脱度(尤度)をベースに計算するか,モデルで推定された”ばらつき”をもとに計算するか*2で,変わってくるってことなのでしょうか?それとも単なるボクがなんらかの凡ミスをしているということなんでしょうか?

と思って一応簡単な疑似データで試してみた。やはり,逆転が起こるみたい。

library(MuMIn)
options(na.action = "na.fail") # Necessary to run for using "dredge"

set.seed(99)
n <- 100

x1 <- rnorm(n, 0, 1)
x2 <- rnorm(n, 0, 1)
x3 <- rnorm(n, 0, 1)
x4 <- rnorm(n, 0, 1)
x5 <- rnorm(n, 0, 1)

yy <- rbinom(n, 50, 1/(1+exp(-8+x1*5+rnorm(n, 0, 0.5))))
d0 <- data.frame(yy, x1, x2, x3, x4, x5)
full.m <- glm(cbind(yy, 50-yy) ~ x1 + x2 + x3 + x4 + x5, binomial, data=d0)
# summary(full.m)

dredge (full.m, rank="AIC")

# AIC best model
m1<- glm(cbind(yy, 50-yy) ~ x1 + x2 , binomial)
r.squaredGLMM(m1) # 0.8196652

# AIC 2nd ranked
m2<- glm(cbind(yy, 50-yy) ~ x1 + x5, binomial)
r.squaredGLMM(m2) # 0.8256106

# AIC 3rd ranked
m3<- glm(cbind(yy, 50-yy) ~ x1 + x2 + x4, binomial)
r.squaredGLMM(m3) # 0.8280399

# AIC 4th ranked
m4<- glm(cbind(yy, 50-yy) ~ x1 + x4 + x5, binomial)
r.squaredGLMM(m4) # 0.8334658

# AIC 5th ranked
m5<- glm(cbind(yy, 50-yy) ~ x1 + x2 + x3, binomial)
r.squaredGLMM(m5) # 0.8214669

*1:正確にはMuMInにあるRの関数を使って

*2:この理解が正しいか,あまり自信はないです