轉載：http://blog.fens.me/r-multi-linear-regression/
前言
本文接上一篇R語言解讀一元線性回歸模型。在許多生活和工作的實際問題中，影響因變量的因素可能不止一個，比如對於知識水平越高的人，收入水平也越高，這樣的一個結論。這其中可能包括了因為更好的家庭條件，所以有了更好的教育；因為在一線城市發展，所以有了更好的工作機會；所處的行業趕上了大的經濟上行周期等。要想解讀這些規律，是復雜的、多維度的，多元回歸分析方法更適合解讀生活的規律。
由於本文為非統計的專業文章，所以當出現與教課書不符的描述，請以教課書為准。本文力求用簡化的語言，來介紹多元線性回歸的知識，同時配合R語言的實現。
目錄
- 多元線性回歸介紹
- 元線性回歸建模
- 模型優化
- 案例：黑色系期貨日K線數據驗證
1. 多元線性回歸介紹
對比一元線性回歸，多元線性回歸是用來確定2個或2個以上變量間關系的統計分析方法。多元線性回歸的基本的分析方法與一元線性回歸方法是類似的，我們首先需要對選取多元數據集並定義數學模型，然后進行參數估計，對估計出來的參數進行顯著性檢驗，殘差分析，異常點檢測，最后確定回歸方程進行模型預測。
由於多元回歸方程有多個自變量，區別於一元回歸方程，有一項很重要的操作就是自變量的優化，挑選出相關性最顯著的自變量，同時去除不顯著的自變量。在R語言中，有很方便地用於優化函數，可以很好的幫助我們來改進回歸模型。
下面就開始多元線性回歸的建模過程。
2. 多元線性回歸建模
做過商品期貨研究的人，都知道黑色系品種是具有產業鏈上下游的關系。鐵礦石是煉鋼的原材料，焦煤和焦炭是煉鋼的能源資源，熱卷即熱軋卷板是以板坯為原料經加熱后制成的鋼板，螺紋鋼是表面帶肋的鋼筋。
由於有產業鏈的關系，假設我們想要預測螺紋鋼的價格，那么影響螺紋鋼價格的因素可以會涉及到原材料，能源資源和同類材料等。比如，鐵礦石價格如果上漲，螺紋鋼就應該要漲價了。
2.1 數據集和數學模型
先從數據開始介紹，這次的數據集，我選擇的期貨黑色系的品種的商品期貨，包括了大連期貨交易所的 焦煤(JM)，焦炭(J)，鐵礦石(I)，上海期貨交易所的 螺紋鋼(RU) 和 熱卷(HC)。
數據集為2016年3月15日，當日白天開盤的交易數據，為黑色系的5個期貨合約的分鍾線的價格數據。
# 數據集已存在df變量中
> head(df,20)
x1 x2 x3 x4 y
2016-03-15 09:01:00 754.5 616.5 426.5 2215 2055
2016-03-15 09:02:00 752.5 614.5 423.5 2206 2048
2016-03-15 09:03:00 753.0 614.0 423.0 2199 2044
2016-03-15 09:04:00 752.5 613.0 422.5 2197 2040
2016-03-15 09:05:00 753.0 615.5 424.0 2198 2043
2016-03-15 09:06:00 752.5 614.5 422.0 2195 2040
2016-03-15 09:07:00 752.0 614.0 421.5 2193 2036
2016-03-15 09:08:00 753.0 615.0 422.5 2197 2043
2016-03-15 09:09:00 754.0 615.5 422.5 2197 2041
2016-03-15 09:10:00 754.5 615.5 423.0 2200 2044
2016-03-15 09:11:00 757.0 616.5 423.0 2201 2045
2016-03-15 09:12:00 756.0 615.5 423.0 2200 2044
2016-03-15 09:13:00 755.5 615.0 423.0 2197 2042
2016-03-15 09:14:00 755.5 615.0 423.0 2196 2042
2016-03-15 09:15:00 756.0 616.0 423.5 2200 2045
2016-03-15 09:16:00 757.5 616.0 424.0 2205 2052
2016-03-15 09:17:00 758.5 618.0 424.0 2204 2051
2016-03-15 09:18:00 759.5 618.5 424.0 2205 2053
2016-03-15 09:19:00 759.5 617.5 424.5 2206 2053
2016-03-15 09:20:00 758.5 617.5 423.5 2201 2050
數據集包括有6列:
- 索引, 為時間
- x1, 為焦炭(j1605)合約的1分鍾線的報價數據
- x2, 為焦煤(jm1605)合約的1分鍾線的報價數據
- x3, 為鐵礦石(i1605)合約的1分鍾線的報價數據
- x4, 為熱卷(hc1605)合約的1分鍾線的報價數據
- y, 為螺紋鋼(rb1605)合約的1分鍾線的報價數據
假設螺紋鋼的價格與其他4個商品的價格有線性關系，那么我們建立以螺紋鋼為因變量，以焦煤、焦炭、鐵礦石和熱卷的為自變量的多元線性回歸模型。用公式表示為：
y = a + b * x1 + c * x2 + d * x3 + e * x4 + ε- y，為因變量，螺紋鋼
- x1，為自變量，焦煤
- x2，為自變量，焦炭
- x3，為自變量，鐵礦石
- x4，為自變量，熱卷
- a，為截距
- b,c,d,e，為自變量系數
- ε, 為殘差，是其他一切不確定因素影響的總和，其值不可觀測。假定ε服從正態分布N(0,σ^2)。
通過對多元線性回歸模型的數學定義，接下來讓我們利用數據集做多元回歸模型的參數估計。
2.2. 回歸參數估計
上面公式中，回歸參數 a, b, c, d,e都是我們不知道的，參數估計就是通過數據來估計出這些參數，從而確定自變量和因變量之前的關系。我們的目標是要計算出一條直線，使直線上每個點的Y值和實際數據的Y值之差的平方和最小，即（Y1實際-Y1預測）^2+（Y2實際-Y2預測）^2+ …… +（Yn實際-Yn預測）^2 的值最小。參數估計時，我們只考慮Y隨X自變量的線性變化的部分，而殘差ε是不可觀測的，參數估計法並不需要考慮殘差。
類似於一元線性回歸，我們用R語言來實現對數據的回歸模型的參數估計，用lm()函數來實現多元線性回歸的建模過程。
# 建立多元線性回歸模型
> lm1<-lm(y~x1+x2+x3+x4,data=df)
# 打印參數估計的結果
> lm1
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Coefficients:
(Intercept) x1 x2 x3 x4
212.8780 0.8542 0.6672 -0.6674 0.4821
這樣我們就得到了y和x關系的方程。
y = 212.8780 + 0.8542 * x1 + 0.6672 * x2 - 0.6674 * x3 + 0.4821 * x42.3. 回歸方程的顯著性檢驗
參考一元線性回歸的顯著性檢驗，多元線性回歸的顯著性檢驗，同樣是需要經過 T檢驗，F檢驗，和R^2(R平方)相關系統檢驗。在R語言中這三種檢驗的方法都已被實現，我們只需要把結果解讀，我們可以summary()函數來提取模型的計算結果。
> summary(lm1)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.9648 -1.3241 -0.0319 1.2403 5.4194
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 212.87796 58.26788 3.653 0.000323 ***
x1 0.85423 0.10958 7.795 2.50e-13 ***
x2 0.66724 0.12938 5.157 5.57e-07 ***
x3 -0.66741 0.15421 -4.328 2.28e-05 ***
x4 0.48214 0.01959 24.609 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.028 on 221 degrees of freedom
Multiple R-squared: 0.9725, Adjusted R-squared: 0.972
F-statistic: 1956 on 4 and 221 DF, p-value: < 2.2e-16
- T檢驗：所自變量都是非常顯著***
- F檢驗：同樣是非常顯著，p-value < 2.2e-16
- 調整后的R^2：相關性非常強為0.972
最后，我們通過的回歸參數的檢驗與回歸方程的檢驗，得到最后多元線性回歸方程為：
y = 212.87796 + 0.85423 * x1 + 0.66724 * x2 - 0.66741 * x3 + 0.48214 * x4
即
螺紋鋼 = 212.87796 + 0.85423 * 焦炭 + 0.66724 * 焦煤 - 0.66741 * 鐵礦石 + 0.48214 * 熱卷
2.4 殘差分析和異常點檢測
在得到的回歸模型進行顯著性檢驗后，還要在做殘差分析(預測值和實際值之間的差)，檢驗模型的正確性，殘差必須服從正態分布N(0,σ^2)。直接用plot()函數生成4種用於模型診斷的圖形，進行直觀地分析。
> par(mfrow=c(2,2))
> plot(lm1)
- 殘差和擬合值(左上)，殘差和擬合值之間數據點均勻分布在y=0兩側，呈現出隨機的分布，紅色線呈現出一條平穩的曲線並沒有明顯的形狀特征。
- 殘差QQ圖(右上)，數據點按對角直線排列，趨於一條直線，並被對角直接穿過，直觀上符合正態分布。
- 標准化殘差平方根和擬合值(左下)，數據點均勻分布在y=0兩側，呈現出隨機的分布，紅色線呈現出一條平穩的曲線並沒有明顯的形狀特征。
- 標准化殘差和杠桿值(右下)，沒有出現紅色的等高線，則說明數據中沒有特別影響回歸結果的異常點。
結論，沒有明顯的異常點，殘差符合假設條件。
2.5. 模型預測
我們得到了多元線性回歸方程的公式，就可以對數據進行預測了。我們可以用R語言的predict()函數來計算預測值y0和相應的預測區間，並把實際值和預測值一起可視化化展示。
> par(mfrow=c(1,1)) #設置畫面布局
# 預測計算
> dfp<-predict(lm1,interval="prediction")
# 打印預測時
> head(dfp,10)
fit lwr upr
2014-03-21 3160.526 3046.425 3274.626
2014-03-24 3193.253 3078.868 3307.637
2014-03-25 3240.389 3126.171 3354.607
2014-03-26 3228.565 3114.420 3342.710
2014-03-27 3222.528 3108.342 3336.713
2014-03-28 3262.399 3148.132 3376.666
2014-03-31 3291.996 3177.648 3406.344
2014-04-01 3305.870 3191.447 3420.294
2014-04-02 3275.370 3161.018 3389.723
2014-04-03 3297.358 3182.960 3411.755
# 合並數據
> mdf<-merge(df$y,dfp)
# 畫圖
> draw(mdf)
圖例說明
- y, 實際價格，紅色線
- fit, 預測價格，綠色線
- lwr，預測最低價，藍色線
- upr，預測最高價，紫色線
從圖中看出，實際價格y和預測價格fit，在大多數的時候都是很貼近的。我們的一個模型就訓練好了！
3. 模型優化
上文中，我們已經很順利的發現了一個非常不錯的模型。如果要進行模型優化，可以用R語言中update()函數進行模型的調整。我們首先檢查一下每個自變量x1,x2,x3,x4和因變量y之間的關系。
pairs(as.data.frame(df))
從圖中，我們可以發現x2與Y的關系，可能是最偏離線性的。那么，我們嘗試對多元線性回歸模型進行調整，從原模型中去掉x2變量。
# 模型調整
> lm2<-update(lm1, .~. -x2)
> summary(lm2)
Call:
lm(formula = y ~ x1 + x3 + x4, data = df)
Residuals:
Min 1Q Median 3Q Max
-6.0039 -1.3842 0.0177 1.3513 4.8028
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 462.47104 34.26636 13.50 < 2e-16 ***
x1 1.08728 0.10543 10.31 < 2e-16 ***
x3 -0.40788 0.15394 -2.65 0.00864 **
x4 0.42582 0.01718 24.79 < 2e-16 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.142 on 222 degrees of freedom
Multiple R-squared: 0.9692, Adjusted R-squared: 0.9688
F-statistic: 2330 on 3 and 222 DF, p-value: < 2.2e-16
當把自變量x2去掉后，自變量x3的T檢驗反而變大了，同時Adjusted R-squared變小了，所以我們這次調整是有問題的。
如果通過生產和原材料的內在邏輯分析，焦煤與焦炭屬於上下游關系。焦煤是生產焦炭的一種原材料，焦炭是焦煤與其他煉焦煤經過配煤焦化形成的產品，一般生產 1 噸焦炭需要1.33 噸煉焦煤，其中焦煤至少占 30% 。
我們把焦煤 和 焦炭的關系改變一下，增加x1*x2的關系匹配到模型，看看效果。
# 模型調整
> lm3<-update(lm1, .~. + x1*x2)
> summary(lm3)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.8110 -1.3501 -0.0595 1.2019 5.3884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7160.32231 7814.50048 0.916 0.361
x1 -8.45530 10.47167 -0.807 0.420
x2 -10.58406 12.65579 -0.836 0.404
x3 -0.64344 0.15662 -4.108 5.63e-05 ***
x4 0.48363 0.01967 24.584 < 2e-16 ***
x1:x2 0.01505 0.01693 0.889 0.375
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.029 on 220 degrees of freedom
Multiple R-squared: 0.9726, Adjusted R-squared: 0.972
F-statistic: 1563 on 5 and 220 DF, p-value: < 2.2e-16
從結果中發現，增加了x1*x2列后，原來的x1,x2和Intercept的T檢驗都不顯著。繼續調整模型，從模型中去掉x1,x2兩個自變量。
# 模型調整
> lm4<-update(lm3, .~. -x1-x2)
> summary(lm4)
Call:
lm(formula = y ~ x3 + x4 + x1:x2, data = df)
Residuals:
Min 1Q Median 3Q Max
-4.9027 -1.2516 -0.0167 1.2748 5.8683
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.950e+02 1.609e+01 43.183 < 2e-16 ***
x3 -6.284e-01 1.530e-01 -4.108 5.61e-05 ***
x4 4.959e-01 1.785e-02 27.783 < 2e-16 ***
x1:x2 1.133e-03 9.524e-05 11.897 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.035 on 222 degrees of freedom
Multiple R-squared: 0.9722, Adjusted R-squared: 0.9718
F-statistic: 2588 on 3 and 222 DF, p-value: < 2.2e-16
從調整后的結果來看，效果還不錯。不過，也並沒有比最初的模型有所提高。
對於模型調整的過程，如果我們手動調整測試時，一般都會基於業務知識來操作。如果是按照數據指標來計算，我們可以用R語言中提供的逐步回歸的優化方法，通過AIC指標來判斷是否需要參數優化。
#對lm1模型做逐步回歸
> step(lm1)
Start: AIC=324.51
y ~ x1 + x2 + x3 + x4
Df Sum of Sq RSS AIC
908.8 324.51
- x3 1 77.03 985.9 340.90
- x2 1 109.37 1018.2 348.19
- x1 1 249.90 1158.8 377.41
- x4 1 2490.56 3399.4 620.65
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Coefficients:
(Intercept) x1 x2 x3 x4
212.8780 0.8542 0.6672 -0.6674 0.4821
通過計算AIC指標，lm1的模型AIC最小時為324.51，每次去掉一個自變量都會讓AIC的值變大，所以我們還是不調整比較好。
對剛才的lm3模型做逐步回歸的模型調整。
#對lm3模型做逐步回歸
> step(lm3)
Start: AIC=325.7 #當前AIC
y ~ x1 + x2 + x3 + x4 + x1:x2
Df Sum of Sq RSS AIC
- x1:x2 1 3.25 908.8 324.51
905.6 325.70
- x3 1 69.47 975.1 340.41
- x4 1 2487.86 3393.5 622.25
Step: AIC=324.51 #去掉x1*x2項的AIC
y ~ x1 + x2 + x3 + x4
Df Sum of Sq RSS AIC
908.8 324.51
- x3 1 77.03 985.9 340.90
- x2 1 109.37 1018.2 348.19
- x1 1 249.90 1158.8 377.41
- x4 1 2490.56 3399.4 620.65
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Coefficients:
(Intercept) x1 x2 x3 x4
212.8780 0.8542 0.6672 -0.6674 0.4821
通過AIC的判斷，去掉X1*X2項后AIC最小，最后的檢驗結果告訴我們，還是原初的模型是最好的。
4. 案例：黑色系期貨日K線數據驗證
最后，我們用上面5個期貨合約的日K線數據測試一下，找到多元回歸關系。
> lm9<-lm(y~x1+x2+x3+x4,data=df) # 日K線數據
> summary(lm9)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)
Residuals:
Min 1Q Median 3Q Max
-173.338 -37.470 3.465 32.158 178.982
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 386.33482 31.07729 12.431 < 2e-16 ***
x1 0.75871 0.07554 10.045 < 2e-16 ***
x2 -0.62907 0.14715 -4.275 2.24e-05 ***
x3 1.16070 0.05224 22.219 < 2e-16 ***
x4 0.46461 0.02168 21.427 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 57.78 on 565 degrees of freedom
Multiple R-squared: 0.9844, Adjusted R-squared: 0.9843
F-statistic: 8906 on 4 and 565 DF, p-value: < 2.2e-16
數據集的基本統計信息。
> summary(df)
Index x1 x2
Min. :2014-03-21 00:00:00 Min. : 606.5 Min. :494.0
1st Qu.:2014-10-21 06:00:00 1st Qu.: 803.5 1st Qu.:613.1
Median :2015-05-20 12:00:00 Median : 939.0 Median :705.8
Mean :2015-05-21 08:02:31 Mean : 936.1 Mean :695.3
3rd Qu.:2015-12-16 18:00:00 3rd Qu.:1075.0 3rd Qu.:773.0
Max. :2016-07-25 00:00:00 Max. :1280.0 Max. :898.0
x3 x4 y
Min. :284.0 Min. :1691 Min. :1626
1st Qu.:374.1 1st Qu.:2084 1st Qu.:2012
Median :434.0 Median :2503 Median :2378
Mean :476.5 Mean :2545 Mean :2395
3rd Qu.:545.8 3rd Qu.:2916 3rd Qu.:2592
Max. :825.0 Max. :3480 Max. :3414
對於日K線數據，黑色系的5個品種，同樣具有非常強的相關關系，那么我們就可以把這個結論應用到實際的交易中了。