統(tǒng)計學(xué)11-多元線性回歸

在統(tǒng)計學(xué)10-回歸一文中介紹了一元線性回歸的概念。假設(shè)我們現(xiàn)在有多個解釋變量，如何構(gòu)造多元線性回歸模型呢？

第一個實例

在這個例子里，我們不僅考慮房子的大小，還考慮其他的解釋變量：
neighborhood：分類變量
bedrooms：數(shù)值變量
bathrooms：數(shù)值變量
style：分類變量

  • 導(dǎo)入庫和數(shù)據(jù)
import numpy as np
import pandas as pd
import statsmodels.api as sms;

df = pd.read_csv('./house_prices.csv')
  • 添加截距
df['intercept'] = 1
  • 分別用簡單線性回歸擬合回歸模型
lm = sms.OLS(df['price'], df[['intercept','area']])
results = lm.fit()
results.summary()

lm = sms.OLS(df['price'], df[['intercept','bedrooms']])
results = lm.fit()
results.summary()

lm = sms.OLS(df['price'], df[['intercept','bathrooms']])
results = lm.fit()
results.summary()
  • 分析結(jié)果
  1. 都具有統(tǒng)計顯著性。
  2. 如果把三個簡單線性回歸模型的決定系數(shù)加起來，值大于1。
    在多元線性回歸里，我們要如何找到 “正確” 的系數(shù)
    在簡單線性回歸那一節(jié)中，你知道了我們?yōu)槭裁匆钚』總€實際數(shù)據(jù)點與模型預(yù)測值之間的平方距離。
  • 用所有變量建立多元線性回歸模型
lm = sms.OLS(df['price'], df[['intercept','area','bedrooms','bathrooms']])
results = lm.fit()
results.summary()
  • 再次分析結(jié)果
  1. 臥室數(shù)量：無統(tǒng)計顯著性。
  2. 浴室數(shù)量：無統(tǒng)計顯著性。
  3. 面積：具有統(tǒng)計顯著性。
  4. 決定系數(shù)小于1。
  5. 不能擬合像“風(fēng)格”之類的分類變量。

另一種求解方式

但在多元線性回歸中，我們要找的數(shù)據(jù)點所處的空間實際上不僅僅是二維的。
可以以線性代數(shù)的方式寫出求解公式：
\bar{β} = (X'X)-X'y

我們將使用 statsmodel 來求得系數(shù)，這些系數(shù)會與上一概念涉及的系數(shù)相似，但你還要用上述等式來解出該系數(shù)，證明所得系數(shù)不是憑空捏造出來的。

X = df[['intercept', 'bathrooms', 'bedrooms', 'area']]
y = df['price']

np.dot(np.dot(np.linalg(np.dot(X.transpose(), X)), X.transpose()), y)

numpy.matrix.transpose：矩陣轉(zhuǎn)置
np.linalg: 矩陣求逆
np.dot：矩陣相乘

啟示

  1. 我們可以將各系數(shù)解釋為： 在模型其它變量不變的情況下，解釋變量每增加一個單位，反應(yīng)變量會隨之增加的預(yù)測幅度。
  2. 當(dāng)分別單獨處理的時候，有些解釋變量是統(tǒng)計顯著的，而當(dāng)將他們一起擬合的時候，一些解釋變量反而不是統(tǒng)計顯著的了。
  3. 根據(jù)我們的現(xiàn)有知識，無法處理分類變量。

虛擬變量

要往線性模型里添加分類變量，就需要把分類變量轉(zhuǎn)變?yōu)?虛擬變量。即為每一個值創(chuàng)建一列：

vv.png

轉(zhuǎn)化后，你需要舍棄一個 虛擬列，才能得到 滿秩 矩陣【1】。

回想一下回歸系數(shù)閉式解那節(jié)內(nèi)容，我們得用 (X'X)-X'y來預(yù)測 β 值。

為了逆轉(zhuǎn)(X'X)，矩陣X 一定要是滿秩的，也就是說，所有X 的列都必須線性獨立。

創(chuàng)建虛擬變量時，如果你不舍棄掉一列，那你就得不到穩(wěn)定解，從 python 里得到的結(jié)果也不會可靠到哪去。

本文重點在于 如果你要用 0 、1 編碼來創(chuàng)建虛擬變量，你就得舍棄一個虛擬列，確保所得矩陣是滿秩的（這樣你從 python 里得到的解才會是可靠的。）

之所以要這么做，原因就在于線性代數(shù)的本質(zhì)，更具體地說，要逆轉(zhuǎn)矩陣，你手里的矩陣必須是滿秩的 （也就是所有列都得線性獨立），因此，你得舍棄掉一個虛擬列，方能得到線性獨立的各列 （和一個滿秩矩陣）。

用python演示創(chuàng)建虛擬變量

pd_dummies = pd.get_dummies(df['neighborhood'])
df_new = df.join(pd_dummies)
lm = sms.OLS(df_new['price'], df_new[['intercept','A']])
results = lm.fit()
results.summary()

如果刪除的虛擬變量為A，剩下的為‘B’和‘C‘，則可以解釋為當(dāng)neighborhood為“A”時，預(yù)測價格為intercept。B和C的斜率解釋為和A相比，B和C對反應(yīng)變量的變化量的影響。
還可以對數(shù)據(jù)進(jìn)行可視化以直觀的觀察分類變量對反應(yīng)變量的影響：

plt.hist(df_new.query("C == 1")['price'], alpha = 0.3, label = 'C');
plt.hist(df_new.query("A == 1")['price'], alpha = 0.3, label = 'A');
plt.hist(df_new.query("B == 1")['price'], alpha = 0.3, label = 'B');
plt.legend();

其他引入虛擬變量的方法

要用 1、0 編碼預(yù)測基準(zhǔn)類別，你要用到截距，此時斜率是與基準(zhǔn)分類比較下的變化；在SAS中，用 1、0、-1 編碼，你則得讓各分類系數(shù)分別乘以 -1，從而填補(bǔ)缺失的分類。

因為如果以 1、0、-1 編碼，對于類別為基準(zhǔn)類別的數(shù)據(jù)行來說，所有的虛擬變量都為-1，因此如果假設(shè)類別列表為A，B，C，且基準(zhǔn)類別為C：
C的平均值 = intercept - B的平均值 - A的平均值。

  • 導(dǎo)入數(shù)據(jù)和庫
import numpy as np
import pandas as pd
import statsmodels.api as sm

df = pd.read_csv('./house_prices.csv')
df2 = df.copy()
df.head(10)
  • 定義創(chuàng)建1，0，-1編碼的虛擬變量。
## The below function creates 1, 0, -1 coded dummy variables.

def dummy_cat(df, col):
    '''
    INPUT:
    df - the dataframe where col is stored
    col - the categorical column you want to dummy (as a string)
    OUTPUT:
    df - the dataframe with the added columns
         for dummy variables using 1, 0, -1 coding
    '''
    for idx, val_0 in enumerate(df[col].unique()):
        if idx + 1 < df[col].nunique():            
            df[val_0] = df[col].apply(lambda x: 1 if x == val_0 else 0)
        else:    
            df[val_0] = df[col].apply(lambda x: -1 if x == val_0 else 0)
            for idx, val_1 in enumerate(df[col].unique()):
                if idx + 1 < df[col].nunique():
                    df[val_1] = df[val_0] + df[val_1]
                else:
                    del df[val_1]
    return df
  • 為分類變量調(diào)用方法
new_df = dummy_cat(df, 'style') # Use on style
new_df.head(10)
  • 線性回歸
new_df['intercept'] = 1

lm = sm.OLS(new_df['price'], new_df[['intercept', 'ranch', 'victorian']])
results = lm.fit()
results.summary()

置信區(qū)間不重疊->顯著不同？

模型假設(shè)及相應(yīng)的解決方法

之前的視頻曾提到 統(tǒng)計學(xué)習(xí)簡介 一書，其中就提到了如下五個假設(shè)：

  1. 因變量-自變量關(guān)系的非線性
  2. 誤差項的相關(guān)性
  3. 非恒定方差和正態(tài)分布誤差
  4. 異常值/高杠桿點
  5. 共線性

這里總結(jié)了判斷上述問題是否存在的方法以及相應(yīng)的解決辦法。這是統(tǒng)計學(xué)家面試時經(jīng)常提出的問題，但該問題是否有實際意義取決于你創(chuàng)建模型的目的。在接下來的概念中，我們會更仔細(xì)地研究某些相關(guān)知識點。下方列出了各知識點的詳盡介紹，我們先來仔細(xì)看看下文涉及的每一項。

線性

線性是假設(shè)因變量和自變量之間真的存在可用線性模型解釋的關(guān)系。如果線性假設(shè)不為真，那你的預(yù)測結(jié)果就不會很準(zhǔn)確，此外，與系數(shù)有關(guān)的線性關(guān)系也就沒什么用了。

為了評估某段線性關(guān)系是否合理，一個很實用的方法是做預(yù)測值 \hat{y} 的殘差(y - \hat{y})圖。如果圖中出現(xiàn)多個曲線部分，那就意味著線性模型實際上可能并不擬合數(shù)據(jù)，自變量和因變量存在其它關(guān)系。創(chuàng)建非線性模型的辦法有很多（甚至可以線性模型的形式來創(chuàng)建）。

在本頁底部的圖片里，這些稱為 偏差 模型。理想來說，我們想要的是像圖片左上角殘差圖那樣的隨機(jī)散點圖。

相關(guān)誤差

如果我們是隨時間變化來收集的數(shù)據(jù)（比如預(yù)測未來股價或利率），或數(shù)據(jù)與空間有關(guān)（如預(yù)測洪澇或干旱地區(qū)），那就很容易出現(xiàn)相關(guān)誤差。通常，我們可以用過去數(shù)據(jù)點提供的信息（針對與時間有關(guān)的數(shù)據(jù)）或用相鄰數(shù)據(jù)點提供的信息（針對與空間有關(guān)的數(shù)據(jù)）來提高預(yù)測結(jié)果。

不考慮相關(guān)誤差的主要問題在于：往往你會利用這一相關(guān)性，得到更好的未來事件預(yù)測數(shù)據(jù)或空間關(guān)聯(lián)事件預(yù)測數(shù)據(jù)。

要判斷是否有相關(guān)誤差，最常用的方法是觀察收集數(shù)據(jù)的域。要是你不確定的話，你可以試試一個叫 Durbin-Watson 的檢驗方法，人們常用該測試來評估誤差相關(guān)性是否造成問題。還有 ARIMA 或 ARMA 模型，人們常用這兩個模型來利用誤差相關(guān)性，以便做出更佳預(yù)測。

非恒定方差和正態(tài)分布誤差

你預(yù)測的值不同，得到的預(yù)測值范圍也不同，那就意味著方差不恒定。非恒定方差對預(yù)測好壞影響不大，但會導(dǎo)致置信區(qū)間和 p 值不準(zhǔn)確，這種時候，在預(yù)測值接近實際值的那部分區(qū)域，系數(shù)的置信區(qū)間會太泛，而在預(yù)測值較遠(yuǎn)離實際值的區(qū)域則會太窄。

通常來說，對數(shù)函數(shù)（或使用其它反應(yīng)變量的變換方式）能夠 “擺脫” 非恒定方差，而要選擇合適的變換方式，我們一般會用 Box-Cox。

用預(yù)測值的殘差圖也可以評估非恒定方差。在本頁底部的圖片中，非恒定方差的標(biāo)簽為 異方差。理想來說，我們要的是一個有異方差殘差的無偏模型（其異方差殘差在一定數(shù)值范圍內(nèi)保持不變）。

雖然本文并不探討殘差的正態(tài)性，如果你想創(chuàng)建可靠的置信區(qū)間，正態(tài)性回歸假設(shè)就十分重要了，更多相關(guān)信息詳見 這里。

異常值/杠桿點

異常值和杠桿點是遠(yuǎn)離數(shù)據(jù)正常趨勢的點。這些點會對你的解造成很大的影響，在現(xiàn)實中，這些點甚至可能是錯誤的。如果從不同來源收集數(shù)據(jù)，你就可能在記錄或收集過程中造成某些數(shù)據(jù)值出錯。

異常值也可能是準(zhǔn)確真實的數(shù)據(jù)點，而不一定是測量或數(shù)據(jù)輸入錯誤。在這種情況下，'修復(fù)'就會變得更為主觀。要如何處理這些異常值往往取決于你的分析目的。線性模型，特別是使用最小二乘法的線性模型，比較容易受到影響，也就是說，大異常值可能會大幅度地左右我們的結(jié)果。當(dāng)然，異常值也有一些解決技巧，也就是我們常說的 正則化。本課不會談及這些技巧，但在 機(jī)器學(xué)習(xí)納米學(xué)位免費課程中，我們對這些技巧做了粗略的介紹。

而在賓夕法尼亞州立大學(xué)提供的完整回歸課程里，就有特別長的篇幅在探討杠桿點的問題，詳見 這里。

共線性（多重共線性）

如果我們的自變量彼此相關(guān)，就會出現(xiàn)多重共線性。多重共線性的一個主要問題在于：它會導(dǎo)致簡單線性回歸系數(shù)偏離我們想要的方向。

要判斷是否有多重共線性，最常見的辦法是借助二變量圖或 方差膨脹因子 (即 VIFs)。

resid-plots.gif

多重共線性與VIF

多元線性回歸模型的重要假設(shè)之一是解釋變量與反應(yīng)變量相關(guān)，而不是彼此相關(guān)。例如前面的例子，bedrooms、bathrooms和房子大小是相關(guān)的。

可以用可視化方法查看解釋變量彼此之間的相關(guān)性：

sb.pairplot(df[['area','bedrooms','bathrooms']]);

除了作圖之外，還可以通過方差膨脹因子VIFs 檢驗解釋變量之間的相關(guān)性。

如果一個屬性的VIF過大，那么可以刪除和這個屬性corr值較大的屬性。

一般情況下，只需了解到多重共線性的原理，以及解決策略即可。其他內(nèi)容可以作為選修內(nèi)容。

計算VIF

y, X = dmatrices('price ~ bedrooms + bathrooms + area' , df, return_type='dataframe')

vif = pd.DataFrame()
vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif["features"] = X.columns

寫成函數(shù)：

def vif_calculator(df, response):
    '''
    INPUT:
    df - a dataframe holding the x and y-variables
    response - the column name of the response as a string
    OUTPUT:
    vif - a dataframe of the vifs
    '''
    df2 = df.drop(response, axis = 1, inplace=False)
    features = "+".join(df2.columns)
    y, X = dmatrices(response + ' ~' + features, df, return_type='dataframe')
    vif = pd.DataFrame()
    vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    vif["features"] = X.columns
    vif = vif.round(1)
    return vif

衡量VIF

一般規(guī)則是，如果VIF大于10，那么模式就會出現(xiàn)共線性。VIFi = 1 / (1-R_i2)

這里R_i2是通過擬合以下線性模型得到的：
xi = b0 + b1x1 + b2x2 + ... + bnxn

也就是說如果R_i2很大，即x_i和其他變量的相關(guān)性具有統(tǒng)計顯著性，那么VIF的值就會很大。

通常，不會發(fā)現(xiàn)只有某一個變量的VIF較大，因為如果VIF較大，意味著有兩個或以上的變量具有相關(guān)性，這些變量都會有較高的VIF。一般做法是刪除最小興趣的變量。

高階項

有時候可能需要擬合類似于非線性模型的模型?？梢酝ㄟ^在模型中增加高階項，如交叉項（x1* x2），二次式，三次式和更高階的項。雖然這些高階項能讓我們較好的預(yù)測反應(yīng)變量，但是應(yīng)該盡量避免使用他們，因為：

  1. 高階項增大了模型的復(fù)雜度。
  2. 使得模型更加難于解釋。（解釋變量單位的平方往往沒有現(xiàn)實意義）

python 演示

# 增加二次式
# 注意，一定也要把各自的低階項加入模型
df['bedrooms_squared'] = df['bedrooms'] * df['bed_rooms']
lm = sm.OLS(new_df['price'], new_df[['intercept', 'bedrooms', 'bedrooms_squared']])
results = lm.fit()
results.summary()

此時問題變得復(fù)雜：如何解釋”臥室數(shù)量的平方“？但是決定系數(shù)并沒有明顯的變化！

要怎么確定高階項？

在創(chuàng)建具備二階、三階甚至更高階變量的模型時，基本上我們要看的是解釋變量和反應(yīng)變量的關(guān)系里有多少條曲線。

如果像下圖一樣，圖中曲線折了一個彎，那你會想添加一個二階項到模型里，因為我們可以明顯地發(fā)現(xiàn)，只用一條線來擬合這個關(guān)系并不能得到最佳效果。


quadraticlinearregression.png

當(dāng)解釋變量和反應(yīng)變量的關(guān)系折了兩個彎，我們就會想添加一個三階的關(guān)系，如下圖所示。


resid2.jpg

詳細(xì)的討論可以參考： https://tamino.wordpress.com/2011/03/31/so-what

解釋交叉項

通過上一節(jié)的討論，我們了解了如何通過觀察散點圖以確定是否采用二次式或者三次式來擬合模型。本節(jié)介紹如何確定應(yīng)該添加交叉項。

當(dāng)添加交叉項時，我們會認(rèn)為變量X1與Y之間的相關(guān)性，依賴于另一個變量X2。

  • 沒有交叉項的情況


    MR.png

這里我們可以看到，對于街區(qū)A和街區(qū)B：

  1. 得到的直線都能很好的擬合這些點。
  2. 得到的直線的斜率是一樣的。也就是說，二者是平行的，街區(qū)和面積不具有相關(guān)性。

  • 有交叉項的情況


    MR2.png

這里我們可以看到，對于街區(qū)A和街區(qū)B：

  1. 得到的直線都能很好的擬合這些點。
  2. 得到的直線的斜率不同。也就是說，價格和面積的關(guān)系取決于房子位于哪個社區(qū)，即變量Area與Price之間的相關(guān)性，依賴于另一個變量Neighborhood。

特征工程與特征選擇

當(dāng)我們進(jìn)行線性回歸或者機(jī)器學(xué)習(xí)時，模型能給我們帶來多大收獲取決于特征工程和特征選擇。特征選擇是指選擇應(yīng)該在模型中使用的解釋變量。
本文已經(jīng)介紹了一些選擇變量的方法，如顯著性水平（p值）和VIF，另一類方法是交叉驗證或正則化。
特征工程包括使用缺失值或刪除行，或者用初始特征的不同比例創(chuàng)建新的一行（稱為scaling：log，exp，standardizing，quadratics），將文本和圖片轉(zhuǎn)化為數(shù)字，或者創(chuàng)建虛擬變量。
用初始特征的不同比例創(chuàng)建新的一行有助于對變量進(jìn)行預(yù)測，但是可能會使得回歸系數(shù)無法解釋，應(yīng)該慎用。

Sebastian 和 Katie 講授的 這個 課程是個很好的學(xué)習(xí)途徑，能幫我們更好地理解本文涉及的許多概念，也對我們接下來更深入地理解機(jī)器學(xué)習(xí)很有幫助。

說到這里，我們已經(jīng)了解了特征工程的多種方式，如虛擬變量和缺失值填充。特別的，如果需要scaling數(shù)據(jù)，有兩個非常著名的庫：

  1. Sklearn，若要查看有關(guān)特征工程的 Sklearn 文獻(xiàn)，請點 這里。從scikit-learn中導(dǎo)入pre-processing，可以看到scale方法。然后還可以看到一個叫做“StandardScaler”的東東，他們可以處理不怎么像標(biāo)準(zhǔn)正態(tài)分布的數(shù)據(jù)集，將他們減去mean再除以標(biāo)準(zhǔn)差。此外還有許多其他用途的scaler，MinMax scaler，MaxAbs Scaler，以及其他處理數(shù)據(jù)的方式quantile transformer，normalization，binarization，encoding categorical features，imputing missing values，custom transformations。
    注意到這些技巧大部分都包括擬合轉(zhuǎn)換，即先擬合數(shù)據(jù)然后按照一定的方式對數(shù)據(jù)進(jìn)行轉(zhuǎn)換，這是scikit-learn函數(shù)的常見用法。具體來說，preprocessing常依序調(diào)用fit、transform對數(shù)據(jù)進(jìn)行擬合和轉(zhuǎn)化，或者直接調(diào)用fit_transform方法。這樣做的目的是使得模型更加健壯，而不會收到scaling的影響?；貧w是一種非常容易受到離散點影響的一種技巧。這樣做可以提高預(yù)測的準(zhǔn)確度，但是確實會改變解釋，
import numpy as np
import pandas as pd
import statsmodels.api as sm;
import sklearn.preprocessing as p

df = pd.DataFrame({'response': [2.4, 3.3, -4.2, 5.6, 1.5, 8.7], 
                         'x1': ['yes','no','yes','maybe','no','yes'],
                         'x2': [-1,-3,np.nan, 0, np.nan, 1],
                         'x3': [2.4, 15, 3.3, 2.4, 1.8, 0.4],
                         'x4': [np.nan, np.nan, 1, 1, 1, 1],
                         'x5': ['A', 'B', np.nan, 'A', 'A', 'A']})

p.scale(df[['x2','x3','x4']])

min_max_scaler = p.MinMaxScaler()
min_max_scaler .fit_transform(df[['x2','x3','x4']])

# 擬合失敗
df['intercept'] = 1
lm = sm.OLS(df['response'], df[['intercept','x2','x3','x4']])
results = lm.fit()
results.summary()

# 另一種常用的縮放特征的方法是減去均值并除以標(biāo)準(zhǔn)偏差
norm = p.StandardScaler()
norm.fit(df[['x2','x3','x4']])
norm.transform(df[['x2','x3','x4']])

sklearn.preprocessing包提供了一些通用的工具方法和transformer類，用于將原始特征向量轉(zhuǎn)換成易于評估的形式。

  1. pandas，有關(guān)在 pandas 里處理缺失值的其它文獻(xiàn)，請點 這里。
df.fillna(df.mean(), inplace = True)
lm = sm.OLS(df['response'], df[['intercept','x2','x3','x4']])
results = lm.fit()
results.summary()

模型評估的測量標(biāo)準(zhǔn)

模型評估的測量標(biāo)準(zhǔn)有很多：R-squared，MSE，AIC，BIC和馬洛斯的CP值等等。這些標(biāo)準(zhǔn)的值取決于數(shù)據(jù)的縮放比例和用于預(yù)測反應(yīng)變量的變量。有些時候，人們傾向于在實踐中采用R-squared和MSE，但是實際上這些方法在模型擬合中存在誤導(dǎo)性，如果我們只在同一個數(shù)據(jù)集上進(jìn)行驗證，無論我們像模型中添加什么變量，都會改善模型，也就是說R-suqared值會上升，MSE會下降。那么，我們怎么知道增加變量是否能真的改善了模型與數(shù)據(jù)的擬合度呢？有一個非常有效的技巧，稱為交叉驗證。

交叉驗證（模型評估）

在交叉驗證中，我們對數(shù)據(jù)子集訓(xùn)練回歸模型，稱為訓(xùn)練數(shù)據(jù)。我們用稱為測試數(shù)據(jù)的數(shù)據(jù)子集測試模型如何較好的運行。在交叉驗證中，我們通過多次切分?jǐn)?shù)據(jù)集，以確保我們的模型能夠通用化（generalize）。

k 折交叉驗證由 Sebastian 的視頻講解，屬于優(yōu)達(dá)學(xué)城的免費課程，很好地解釋了交叉驗證的工作原理。

在交叉驗證中，我們希望訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)的兩個集合都做到最大化。訓(xùn)練數(shù)據(jù)集越大，學(xué)習(xí)效果就越好；測試數(shù)據(jù)集越大，驗證效果越好。顯然這里需要折衷。
折衷的要點是將訓(xùn)練數(shù)據(jù)平分到相同大小的k個容器中，比方說有200個數(shù)據(jù)，如果k=10，那么每個容器中有20條數(shù)據(jù)。在k階交叉驗證中，將運行k次單獨的學(xué)習(xí)實驗，在每次實驗中，從k個子集中挑選一個作為驗證集，剩下k-1容器放在一起作為訓(xùn)練集，并在驗證集上驗證性能。交叉驗證的要點是這個操作會運行多次。在這個例子中，k=10，因此運行10次，然后對十個不同的運行結(jié)果的表現(xiàn)進(jìn)行平均，也就是說將這k次實驗的測試結(jié)果取平均值。這樣，學(xué)習(xí)算法的評估將更加準(zhǔn)確，而且可以對全部數(shù)據(jù)進(jìn)行驗證。

可以在這里找到參考文檔。

通過python演示數(shù)據(jù)集的切分

  • 引入庫和數(shù)據(jù)
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
df = pd.read_csv('./data/house_prices.csv')
  • 切分?jǐn)?shù)據(jù)集
    這一步必須在特征工程之前做。因為如果你在切分前得到數(shù)據(jù)的平均數(shù)，那么這就意味這測試集的數(shù)據(jù)包含了訓(xùn)練集的數(shù)據(jù)，而這不是我們期待的結(jié)果。
X = df[['neighborhood','area','bedroom','bathroom']]
y = df[['price']]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)
  • 初始化線性回歸模型
lm = LinearRegression()
  • 對訓(xùn)練數(shù)據(jù)集進(jìn)行擬合
lm.fit(X_train, y_train)
  • 對測試數(shù)據(jù)集進(jìn)行預(yù)測
y_preds = lm.predict(X_tests)
  • 對預(yù)測結(jié)果進(jìn)行打分
    應(yīng)該用測試數(shù)據(jù)集而不是訓(xùn)練數(shù)據(jù)集計算分?jǐn)?shù)，因為測試數(shù)據(jù)集體現(xiàn)了模型在沒有遇到過的數(shù)據(jù)集上的表現(xiàn)。
    根據(jù)分?jǐn)?shù)，可以選擇在模型中應(yīng)該保留哪些屬性。
# 默認(rèn)的打分方法為R-squared
lm.score(X_test, y_test)

案例研究

-導(dǎo)入數(shù)據(jù)集和庫

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from patsy import dmatrices
import matplotlib.pyplot as plt
%matplotlib inline

np.random.seed(42)
boston_data = load_boston()
df = pd.DataFrame()
df['MedianHomePrice'] = boston_data.target
df2 = pd.DataFrame(boston_data.data)
df2.columns = boston_data.feature_names
df = df.join(df2)
df.head()
  • 準(zhǔn)備X，Y
y = df['MedianHomePrice']
X = df.drop(['MedianHomePrice'], axis=1)
X.head(10)
  • split 數(shù)據(jù)
X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.20, random_state=0)
X_train.head()
  • 獲取數(shù)據(jù)集中每個特征的匯總
X_train.isnull().sum()
X_train.info()
X_train['AGE'].median()
X_train['CHAS'].mean()

一般情況下，可以查看：

  1. 數(shù)據(jù)集的缺失值數(shù)
  2. 數(shù)據(jù)集的分類變量數(shù)
  3. 數(shù)據(jù)集某個變量的中位數(shù)
    。。。
  • 使用 corr方法對每個變量進(jìn)行相互比較
df2 = pd.DataFrame()
df2['MedianHomePrice'] = y_train
df2[list(df.columns)[1:]] = X_train
df2.corr()

通過corr的結(jié)果可以知道兩方面的信息：

  1. 和反應(yīng)變量相關(guān)性較強(qiáng)的解釋變量。
  2. 解釋變量之間是否有較強(qiáng)的相關(guān)性。

corr的使用可以參照這里。

  • 縮放數(shù)據(jù)集中的所有 x 變量
# 使用StandardScaler來縮放數(shù)據(jù)集中的所有 x 變量
norm = StandardScaler()
norm.fit(X_train)
# 將結(jié)果存儲在 X_scaled_train中
X_scaled_train = norm.transform(X_train)

StandardScaler的使用方法參照這里。

  • 創(chuàng)建一個 pandas 數(shù)據(jù)幀存儲縮放的 x 變量以及training響應(yīng)變量
X_scaled_train = pd.DataFrame(X_scaled_train)
X_scaled_train.columns = list(df.columns)[1:]
df3 = X_scaled_train.join(y_train.reset_index())
  • 用所有的縮放特征來擬合線性模型，以預(yù)測此響應(yīng)（平均房價)
df3['intercept'] = 1
X_vars_full = df3.drop(['MedianHomePrice','index'] , axis=1, inplace=False)
lm = sm.OLS(df3['MedianHomePrice'], X_vars_full)
results = lm.fit()
results.summary()

觀察由線性模型得出的p值，標(biāo)注變量為具有統(tǒng)計顯著性或無顯著性。

  • 用于計算vif的函數(shù)
def vif_calculator(df, response):
    '''
    INPUT:
    df - a dataframe holding the x and y-variables
    response - the column name of the response as a string
    OUTPUT:
    vif - a dataframe of the vifs
    '''
    df2 = df.drop(response, axis = 1, inplace=False)
    features = "+".join(df2.columns)
    y, X = dmatrices(response + ' ~' + features, df, return_type='dataframe')
    vif = pd.DataFrame()
    vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    vif["features"] = X.columns
    vif = vif.round(1)
    return vif
  • 計算vif
vif_calculator(df3, 'MedianHomePrice')

可以根據(jù)經(jīng)驗設(shè)置VIF值的閾值，刪除VIF高于該閾值的變量。如果需要的或者在乎的變量的VIF較高，那么就需要考慮刪除和VIF較高的變量相關(guān)性較強(qiáng)的變量。

  • 為模型打分
X_data = df.drop('MedianHomePrice', axis=1, inplace=False)
X_train, X_test, y_train, y_test = train_test_split(
                 X_data, df['MedianHomePrice'], test_size=0.2, random_state=0)


lm_full = LinearRegression()
lm_full.fit(X_train, y_train)
# 返回預(yù)測結(jié)果的R-sqaured值
lm_full.score(X_test, y_test)
  • 還可以對不同的模型進(jìn)行比較(假設(shè)RAD的p值較大)
X_train_red = X_train.drop(['AGE','NOX','TAX'] , axis=1, inplace=False)
X_test_red = X_test.drop(['AGE','NOX','TAX'] , axis=1, inplace=False)

X_train_red2 = X_train.drop(['AGE','NOX','TAX','RAD'] , axis=1, inplace=False)
X_test_red2 = X_test.drop(['AGE','NOX','TAX','RAD'] , axis=1, inplace=False)

lm_red = LinearRegression()
lm_red.fit(X_train_red, y_train)
print(lm_red.score(X_test_red, y_test))

lm_red2 = LinearRegression()
lm_red2.fit(X_train_red2, y_train)
print(lm_red2.score(X_test_red2, y_test))

LinearRegression的使用可以參照文檔。

決定系數(shù)高的模型可能只是因為變量較多，而不意味著這個模型真的能在新 （測試）數(shù)據(jù)集上有更好的效果。

參考書

  1. 多元線性回歸涉及很多概念和知識，統(tǒng)計學(xué)習(xí)簡介尤其是該書的第三章是很好的資源。
  2. 對于線性代數(shù)，可觀看可汗學(xué)院線性代數(shù)免費課程。

【1】滿秩矩陣是一個很重要的概念, 它是判斷一個矩陣是否可逆。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容