R語言解讀多元線性回歸模型

在許多生活和工作的實(shí)際問題中，影響因變量的因素可能不止一個(gè)，比如對(duì)于知識(shí)水平越高的人，收入水平也越高，這樣的一個(gè)結(jié)論。這其中可能包括了因?yàn)楦玫募彝l件，所以有了更好的教育;因?yàn)樵谝痪€城市發(fā)展，所以有了更好的工作機(jī)會(huì);所處的行業(yè)趕上了大的經(jīng)濟(jì)上行周期等。要想解讀這些規(guī)律，是復(fù)雜的、多維度的，多元回歸分析方法更適合解讀生活的規(guī)律。

由于本文為非統(tǒng)計(jì)的專業(yè)文章，所以當(dāng)出現(xiàn)與教課書不符的描述，請(qǐng)以教課書為準(zhǔn)。本文力求用簡(jiǎn)化的語言，來介紹多元線性回歸的知識(shí)，同時(shí)配合R語言的實(shí)現(xiàn)。

1. 多元線性回歸介紹

對(duì)比一元線性回歸，多元線性回歸是用來確定2個(gè)或2個(gè)以上變量間關(guān)系的統(tǒng)計(jì)分析方法。多元線性回歸的基本的分析方法與一元線性回歸方法是類似的，我們首先需要對(duì)選取多元數(shù)據(jù)集并定義數(shù)學(xué)模型，然后進(jìn)行參數(shù)估計(jì)，對(duì)估計(jì)出來的參數(shù)進(jìn)行顯著性檢驗(yàn)，殘差分析，異常點(diǎn)檢測(cè)，最后確定回歸方程進(jìn)行模型預(yù)測(cè)。

由于多元回歸方程有多個(gè)自變量，區(qū)別于一元回歸方程，有一項(xiàng)很重要的操作就是自變量的優(yōu)化，挑選出相關(guān)性最顯著的自變量，同時(shí)去除不顯著的自變量。在R語言中，有很方便地用于優(yōu)化函數(shù)，可以很好的幫助我們來改進(jìn)回歸模型。

下面就開始多元線性回歸的建模過程。

2. 多元線性回歸建模

做過商品期貨研究的人，都知道黑色系品種是具有產(chǎn)業(yè)鏈上下游的關(guān)系。鐵礦石是煉鋼的原材料，焦煤和焦炭是煉鋼的能源資源，熱卷即熱軋卷板是以板坯為原料經(jīng)加熱后制成的鋼板，螺紋鋼是表面帶肋的鋼筋。

由于有產(chǎn)業(yè)鏈的關(guān)系，假設(shè)我們想要預(yù)測(cè)螺紋鋼的價(jià)格，那么影響螺紋鋼價(jià)格的因素可以會(huì)涉及到原材料，能源資源和同類材料等。比如，鐵礦石價(jià)格如果上漲，螺紋鋼就應(yīng)該要漲價(jià)了。

2.1 數(shù)據(jù)集和數(shù)學(xué)模型

先從數(shù)據(jù)開始介紹，這次的數(shù)據(jù)集，我選擇的期貨黑色系的品種的商品期貨，包括了大連期貨交易所的 焦煤(JM)，焦炭(J)，鐵礦石(I)，上海期貨交易所的 螺紋鋼(RU) 和 熱卷(HC)。

數(shù)據(jù)集為2016年3月15日，當(dāng)日白天開盤的交易數(shù)據(jù)，為黑色系的5個(gè)期貨合約的分鐘線的價(jià)格數(shù)據(jù)。

數(shù)據(jù)集包括有6列:索引, 為時(shí)間

x1, 為焦炭(j1605)合約的1分鐘線的報(bào)價(jià)數(shù)據(jù)

x2, 為焦煤(jm1605)合約的1分鐘線的報(bào)價(jià)數(shù)據(jù)

x3, 為鐵礦石(i1605)合約的1分鐘線的報(bào)價(jià)數(shù)

x4, 為熱卷(hc1605)合約的1分鐘線的報(bào)價(jià)數(shù)據(jù)

y, 為螺紋鋼(rb1605)合約的1分鐘線的報(bào)價(jià)數(shù)據(jù)

假設(shè)螺紋鋼的價(jià)格與其他4個(gè)商品的價(jià)格有線性關(guān)系，那么我們建立以螺紋鋼為因變量，以焦煤、焦炭、鐵礦石和熱卷的為自變量的多元線性回歸模型。用公式表示為：

y，為因變量，螺紋鋼

x1，為自變量，焦煤

x2，為自變量，焦炭x3，為自變量，鐵礦石

x4，為自變量，熱卷

a，為截距

b,c,d,e，為自變量系數(shù)

ε, 為殘差，是其他一切不確定因素影響的總和，其值不可觀測(cè)。假定ε服從正態(tài)分布N(0,σ^2)。

通過對(duì)多元線性回歸模型的數(shù)學(xué)定義，接下來讓我們利用數(shù)據(jù)集做多元回歸模型的參數(shù)估計(jì)。

2.2. 回歸參數(shù)估計(jì)

上面公式中，回歸參數(shù) a, b, c, d,e都是我們不知道的，參數(shù)估計(jì)就是通過數(shù)據(jù)來估計(jì)出這些參數(shù)，從而確定自變量和因變量之前的關(guān)系。我們的目標(biāo)是要計(jì)算出一條直線，使直線上每個(gè)點(diǎn)的Y值和實(shí)際數(shù)據(jù)的Y值之差的平方和最小，即(Y1實(shí)際-Y1預(yù)測(cè))^2+(Y2實(shí)際-Y2預(yù)測(cè))^2+ …… +(Yn實(shí)際-Yn預(yù)測(cè))^2 的值最小。參數(shù)估計(jì)時(shí)，我們只考慮Y隨X自變量的線性變化的部分，而殘差ε是不可觀測(cè)的，參數(shù)估計(jì)法并不需要考慮殘差。

類似于一元線性回歸，我們用R語言來實(shí)現(xiàn)對(duì)數(shù)據(jù)的回歸模型的參數(shù)估計(jì)，用lm()函數(shù)來實(shí)現(xiàn)多元線性回歸的建模過程。

這樣我們就得到了y和x關(guān)系的方程。

2.3. 回歸方程的顯著性檢驗(yàn)

參考一元線性回歸的顯著性檢驗(yàn)，多元線性回歸的顯著性檢驗(yàn)，同樣是需要經(jīng)過 T檢驗(yàn)，F(xiàn)檢驗(yàn)，和R^2(R平方)相關(guān)系統(tǒng)檢驗(yàn)。在R語言中這三種檢驗(yàn)的方法都已被實(shí)現(xiàn)，我們只需要把結(jié)果解讀，我們可以summary()函數(shù)來提取模型的計(jì)算結(jié)果。

T檢驗(yàn)：所自變量都是非常顯著***

F檢驗(yàn)：同樣是非常顯著，p-value < 2.2e-16

調(diào)整后的R^2：相關(guān)性非常強(qiáng)為0.972

最后，我們通過的回歸參數(shù)的檢驗(yàn)與回歸方程的檢驗(yàn)，得到最后多元線性回歸方程為：

2.4 殘差分析和異常點(diǎn)檢測(cè)

在得到的回歸模型進(jìn)行顯著性檢驗(yàn)后，還要在做殘差分析(預(yù)測(cè)值和實(shí)際值之間的差)，檢驗(yàn)?zāi)Ｐ偷恼_性，殘差必須服從正態(tài)分布N(0,σ^2)。直接用plot()函數(shù)生成4種用于模型診斷的圖形，進(jìn)行直觀地分析。

殘差和擬合值(左上)，殘差和擬合值之間數(shù)據(jù)點(diǎn)均勻分布在y=0兩側(cè)，呈現(xiàn)出隨機(jī)的分布，紅色線呈現(xiàn)出一條平穩(wěn)的曲線并沒有明顯的形狀特征。

殘差QQ圖(右上)，數(shù)據(jù)點(diǎn)按對(duì)角直線排列，趨于一條直線，并被對(duì)角直接穿過，直觀上符合正態(tài)分布。

標(biāo)準(zhǔn)化殘差平方根和擬合值(左下)，數(shù)據(jù)點(diǎn)均勻分布在y=0兩側(cè)，呈現(xiàn)出隨機(jī)的分布，紅色線呈現(xiàn)出一條平穩(wěn)的曲線并沒有明顯的形狀特征。

標(biāo)準(zhǔn)化殘差和杠桿值(右下)，沒有出現(xiàn)紅色的等高線，則說明數(shù)據(jù)中沒有特別影響回歸結(jié)果的異常點(diǎn)。

結(jié)論，沒有明顯的異常點(diǎn)，殘差符合假設(shè)條件。

2.5. 模型預(yù)測(cè)

我們得到了多元線性回歸方程的公式，就可以對(duì)數(shù)據(jù)進(jìn)行預(yù)測(cè)了。我們可以用R語言的predict()函數(shù)來計(jì)算預(yù)測(cè)值y0和相應(yīng)的預(yù)測(cè)區(qū)間，并把實(shí)際值和預(yù)測(cè)值一起可視化化展示。

圖例說明：

y, 實(shí)際價(jià)格，紅色線

fit, 預(yù)測(cè)價(jià)格，綠色線

lwr，預(yù)測(cè)最低價(jià)，藍(lán)色線

upr，預(yù)測(cè)最高價(jià)，紫色線

從圖中看出，實(shí)際價(jià)格y和預(yù)測(cè)價(jià)格fit，在大多數(shù)的時(shí)候都是很貼近的。我們的一個(gè)模型就訓(xùn)練好了!

3. 模型優(yōu)化

上文中，我們已經(jīng)很順利的發(fā)現(xiàn)了一個(gè)非常不錯(cuò)的模型。如果要進(jìn)行模型優(yōu)化，可以用R語言中update()函數(shù)進(jìn)行模型的調(diào)整。我們首先檢查一下每個(gè)自變量x1,x2,x3,x4和因變量y之間的關(guān)系。

從圖中，我們可以發(fā)現(xiàn)x2與Y的關(guān)系，可能是最偏離線性的。那么，我們嘗試對(duì)多元線性回歸模型進(jìn)行調(diào)整，從原模型中去掉x2變量。

當(dāng)把自變量x2去掉后，自變量x3的T檢驗(yàn)反而變大了，同時(shí)Adjusted R-squared變小了，所以我們這次調(diào)整是有問題的。

如果通過生產(chǎn)和原材料的內(nèi)在邏輯分析，焦煤與焦炭屬于上下游關(guān)系。焦煤是生產(chǎn)焦炭的一種原材料，焦炭是焦煤與其他煉焦煤經(jīng)過配煤焦化形成的產(chǎn)品，一般生產(chǎn) 1 噸焦炭需要1.33 噸煉焦煤，其中焦煤至少占 30% 。

我們把焦煤 和 焦炭的關(guān)系改變一下，增加x1*x2的關(guān)系匹配到模型，看看效果。

從結(jié)果中發(fā)現(xiàn)，增加了x1*x2列后，原來的x1,x2和Intercept的T檢驗(yàn)都不顯著。繼續(xù)調(diào)整模型，從模型中去掉x1,x2兩個(gè)自變量。

從調(diào)整后的結(jié)果來看，效果還不錯(cuò)。不過，也并沒有比最初的模型有所提高。

對(duì)于模型調(diào)整的過程，如果我們手動(dòng)調(diào)整測(cè)試時(shí)，一般都會(huì)基于業(yè)務(wù)知識(shí)來操作。如果是按照數(shù)據(jù)指標(biāo)來計(jì)算，我們可以用R語言中提供的逐步回歸的優(yōu)化方法，通過AIC指標(biāo)來判斷是否需要參數(shù)優(yōu)化。

通過計(jì)算AIC指標(biāo)，lm1的模型AIC最小時(shí)為324.51，每次去掉一個(gè)自變量都會(huì)讓AIC的值變大，所以我們還是不調(diào)整比較好。

對(duì)剛才的lm3模型做逐步回歸的模型調(diào)整。

通過AIC的判斷，去掉X1*X2項(xiàng)后AIC最小，最后的檢驗(yàn)結(jié)果告訴我們，還是原初的模型是最好的。

4. 案例：黑色系期貨日K線數(shù)據(jù)驗(yàn)證

最后，我們用上面5個(gè)期貨合約的日K線數(shù)據(jù)測(cè)試一下，找到多元回歸關(guān)系。

數(shù)據(jù)集的基本統(tǒng)計(jì)信息。

對(duì)于日K線數(shù)據(jù)，黑色系的5個(gè)品種，同樣具有非常強(qiáng)的相關(guān)關(guān)系，那么我們就可以把這個(gè)結(jié)論應(yīng)用到實(shí)際的交易中了。

本文通過多元回歸的統(tǒng)計(jì)分析方法，介紹多元回歸在金融市場(chǎng)的基本應(yīng)用。我們通過建立因變量和多個(gè)自變量的模型，從而發(fā)現(xiàn)生活中更復(fù)雜的規(guī)律，并建立有效的驗(yàn)證指標(biāo)。讓我們的技術(shù)優(yōu)勢(shì)，去金融市場(chǎng)搶錢吧。