免責(zé)聲明：本文是通過(guò)網(wǎng)絡(luò)收集并結(jié)合自身學(xué)習(xí)等途徑合法獲取，僅作為學(xué)習(xí)交流使用，其版權(quán)歸出版社或者原創(chuàng)作者所有，并不對(duì)涉及的版權(quán)問(wèn)題負(fù)責(zé)。若原創(chuàng)作者或者出版社認(rèn)為侵權(quán)，請(qǐng)聯(lián)系及時(shí)聯(lián)系，我將立即刪除文章，十分感謝！
注：來(lái)源劉順祥《從零開(kāi)始學(xué)Python數(shù)據(jù)分析與挖掘》，版權(quán)歸原作者所有，僅供學(xué)習(xí)使用，不用于商業(yè)用途，如有侵權(quán)請(qǐng)留言聯(lián)系刪除，感謝合作。
前言線性回歸模型屬于經(jīng)典的統(tǒng)計(jì)學(xué)模型，該模型的應(yīng)用場(chǎng)景是根據(jù)已知的變量（自變量）來(lái)預(yù)測(cè)某個(gè)連續(xù)的數(shù)值變量（因變量）。例如，餐廳根據(jù)每天的營(yíng)業(yè)數(shù)據(jù)（包括菜譜價(jià)格、就餐人數(shù)、預(yù)定人數(shù)、特價(jià)菜折扣等）預(yù)測(cè)就餐規(guī)?；驙I(yíng)業(yè)額；網(wǎng)站根據(jù)訪問(wèn)的歷史數(shù)據(jù)（包括新用戶的注冊(cè)量、老用戶的活躍度、網(wǎng)頁(yè)內(nèi)容的更新頻率等）預(yù)測(cè)用戶的支付轉(zhuǎn)化率；醫(yī)院根據(jù)患者的病歷數(shù)據(jù)（如體檢指標(biāo)、藥物服用情況、平時(shí)的飲食習(xí)慣等）預(yù)測(cè)某種疾病發(fā)生的概率。站在數(shù)據(jù)挖掘的角度看待線性回歸模型，它屬于一種有監(jiān)督的學(xué)習(xí)算法，即在建模過(guò)程中必須同時(shí)具備自變量x和因變量y。本章的重點(diǎn)就是介紹有關(guān)線性回歸模型背后的數(shù)學(xué)原理和應(yīng)用實(shí)戰(zhàn)，通過(guò)本章內(nèi)容的學(xué)習(xí)，讀者將會(huì)掌握如下內(nèi)容：
一元線性回歸模型也被稱為簡(jiǎn)單線性回歸模型，是指模型中只含有一個(gè)自變量和一個(gè)因變量，用來(lái)建模的數(shù)據(jù)集可以表示成{(x1,y1),(x2,y2),…,(xn,yn)}。其中，xi表示自變量x的第i個(gè)值，yi表示因變量y的第i個(gè)值，n表示數(shù)據(jù)集的樣本量。當(dāng)模型構(gòu)建好之后，就可以根據(jù)其他自變量x的值，預(yù)測(cè)因變量y的值，該模型的數(shù)學(xué)公式可以表示成：
y=a+bx+ε 如上公式所示，該模型特別像初中所學(xué)的一次函數(shù)。其中，a為模型的截距項(xiàng)，b為模型的斜率項(xiàng)，ε為模型的誤差項(xiàng)。模型中的a和b統(tǒng)稱為回歸系數(shù)，誤差項(xiàng)ε的存在主要是為了平衡等號(hào)兩邊的值，通常被稱為模型無(wú)法解釋的部分。 為了使讀者理解簡(jiǎn)單線性回歸模型的數(shù)學(xué)公式，這里不妨以收入數(shù)據(jù)集為例，探究工作年限與收入之間的關(guān)系。在第6章的數(shù)據(jù)可視化部分已經(jīng)介紹了有關(guān)散點(diǎn)圖的繪制，下面將繪制工作年限與收入的散點(diǎn)圖，并根據(jù)散點(diǎn)圖添加一條擬合線：
# 工作年限與收入之間的散點(diǎn)圖# 導(dǎo)入第三方模塊import pandas as pdimport matplotlib.pyplot as pltimport seaborn as sns# 導(dǎo)入數(shù)據(jù)集income = pd.read_csv(r'D:\projects\notebook\data\Salary_Data.csv')# 繪制散點(diǎn)圖sns.lmplot(x = 'YearsExperience', y = 'Salary', data = income, ci = None)# 顯示圖形plt.show()結(jié)果:
圖7-1反映的就是自變量YearsExperience與因變量Salary之間的散點(diǎn)圖，從散點(diǎn)圖的趨勢(shì)來(lái)看，工作年限與收入之間存在明顯的正相關(guān)關(guān)系，即工作年限越長(zhǎng)，收入水平越高。圖中的直線就是關(guān)于散點(diǎn)的線性回歸擬合線，從圖中可知，每個(gè)散點(diǎn)基本上都是圍繞在擬合線附近。雖然通過(guò)可視化的方法可以得知散點(diǎn)間的關(guān)系和擬合線，但如何得到這條擬合線的數(shù)學(xué)表達(dá)式呢？
擬合線的求解本節(jié)的內(nèi)容就是關(guān)于簡(jiǎn)單線性回歸模型的求解，即如何根據(jù)自變量x和因變量y，求解回歸系數(shù)a和b。前面已經(jīng)提到，誤差項(xiàng)ε是為了平衡等號(hào)兩邊的值，如果擬合線能夠精確地捕捉到每一個(gè)點(diǎn)（所有的散點(diǎn)全部落在擬合線上），那么對(duì)應(yīng)的誤差項(xiàng)ε應(yīng)該為0。按照這個(gè)思路來(lái)看，要想得到理想的擬合線，就必須使誤差項(xiàng)ε達(dá)到最小。由于誤差項(xiàng)是y與a+bx的差，結(jié)果可能為正值或負(fù)值，因此誤差項(xiàng)ε達(dá)到最小的問(wèn)題需轉(zhuǎn)換為誤差平方和最小的問(wèn)題（最小二乘法的思路）。誤差平方和的公式可以表示為：
由于建模時(shí)的自變量值和因變量值都是已知的，因此求解誤差平方和最小值的問(wèn)題就是求解函數(shù)J(a,b)的最小值，而該函數(shù)的參數(shù)就是回歸系數(shù)a和b。該目標(biāo)函數(shù)其實(shí)就是一個(gè)二元二次函數(shù)，如需使得目標(biāo)函數(shù)J(a,b)達(dá)到最小，可以使用偏導(dǎo)數(shù)的方法求解出參數(shù)a和b，進(jìn)而得到目標(biāo)函數(shù)的最小值。關(guān)于目標(biāo)函數(shù)的求導(dǎo)過(guò)程如下：
第一步：展開(kāi)平方項(xiàng)
第二步：設(shè)偏導(dǎo)數(shù)為0
第三步：和公式轉(zhuǎn)換
第四步：化解
第五步：將參數(shù)a帶入，求解b
如上推導(dǎo)結(jié)果所示，參數(shù)a和b的值都是關(guān)于自變量x和因變量y的公式。接下來(lái)，根據(jù)該公式，利用Pyhton計(jì)算出回歸模型的參數(shù)值a和b。
# 簡(jiǎn)單線性回歸模型的參數(shù)求解# 樣本量n = income.shape[0]# 計(jì)算自變量、因變量、自變量平方、自變量與因變量乘積的和sum_x = income.YearsExperience.sum()sum_y = income.Salary.sum()sum_x2 = income.YearsExperience.pow(2).sum()xy = income.YearsExperience * income.Salarysum_xy = xy.sum()# 根據(jù)公式計(jì)算回歸模型的參數(shù)b = (sum_xy-sum_x*sum_y/n)/(sum_x2-sum_x**2/n)a = income.Salary.mean()-b*income.YearsExperience.mean()# 打印出計(jì)算結(jié)果print('回歸參數(shù)a的值：',a)print('回歸參數(shù)b的值：',b)結(jié)果：
回歸參數(shù)a的值： 25792.200198668666回歸參數(shù)b的值： 9449.962321455081如上所示，利用Python的計(jì)算功能，最終得到模型的回歸參數(shù)值。你可能會(huì)覺(jué)得麻煩，為了計(jì)算回歸模型的參數(shù)還得人工寫(xiě)代碼，是否有現(xiàn)成的第三方模塊可以直接調(diào)用呢？答案是肯定的，這個(gè)模塊就是statsmodels，它是專門(mén)用于統(tǒng)計(jì)建模的第三方模塊，如需實(shí)現(xiàn)線性回歸模型的參數(shù)求解，可以調(diào)用子模塊中的ols函數(shù)。有關(guān)該函數(shù)的語(yǔ)法及參數(shù)含義可見(jiàn)下方：
ols(formula, data, subset=None, drop_cols=None) 這是一個(gè)語(yǔ)法非常簡(jiǎn)單的函數(shù)，而且參數(shù)也通俗易懂，但該函數(shù)的功能卻很強(qiáng)大，不僅可以計(jì)算模型的參數(shù)，還可以對(duì)模型的參數(shù)和模型本身做顯著性檢驗(yàn)、計(jì)算模型的決定系數(shù)等。接下來(lái)，利用該函數(shù)計(jì)算模型的參數(shù)值，進(jìn)而驗(yàn)證手工方式計(jì)算的參數(shù)是否正確：
# 導(dǎo)入第三方模塊import statsmodels.api as sm# 利用收入數(shù)據(jù)集，構(gòu)建回歸模型fit = sm.formula.ols('Salary ~ YearsExperience', data = income).fit()# 返回模型的參數(shù)值fit.params結(jié)果：
Intercept 25792.200199YearsExperience 9449.962321dtype: float64如上結(jié)果所示，Intercept表示截距項(xiàng)對(duì)應(yīng)的參數(shù)值，YearsExperience表示自變量工作年限對(duì)應(yīng)的參數(shù)值。對(duì)比發(fā)現(xiàn)，函數(shù)計(jì)算出來(lái)的參數(shù)值與手工計(jì)算的結(jié)果完全一致，所以，關(guān)于收入的簡(jiǎn)單線性回歸模型可以表示成：
Salary=25792.20+9449.96YearsExperience讀者會(huì)不會(huì)覺(jué)得一元線性回歸模型比較簡(jiǎn)單呢？它反映的是單個(gè)自變量對(duì)因變量的影響，然而實(shí)際情況中，影響因變量的自變量往往不止一個(gè)，從而需要將一元線性回歸模型擴(kuò)展到多元線性回歸模型。如果構(gòu)建多元線性回歸模型的數(shù)據(jù)集包含n個(gè)觀測(cè)、p+1個(gè)變量（其中p個(gè)自變量和1個(gè)因變量），則這些數(shù)據(jù)可以寫(xiě)成下方的矩陣形式：
其中，xij代表第個(gè)i行的第j個(gè)變量值。如果按照一元線性回歸模型的邏輯，那么多元線性回歸模型應(yīng)該就是因變量y與自變量X的線性組合，即可以將多元線性回歸模型表示成：
y=β0+β1x1+β2x2+…+βpxn+ε 1根據(jù)線性代數(shù)的知識(shí)，可以將上式表示成y=Xβ+ε。其中，β為p×1的一維向量，代表了多元線性回歸模型的偏回歸系數(shù)；ε為n×1的一維向量，代表了模型擬合后每一個(gè)樣本的誤差項(xiàng)。
7.2.1 回歸模型的參數(shù)求解在多元線性回歸模型所涉及的數(shù)據(jù)中，因變量y是一維向量，而自變量X為二維矩陣，所以對(duì)于參數(shù)的求解不像一元線性回歸模型那樣簡(jiǎn)單，但求解的思路是完全一致的。為了使讀者掌握多元線性回歸模型參數(shù)的求解過(guò)程，這里把詳細(xì)的推導(dǎo)步驟羅列到下方：
第一步：構(gòu)建目標(biāo)函數(shù)
根據(jù)線性代數(shù)的知識(shí)，可以將向量的平方和公式轉(zhuǎn)換為向量的內(nèi)積，接下來(lái)需要對(duì)該式進(jìn)行平方項(xiàng)的展現(xiàn)。
第二步：展開(kāi)平方項(xiàng)
J(β)=(y-Xβ)'(y-Xβ) =(y'-β'X')(y-Xβ) =(y'y-y'Xβ-β'X'y+β'X'Xβ)由于上式中的y’Xβ和β’X’y均為常數(shù)，并且常數(shù)的轉(zhuǎn)置就是其本身，因此y’Xβ和β’X’y是相等的。接下來(lái)，對(duì)目標(biāo)函數(shù)求參數(shù)β的偏導(dǎo)數(shù)。
第三步：求偏導(dǎo)
如上式所示，根據(jù)高等數(shù)學(xué)的知識(shí)可知，欲求目標(biāo)函數(shù)的極值，一般都需要對(duì)目標(biāo)函數(shù)求導(dǎo)數(shù)，再令導(dǎo)函數(shù)為0，進(jìn)而根據(jù)等式求得到函數(shù)中的參數(shù)值。
第四步：計(jì)算偏回歸系數(shù)的值
X'Xβ=X'y β=(X'X)-1X'y經(jīng)過(guò)如上四步的推導(dǎo)，最終可以得到偏回歸系數(shù)β與自變量X、因變量y的數(shù)學(xué)關(guān)系。這個(gè)求解過(guò)程也被成為“最小二乘法”?；谝阎钠貧w系數(shù)β就可以構(gòu)造多元線性回歸模型。前文也提到，構(gòu)建模型的最終目的是為了預(yù)測(cè)，即根據(jù)其他已知的自變量X的值預(yù)測(cè)未知的因變量y的值。
7.2.2 回歸模型的預(yù)測(cè)如果已經(jīng)得知某個(gè)多元線性回歸模型y=β0+β1x1+β2x2+…+βpxn，當(dāng)有其他新的自變量值時(shí)，就可以將這些值帶入如上的公式中，最終得到未知的y值。在Python中，實(shí)現(xiàn)線性回歸模型的預(yù)測(cè)可以使用predict“方法”，關(guān)于該“方法”的參數(shù)含義如下：
predict(exog=None, transform=True)接下來(lái)將基于statsmodels模塊對(duì)多元線性回歸模型的參數(shù)進(jìn)行求解，進(jìn)而依據(jù)其他新的自變量值實(shí)現(xiàn)模型的預(yù)測(cè)功能。這里不妨以某產(chǎn)品的利潤(rùn)數(shù)據(jù)集為例，該數(shù)據(jù)集包含5個(gè)變量，分別是產(chǎn)品的研發(fā)成本、管理成本、市場(chǎng)營(yíng)銷成本、銷售市場(chǎng)和銷售利潤(rùn)，數(shù)據(jù)集的部分截圖如表7-1所示。
表7-1中數(shù)據(jù)集中的Profit變量為因變量，其他變量將作為模型的自變量。需要注意的是，數(shù)據(jù)集中的State變量為字符型的離散變量，是無(wú)法直接帶入模型進(jìn)行計(jì)算的，所以建模時(shí)需要對(duì)該變量進(jìn)行特殊處理。有關(guān)產(chǎn)品利潤(rùn)的建模和預(yù)測(cè)過(guò)程如下代碼所示：
# 多元線性回歸模型的構(gòu)建和預(yù)測(cè)# 導(dǎo)入模塊from sklearn import model_selection# 導(dǎo)入數(shù)據(jù)Profit = pd.read_excel(r'D:\projects\notebook\data\Predict to Profit.xlsx')# 將數(shù)據(jù)集拆分為訓(xùn)練集和測(cè)試集train, test = model_selection.train_test_split(Profit, test_size = 0.2, random_state=1234)# 根據(jù)train數(shù)據(jù)集建模model = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + C(State)', data = train).fit()print('模型的偏回歸系數(shù)分別為：\n', model.params)# 刪除test數(shù)據(jù)集中的Profit變量，用剩下的自變量進(jìn)行預(yù)測(cè)test_X = test.drop(labels = 'Profit', axis = 1)pred = model.predict(exog = test_X)print('對(duì)比預(yù)測(cè)值和實(shí)際值的差異：\n',pd.DataFrame({'Prediction':pred,'Real':test.Profit}))結(jié)果：
模型的偏回歸系數(shù)分別為： Intercept 58581.516503C(State)[T.Florida] 927.394424C(State)[T.New York] -513.468310RD_Spend 0.803487Administration -0.057792Marketing_Spend 0.013779dtype: float64對(duì)比預(yù)測(cè)值和實(shí)際值的差異： Prediction Real8 150621.345801 152211.7748 55513.218079 35673.4114 150369.022458 132602.6542 74057.015562 71498.4929 103413.378282 101004.6444 67844.850378 65200.334 173454.059691 166187.9431 99580.888894 97483.5613 128147.138396 134307.3518 130693.433835 124266.90如上結(jié)果所示，得到多元線性回歸模型的回歸系數(shù)及測(cè)試集上的預(yù)測(cè)值，為了比較，將預(yù)測(cè)值和測(cè)試集中的真實(shí)Profit值羅列在一起。針對(duì)如上代碼需要說(shuō)明三點(diǎn)：
對(duì)于輸出的回歸系數(shù)結(jié)果，讀者可能會(huì)感到疑惑，為什么字符型變量State對(duì)應(yīng)兩個(gè)回歸系數(shù)，而且標(biāo)注了Florida和New York。那是因?yàn)樽址妥兞縎tate含有三種不同的值，分別是California、Florida和New York，在建模時(shí)將該變量當(dāng)作啞變量處理，所以三種不同的值就會(huì)衍生出兩個(gè)變量，分別是State[Florida]和State[New York]，而另一個(gè)變量State[California]就成了對(duì)照組。正如建模中的代碼所示，將State變量套在C()中，就表示State變量需要進(jìn)行啞變量處理。但是這樣做會(huì)存在一個(gè)缺陷，那就是無(wú)法指定變量中的某個(gè)值作為對(duì)照組，正如模型結(jié)果中默認(rèn)將State變量的California值作為對(duì)照組（因?yàn)樵撝翟谌齻€(gè)值中的字母順序是第一個(gè)）。如需解決這個(gè)缺陷，就要通過(guò)pandas模塊中的get_dummies函數(shù)生成啞變量，然后將所需的對(duì)照組對(duì)應(yīng)的啞變量刪除即可。為了使讀者明白該解決方案，這里不妨重新建模，并以State變量中的New York值作為對(duì)照組，代碼如下：
# 生成由State變量衍生的啞變量dummies = pd.get_dummies(Profit.State)# 將啞變量與原始數(shù)據(jù)集水平合并Profit_New = pd.concat([Profit,dummies], axis = 1)# 刪除State變量和California變量（因?yàn)镾tate變量已被分解為啞變量，New York變量需要作為參照組）Profit_New.drop(labels = ['State','New York'], axis = 1, inplace = True)# 拆分?jǐn)?shù)據(jù)集Profit_Newtrain, test = model_selection.train_test_split(Profit_New, test_size = 0.2, random_state=1234)# 建模model2 = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + Florida + California', data = train).fit()print('模型的偏回歸系數(shù)分別為：\n', model2.params)結(jié)果：
模型的偏回歸系數(shù)分別為： Intercept 58068.048193RD_Spend 0.803487Administration -0.057792Marketing_Spend 0.013779Florida 1440.862734California 513.468310dtype: float64如上結(jié)果所示，從離散變量State中衍生出來(lái)的啞變量在回歸系數(shù)的結(jié)果里只保留了Florida和California，而New York變量則作為了參照組。以該模型結(jié)果為例，得到的模型公式可以表達(dá)為：
Profit=58068.05+0.80RD_Spend-0.06Administation+0.01Marketing_Spend+1440.86Florida+513.47California 雖然模型的回歸系數(shù)求解出來(lái)了，但從統(tǒng)計(jì)學(xué)的角度該如何解釋模型中的每個(gè)回歸系數(shù)呢？下面分別以研發(fā)成本RD_Spend變量和啞變量Florida為例，解釋這兩個(gè)變量對(duì)模型的作用：在其他變量不變的情況下，研發(fā)成本每增加1美元，利潤(rùn)會(huì)增加0.80美元；在其他變量不變的情況下，以New York為基準(zhǔn)線，如果在Florida銷售產(chǎn)品，利潤(rùn)會(huì)增加1440.86美元。關(guān)于產(chǎn)品利潤(rùn)的多元線性回歸模型已經(jīng)構(gòu)建完成，但是該模型的好與壞并沒(méi)有相應(yīng)的結(jié)論，還需要進(jìn)行模型的顯著性檢驗(yàn)和回歸系數(shù)的顯著性檢驗(yàn)。在下一節(jié)，將重點(diǎn)介紹有關(guān)線性回歸模型中的幾點(diǎn)重要假設(shè)檢驗(yàn)。
7.3 回歸模型的假設(shè)檢驗(yàn)模型的顯著性檢驗(yàn)是指構(gòu)成因變量的線性組合是否有效，即整個(gè)模型中是否至少存在一個(gè)自變量能夠真正影響到因變量的波動(dòng)。該檢驗(yàn)是用來(lái)衡量模型的整體效應(yīng)?；貧w系數(shù)的顯著性檢驗(yàn)是為了說(shuō)明單個(gè)自變量在模型中是否有效，即自變量對(duì)因變量是否具有重要意義。這種檢驗(yàn)則是出于對(duì)單個(gè)變量的肯定與否。模型的顯著性檢驗(yàn)和回歸系數(shù)的顯著性檢驗(yàn)分別使用統(tǒng)計(jì)學(xué)中的F檢驗(yàn)法和t檢驗(yàn)法，接下來(lái)將介紹有關(guān)F檢驗(yàn)和t檢驗(yàn)的理論知識(shí)和實(shí)踐操作。
7.3.1 模型的顯著性檢驗(yàn)——F檢驗(yàn)在統(tǒng)計(jì)學(xué)中，有關(guān)假設(shè)檢驗(yàn)的問(wèn)題，都有一套成熟的步驟。首先來(lái)看一下如何應(yīng)用F檢驗(yàn)法完成模型的顯著性檢驗(yàn)，具體的檢驗(yàn)步驟如下：（1）提出問(wèn)題的原假設(shè)和備擇假設(shè)。（2）在原假設(shè)的條件下，構(gòu)造統(tǒng)計(jì)量F。（3）根據(jù)樣本信息，計(jì)算統(tǒng)計(jì)量的值。（4）對(duì)比統(tǒng)計(jì)量的值和理論F分布的值，如果計(jì)算的統(tǒng)計(jì)量值超過(guò)理論的值，則拒絕原假設(shè)，否則需接受原假設(shè)。
下面將按照上述四個(gè)步驟對(duì)構(gòu)造的多元線性回歸模型進(jìn)行F檢驗(yàn)，進(jìn)一步確定該模型是否可用，詳細(xì)操作步驟如下：步驟一：提出假設(shè)
H0為原假設(shè)，該假設(shè)認(rèn)為模型的所有偏回歸系數(shù)全為0，即認(rèn)為沒(méi)有一個(gè)自變量可以構(gòu)成因變量的線性組合； H1為備擇假設(shè)，正好是原假設(shè)的對(duì)立面，即p個(gè)自變量中，至少有一個(gè)變量可以構(gòu)成因變量的線性組合。就F檢驗(yàn)而言，研究者往往是更加希望通過(guò)數(shù)據(jù)來(lái)推翻原假設(shè)H0，而接受備擇假設(shè)H1的結(jié)論。
步驟二：構(gòu)造統(tǒng)計(jì)量為了使讀者理解F統(tǒng)計(jì)量的構(gòu)造過(guò)程，可以先觀看圖7-2，然后掌握總的離差平方和、回歸離差平方和與誤差平方和的概念與差異。
假設(shè)圖中的斜線代表某條線性擬合線，點(diǎn)p(x,y)代表數(shù)據(jù)集中的某個(gè)點(diǎn)，則為點(diǎn)x處的預(yù)測(cè)值；為真實(shí)值與預(yù)測(cè)值之間的差異；為預(yù)測(cè)值與總體平均值之間的差異；為真實(shí)值與總體平均值之間的差異。如果將這些差異向量做平方和的計(jì)算，則會(huì)得到：
如上公式所示，公式中的ESS稱為誤差平方和，衡量的是因變量的實(shí)際值與預(yù)測(cè)值之間的離差平方和，會(huì)隨著模型的變化而變動(dòng)（因?yàn)槟Ｐ偷淖兓瘯?huì)導(dǎo)致預(yù)測(cè)值的變動(dòng)）；RSS為回歸離差平方和，衡量的是因變量的預(yù)測(cè)值與實(shí)際均值之間的離差平方和，同樣會(huì)隨著模型的變化而變動(dòng)；TSS為總的離差平方和，衡量的是因變量的值與其均值之間的離差平方和，而其值并不會(huì)隨模型的變化而變動(dòng)，即它是一個(gè)固定值。根據(jù)統(tǒng)計(jì)計(jì)算，這三個(gè)離差平方和之間存在這樣的等式關(guān)系：TSS=ESS+RSS。由于TSS的值不會(huì)隨模型的變化而變動(dòng)，因此ESS與RSS之間存在嚴(yán)格的負(fù)向關(guān)系，即ESS的降低會(huì)導(dǎo)致RSS的增加。正如7.1.1節(jié)所介紹的內(nèi)容，線性回歸模型的參數(shù)求解是依據(jù)誤差平方和最小的理論，如果根據(jù)線性回歸模型得到的ESS值達(dá)到最小，那么對(duì)應(yīng)的RSS值就會(huì)達(dá)到最大，進(jìn)而RSS與ESS的商也會(huì)達(dá)到最大。按照這個(gè)邏輯，便可以構(gòu)造F統(tǒng)計(jì)量，該統(tǒng)計(jì)量可以表示成回歸離差平方和RSS與誤差平方和ESS的公式：
其中，p和n-p-1分別為RSS和ESS的自由度。模型擬合得越好， ESS就會(huì)越小， RSS則會(huì)越大，得到的F統(tǒng)計(jì)量也就越大。步驟三：計(jì)算統(tǒng)計(jì)量下面按照F統(tǒng)計(jì)量的公式，運(yùn)用Python計(jì)算該統(tǒng)計(jì)量的值，詳細(xì)的計(jì)算過(guò)程可見(jiàn)下方代碼：
# 導(dǎo)入第三方模塊import numpy as np# 計(jì)算建模數(shù)據(jù)中，因變量的均值ybar = train.Profit.mean()# 統(tǒng)計(jì)變量個(gè)數(shù)和觀測(cè)個(gè)數(shù)p = model2.df_modeln = train.shape[0]# 計(jì)算回歸離差平方和RSS = np.sum((model2.fittedvalues-ybar) ** 2)# 計(jì)算誤差平方和ESS = np.sum(model2.resid ** 2)# 計(jì)算F統(tǒng)計(jì)量的值F = (RSS/p)/(ESS/(n-p-1))print('F統(tǒng)計(jì)量的值：',F)結(jié)果:
F統(tǒng)計(jì)量的值： 174.63721716844674為了驗(yàn)證手工計(jì)算的結(jié)果是否正確，可以通過(guò)fvalue“方法”直接獲得模型的F統(tǒng)計(jì)量值，如下結(jié)果所示，經(jīng)過(guò)對(duì)比發(fā)現(xiàn)，手工計(jì)算的結(jié)果與模型自帶的F統(tǒng)計(jì)量值完全一致：
model2.fvalue結(jié)果：
174.63721715703542步驟四：對(duì)比結(jié)果下結(jié)論最后一步所要做的是對(duì)比F統(tǒng)計(jì)量的值與理論F分布的值，如果讀者手中有F分布表，可以根據(jù)置信水平（0.05）和自由度（5,34）查看對(duì)應(yīng)的分布值。為了簡(jiǎn)單起見(jiàn)，這里直接調(diào)用Python函數(shù)計(jì)算理論分布值：
# 導(dǎo)入模塊from scipy.stats import f# 計(jì)算F分布的理論值F_Theroy = f.ppf(q=0.95, dfn = p,dfd = n-p-1)結(jié)果：
F分布的理論值為： 2.502635007415366如上結(jié)果所示，在原假設(shè)的前提下，計(jì)算出來(lái)的F統(tǒng)計(jì)量值174.64遠(yuǎn)遠(yuǎn)大于F分布的理論值2.50，所以應(yīng)當(dāng)拒絕原假設(shè)，即認(rèn)為多元線性回歸模型是顯著的，也就是說(shuō)回歸模型的偏回歸系數(shù)都不全為0。
7.3.2 回歸系數(shù)的顯著性檢驗(yàn)——t檢驗(yàn)?zāi)Ｐ屯ㄟ^(guò)了顯著性檢驗(yàn)，只能說(shuō)明關(guān)于因變量的線性組合是合理的，但并不能說(shuō)明每個(gè)自變量對(duì)因變量都具有顯著意義，所以還需要對(duì)模型的回歸系數(shù)做顯著性檢驗(yàn)。關(guān)于系數(shù)的顯著性檢驗(yàn)，需要使用t檢驗(yàn)法，構(gòu)造t統(tǒng)計(jì)量。接下來(lái)按照模型顯著性檢驗(yàn)的四個(gè)步驟，對(duì)偏回歸系數(shù)進(jìn)行顯著性檢驗(yàn)。步驟一：提出假設(shè)
如前文所提，t檢驗(yàn)的出發(fā)點(diǎn)就是驗(yàn)證每一個(gè)自變量是否能夠成為影響因變量的重要因素。t檢驗(yàn)的原假設(shè)是假定第j變量的偏回歸系數(shù)為0，即認(rèn)為該變量不是因變量的影響因素；而備擇假設(shè)則是相反的假定，認(rèn)為第j變量是影響因變量的重要因素。
步驟二：構(gòu)造統(tǒng)計(jì)量
其中，為線性回歸模型的第j個(gè)系數(shù)估計(jì)值；βj為原假設(shè)中的假定值，即0；為回歸系數(shù)的標(biāo)準(zhǔn)誤，對(duì)應(yīng)的計(jì)算公式如下：
其中，為誤差平方和，cjj為矩陣(X’X)-1主對(duì)角線上第j個(gè)元素。
步驟三：計(jì)算統(tǒng)計(jì)量如果讀者對(duì)t統(tǒng)計(jì)量值的計(jì)算比較感興趣，可以使用如上公式完成統(tǒng)計(jì)量的計(jì)算，這里就不手工計(jì)算了。為了方便起見(jiàn)，可以直接調(diào)用summary“方法”，輸出線性回歸模型的各項(xiàng)指標(biāo)值：
# 模型的概覽信息model2.summary()結(jié)果：
如上結(jié)果所示，模型的概覽信息包含三個(gè)部分，第一部分主要是有關(guān)模型的信息，例如模型的判決系數(shù)R2，用來(lái)衡量自變量對(duì)因變量的解釋程度、模型的F統(tǒng)計(jì)量值，用來(lái)檢驗(yàn)?zāi)Ｐ偷娘@著性、模型的信息準(zhǔn)則AIC或BIC，用來(lái)對(duì)比模型擬合效果的好壞等；第二部分主要包含偏回歸系數(shù)的信息，例如回歸系數(shù)的估計(jì)值Coef、t統(tǒng)計(jì)量值、回歸系數(shù)的置信區(qū)間等；第三部分主要涉及模型誤差項(xiàng)ε的有關(guān)信息，例如用于檢驗(yàn)誤差項(xiàng)獨(dú)立性的杜賓－瓦特森統(tǒng)計(jì)量Durbin-Watson、用于衡量誤差項(xiàng)是否服從正態(tài)分布的JB統(tǒng)計(jì)量以及有關(guān)誤差項(xiàng)偏度Skew和峰度Kurtosis的計(jì)算值等。
步驟四：對(duì)比結(jié)果下結(jié)論在第二部分的內(nèi)容中，含有每個(gè)偏回歸系數(shù)的t統(tǒng)計(jì)量值，它的計(jì)算就是由估計(jì)值coef和標(biāo)準(zhǔn)誤std err的商所得的。同時(shí)，每個(gè)t統(tǒng)計(jì)量值都對(duì)應(yīng)了概率值p，用來(lái)判別統(tǒng)計(jì)量是否顯著的直接辦法，通常概率值p小于0.05時(shí)表示拒絕原假設(shè)。從返回的結(jié)果可知，只有截距項(xiàng)Intercept和研發(fā)成本RD_Spend對(duì)應(yīng)的p值小于0.05，才說(shuō)明其余變量都沒(méi)有通過(guò)系數(shù)的顯著性檢驗(yàn)，即在模型中這些變量不是影響利潤(rùn)的重要因素。
7.4 回歸模型的診斷當(dāng)回歸模型構(gòu)建好之后，并不意味著建模過(guò)程的結(jié)束，還需要進(jìn)一步對(duì)模型進(jìn)行診斷，目的就是使診斷后的模型更加健壯。統(tǒng)計(jì)學(xué)家在發(fā)明線性回歸模型的時(shí)候就提出了一些假設(shè)前提，只有在滿足這些假設(shè)前提的情況下，所得的模型才是合理的。本節(jié)的主要內(nèi)容就是針對(duì)如下幾點(diǎn)假設(shè)，完成模型的診斷工作：
除了上面提到的五點(diǎn)假設(shè)之外，還需要注意的是，線性回歸模型對(duì)異常值是非常敏感的，即模型的構(gòu)建過(guò)程非常容易受到異常值的影響，所以診斷過(guò)程中還需要對(duì)原始數(shù)據(jù)的觀測(cè)進(jìn)行異常點(diǎn)識(shí)別和處理。接下來(lái)，結(jié)合理論知識(shí)與Python代碼逐一展開(kāi)模型的診斷過(guò)程。
7.4.1 正態(tài)性檢驗(yàn)雖然模型的前提假設(shè)是對(duì)殘差項(xiàng)要求服從正態(tài)分布，但是其實(shí)質(zhì)就是要求因變量服從正態(tài)分布。對(duì)于多元線性回歸模型y=Xβ+ε來(lái)說(shuō)，等式右邊的自變量屬于已知變量，而等式左邊的因變量為未知變量（故需要通過(guò)建模進(jìn)行預(yù)測(cè)）。所以，要求誤差項(xiàng)服從正態(tài)分布，就是要求因變量服從正態(tài)分布，關(guān)于正態(tài)性檢驗(yàn)通常運(yùn)用兩類方法，分別是定性的圖形法（直方圖、PP圖或QQ圖）和定量的非參數(shù)法（Shapiro檢驗(yàn)和K-S檢驗(yàn)），接下來(lái)通過(guò)具體的代碼對(duì)原數(shù)據(jù)集中的利潤(rùn)變量進(jìn)行正態(tài)性檢驗(yàn)。
1．直方圖法
# 正態(tài)性檢驗(yàn)# 直方圖法# 導(dǎo)入第三方模塊import scipy.stats as stats# 中文和負(fù)號(hào)的正常顯示plt.rcParams['font.sans-serif'] = ['SimHei'] # 用來(lái)正常顯示中文標(biāo)簽plt.rcParams['axes.unicode_minus'] = False # 用來(lái)正常顯示負(fù)號(hào)# 繪制直方圖sns.distplot(a = Profit_New.Profit, bins = 10, fit = stats.norm, norm_hist = True, hist_kws = {'color':'steelblue', 'edgecolor':'black'}, kde_kws = {'color':'black', 'linestyle':'--', 'label':'核密度曲線'}, fit_kws = {'color':'red', 'linestyle':':', 'label':'正態(tài)密度曲線'})# 顯示圖例plt.legend()# 顯示圖形plt.show()結(jié)果：
圖7-4中繪制了因變量Profit的直方圖、核密度曲線和理論正態(tài)分布的密度曲線，添加兩條曲線的目的就是比對(duì)數(shù)據(jù)的實(shí)際分布與理論分布之間的差異。如果兩條曲線近似或吻合，就說(shuō)明該變量近似服從正態(tài)分布。從圖中看，核密度曲線與正態(tài)密度曲線的趨勢(shì)比較吻合，故直觀上可以認(rèn)為利潤(rùn)變量服從正態(tài)分布。
2．PP圖與QQ圖
# 殘差的正態(tài)性檢驗(yàn)（PP圖和QQ圖法）pp_qq_plot = sm.ProbPlot(Profit_New.Profit)# 繪制PP圖pp_qq_plot.ppplot(line = '45')plt.title('P-P圖')# 繪制QQ圖pp_qq_plot.qqplot(line = 'q')plt.title('Q-Q圖')# 顯示圖形plt.show()結(jié)果：
PP圖的思想是比對(duì)正態(tài)分布的累計(jì)概率值和實(shí)際分布的累計(jì)概率值，而QQ圖則比對(duì)正態(tài)分布的分位數(shù)和實(shí)際分布的分位數(shù)。判斷變量是否近似服從正態(tài)分布的標(biāo)準(zhǔn)是：如果散點(diǎn)都比較均勻地散落在直線上，就說(shuō)明變量近似服從正態(tài)分布，否則就認(rèn)為數(shù)據(jù)不服從正態(tài)分布。從圖7-5可知，不管是PP圖還是QQ圖，繪制的散點(diǎn)均落在直線的附近，沒(méi)有較大的偏離，故認(rèn)為利潤(rùn)變量近似服從正態(tài)分布。
3．Shapiro檢驗(yàn)和K-S檢驗(yàn)這兩種檢驗(yàn)方法均屬于非參數(shù)方法，它們的原假設(shè)被設(shè)定為變量服從正態(tài)分布，兩者的最大區(qū)別在于使用的數(shù)據(jù)量不一樣，若數(shù)據(jù)量低于5000，則使用shapiro檢驗(yàn)法比較合理，否則使用K-S檢驗(yàn)法。scipy的子模塊stats提供了專門(mén)的檢驗(yàn)函數(shù)，分別是shapiro函數(shù)和kstest函數(shù)，由于利潤(rùn)數(shù)據(jù)集的樣本量小于5000，故下面運(yùn)用shapiro函數(shù)對(duì)利潤(rùn)做定量的正態(tài)性檢驗(yàn)：
# 導(dǎo)入模塊import scipy.stats as statsstats.shapiro(Profit_New.Profit)結(jié)果：
(0.9793398380279541, 0.537902295589447)如上結(jié)果所示，元組中的第一個(gè)元素是shapiro檢驗(yàn)的統(tǒng)計(jì)量值，第二個(gè)元素是對(duì)應(yīng)的概率值p。由于p值大于置信水平0.05，故接受利潤(rùn)變量服從正態(tài)分布的原假設(shè)。為了應(yīng)用K-S檢驗(yàn)的函數(shù)kstest，這里隨機(jī)生成正態(tài)分布變量x1和均勻分布變量x2，具體操作代碼如下：
# 生成正態(tài)分布和均勻分布隨機(jī)數(shù)rnorm = np.random.normal(loc = 5, scale=2, size = 10000)runif = np.random.uniform(low = 1, high = 100, size = 10000)# 正態(tài)性檢驗(yàn)KS_Test1 = stats.kstest(rvs = rnorm, args = (rnorm.mean(), rnorm.std()), cdf = 'norm')KS_Test2 = stats.kstest(rvs = runif, args = (runif.mean(), runif.std()), cdf = 'norm')print(KS_Test1)print(KS_Test2)結(jié)果：
KstestResult(statistic=0.0054919289086649, pvalue=0.9236250840080669)KstestResult(statistic=0.05906629005008068, pvalue=9.941845404073026e-31)如上結(jié)果所示，正態(tài)分布隨機(jī)數(shù)的檢驗(yàn)p值大于置信水平0.05，則需接受原假設(shè)；均勻分布隨機(jī)數(shù)的檢驗(yàn)p值遠(yuǎn)遠(yuǎn)小于0.05，則需拒絕原假設(shè)。需要說(shuō)明的是，如果使用kstest函數(shù)對(duì)變量進(jìn)行正態(tài)性檢驗(yàn)，必須指定args參數(shù)，它用于傳遞被檢驗(yàn)變量的均值和標(biāo)準(zhǔn)差。 如果因變量的檢驗(yàn)結(jié)果不滿足正態(tài)分布時(shí)，需要對(duì)因變量做某種數(shù)學(xué)轉(zhuǎn)換，使用比較多的轉(zhuǎn)換方法有:
7.4.2 多重共線性檢驗(yàn)多重共線性是指模型中的自變量之間存在較高的線性相關(guān)關(guān)系，它的存在會(huì)給模型帶來(lái)嚴(yán)重的后果，例如由“最小二乘法”得到的偏回歸系數(shù)無(wú)效、增大偏回歸系數(shù)的方差、模型缺乏穩(wěn)定性等，所以，對(duì)模型的多重共線性檢驗(yàn)就顯得尤其重要了。
關(guān)于多重共線性的檢驗(yàn)可以使用方差膨脹因子VIF來(lái)鑒定，如果VIF大于10，則說(shuō)明變量間存在多重共線性；如果VIF大于100，則表名變量間存在嚴(yán)重的多重共線性。方差膨脹因子VIF的計(jì)算步驟如下：（1）構(gòu)造每一個(gè)自變量與其余自變量的線性回歸模型，例如，數(shù)據(jù)集中含有p個(gè)自變量，則第一個(gè)自變量與其余自變量的線性組合可以表示為：
x1=c0+α2x1+…+αpxp+ε（2）根據(jù)如上線性回歸模型得到相應(yīng)的判決系數(shù)R2，進(jìn)而計(jì)算第一個(gè)自變量的方差膨脹因子VIF：
Python中的statsmodels模塊提供了計(jì)算方差膨脹因子VIF的函數(shù)，下面利用該函數(shù)計(jì)算兩個(gè)自變量的方差膨脹因子：
# 導(dǎo)入statsmodels模塊中的函數(shù)from statsmodels.stats.outliers_influence import variance_inflation_factor# 自變量X(包含RD_Spend、Marketing_Spend和常數(shù)列1)X = sm.add_constant(Profit_New.ix[:,['RD_Spend','Marketing_Spend']])# 構(gòu)造空的數(shù)據(jù)框，用于存儲(chǔ)VIF值vif = pd.DataFrame()vif['features'] = X.columnsvif['VIF Factor'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]# 返回VIF值vif結(jié)果：
如上結(jié)果所示，兩個(gè)自變量對(duì)應(yīng)的方差膨脹因子均低于10，說(shuō)明構(gòu)建模型的數(shù)據(jù)并不存在多重共線性。如果發(fā)現(xiàn)變量之間存在多重共線性的話，可以考慮刪除變量或者重新選擇模型（如嶺回歸模型或LASSO模型）。
7.4.3 線性相關(guān)性檢驗(yàn)線性相關(guān)性檢驗(yàn)，顧名思義，就是確保用于建模的自變量和因變量之間存在線性關(guān)系。關(guān)于線性關(guān)系的判斷，可以使用Pearson相關(guān)系數(shù)和可視化方法進(jìn)行識(shí)別，有關(guān)Pearson相關(guān)系數(shù)的計(jì)算公式如下：
其中，COV(x,y)為自變量x與因變量y之間的協(xié)方差，D(x)和D(y)分別為自變量x和因變量y的方差。Pearson相關(guān)系數(shù)的計(jì)算可以直接使用數(shù)據(jù)框的corrwith“方法”，該方法最大的好處是可以計(jì)算任意指定變量間的相關(guān)系數(shù)。下面使用該方法計(jì)算因變量與每個(gè)自變量之間的相關(guān)系數(shù)，具體代碼如下：
# 計(jì)算數(shù)據(jù)集Profit_New中每個(gè)自變量與因變量利潤(rùn)之間的相關(guān)系數(shù)Profit_New.drop('Profit', axis = 1).corrwith(Profit_New.Profit)結(jié)果：
RD_Spend 0.978437Administration 0.205841Marketing_Spend 0.739307California -0.083258Florida 0.088008如上結(jié)果所示，自變量中只有研發(fā)成本和市場(chǎng)營(yíng)銷成本與利潤(rùn)之間存在較高的相關(guān)系數(shù)，相關(guān)性分別達(dá)到0.978和0.739，而其他變量與利潤(rùn)之間幾乎沒(méi)有線性相關(guān)性可言。通常情況下，可以參考表7-3判斷相關(guān)系數(shù)對(duì)應(yīng)的相關(guān)程度：
以管理成本Administration為例，與利潤(rùn)之間的相關(guān)系數(shù)只有0.2，被認(rèn)定為不相關(guān)，這里的不相關(guān)只能說(shuō)明兩者之間不存在線性關(guān)系。如果利潤(rùn)和管理成本之間存在非線性關(guān)系時(shí)，Pearson相關(guān)系數(shù)也同樣會(huì)很小，所以還需要通過(guò)可視化的方法，觀察自變量與因變量之間的散點(diǎn)關(guān)系。讀者可以應(yīng)用matplotlib模塊中的scatter函數(shù)繪制五個(gè)自變量與因變量之間的散點(diǎn)圖，那樣做可能會(huì)使代碼顯得冗長(zhǎng)。這里介紹另一個(gè)繪制散點(diǎn)圖的函數(shù)，那就是seaborn模塊中的pairplot函數(shù)，它可以繪制多個(gè)變量間的散點(diǎn)圖矩陣。
# 散點(diǎn)圖矩陣# 導(dǎo)入模塊import matplotlib.pyplot as pltimport seaborn# 繪制散點(diǎn)圖矩陣seaborn.pairplot(Profit_New.ix[:,['RD_Spend','Administration','Marketing_Spend','Profit']])# 顯示圖形結(jié)果：
如圖7-6所示，由于California與Florida都是啞變量，故沒(méi)有將其放入散點(diǎn)圖矩陣中。從圖中結(jié)果可知，研發(fā)成本與利潤(rùn)之間的散點(diǎn)圖幾乎為一條向上傾斜的直線（見(jiàn)左下角的散點(diǎn)圖），說(shuō)明兩種變量確實(shí)存在很強(qiáng)的線性關(guān)系；市場(chǎng)營(yíng)銷成本與利潤(rùn)的散點(diǎn)圖同樣向上傾斜，但很多點(diǎn)的分布還是比較分散的（見(jiàn)第一列第三行的散點(diǎn)圖）；管理成本與利潤(rùn)之間的散點(diǎn)圖呈水平趨勢(shì)，而且分布也比較寬，說(shuō)明兩者之間確實(shí)沒(méi)有任何關(guān)系（見(jiàn)第一列第二行的散點(diǎn)圖）。以7.2.2節(jié)中重構(gòu)的模型model2為例，綜合考慮相關(guān)系數(shù)、散點(diǎn)圖矩陣和t檢驗(yàn)的結(jié)果，最終確定只保留模型model2中的RD_Spend和Marketing_Spend兩個(gè)自變量，下面重新對(duì)該模型做修正：
# 模型修正model3 = sm.formula.ols('Profit ~ RD_Spend + Marketing_Spend', data = train).fit()# 模型回歸系數(shù)的估計(jì)值model3.params結(jié)果：
Intercept 51902.112471RD_Spend 0.785116Marketing_Spend 0.019402dtype: float64如上結(jié)果所示，返回的是模型兩個(gè)自變量的系數(shù)估計(jì)值，可以將多元線性回歸模型表示成：
Profit=51902.11+0.79RD_Spend+0.02Marketing_Spend7.4.4 異常值檢驗(yàn)
由于多元線性回歸模型容易受到極端值的影響，故需要利用統(tǒng)計(jì)方法對(duì)觀測(cè)樣本進(jìn)行異常點(diǎn)檢測(cè)。如果在建模過(guò)程中發(fā)現(xiàn)異常數(shù)據(jù)，需要對(duì)數(shù)據(jù)集進(jìn)行整改，如刪除異常值或衍生出是否為異常值的啞變量。對(duì)于線性回歸模型來(lái)說(shuō)，通常利用帽子矩陣、DFFITS準(zhǔn)則、學(xué)生化殘差或Cook距離進(jìn)行異常點(diǎn)檢測(cè)。接下來(lái)，分別對(duì)這四種檢測(cè)方法做簡(jiǎn)單介紹。
1．帽子矩陣帽子矩陣的設(shè)計(jì)思路就是考察第i個(gè)樣本對(duì)預(yù)測(cè)值:
的影響大小，根據(jù)7.2.1節(jié)中推導(dǎo)得到的回歸系數(shù)求解公式，可以將多元線性回歸模型表示成：
其中，H=X(X’X)-1X’，H就稱為帽子矩陣，全都是關(guān)于自變量X的計(jì)算。判斷樣本是否為異常點(diǎn)的方法，可以使用如下公式：
其中，hii為帽子矩陣H的第i個(gè)主對(duì)角線元素，p為自變量個(gè)數(shù)，n為用于建模數(shù)據(jù)集的樣本量。如果對(duì)角線元素滿足上面的公式，則代表第i個(gè)樣本為異常觀測(cè)。
2．DFFITS準(zhǔn)則DFFITS準(zhǔn)則借助于帽子矩陣，構(gòu)造了另一個(gè)判斷異常點(diǎn)的統(tǒng)計(jì)量，該統(tǒng)計(jì)量可以表示成如下公式：
其中，hii為帽子矩陣H的第i個(gè)主對(duì)角線元素，εi為第i個(gè)樣本點(diǎn)的預(yù)測(cè)誤差，σ為誤差項(xiàng)的標(biāo)準(zhǔn)差，判斷樣本為異常點(diǎn)的方法，可以使用如下規(guī)則：
需要注意的是，在DFFITS準(zhǔn)則的公式中，乘積的第二項(xiàng)實(shí)際上是學(xué)生化殘差，也可以用來(lái)判定第i個(gè)樣本是否為異常點(diǎn)，判斷標(biāo)準(zhǔn)如下：
3．Cook距離Cook距離是一種相對(duì)抽象的判斷準(zhǔn)則，無(wú)法通過(guò)具體的臨界值判斷樣本是否為異常點(diǎn)，對(duì)于該距離，Cook統(tǒng)計(jì)量越大的點(diǎn)，其成為異常點(diǎn)的可能性越大。Cook統(tǒng)計(jì)量可以表示為如下公式：
其中，ri為學(xué)生化殘差。如果使用如上四種方法判別數(shù)據(jù)集的第i個(gè)樣本是否為異常點(diǎn)，前提是已經(jīng)構(gòu)造好一個(gè)線性回歸模型，然后基于get_influence“方法”獲得四種統(tǒng)計(jì)量的值。為了檢驗(yàn)?zāi)Ｐ椭袛?shù)據(jù)集的樣本是否存在異常，這里沿用上節(jié)中構(gòu)造的模型model3，具體代碼如下：
# 異常值檢驗(yàn)outliers = model3.get_influence()# 高杠桿值點(diǎn)（帽子矩陣）leverage = outliers.hat_matrix_diag# dffits值dffits = outliers.dffits[0]# 學(xué)生化殘差resid_stu = outliers.resid_studentized_external# cook距離cook = outliers.cooks_distance[0]# 合并各種異常值檢驗(yàn)的統(tǒng)計(jì)量值contat1 = pd.concat([pd.Series(leverage, name = 'leverage'),pd.Series(dffits, name = 'dffits'), pd.Series(resid_stu,name = 'resid_stu'),pd.Series(cook, name = 'cook')],axis = 1)# 重設(shè)train數(shù)據(jù)的行索引train.index = range(train.shape[0])# 將上面的統(tǒng)計(jì)量與train數(shù)據(jù)集合并profit_outliers = pd.concat([train,contat1], axis = 1)profit_outliers.head()結(jié)果:
如表7-4所示，合并了train數(shù)據(jù)集和四種統(tǒng)計(jì)量的值，接下來(lái)要做的就是選擇一種或多種判斷方法，將異常點(diǎn)查詢出來(lái)。為了簡(jiǎn)單起見(jiàn)，這里使用標(biāo)準(zhǔn)化殘差，當(dāng)標(biāo)準(zhǔn)化殘差大于2時(shí)，即認(rèn)為對(duì)應(yīng)的數(shù)據(jù)點(diǎn)為異常值。
# 計(jì)算異常值數(shù)量的比例outliers_ratio = sum(np.where((np.abs(profit_outliers.resid_stu)>2),1,0))/profit_outliers.shape[0]outliers_ratio結(jié)果：
0.02564102564102564如上結(jié)果所示，通過(guò)標(biāo)準(zhǔn)化殘差監(jiān)控到了異常值，并且異常比例為2.5%。對(duì)于異常值的處理辦法，可以使用兩種策略，如果異常樣本的比例不高（如小于等于5%），可以考慮將異常點(diǎn)刪除；如果異常樣本的比例比較高，選擇刪除會(huì)丟失一些重要信息，所以需要衍生啞變量，即對(duì)于異常點(diǎn)，設(shè)置啞變量的值為1，否則為0。如上可知，建模數(shù)據(jù)的異常比例只有2.5%，故考慮將其刪除。
# 挑選出非異常的觀測(cè)點(diǎn)none_outliers = profit_outliers.ix[np.abs(profit_outliers.resid_stu)<=2,]# 應(yīng)用無(wú)異常值的數(shù)據(jù)集重新建模model4 = sm.formula.ols('Profit ~ RD_Spend + Marketing_Spend', data = none_outliers).fit()model4.params結(jié)果：
Intercept 51827.416821RD_Spend 0.797038Marketing_Spend 0.017740dtype: float64如上結(jié)果所示，經(jīng)過(guò)異常點(diǎn)的排除，重構(gòu)模型的偏回歸系數(shù)發(fā)生了變動(dòng)，故可以將模型寫(xiě)成如下公式：
Profit=51827.42+0.80RD_Spend+0.02Marketing_Spend7.4.5 獨(dú)立性檢驗(yàn)殘差的獨(dú)立性檢驗(yàn)，說(shuō)白了也是對(duì)因變量y的獨(dú)立性檢驗(yàn)，因?yàn)樵诰€性回歸模型的等式左右只有y和殘差項(xiàng)ε屬于隨機(jī)變量，如果再加上正態(tài)分布，就構(gòu)成了殘差項(xiàng)獨(dú)立同分布于正態(tài)分布的假設(shè)。關(guān)于殘差的獨(dú)立性檢驗(yàn)通常使用Durbin-Watson統(tǒng)計(jì)量值來(lái)測(cè)試，如果DW值在2左右，則表明殘差項(xiàng)之間是不相關(guān)的；如果與2偏離的較遠(yuǎn)，則說(shuō)明不滿足殘差的獨(dú)立性假設(shè)。對(duì)于DW統(tǒng)計(jì)量的值，其實(shí)都不需要另行計(jì)算，因?yàn)樗谀Ｐ偷母庞[信息中，以模型model4為例：
# Durbin-Watson統(tǒng)計(jì)量# 模型概覽model4.summary()結(jié)果:
如表7-5所示，殘差項(xiàng)對(duì)應(yīng)的DW統(tǒng)計(jì)量值為2.065，比較接近于2，故可以認(rèn)為模型的殘差項(xiàng)之間是滿足獨(dú)立性這個(gè)假設(shè)前提的。
7.4.6 方差齊性檢驗(yàn)方差齊性是要求模型殘差項(xiàng)的方差不隨自變量的變動(dòng)而呈現(xiàn)某種趨勢(shì)，否則，殘差的趨勢(shì)就可以被自變量刻畫(huà)。如果殘差項(xiàng)不滿足方差齊性（方差為一個(gè)常數(shù)），就會(huì)導(dǎo)致偏回歸系數(shù)不具備有效性，甚至導(dǎo)致模型的預(yù)測(cè)也不準(zhǔn)確。所以，建模后需要驗(yàn)證殘差項(xiàng)是否滿足方差齊性。關(guān)于方差齊性的檢驗(yàn)，一般可以使用兩種方法，即圖形法（散點(diǎn)圖）和統(tǒng)計(jì)檢驗(yàn)法（BP檢驗(yàn)）。
1．圖形法如上所說(shuō)，方差齊性是指殘差項(xiàng)的方差不隨自變量的變動(dòng)而變動(dòng)，所以只需要繪制殘差與自變量之間的散點(diǎn)圖，就可以發(fā)現(xiàn)兩者之間是否存在某種趨勢(shì)：
# 方差齊性檢驗(yàn)# 設(shè)置第一張子圖的位置ax1 = plt.subplot2grid(shape = (2,1), loc = (0,0))# 繪制散點(diǎn)圖ax1.scatter(none_outliers.RD_Spend, (model4.resid-model4.resid.mean())/model4.resid.std())# 添加水平參考線ax1.hlines(y = 0 ,xmin = none_outliers.RD_Spend.min(),xmax = none_outliers.RD_Spend.max(), color = 'red', linestyles = '--')# 添加x軸和y軸標(biāo)簽ax1.set_xlabel('RD_Spend')ax1.set_ylabel('Std_Residual')# 設(shè)置第二張子圖的位置ax2 = plt.subplot2grid(shape = (2,1), loc = (1,0))# 繪制散點(diǎn)圖ax2.scatter(none_outliers.Marketing_Spend, (model4.resid-model4.resid.mean())/model4.resid.std())# 添加水平參考線ax2.hlines(y = 0 ,xmin = none_outliers.Marketing_Spend.min(),xmax = none_outliers.Marketing_Spend.max(), color = 'red', linestyles = '--')# 添加x軸和y軸標(biāo)簽ax2.set_xlabel('Marketing_Spend')ax2.set_ylabel('Std_Residual')# 調(diào)整子圖之間的水平間距和高度間距plt.subplots_adjust(hspace=0.6, wspace=0.3)# 顯示圖形plt.show()結(jié)果：
如圖7-7所示，標(biāo)準(zhǔn)化殘差并沒(méi)有隨自變量的變動(dòng)而呈現(xiàn)喇叭形，所有的散點(diǎn)幾乎均勻地分布在參考線y=0的附近。所以，可以說(shuō)明模型的殘差項(xiàng)滿足方差齊性的前提假設(shè)。
2．BP檢驗(yàn)方差齊性檢驗(yàn)的另一個(gè)統(tǒng)計(jì)方法是BP檢驗(yàn)，它的原假設(shè)是殘差的方差為一個(gè)常數(shù)，通過(guò)構(gòu)造拉格朗日乘子LM統(tǒng)計(jì)量，實(shí)現(xiàn)方差齊性的檢驗(yàn)。該檢驗(yàn)可以借助于statsmodels模塊中的het_breushpagan函數(shù)完成，具體代碼如下：
# BP檢驗(yàn)sm.stats.diagnostic.het_breushpagan(model4.resid, exog_het = model4.model.exog)結(jié)果：
(1.467510366830809, 0.48010272699006995, 0.7029751237162342, 0.5019659740962923)如上結(jié)果所示，元組中一共包含四個(gè)值，第一個(gè)值1.468為L(zhǎng)M統(tǒng)計(jì)量；第二個(gè)值是統(tǒng)計(jì)量對(duì)應(yīng)的概率p值，該值大于0.05，說(shuō)明接受殘差方差為常數(shù)的原假設(shè)；第三個(gè)值為F統(tǒng)計(jì)量，用于檢驗(yàn)殘差平方項(xiàng)與自變量之間是否獨(dú)立，如果獨(dú)立則表明殘差方差齊性；第四個(gè)值則為F統(tǒng)計(jì)量的概率p值，同樣大于0.05，則進(jìn)一步表示殘差項(xiàng)滿足方差齊性的假設(shè)。如果模型的殘差不滿足齊性的話，可以考慮兩類方法來(lái)解決，一類是模型變換法，另一類是“加權(quán)最小二乘法”（可以使用statsmodels模塊中的wls函數(shù)）。對(duì)于模型變換法來(lái)說(shuō)，主要考慮殘差與自變量之間的關(guān)系，如果殘差與某個(gè)自變量x成正比，則需將原模型的兩邊同除以；如果殘差與某個(gè)自變量x的平方成正比，則需將原始模型的兩邊同除以x；對(duì)于加權(quán)最小二乘法來(lái)說(shuō)，關(guān)鍵是如何確定權(quán)重，根據(jù)多方資料的搜索和驗(yàn)證，一般選擇如下三種權(quán)重來(lái)進(jìn)行對(duì)比測(cè)試：
3．回歸模型的預(yù)測(cè)經(jīng)過(guò)前文的模型構(gòu)造、假設(shè)檢驗(yàn)和模型診斷，最終確定合理的模型model4。接下來(lái)要做的就是利用該模型完成測(cè)試集上的預(yù)測(cè)，具體代碼如下：
# 模型預(yù)測(cè)# model4對(duì)測(cè)試集的預(yù)測(cè)pred4 = model4.predict(exog = test.ix[:,['RD_Spend','Marketing_Spend']])# 對(duì)比預(yù)測(cè)值與實(shí)際值pd.DataFrame({'Prediction':pred4,'Real':test.Profit})# 繪制預(yù)測(cè)值與實(shí)際值的散點(diǎn)圖plt.scatter(x = test.Profit, y = pred4)# 添加斜率為1，截距項(xiàng)為0的參考線plt.plot([test.Profit.min(),test.Profit.max()],[test.Profit.min(),test.Profit.max()], color = 'red', linestyle = '--')# 添加軸標(biāo)簽plt.xlabel('實(shí)際值')plt.ylabel('預(yù)測(cè)值')# 顯示圖形plt.show()結(jié)果：
如圖7-8所示，繪制了有關(guān)模型在測(cè)試集上的預(yù)測(cè)值和實(shí)際值的散點(diǎn)圖，該散點(diǎn)圖可以用來(lái)衡量預(yù)測(cè)值與實(shí)際值之間的距離差異。如果兩者非常接近，那么得到的散點(diǎn)圖一定會(huì)在對(duì)角線附近微微波動(dòng)。從圖7-8的結(jié)果來(lái)看，大部分的散點(diǎn)都落在對(duì)角線附近，說(shuō)明模型的預(yù)測(cè)效果還是不錯(cuò)的。
7.5 本章小結(jié)本章重點(diǎn)介紹了有關(guān)線性回歸模型的理論知識(shí)與應(yīng)用實(shí)戰(zhàn)，內(nèi)容包含模型回歸系數(shù)的求解、模型及系數(shù)的顯著性檢驗(yàn)、模型的假設(shè)診斷及預(yù)測(cè)。在實(shí)際應(yīng)用中，如果因變量為數(shù)值型變量，可以考慮使用線性回歸模型，但是前提得滿足幾點(diǎn)假設(shè)，如因變量服從正態(tài)分布、自變量間不存在多重共線性、自變量與因變量之間存在線性關(guān)系、用于建模的數(shù)據(jù)集不存在異常點(diǎn)、殘差項(xiàng)滿足方差異性和獨(dú)立性。 為了使讀者掌握有關(guān)線性回歸模型的函數(shù)和“方法”，這里將其重新梳理一下，以便讀者查閱和記憶：
聯(lián)系客服
