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