很多人在做數(shù)據(jù)分析時(shí)會(huì)經(jīng)常用到一元線性回歸，這是描述兩個(gè)變量間統(tǒng)計(jì)關(guān)系的最簡(jiǎn)單的回歸模型。但現(xiàn)實(shí)問(wèn)題中，我們往往會(huì)碰到多個(gè)變量間的線性關(guān)系的問(wèn)題，這時(shí)就要用到多元線性回歸，多元線性回歸是一元回歸的一種推廣，其在實(shí)際應(yīng)用中非常廣泛，本文就用python代碼來(lái)展示一下如何用多元線性回歸來(lái)解決實(shí)際問(wèn)題。

圖1. 多元回歸模型中要用到的公式

如圖1所示，我們假設(shè)隨機(jī)變量y與一般變量x1、x2、...、xp之間線性回歸模型為（1）式，式中y為因變量，x1、x2、...、xp是自變量，β1、β2、...、βp是回歸系數(shù)，β0是回歸常數(shù)。對(duì)于一個(gè)實(shí)際問(wèn)題，如果我們獲得n組觀測(cè)數(shù)據(jù)（xi1，xi2，...，xip；y）（i = 1，2，...，n），則我們可以把這n組觀測(cè)數(shù)據(jù)寫(xiě)成矩陣形式y=Xβ+ε。

在求出了回歸方程之后，我們往往還要對(duì)回歸方程進(jìn)行顯著性檢驗(yàn)。這里的顯著性檢驗(yàn)主要包括三部分。第一個(gè)是F檢驗(yàn)，也就是檢驗(yàn)自變量x1、x2、...、xp從整體上對(duì)y是否有明顯的影響，主要用到（2）、（3）、（4）式，其中（2）和（3）式是一個(gè)式子，不過(guò)是用不同符號(hào)表示；第二個(gè)是t檢驗(yàn)，是對(duì)每個(gè)自變量進(jìn)行顯著性檢驗(yàn)，就是看每個(gè)自變量是否對(duì)y有顯著性影響，這和前面從整體上檢驗(yàn)還是有區(qū)別的；第三個(gè)是擬合優(yōu)度，也就是R2，其取值在0到1之間，越接近1，表明回歸擬合的效果越好，越接近于0，則效果越差，但R只能直觀反映擬合的效果，不能代替F檢驗(yàn)作為嚴(yán)格的顯著性檢驗(yàn)。

上面是多元線性回歸的一個(gè)簡(jiǎn)單介紹，其詳細(xì)原理內(nèi)容較多，有興趣的讀者可以去查閱一下相關(guān)文獻(xiàn)，這里不再贅述，只重點(diǎn)講解如何用python進(jìn)行分析。下面我們還是用代碼來(lái)展示一下多元線性回歸的分析過(guò)程。

這里我們用到的數(shù)據(jù)來(lái)源于2013年《中國(guó)統(tǒng)計(jì)年鑒》，數(shù)據(jù)以居民的消費(fèi)性支出為因變量y，其他9個(gè)變量為自變量，其中x1是居民的食品花費(fèi)，x2是衣著花費(fèi)，x3是居住花費(fèi)，x4是醫(yī)療保健花費(fèi)，x5是文教娛樂(lè)花費(fèi)，x6是職工平均工資，x7是地區(qū)的人均GDP，x8是地區(qū)的消費(fèi)價(jià)格指數(shù)，x9是地區(qū)的失業(yè)率。在這所有變量里面，x1至x7以及y的單位是元，x9是百分?jǐn)?shù)，x8沒(méi)有單位，因?yàn)槠涫窍M(fèi)價(jià)格指數(shù)。數(shù)據(jù)的總體大小為31x10，即31行、10列，大體內(nèi)容如圖2所示。

圖2. 數(shù)據(jù)集部分內(nèi)容

首先還是導(dǎo)入需要的庫(kù)。

import?numpy?as?np
import?pandas?as?pd
import?statsmodels.api?as?sm

接下來(lái)是數(shù)據(jù)預(yù)處理，因?yàn)樵瓟?shù)據(jù)的列標(biāo)太長(zhǎng)，我們要處理一下，去除其中的中文，只留下英文名稱(chēng)。

file?=?r'C:\Users\data.xlsx'
data?=?pd.read_excel(file)
data.columns?=?['y',?'x1',?'x2',?'x3',?'x4',?'x5',?'x6',?'x7',?'x8',?'x9']

然后我們就開(kāi)始生成多元線性模型，代碼如下。

x?=?sm.add_constant(data.iloc[:,1:])?#生成自變量
y?=?data['y']?#生成因變量
model?=?sm.OLS(y,?x)?#生成模型
result?=?model.fit()?#模型擬合
result.summary()?#模型描述

很明顯，這里的自變量是指x1到x9這9個(gè)自變量，代碼data.iloc[:,1:]就是去掉原數(shù)據(jù)中第一列，也就是y那一列的數(shù)據(jù)，result.summary()則是生成一份結(jié)果描述，其內(nèi)容如圖3所示。

圖3. 包含所有自變量的回歸結(jié)果

在這個(gè)結(jié)果中，我們主要看“coef”、“t”和“P>|t|”這三列。coef就是前面說(shuō)過(guò)的回歸系數(shù)，const這個(gè)值就是回歸常數(shù)，所以我們得到的這個(gè)回歸模型就是y = 320.640948 + 1.316588 x1 + 1.649859 x2 + 2.17866 x3 - 0.005609 x4 + 1.684283 x5 + 0.01032 x6 + 0.003655 x7 -19.130576 x8 + 50.515575 x9。而“t”和“P>|t|”這兩列是等價(jià)的，使用時(shí)選擇其中一個(gè)就行，其主要用來(lái)判斷每個(gè)自變量和y的線性顯著關(guān)系，后面我們會(huì)講到。從圖中還可以看出，Prob (F-statistic)為4.21e-20，這個(gè)值就是我們常用的P值，其接近于零，說(shuō)明我們的多元線性方程是顯著的，也就是y與x1、x2、...、x9有著顯著的線性關(guān)系，而R-squared是0.992，也說(shuō)明這個(gè)線性關(guān)系比較顯著。理論上，這個(gè)多元線性方程已經(jīng)求出來(lái)了，而且效果還不錯(cuò)，我們就可以用其進(jìn)行預(yù)測(cè)了，但這里我們還是要進(jìn)行更深一步的探討。前面說(shuō)過(guò)，y與x1、x2、...、x9有著顯著的線性關(guān)系，這里要注意x1到x9這9個(gè)變量被看作是一個(gè)整體，y與這個(gè)整體有顯著的線性關(guān)系，但不代表y與其中的每個(gè)自變量都有顯著的線性關(guān)系，我們?cè)谶@里要找出那些與y的線性關(guān)系不顯著的自變量，然后把它們剔除，只留下關(guān)系顯著的，這就是前面說(shuō)過(guò)的t檢驗(yàn)，t檢驗(yàn)的原理內(nèi)容有些復(fù)雜，有興趣的讀者可以自行查閱資料，這里不再贅述。我們可以通過(guò)圖3中“P>|t|”這一列來(lái)判斷，這一列中我們可以選定一個(gè)閾值，比如統(tǒng)計(jì)學(xué)常用的就是0.05、0.02或0.01，這里我們就用0.05，凡是P>|t|這列中數(shù)值大于0.05的自變量，我們都把它剔除掉，這些就是和y線性關(guān)系不顯著的自變量，所以都舍去，請(qǐng)注意這里指的自變量是x1到x9，不包括圖3中const這個(gè)值。但是這里有一個(gè)原則，就是一次只能剔除一個(gè)，剔除的這個(gè)往往是P值最大的那個(gè)，比如圖3中P值最大的是x4，那么就把它剔除掉，然后再用剩下的x1、x2、x3、x5、x6、x7、x8、x9來(lái)重復(fù)上述建模過(guò)程，再找出P值最大的那個(gè)自變量，把它剔除，如此重復(fù)這個(gè)過(guò)程，直到所有P值都小于等于0.05，剩下的這些自變量就是我們需要的自變量，這些自變量和y的線性關(guān)系都比較顯著，我們要用這些自變量來(lái)進(jìn)行建模。

我們可以將上述過(guò)程寫(xiě)成一個(gè)函數(shù)，命名為looper，代碼如下。

def?looper(limit):
????cols?=?['x1',?'x2',?'x3',?'x5',?'x6',?'x7',?'x8',?'x9']
????for?i?in?range(len(cols)):
????????data1?=?data[cols]
????????x?=?sm.add_constant(data1)?#生成自變量
????????y?=?data['y']?#生成因變量
????????model?=?sm.OLS(y,?x)?#生成模型
????????result?=?model.fit()?#模型擬合
????????pvalues?=?result.pvalues?#得到結(jié)果中所有P值
????????pvalues.drop('const',inplace=True)?#把const取得
????????pmax?=?max(pvalues)?#選出最大的P值
????????if?pmax>limit:
????????????ind?=?pvalues.idxmax()?#找出最大P值的index
????????????cols.remove(ind)?#把這個(gè)index從cols中刪除
????????else:
????????????return?result

result?=?looper(0.05)
result.summary()

其結(jié)果如圖4所示。從結(jié)果中可以看到最后剩下的有效變量為x1、x2、x3和x5，我們得到的多元線性模型為y = -1694.6269 + 1.3642 x1 + 1.7679 x2 + 2.2894 x3 + 1.7424 x5，這個(gè)就是我們最終要用到的有效的多元線性模型。

圖4. 剔除無(wú)效變量后的回歸模型

那么問(wèn)題來(lái)了，前面我們得到的包含所有自變量的多元線性模型和這個(gè)剔除部分變量的模型，我們要選擇哪一個(gè)，畢竟第一個(gè)模型的整體線性效果也挺顯著，依據(jù)筆者的經(jīng)驗(yàn)，這個(gè)還是要看具體的項(xiàng)目要求。因?yàn)槲覀儗?shí)際項(xiàng)目中遇到的問(wèn)題都是現(xiàn)實(shí)生活中真實(shí)存在的例子，不再是單純的數(shù)學(xué)題了，比如本例中的x8消費(fèi)價(jià)格指數(shù)和x9地區(qū)的失業(yè)率，這兩個(gè)肯定對(duì)y是有一定影響的，如果盲目剔除，可能會(huì)對(duì)最終的結(jié)果產(chǎn)生不良影響，所以我們還是要根據(jù)實(shí)際需求來(lái)做決定。

最后還有一個(gè)問(wèn)題要討論一下，就是本例中沒(méi)有對(duì)原始數(shù)據(jù)進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化。那么我們?cè)跀?shù)據(jù)分析中是否要對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化？

這個(gè)也是要視情況而定。像本例中這些數(shù)據(jù)都是帶有具體的量綱和單位，那么就不要對(duì)其進(jìn)行標(biāo)準(zhǔn)化，我們得到的這個(gè)線性回歸模型是在原始變量基礎(chǔ)上進(jìn)行擬合所得的結(jié)果，這個(gè)式子是包含物理單位的，說(shuō)白了它們都是有一定實(shí)際意義的。在這種情況下，我們輸入特定的自變量的值，即可得到相應(yīng)的y值，預(yù)測(cè)效果直截了當(dāng)，這是采取原始數(shù)據(jù)進(jìn)行線性擬合的好處。

如果我們對(duì)原始數(shù)據(jù)采取了標(biāo)準(zhǔn)化處理，情況就不同了。標(biāo)準(zhǔn)化處理后自變量、因變量的物理單位沒(méi)有了，我們拿此時(shí)的模型做預(yù)測(cè)時(shí)就會(huì)十分麻煩，要對(duì)新的自變量取值進(jìn)行標(biāo)準(zhǔn)化，得到的y還是一個(gè)標(biāo)準(zhǔn)化后的數(shù)據(jù)，一眼看不到它的實(shí)際大小和物理意義。當(dāng)然有些純數(shù)學(xué)問(wèn)題，其變量沒(méi)有單位，這時(shí)候可以對(duì)其進(jìn)行標(biāo)準(zhǔn)化，這會(huì)有利于對(duì)問(wèn)題的分析。所以這個(gè)還是要視情況而定。

贊 賞 作 者

▼點(diǎn)擊成為社區(qū)會(huì)員? ?喜歡就點(diǎn)個(gè)在看吧