数学いらずの医科統計学PART7 CHAPTER38で、重回帰分析における過剰適合について説明している。
パラメータ取ったんならたくさん使って推定しようぜ!!というのはよくわかる気持ちだが、使い過ぎるとパラメタ間の相互作用とかなんかいろいろな都合でダメになることが多々あるらしい。
パラメタがあるなら頑張って
の総当りを考えてもいいけど、組み合わせ爆発が起きてｱﾎﾞｰﾝ…
いい感じに計算量を減らす方法として、ステップワイズ法というものがあり、1回ごとに説明能力の高いパラメタを足していく変数増加法と、説明能力の低いパラメタを減らしていく変数減少法がある。
過剰適合の例として、A Note on Screening Regression Equationsのシミュレーションが載っている。
いま、適当な変数を50個、100サンプルについてランダムサンプリングで作成する。そして、これまた適当に推定する変数
をランダムサンプリングで作成する。こうすると、
と
の関係はまったくのランダムである。
この状態で、50個のパラメタから
を予測する。すると、重回帰のp値は、完全に一様分布になる。
1回重回帰を行うと、50個のパラメタについてそれぞれp値が出る。この中で、小さいp値15個を採用して、これらだけでまた重回帰を再度行う。
するとあら不思議、p値は小さくなってしまうのでした。
n <- 100 # サンプル数 X <- 50 # パラメタの数 y <- rnorm(n) # 推定されるパラメタを適当に発生 niter <- 2000 # シミュレーション回数 res1 <- res2 <- numeric(niter) for(i in seq(niter)){ # 全パラメタを使って重回帰を行う x <- replicate(X, runif(n, min=0, max=10)) # パラメタ群の作成 g0 <- glm(y ~ x) sg0 <- summary(aov(g0)) res1[i] <- sg0[[1]][1, 5] # ロジスティック回帰全体のp値。多分。 # 選抜された小さいp値を使って重回帰を再度行う。 p0 <- summary(g0)$coefficients[, "Pr(>|t|)"][-1] # Intercept を除いたパラメータ群のp値 idx <- head(order(p0), 15) g1 <- glm(y ~ x[,idx]) p1 <- summary(g1)$coefficients[, "Pr(>|t|)"][-1] res2[i] <- summary(aov(g1))[[1]][1, 5] } yl <- c(0, 1) plot(sort(res1), xlab="niter", ylab="p value of logistic regression", ylim=yl)
50個のパラメータについてp値を見ると、いくつかは0.05を下回る。重回帰全体のp値は一様分布に従っている。今回は0.275だった。
summary(g0) summary(aov(g0))
Call: glm(formula = y ~ x) Deviance Residuals: Min 1Q Median 3Q Max -1.75609 -0.44399 -0.04823 0.46249 1.86142 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.695361 1.520671 0.457 0.64950 x1 0.046170 0.052140 0.885 0.38021 x2 -0.018032 0.047290 -0.381 0.70463 x3 0.086522 0.044042 1.965 0.05515 . x4 0.036942 0.048609 0.760 0.45091 x5 -0.007975 0.049185 -0.162 0.87186 x6 0.024968 0.057504 0.434 0.66606 x7 -0.065707 0.054875 -1.197 0.23691 x8 -0.086664 0.042464 -2.041 0.04667 * x9 -0.005942 0.047892 -0.124 0.90176 x10 -0.015981 0.050971 -0.314 0.75520 x11 -0.017552 0.047747 -0.368 0.71475 x12 0.051322 0.048781 1.052 0.29791 x13 -0.032189 0.047477 -0.678 0.50096 x14 -0.026207 0.049267 -0.532 0.59717 x15 -0.010495 0.051922 -0.202 0.84064 x16 0.060763 0.044115 1.377 0.17466 x17 0.068961 0.044256 1.558 0.12562 x18 -0.034572 0.045655 -0.757 0.45254 x19 0.070264 0.046895 1.498 0.14046 x20 -0.044588 0.051580 -0.864 0.39156 x21 0.024617 0.045483 0.541 0.59079 x22 0.044463 0.042775 1.039 0.30369 x23 0.087321 0.053239 1.640 0.10737 x24 -0.164941 0.049773 -3.314 0.00174 ** x25 0.021186 0.047068 0.450 0.65460 x26 0.053364 0.048168 1.108 0.27332 x27 -0.043749 0.045705 -0.957 0.34316 x28 0.051669 0.045800 1.128 0.26476 x29 -0.011417 0.046352 -0.246 0.80647 x30 0.097146 0.047152 2.060 0.04470 * x31 -0.032972 0.043417 -0.759 0.45123 x32 0.028890 0.045448 0.636 0.52794 x33 -0.022130 0.044069 -0.502 0.61780 x34 -0.098116 0.043913 -2.234 0.03006 * x35 0.010344 0.048201 0.215 0.83097 x36 -0.015528 0.051402 -0.302 0.76386 x37 -0.055377 0.058149 -0.952 0.34561 x38 -0.091507 0.047597 -1.923 0.06036 . x39 -0.041357 0.045481 -0.909 0.36764 x40 0.054348 0.043034 1.263 0.21260 x41 0.011998 0.046589 0.258 0.79785 x42 -0.044907 0.045985 -0.977 0.33359 x43 0.100701 0.062550 1.610 0.11384 x44 -0.095154 0.051791 -1.837 0.07224 . x45 -0.024633 0.046039 -0.535 0.59504 x46 -0.004963 0.043606 -0.114 0.90985 x47 0.016534 0.044107 0.375 0.70938 x48 -0.070204 0.049319 -1.423 0.16094 x49 -0.070327 0.051151 -1.375 0.17542 x50 0.069107 0.041615 1.661 0.10317 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 0.9408158) Null deviance: 101.94 on 99 degrees of freedom Residual deviance: 46.10 on 49 degrees of freedom AIC: 310.35 Number of Fisher Scoring iterations: 2
Df Sum Sq Mean Sq F value Pr(>F) x 50 55.84 1.1167 1.187 0.275 Residuals 49 46.10 0.9408
選抜されたパラメタでは、15個のパラメタでさらに小さいp値が出る傾向にあり、重回帰全体でも0.00033とさらに小さなp値になった。
summary(g1) summary(aov(g1))
Call: glm(formula = y ~ x[, idx]) Deviance Residuals: Min 1Q Median 3Q Max -1.88278 -0.45673 -0.03604 0.44589 2.00477 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.32692 0.58424 0.560 0.57726 x[, idx]1 -0.15765 0.03574 -4.411 3.03e-05 *** x[, idx]2 -0.05452 0.03166 -1.722 0.08868 . x[, idx]3 0.06987 0.03150 2.218 0.02926 * x[, idx]4 -0.06702 0.03255 -2.059 0.04260 * x[, idx]5 0.03643 0.03116 1.169 0.24564 x[, idx]6 -0.06498 0.03613 -1.798 0.07569 . x[, idx]7 -0.08877 0.03401 -2.610 0.01071 * x[, idx]8 0.05185 0.03155 1.644 0.10398 x[, idx]9 0.06842 0.03532 1.937 0.05608 . x[, idx]10 0.08551 0.03707 2.307 0.02353 * x[, idx]11 0.05927 0.03084 1.922 0.05803 . x[, idx]12 0.06942 0.03317 2.093 0.03939 * x[, idx]13 -0.08868 0.03129 -2.834 0.00575 ** x[, idx]14 0.07796 0.03302 2.361 0.02056 * x[, idx]15 -0.03670 0.03352 -1.095 0.27675 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 0.7708218) Null deviance: 101.936 on 99 degrees of freedom Residual deviance: 64.749 on 84 degrees of freedom AIC: 274.32 Number of Fisher Scoring iterations: 2
Df Sum Sq Mean Sq F value Pr(>F) x[, idx] 15 37.19 2.4791 3.216 0.00033 *** Residuals 84 64.75 0.7708 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
plot(res1, res2, xlab="Original", ylab="Selected", xlim=c(0, 1), ylim=c(0, 1))
ググっていたら見つけた epicalc パッケージ。オッズ比にしてくれるっぽい。
library(epicalc) logistic.display(g0)
OR lower95ci upper95ci Pr(>|Z|) x1 1.0472521 0.9455177 1.1599327 0.380214365 x2 0.9821300 0.8951904 1.0775132 0.704632755 x3 1.0903755 1.0002014 1.1886794 0.055152306 x4 1.0376325 0.9433390 1.1413513 0.450906890 x5 0.9920569 0.9008871 1.0924531 0.871862956 x6 1.0252818 0.9160013 1.1475997 0.666055908 x7 0.9364049 0.8409179 1.0427344 0.236914763 x8 0.9169849 0.8437553 0.9965700 0.046669206 ... 以下略
