python多元回歸分析案例(Python數據分析系列 15.多元回歸分析)

數據與智能 本公眾號關注大數據與人工智能技術。由一批具備多年實戰經驗的技術極客參與運營管理，持續輸出大數據、數據分析、推薦系統、機器學習、人工智能等方向的原創文章，每周至少輸出7篇精品原創。同時，我們會關注和分享大數據與人工智能行業動態。歡迎關注。

來源 | Data Science from Scratch， Second Edition

作者 | Joel Grus

譯者 | cloverErna

校對 | gongyouliu

編輯 | auroral-L

全文共4710字，預計閱讀時間40分鐘。

第十五章 多元回歸分析

1. 模型

2. 最小二乘模型的進一步假設

3. 擬合模型

4. 解釋模型

5. 擬合優度

6. 題外話：Bootstrap

7. 回歸系數的標準誤差

8. 正則化

9. 延伸學習

我不會盯著一個問題看，并往其中添加無用的變量。

——比爾· 帕塞爾斯

雖然副總對你的預測模型很滿意，但是他認為你還可以做得更好。為此，你收集了額外的數據：對于每一個用戶，你不僅了解他每天工作多少小時，同時還調查了他是否擁有博士學位。你希望通過這些補充資料來改進模型。

因此，你提出了一個帶有更多自變量的線性模型：

顯然，用戶是否擁有博士學位并非一個數值問題，但如同第 11 章所提到的，我們可以引入一個虛擬變量，當這個變量等于 1 的時候，表示用戶擁有博士學位，反之則表示沒有博士學位，這樣就能像其他變量一樣將其視為一個數值了。

1. 模型

回想一下，我們在第 14 章中所擬合的模型形式如下所示：

現在，如果每個輸入 xi 不再是單個數字，而是一個由 k 個數字 xi 1,…, xi k 組成的向量，那么我們的多元回歸模型則應該為：

對于多元回歸分析來說，參數向量通常稱為 β。我們希望這個向量也包括一個常數項，為此，只要向我們的數據中添加一列即可：

同時：

那么，我們的模型可以用下列函數實現：

在這種特殊情況下，我們的自變量x將是一個向量列表，每個向量都是這樣的：

2. 最小二乘模型的進一步假設

對于我們的模型（以及我們的解決方案）來說，需要添加另外兩個假設，才能夠言之有理。第一個假設是 x 的各個列是線性無關的，即任何一列絕對不會是其他列的加權和。如果這個假設不成立，則無法估計 beta。為了了解極端的情形，我們可以想象數據中有一個額外的字段 num_acquaintances，并且對于每一個用戶來說它都等于 num_friends。

那么，對于任何 beta，如果 num_friends 的系數增加了某個數值，而num_acquaintances 的系數同時減小相同數值的話，那么模型的預測就會保持不變。也就是說，我們根本就沒有辦法確定 num_friends 的系數。（通常來說，對于這個假設的違背情況一般不會這么明顯。）

第二個重要的假設是 x 的各列與誤差 ε 無關。如果這個假設不成立，對于 beta 的估計就會出現系統性的錯誤。

比如，在第 14 章中，我們建立的模型的預測結果為，用戶每增加一個朋友，每天花在網站上的時間就會多出 0.90 分鐘。

想象一下，還有下列情況：

? 工作時間越長的人在網站上花的時間越少。

? 朋友更多的人傾向于工作更長時間。

也就是說，假設“實際的”模型為：

其中β2是負的，而工作時間和朋友則是正相關的。在這種情況下，當我們最小化單變量模型的誤差時：

我們會低估 β1。

考慮一下，如果這個單變量模型已知 β1 的“實際”值，那么這時再用它來預測將會如何。（亦即，這個值來自令誤差最小化的“實際”模型。）這時候，對工作時間比較長的用戶來說，產生的預測值往往太小；對工作時間比較短的用戶來說，產生的預測值往往又過大，這是因為 β2>0，但是我們“忘了”把它考慮進去。由于工作時間與朋友的數量是呈正相關的，這就意味著對于朋友數量較多的用戶來說，模型給出的預測值往往太小；對于朋友數量較少的用戶來說，模型給出的預測值往往太大。

這樣做的結果是，我們可以通過降低 β1的估計值來減少（單變量模型）的誤差，即誤差最小化的 β1 是小于其“實際”值的。也就是說，在這種情況下，單變量的最小二乘解偏向于低估 β1。一般而言，當自變量具有與此類似的誤差時，我們的最小二乘解給出的 β 是有偏估計。

3. 擬合模型

就像對簡單線性模型所做的那樣，我們這里也需要尋找一個能夠最小化誤差的平方和的beta。要想以手動方式找到一個精確的解可不是一件容易的事，因此，我們轉而求助于梯度下降。下面我們首先創建一個待最小化的誤差函數。誤差函數幾乎與我們在第14章中使用的完全相同，除了它不期望參數，而是需要一個任意長度的向量：

如果你熟悉微積分，可以通過下面的方式進行計算：

否則的話，就按照我說的來。

至此，我們就可以利用隨機梯度下降法來尋找最優的 beta 了。讓我們首先寫出一個可以用于任何數據集的least_squares_fit函數：

然后，我們可以將其應用于我們的數據：

在實際實踐中，你不會使用梯度下降來估計線性回歸；你將使用超出本書范圍的線性代數技術得到精確的系數。如果你這樣做了，你就會得到以下公式：

這和我們發現的很接近。

4. 解釋模型

你應該把模型的系數看作在其他條件相同的情況下每個因素的影響力的大小。在其他條件相同的情況下，每增加一個朋友，每天花在網站上的時間就會多出一分鐘。在其他條件相同的情況下，用戶在工作日每多工作一個小時，每天花在網站上的時間就會減少兩分鐘。在其他條件相同的情況下，擁有博士學位的用戶每天用在網站上的時間會多出一分鐘。

但是，它沒有（直接）反映變量之間的任何相互作用。較之于朋友較少的人而言，工作時間對朋友較多的人的影響很可能是不一樣的，而這個模型并沒有捕捉到這一點。要想處理這種情況，一種方法是引入一個新變量，即“朋友數量”與“工作時間”之積。這樣實際上就使得“工作時間”的系數可以隨著朋友數量的增加而增加（或減少）。

還有一種可能，就是朋友越多，花在網站上的時間就越多，但是達到一個上限之后，更多的朋友反而會導致花在網站上的時間變少。（或許是因為朋友太多了反而上網體驗會很糟糕？）我們可以設法讓模型捕獲到這一點，方法是添加另一個變量，即朋友數量的平方。

一旦我們開始添加變量，我們就需要考慮它們的系數“問題”。對于添加的乘積、對數、二次冪以及更高次冪來說，其數量上是沒有限制的。

5. 擬合優度

現在，我們再來看看 R 的平方值。

現在已經增加到0.68個：

但是不要忘了，只要向回歸模型中添加新的變量就必然導致 R 的平方變大。歸根結底，前面的簡單回歸模型只是這里的多元回歸模型的特例而已，即“工作時間”和“博士”這兩列的系數都等于 0。因此，最優的多元回歸模型，其誤差一定不會高于簡單回歸模型。

因此，對于多元回歸分析而言，我們還要考察系數的標準誤差，即衡量每個 βi 的估計值的可靠程度。

總的來說，回歸模型通常能夠很好地擬合我們的數據，但是，如果某些自變量是相關的（或不相關的），那么其系數就未必有多大的意義了。

對于這些誤差，傳統的度量方法通常都帶有一個前提假設，即誤差 εi 是獨立的正態隨機變量，其平均值為 0，標準偏差為σ（未知）。那樣的話，我們（或者說我們的統計軟件）就可以使用線性代數來確定每個系數的標準誤差了。這個誤差越大，說明模型的系數越不靠譜。令人遺憾的是，我們不打算從零開始介紹這類線性代數。

6. 題外話：Bootstrap

假設我們有一個含有 n 個數據點的樣本，并且這些點是按照某種（我們不知道的）概率分布生成的：

在第 5 章中，我們曾經編寫了一個計算觀測數據中位數的函數，現在拿它來估算該分布本身的中位數。

但是，我們該如何了解這些估計值的可靠性呢？如果樣品中所有的數據都非常接近 100，則實際的中位數很可能也非常接近 100。如果樣本中一半左右的數據接近 0，而另一半則接近 200，那么我們就很難確信中位數到底接近多少。

如果我們能夠不斷獲得新的樣本，那么就可以計算出每個新樣本的中位數，并觀察這些中位數的分布情況。但是，一般這是不現實的。相反，我們可以利用 Bootstrap 來獲得新的數據集，即選擇 n 個數據點并用原來的數據將其替換，然后計算合成的數據集的中位數：

例如，考慮下列兩個數據集：

如果你計算每個數據集的中位數，會發現它們都非常接近 100。然而，如果你考察下面的語句：

大部分情況下你看到的數字確實非常接近 100。然而，如果你考察下面的語句：

你會發現，不僅有許多數字接近 0，而且還有許多數字接近 200。

第一組中位數的 standard_deviation 接近 0，而第二組中位數的 standard_deviation 接近100:

（這種極端的情況通過人工檢查數據很容易弄清楚，但一般情況下都不是真的。）

7. 回歸系數的標準誤差

我們可以采用同樣的方法來估計回歸系數的標準誤差。我們可以對數據重復采用bootstrap_sample 樣本，并根據這些樣本估算 beta。如果某個自變量（如 num_friends）的系數在各個樣本上變化不大，那么就可以確信我們的估計是比較嚴密的。如果這個系數隨著樣本的不同而起伏較大，那么我們就不能完全相信我們的估計。

唯一需要說明的是，采樣前，我們需要把數據 X 和數據 Y 放到一起（用 zip），以確保對自變量和因變量一起進行采樣。這就意味著 bootstrap_sample 將返回一個由(x_i, y_i)數據對組成的列表，因此我們需要將其重新組合成一個 x_sample 和一個 y_sample：

之后，我們就可以估算每個系數的標準偏差了：

（如果我們收集了100多個樣本并使用5000多個迭代來估計每個測試版，我們可能會得到更好的估計，但我們沒有一整天時間。）

我們可以使用它們來檢驗諸如“βi 等于 0 嗎？”之類的假設。在滿足 βi =0（以及與 εi 分布有關的其他假設）的條件下，則有：

也就是說，這個統計量等于我們估算的 βj 除以估算的其標準誤差，它符合具有“n-k 個自由度”的學生 t 分布（Student’s t-distribution）。

如果我們有一個 students_t_cdf 函數，那么就可以計算每個最小二乘系數的 p 值，從而指出實際的系數為 0 時觀察到這個值的可能性有多大。令人遺憾的是，實際上我們沒有這樣的函數。（雖然我們不想從頭做起。）

然而，隨著自由度變大，t 分布越接近標準正態分布。在這種情況下，即 n 比 k 大得多的情況下，我們便可以使用 normal_cdf 了，并且我們覺得它效果還不錯：

（在其他情況下，我們很可能會使用一個知道如何計算 t 分布和精確的標準誤差的統計

軟件。）

雖然大多數系數的 p 值都非常小（但非 0 值），但是“博士學位”的系數與零沒有“顯著”區別，也就是說“博士學位”的系數很可能是隨機的，無意義的。

在對回歸分析要求更加精細的情形下，你可能需要對數據的各種假設進行更加細致的測試，比如“至少有一個 βj 是非 0 值”，或者“β1 等于 β2 且 β3 等于 β4”等，以便進行 F 測試，但是這些內容已經超出了本書的討論范圍。

8. 正則化

在實踐中，線性回歸經常需要處理具有很多變量的數據集，這時就需要用到另外兩個技巧。首先，涉及的變量越多，模型越容易對訓練集產生過擬合現象。其次，非零系數越多，越難以搞清楚它們的意義。如果我們的目標是解釋某些現象，一個只考慮三方面因素的稀疏型模型通常要比涉及數百個因素的模型要更好一些。

正則化是指給誤差項添加一個懲罰項，并且該懲罰項會隨著 beta 的增大而增大。然后，我們開始設法將誤差項和懲罰項的組合值最小化。因此，懲罰項權重越大，就越能防止系數過大。

例如，在嶺回歸（ridge regression）中，我們添加了一個與 beta_i 的平方之和成正比的懲罰項。（當然，我們一般不會懲罰 beta_0，因為它是個常數項。）

然后，我們可以用通常的方式將其插入到梯度下降中：

然后我們只需要修改least_squares_fit函數來使用sqerror_ridge_gradient而不是sqerror_gradient。（我不會在這里重復該代碼。）

如果alpha設置為0，根本沒有懲罰，我們得到與以前相同的結果：

隨著 alpha 的增大，擬合優度會變差，但是 beta 會變小：

特別地，隨著懲罰項的增大，“博士學位”的系數會變成 0，這與我們之前的結果是一致的，即它與 0 沒有顯著區別。

注意

在利用這個方法之前，通常需要調整數據的規模。因為即使是同一個模型，如果將幾年數據一下變為幾百年的數據，那么它的最小二乘法系數就會增加上百倍，那樣得到的懲罰肯定也會驟增。

還有一個方法是 lasso 回歸，它用的懲罰方式如下所示：

總的說來，嶺回歸的懲罰項會縮小系數，但是，lasso 的懲罰項卻趨向于迫使系數變為 0值，這使得它更適于學習稀疏模型。令人遺憾的是，它不適用于梯度下降法，這意味著我們將無法從頭開始解決這個問題。

9. 延伸學習

? 回歸分析具有深厚而廣闊的理論背景。要想了解這些背景理論，你需要閱讀相應的教科

書，至少也得閱讀大量的維基百科文章。

? scikit-learn 的 linear_model 模塊提供了一個 LinearRegression 模型，它跟我們的模型頗為相近。此外，它還提供了Ridge 回歸和 Lasso 回歸，以及其他類型的正則化算法。

? 另一個相關的 Python 模塊是 Statsmodels（https://www.statsmodels.org/），它也包含了線性回歸模型及許多其他內容。