R語言——多元線性回歸

1、多元線性回歸模型

1.1多元回歸模型與多元回歸方程

設因變量為y，k個自變量分別為，描述因變量y如何依賴於自變量和誤差項ε的方程稱為多元回歸模型。其一般形式可表示為：

式中，為模型的參數，ε為隨機誤差項。

上式表明，y是的線性函數加上隨機誤差項ε。隨機誤差項的解釋見：隨機誤差項。

與一元線性回歸類似，在多元線性回歸模型中，對誤差項同樣有三個基本假設：

1. 誤差項期望為0；
2. 對於自變量的所有值，ε的值都相等；
3. 誤差項ε是一個服從正態分布的隨機變量，且相互獨立。

根據回歸模型的假定，有

，

上式被稱為多元回歸方程，它描述了因變量y的期望值與自變量的關系。

 1.2估計的多元回歸方程

回歸方程中的參數是未知的，正是我們感興趣的值。因此，當用樣本數據計算出來的來去估計未知參數時，就得到了估計的多元回歸方程，其一般形式為：

式中，是參數的估計值，表示當不變時，每變動一個單位，y的平均變動量。其余偏回歸系數含義類似。

1.3 參數的最小二乘估計

同一元線性回歸的最小二乘法估計，使得殘差平方和達到最小的即為參數的最小二乘估計。

由於多元回歸參數的最小二乘參數估計計算量比較大，手工幾乎無法完成，需要借助於計算機。

本文選取了UCI數據集中一個關於計算機CPU性能的數據，該數據可以從網站http://archive.ics.uci.edu/ml/machine-learning-databases/cpu-performance/上下載。

在R中讀取該數據，先對該數據集的數據結構進行探究。

> setwd("D:/Rdata/complex_data_analysis/")
> cpu<-read.table("cpu.txt",header = TRUE,sep=",")
> str(cpu)
'data.frame':	209 obs. of  10 variables:
 $ name : Factor w/ 30 levels "adviser","amdahl",..: 1 2 2 2 2 2 2 2 2 2 ...
 $ Model: Factor w/ 209 levels "100","1100/61-h1",..: 30 63 64 65 66 67 75 76 77 78 ...
 $ MYCT : int  125 29 29 29 29 26 23 23 23 23 ...
 $ MMIN : int  256 8000 8000 8000 8000 8000 16000 16000 16000 32000 ...
 $ MMAX : int  6000 32000 32000 32000 16000 32000 32000 32000 64000 64000 ...
 $ CACH : int  256 32 32 32 32 64 64 64 64 128 ...
 $ CHMIN: int  16 8 8 8 8 8 16 16 16 32 ...
 $ CHMAX: int  128 32 32 32 16 32 32 32 32 64 ...
 $ PRP  : int  198 269 220 172 132 318 367 489 636 1144 ...
 $ ERP  : int  199 253 253 253 132 290 381 381 749 1238 ...

　　

根據R的輸出結果可以看出，該數據有209個觀測值，共10個變量。其變量名的細節可以從網站中關於該數據的描述中獲得。

Name：供應商名稱，

Model：型號名稱

MYCT：機器周期時間，以納秒為單位

MMIN：最小主內存千字節

MMAX：最大主內存千字節

CACH：高速緩存，以千字節為單位

CHMIN：單位最小信道

CHMAX：單位最大信道

PRP：發表的相關性能

ERP：根據原始文章估計的相對性能

由於供應商和型號名稱屬於分類數據，另外ERP是數據提供者根據他的一篇文章預測出來的，這三個變量都需要移除。

>cpu<-cpu[,c(-1,-2,-10)]
> str(cpu)  #查看需要進行分析的數據的結構
'data.frame':	209 obs. of  7 variables:
 $ MYCT : int  125 29 29 29 29 26 23 23 23 23 ...
 $ MMIN : int  256 8000 8000 8000 8000 8000 16000 16000 16000 32000 ...
 $ MMAX : int  6000 32000 32000 32000 16000 32000 32000 32000 64000 64000 ...
 $ CACH : int  256 32 32 32 32 64 64 64 64 128 ...
 $ CHMIN: int  16 8 8 8 8 8 16 16 16 32 ...
 $ CHMAX: int  128 32 32 32 16 32 32 32 32 64 ...
 $ PRP  : int  198 269 220 172 132 318 367 489 636 1144 ...
> head(cpu)  #i查看前幾行數據特征
  MYCT MMIN  MMAX CACH CHMIN CHMAX PRP
1  125  256  6000  256    16   128 198
2   29 8000 32000   32     8    32 269
3   29 8000 32000   32     8    32 220
4   29 8000 32000   32     8    32 172
5   29 8000 16000   32     8    16 132
6   26 8000 32000   64     8    32 318
> lm.cpu<-lm(PRP~.,data = cpu)
> lm.cpu$coefficients  #輸出估計的回歸系數
 (Intercept)         MYCT         	MMIN         MMAX         CACH        	CHMIN 		CHMAX 
-55.89393361   0.04885490   0.01529257   0.00557139   0.64140143  -0.27035755 	1.48247217

　　

與一元線性回歸類似，R中也采用lm()函數進行多元回歸分析，其中"PRP~."表示PRP為因變量，其他所有變量為自變量，data為所使用的數據框。

輸出所有的系數——參數的最小二乘估計。

2、回歸方程的擬合度

2.1 多重判定系數

與一元線性回歸類似，對多元線性回歸方程，需要用多重判定系數來評價其擬合程度。其公式與一元線性回歸的判定系數公式相同，即：

。

但是對於多重判定系數，有一點需要注意，即自變量個數的增加將影響到因變量中被估計的回歸方程所解釋的變差數量。當增加自變量時，會使預測誤差變得較小，從而減少殘差平方和SSE，因此會增大判定系數R²。如果模型中增加一個自變量，即使這個自變量在統計上並不顯著，R²也會變大。因此，為了避免增加自變量的個數而高估R²，可以使用自變量個數k以及樣本量n去調整R²，得到估計的R²，其公式為：

，調整以后的永遠小於R²，也不會由於模型中自變量的個數的增加越來越接近1.。

R²的平方根稱為多重相關系數，也稱為復相關系數。在上面的結果中，可以輸出R²。

> summary(lm.cpu)

Call:
lm(formula = PRP ~ ., data = cpu)

Residuals:
    Min      1Q  Median      3Q     Max 
-195.82  -25.17    5.40   26.52  385.75 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.589e+01  8.045e+00  -6.948 5.00e-11 ***
MYCT         4.885e-02  1.752e-02   2.789   0.0058 ** 
MMIN         1.529e-02  1.827e-03   8.371 9.42e-15 ***
MMAX         5.571e-03  6.418e-04   8.681 1.32e-15 ***
CACH         6.414e-01  1.396e-01   4.596 7.59e-06 ***
CHMIN       -2.704e-01  8.557e-01  -0.316   0.7524    
CHMAX        1.482e+00  2.200e-01   6.737 1.65e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 59.99 on 202 degrees of freedom
Multiple R-squared:  0.8649,	Adjusted R-squared:  0.8609 
F-statistic: 215.5 on 6 and 202 DF,  p-value: < 2.2e-16

　　

上式輸出了R²為0.8649，調整的R²為0.8609，說明模型擬合地還不錯。

3、顯著性檢驗

在一元線性回歸中，我們已經學會了對回歸模型進行線性關系檢驗與回歸系數檢驗，這里同樣也需要進行這兩步。但是在這里，線性關系的檢驗與回歸系數的檢驗並不等價。在多元線性回歸模型中，只要有一個變量與因變量的線性關系顯著，F檢驗就能通過，但這不意味着所有的變量都與因變量的線性關系顯著。

3.1 線性關系檢驗

第一步：提出假設。

至少有一個不等於0

第二步：計算檢驗統計量F。

第三步：做出統計決策。

根據F的值與查F分布表得出的Fα，做出決策。

根據上述結果，p值為2.2e-16，幾乎為0，遠遠 小於常用的α=0.05，因此線性關系顯著。

3.2 回歸系數的檢驗

第一步：提出假設。

第二步：計算檢驗統計量t。

，式中，是回歸系數的抽樣分布的標准差，即

第三步：做出統計決策。

根據上述輸出結果，只有一個回歸系數——MMIN不顯著，其他都比較顯著。

4、多重共線性

4.1 多重共線性及其判別

當模型中兩個自變量之間高度相關時，則稱模型中存在多重共線性。判斷模型的多重共線性最直觀的方法是求自變量的相關系數矩陣。當然，還可以采取以下幾種方法判斷模型中是否存在多重共線性。

1. 模型的線性關系檢驗（F檢驗）顯著時，幾乎所有的回歸系數檢驗（t檢驗）都不顯著；
2. 回歸系數的正負號與預期的相反；
3. 容忍度與方差擴大因子。

某個自變量的容忍度等於1減去以該變量為因變量，其它變量為自變量所擬合的模型的R²,即1- R²。通常認為容忍度小於0.1時，存在嚴重的多重共線性。方差擴大因子（VIF）為容忍度的倒數，當VIF大於10時，一般可認為存在嚴重的多重共線性。下面求出自變量的相關系數矩陣。

> cor(cpu[,1:6])
            MYCT       MMIN       MMAX       CACH      CHMIN      CHMAX
MYCT   1.0000000 -0.3356422 -0.3785606 -0.3209998 -0.3010897 -0.2505023
MMIN  -0.3356422  1.0000000  0.7581573  0.5347291  0.5171892  0.2669074
MMAX  -0.3785606  0.7581573  1.0000000  0.5379898  0.5605134  0.5272462
CACH  -0.3209998  0.5347291  0.5379898  1.0000000  0.5822455  0.4878458
CHMIN -0.3010897  0.5171892  0.5605134  0.5822455  1.0000000  0.5482812
CHMAX -0.2505023  0.2669074  0.5272462  0.4878458  0.5482812  1.0000000

　　

可以看出，MMIN和MMAX的相關系數為0.76，比較大，另外，CHMIN與其他變量（除了CHMAX）的相關系數絕對值也超過了0.5，初步認為存在多重共線性。

4.2 多重共線性的處理

（1）將一個或多個自變量從模型中剔除，使保留的自變量之間盡可能不相關。

（2）如果非要保留模型中所有的自變量，那就應該：

• 避免根據t統計量對單個參數進行檢驗；
• 對因變量的推斷限定在自變量樣本的范圍內。

5、預測

R中可以根據估計的回歸方程進行預測，期預測結果可以表示為：

> predict.lm(lm.cpu,interval = "prediction")[1:6,]
       fit       lwr      upr
1 337.1856 201.45329 472.9180
2 311.9490 192.27143 431.6266
3 311.9490 192.27143 431.6266
4 311.9490 192.27143 431.6266
5 199.0872  79.60871 318.5657
6 332.3273 212.77407 451.8805

　　

上面求出預測區間以及相應的估計值。

6、變量選擇與逐步回歸

當模型中存在多重共線性時，可以使用以下幾種方法進行處理：

1. 向前選擇
2. 向后剔除
3. 逐步回歸

R中提供了step()函數實現這個功能。

> stpcpu<-step(lm.cpu)
Start:  AIC=1718.24
PRP ~ MYCT + MMIN + MMAX + CACH + CHMIN + CHMAX

        Df Sum of Sq    RSS    AIC
- CHMIN  1       359 727279 1716.3
<none>               726920 1718.2
- MYCT   1     27985 754905 1724.1
- CACH   1     76009 802929 1737.0
- CHMAX  1    163347 890267 1758.6
- MMIN   1    252179 979099 1778.5
- MMAX   1    271177 998097 1782.5

Step:  AIC=1716.34
PRP ~ MYCT + MMIN + MMAX + CACH + CHMAX

        Df Sum of Sq    RSS    AIC
<none>               727279 1716.3
- MYCT   1     28343 755623 1722.3
- CACH   1     78715 805995 1735.8
- CHMAX  1    177114 904393 1759.9
- MMIN   1    258252 985531 1777.8
- MMAX   1    270856 998135 1780.5

　　

該函數根據AIC信息准則對變量進行選擇，取使AIC最小的回歸模型。最終結果保留了5個自變量，CHMIN變量被移除。

> summary(stpcpu)  #輸出變量選擇后的最終結果

Call:
lm(formula = PRP ~ MYCT + MMIN + MMAX + CACH + CHMAX, data = cpu)

Residuals:
    Min      1Q  Median      3Q     Max 
-193.37  -24.95    5.76   26.64  389.66 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.608e+01  8.007e+00  -7.003 3.59e-11 ***
MYCT         4.911e-02  1.746e-02   2.813   0.0054 ** 
MMIN         1.518e-02  1.788e-03   8.490 4.34e-15 ***
MMAX         5.562e-03  6.396e-04   8.695 1.18e-15 ***
CACH         6.298e-01  1.344e-01   4.687 5.07e-06 ***
CHMAX        1.460e+00  2.076e-01   7.031 3.06e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 59.86 on 203 degrees of freedom
Multiple R-squared:  0.8648,	Adjusted R-squared:  0.8615 
F-statistic: 259.7 on 5 and 203 DF,  p-value: < 2.2e-16

　　

處理后的R²幾乎沒有變化，但是調整的R²略微增加。