[[339266]]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

def looper(limit):  
    cols = ['x1', 'x2', 'x3', 'x5', 'x6', 'x7', 'x8', 'x9']  
    for i in range(len(cols)):  
        datadata1 = 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. 剔除無效變量后的回歸模型

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

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

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

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