R語言線性回歸之逐步回歸

“最優“回歸方程的選擇”
當變量中含有對Y影響不大的變量時，可能因為誤差平方和的自由度減小而使方差的估計增大，從而影響回歸預測的精度，適當的選擇一個變量建立一個最優的回歸方程十重要。此處采用逐步回歸法。

逐步回歸法計算
#水泥熱量與四種成分的關系
cement<-data.frame(
    X1=c( 7,  1, 11, 11,  7, 11,  3,  1,  2, 21,  1, 11, 10),
    X2=c(26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68),
    X3=c( 6, 15,  8,  8,  6,  9, 17, 22, 18,  4, 23,  9,  8),
    X4=c(60, 52, 20, 47, 33, 22,  6, 44, 22, 26, 34, 12, 12),
    Y =c(78.5, 74.3, 104.3,  87.6,  95.9, 109.2, 102.7, 72.5,
         93.1,115.9,  83.8, 113.3, 109.4)
)
lm.sol<-lm(Y ~ X1+X2+X3+X4, data=cement)
summary(lm.sol)

Call:
lm(formula = Y ~ X1 + X2 + X3 + X4, data = cement)

Residuals:
    Min      1Q  Median      3Q     Max
-3.1750 -1.6709  0.2508  1.3783  3.9254

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  62.4054    70.0710   0.891   0.3991  
X1            1.5511     0.7448   2.083   0.0708 .
X2            0.5102     0.7238   0.705   0.5009  
X3            0.1019     0.7547   0.135   0.8959  
X4           -0.1441     0.7091  -0.203   0.8441  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.446 on 8 degrees of freedom
Multiple R-squared:  0.9824,    Adjusted R-squared:  0.9736
F-statistic: 111.5 on 4 and 8 DF,  p-value: 4.756e-07
#回歸系數沒有一項通過檢測
#用step( )做回歸分析
lm.ste<-step(lm.sol)
Start:  AIC=26.94
Y ~ X1 + X2 + X3 + X4

       Df Sum of Sq    RSS    AIC
- X3    1    0.1091 47.973 24.974
- X4    1    0.2470 48.111 25.011
- X2    1    2.9725 50.836 25.728
<none>              47.864 26.944
- X1    1   25.9509 73.815 30.576

Step:  AIC=24.97
Y ~ X1 + X2 + X4

       Df Sum of Sq    RSS    AIC
<none>               47.97 24.974
- X4    1      9.93  57.90 25.420
- X2    1     26.79  74.76 28.742
- X1    1    820.91 868.88 60.629
用全部變量做回歸分析時，AIC值為26.94。接下來顯示如果去除X3，則AIC = 24.97,去掉x4則為25.01，去掉X3可以使AIC達到最小，R軟件自動去掉x3.

summary(lm.ste)
Call:
lm(formula = Y ~ X1 + X2 + X4, data = cement)

Residuals:
    Min      1Q  Median      3Q     Max
-3.0919 -1.8016  0.2562  1.2818  3.8982

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  71.6483    14.1424   5.066 0.000675 ***
X1            1.4519     0.1170  12.410 5.78e-07 ***
X2            0.4161     0.1856   2.242 0.051687 .  
X4           -0.2365     0.1733  -1.365 0.205395    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.309 on 9 degrees of freedom
Multiple R-squared:  0.9823,    Adjusted R-squared:  0.9764
F-statistic: 166.8 on 3 and 9 DF,  p-value: 3.323e-08

回歸系數顯著性水平有大提高，但是X2,X2系數檢驗不理想。從step（）可以看出去掉x4，AIC從24.97變為25.42，是增加的最少的，除AIC準則外，殘差平方各也是逐步回歸的重要指標之一，從直觀上看，擬合越好的直線，殘差平方和應該最小，去掉x4后，殘差平方和上升了9.93，也是最少的。從這兩項指標看，應該去掉x4.

lm.opt <- lm(Y ~ X1+X2,data = cement);
summary(lm.opt)

Call:
lm(formula = Y ~ X1 + X2, data = cement)

Residuals:
   Min     1Q Median     3Q    Max
-2.893 -1.574 -1.302  1.363  4.048

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 52.57735    2.28617   23.00 5.46e-10 ***
X1           1.46831    0.12130   12.11 2.69e-07 ***
X2           0.66225    0.04585   14.44 5.03e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.406 on 10 degrees of freedom
Multiple R-squared:  0.9787,    Adjusted R-squared:  0.9744
F-statistic: 229.5 on 2 and 10 DF,  p-value: 4.407e-09

這個結果還算滿意

Y = 52.58 + 1.46831*X1 + 0.66225*X2

改變step( )中的某些參數，可能得到不同的結果。

lm.ste<-step(lm.sol, trace=0, k=3); lm.ste

Call:
lm(formula = Y ~ X1 + X2, data = cement)

Coefficients:
(Intercept)           X1           X2  
    52.5773       1.4683       0.663  
直接去掉X3和X4。
從增加變量的角度考慮逐步回歸

lm0<-lm(Y~1, data=cement)
lm.ste<-step(lm0, scope = ~X1+X2+X3+X4, k=4)
Start:  AIC=73.44
Y ~ 1

       Df Sum of Sq     RSS    AIC
+ X4    1   1831.90  883.87 62.852
+ X2    1   1809.43  906.34 63.178
+ X1    1   1450.08 1265.69 67.519
+ X3    1    776.36 1939.40 73.067
<none>              2715.76 73.444

Step:  AIC=62.85
Y ~ X4

       Df Sum of Sq     RSS    AIC
+ X1    1    809.10   74.76 34.742
+ X3    1    708.13  175.74 45.853
<none>               883.87 62.852
+ X2    1     14.99  868.88 66.629
- X4    1   1831.90 2715.76 73.444

Step:  AIC=34.74
Y ~ X4 + X1

       Df Sum of Sq     RSS    AIC
+ X2    1     26.79   47.97 32.974
+ X3    1     23.93   50.84 33.728
<none>                74.76 34.742
- X1    1    809.10  883.87 62.852
- X4    1   1190.92 1265.69 67.519

Step:  AIC=32.97
Y ~ X4 + X1 + X2

       Df Sum of Sq    RSS    AIC
- X4    1      9.93  57.90 31.420
<none>               47.97 32.974
- X2    1     26.79  74.76 34.742
+ X3    1      0.11  47.86 36.944
- X1    1    820.91 868.88 66.629

Step:  AIC=31.42
Y ~ X1 + X2

       Df Sum of Sq     RSS    AIC
<none>                57.90 31.420
+ X4    1      9.93   47.97 32.974
+ X3    1      9.79   48.11 33.011
- X1    1    848.43  906.34 63.178
- X2    1   1207.78 1265.69 67.519
這里取k4，最后還剩下x1與x2。
在R中，還有兩個函數可以做逐步回歸，一個是add1( )函數，用于增加變量，一個是drop1( )函數，用于減小變量。事實上，step( )就是使用這兩個函數來自動增加和減小變量。

add1(lm0, scope = ~X1+X2+X3+X4, test="F")
Single term additions

Model:
Y ~ 1
       Df Sum of Sq     RSS    AIC F value    Pr(>F)    
<none>              2715.76 71.444                      
X1      1   1450.08 1265.69 63.519 12.6025 0.0045520 **
X2      1   1809.43  906.34 59.178 21.9606 0.0006648 ***
X3      1    776.36 1939.40 69.067  4.4034 0.0597623 .  
X4      1   1831.90  883.87 58.852 22.7985 0.0005762 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
打算減少變量

drop1(lm.sol, test="F")
Single term deletions

Model:
Y ~ X1 + X2 + X3 + X4
       Df Sum of Sq    RSS    AIC F value  Pr(>F)  
<none>              47.864 26.944                  
X1      1   25.9509 73.815 30.576  4.3375 0.07082 .
X2      1    2.9725 50.836 25.728  0.4968 0.50090  
X3      1    0.1091 47.973 24.974  0.0182 0.89592  
X4      1    0.2470 48.111 25.011  0.0413 0.84407  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
根據每步計算的結果情況，人工選擇增加還是去掉某些變量。