モデル選択と最適なパラメータ数

MikuHatsune2013-07-19

数学いらずの医科統計学PART7 CHAPTER35で、モデルの選択について言及している。
パラメータの数が増えると、各パラメータの信頼区間が大きくなると言っているが、正直よくわからなかった。
 
パラメータの最適数の推定法として、AICというものが有名。Rなら回帰オブジェクトに summary 関数を実行した時にときたま AIC が計算されて返ってくるし、求めたかったら AICcmodavg パッケージが使える。
AICcmodavg 内のデータセット cement を使って、最適なパラメータ設定を考える。
いま、cement$y をパラメータ x1, x2, x3, x4 の4つで回帰することを考える。
y \~ a_1x_1+a_2x_2+a_3x_3+a_4x_4+b
という感じで、最大4つのパラメータを使って回帰できるが、ぶっちゃけ4つ使う必要あるの?という問題。
これを AIC 最小となるパラメータの組み合わせを使って回帰する。
いま、4つのパラメータから少なくとも1つ使う組み合わせは、2^4-1=15通りあるから、これを強引にループ回して AIC 最小となる組み合わせを探す。
結果としては{x1, x2}のパラメータを使って説明するときが AIC 最小のようだ。

library(gplots)
library(AICcmodavg)
data(cement)

cmb <- unique(apply(combinations(4, 4, repeats.allowed=TRUE), 1, unique)) # 4つの要素の組み合わせでグループをつくる。
res <- rep(0, length(cmb)) # AIC を納める。
lms <- vector("list", length(res)) # lm の結果を納める。
mat <- matrix(0, nrow(cement), length(cmb)) # 回帰した y の値を納める。
for(i in seq(cmb)){
	x <- cmb[[i]]
	tmp <- cement[, c(x, 5)]
	lm0 <- lm(y ~ ., data = tmp)
	res[i] <- AICc(lm0, return.K = FALSE)
	lms[[i]] <- lm0
	cof <- lm0$coefficients
	mat[, i] <- rowSums(sweep(cbind(1, tmp[,-ncol(tmp)]), 2, cof, "*"))
}

# AIC の順位をプロットする。
group <- mapply(function(x) paste("{", x, "}", sep=""), lapply(cmb, paste, collapse=", "))
idx <- order(res)
plot(res[idx], xlab="", ylab="AIC", xaxt="n")
axis(1, seq(cmb), label=NA)
text(seq(cmb), par()$usr[3], group[idx], xpd=TRUE, srt=45, adj=c(1.3, 1.3))

# 回帰した値が真の y の値とどれくらい近そうかプロット。
cols <- rainbow(length(cmb))
matplot(cement$y, mat, xlim=c(min(mat), max(mat)+10), pch=seq(cmb), ylab="Regression", col=cols)
abline(0, 1, lty=2)
legend("bottomright", legend=paste("AIC =", round(res, 2)), bty="n", col=cols, pch=seq(cmb))