重回帰法の落とし穴

MikuHatsune2013-08-01

数学いらずの医科統計学PART7 CHAPTER38で、重回帰分析における過剰適合について説明している。
パラメータ取ったんならたくさん使って推定しようぜ!!というのはよくわかる気持ちだが、使い過ぎるとパラメタ間の相互作用とかなんかいろいろな都合でダメになることが多々あるらしい。
パラメタがkあるなら頑張って2^k総当りを考えてもいいけど、組み合わせ爆発が起きてアボーン…
いい感じに計算量を減らす方法として、ステップワイズ法というものがあり、1回ごとに説明能力の高いパラメタを足していく変数増加法と、説明能力の低いパラメタを減らしていく変数減少法がある。
 
過剰適合の例として、A Note on Screening Regression Equationsのシミュレーションが載っている。
いま、適当な変数xを50個、100サンプルについてランダムサンプリングで作成する。そして、これまた適当に推定する変数yをランダムサンプリングで作成する。こうすると、xyの関係はまったくのランダムである。
この状態で、50個のパラメタxからyを予測する。すると、重回帰のp値は、完全に一様分布になる。
1回重回帰を行うと、50個のパラメタxについてそれぞれ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
... 以下略