python實現邏輯回歸

本文基於yhat上Logistic Regression in Python，作了中文翻譯，並相應補充了一些內容。本文並不研究邏輯回歸具體算法實現，而是使用了一些算法庫，旨在幫助需要用Python來做邏輯回歸的訓練和預測的讀者快速上手。

邏輯回歸是一項可用於預測二分類結果(binary outcome)的統計技術，廣泛應用於金融、醫學、犯罪學和其他社會科學中。邏輯回歸使用簡單且非常有效，你可以在許多機器學習、應用統計的書中的前幾章中找到個關於邏輯回歸的介紹。邏輯回歸在許多統計課程中都會用到。

我們不難找到使用R語言的高質量的邏輯回歸實例，如UCLA的教程R Data Analysis Examples: Logit Regression就是一個很好的資源。Python是機器學習領域最流行的語言之一，並且已有許多Python的資源涵蓋了支持向量積和文本分類等話題，但少有關於邏輯回歸的資料。

本文介紹了如何使用Python來完成邏輯回歸。

 

簡介

示例代碼中使用了一些算法包，請確保在運行這些代碼前，你的電腦已經安裝了如下包:

• numpy: Python的語言擴展，定義了數字的數組和矩陣
• pandas: 直接處理和操作數據的主要package
• statsmodels: 統計和計量經濟學的package，包含了用於參數評估和統計測試的實用工具
• pylab: 用於生成統計圖

可參考 Windows安裝Python機器學習包 或 Ubuntu/CentOS安裝Python機器學習包 來搭建所需要的環境。

邏輯回歸的實例

在此使用與Logit Regression in R相同的數據集來研究Python中的邏輯回歸，目的是要辨別不同的因素對研究生錄取的影響。

數據集中的前三列可作為預測變量(predictor variables)：

• gpa
• gre分數
• rank表示本科生母校的聲望

第四列admit則是二分類目標變量(binary target variable)，它表明考生最終是否被錄用。

加載數據

使用 pandas.read_csv加載數據，這樣我們就有了可用於探索數據的DataFrame。

1. import pandas as pd
2. import statsmodels.api as sm
3. import pylab as pl
4. import numpy as np
5.  
6. # 加載數據
7. # 備用地址: http://cdn.powerxing.com/files/lr-binary.csv
8. df = pd.read_csv("http://www.ats.ucla.edu/stat/data/binary.csv")
9.  
10. # 瀏覽數據集
11. print df.head()
12. # admit gre gpa rank
13. # 0 0 380 3.61 3
14. # 1 1 660 3.67 3
15. # 2 1 800 4.00 1
16. # 3 1 640 3.19 4
17. # 4 0 520 2.93 4
18.  
19. # 重命名'rank'列，因為dataframe中有個方法名也為'rank'
20. df.columns = ["admit", "gre", "gpa", "prestige"]
21. print df.columns
22. # array([admit, gre, gpa, prestige], dtype=object)

Python

注意到有一列屬性名為rank，但因為rank也是pandas dataframe中一個方法的名字，因此需要將該列重命名為”prestige”.

統計摘要(Summary Statistics) 以及 查看數據

現在我們就將需要的數據正確載入到Python中了，現在來看下數據。我們可以使用pandas的函數describe來給出數據的摘要–describe與R語言中的summay類似。這里也有一個用於計算標准差的函數std，但在describe中已包括了計算標准差。

我特別喜歡pandas的pivot_table/crosstab聚合功能。crosstab可方便的實現多維頻率表(frequency tables)(有點像R語言中的table)。你可以用它來查看不同數據所占的比例。

1. # summarize the data
2. print df.describe()
3. # admit gre gpa prestige
4. # count 400.000000 400.000000 400.000000 400.00000
5. # mean 0.317500 587.700000 3.389900 2.48500
6. # std 0.466087 115.516536 0.380567 0.94446
7. # min 0.000000 220.000000 2.260000 1.00000
8. # 25% 0.000000 520.000000 3.130000 2.00000
9. # 50% 0.000000 580.000000 3.395000 2.00000
10. # 75% 1.000000 660.000000 3.670000 3.00000
11. # max 1.000000 800.000000 4.000000 4.00000
12.  
13. # 查看每一列的標准差
14. print df.std()
15. # admit 0.466087
16. # gre 115.516536
17. # gpa 0.380567
18. # prestige 0.944460
19.  
20. # 頻率表，表示prestige與admin的值相應的數量關系
21. print pd.crosstab(df['admit'], df['prestige'], rownames=['admit'])
22. # prestige 1 2 3 4
23. # admit
24. # 0 28 97 93 55
25. # 1 33 54 28 12
26.  
27. # plot all of the columns
28. df.hist()
29. pl.show()

Python

運行代碼后，繪制的柱狀統計圖如下所示：

使用pylab繪制的柱狀統計圖

虛擬變量(dummy variables)

虛擬變量，也叫啞變量，可用來表示分類變量、非數量因素可能產生的影響。在計量經濟學模型，需要經常考慮屬性因素的影響。例如，職業、文化程度、季節等屬性因素往往很難直接度量它們的大小。只能給出它們的“Yes—D=1”或”No—D=0”，或者它們的程度或等級。為了反映屬性因素和提高模型的精度，必須將屬性因素“量化”。通過構造0-1型的人工變量來量化屬性因素。

pandas提供了一系列分類變量的控制。我們可以用get_dummies來將”prestige”一列虛擬化。

get_dummies為每個指定的列創建了新的帶二分類預測變量的DataFrame，在本例中，prestige有四個級別：1，2，3以及4（1代表最有聲望），prestige作為分類變量更加合適。當調用get_dummies時，會產生四列的dataframe，每一列表示四個級別中的一個。

1. # 將prestige設為虛擬變量
2. dummy_ranks = pd.get_dummies(df['prestige'], prefix='prestige')
3. print dummy_ranks.head()
4. # prestige_1 prestige_2 prestige_3 prestige_4
5. # 0 0 0 1 0
6. # 1 0 0 1 0
7. # 2 1 0 0 0
8. # 3 0 0 0 1
9. # 4 0 0 0 1
10.  
11. # 為邏輯回歸創建所需的data frame
12. # 除admit、gre、gpa外，加入了上面常見的虛擬變量（注意，引入的虛擬變量列數應為虛擬變量總列數減1，減去的1列作為基准）
13. cols_to_keep = ['admit', 'gre', 'gpa']
14. data = df[cols_to_keep].join(dummy_ranks.ix[:, 'prestige_2':])
15. print data.head()
16. # admit gre gpa prestige_2 prestige_3 prestige_4
17. # 0 0 380 3.61 0 1 0
18. # 1 1 660 3.67 0 1 0
19. # 2 1 800 4.00 0 0 0
20. # 3 1 640 3.19 0 0 1
21. # 4 0 520 2.93 0 0 1
22.  
23. # 需要自行添加邏輯回歸所需的intercept變量
24. data['intercept'] = 1.0

Python

這樣，數據原本的prestige屬性就被prestige_x代替了，例如原本的數值為2，則prestige_2為1，prestige_1、prestige_3、prestige_4都為0。

將新的虛擬變量加入到了原始的數據集中后，就不再需要原來的prestige列了。在此要強調一點，生成m個虛擬變量后，只要引入m-1個虛擬變量到數據集中，未引入的一個是作為基准對比的。

最后，還需加上常數intercept，statemodels實現的邏輯回歸需要顯式指定。

執行邏輯回歸

實際上完成邏輯回歸是相當簡單的，首先指定要預測變量的列，接着指定模型用於做預測的列，剩下的就由算法包去完成了。

本例中要預測的是admin列，使用到gre、gpa和虛擬變量prestige_2、prestige_3、prestige_4。prestige_1作為基准，所以排除掉，以防止多元共線性(multicollinearity)和引入分類變量的所有虛擬變量值所導致的陷阱(dummy variable trap)。

1. # 指定作為訓練變量的列，不含目標列`admit`
2. train_cols = data.columns[1:]
3. # Index([gre, gpa, prestige_2, prestige_3, prestige_4], dtype=object)
4.  
5. logit = sm.Logit(data['admit'], data[train_cols])
6.  
7. # 擬合模型
8. result = logit.fit()

Python

在這里是使用了statesmodels的Logit函數，更多的模型細節可以查閱statesmodels的文檔

使用訓練模型預測數據

（本小節是博主補充的）通過上述步驟，我們就得到了訓練后的模型。基於這個模型，我們就可以用來預測數據，代碼如下：

1. # 構建預測集
2. # 與訓練集相似，一般也是通過 pd.read_csv() 讀入
3. # 在這邊為方便，我們將訓練集拷貝一份作為預測集（不包括 admin 列）
4. import copy
5. combos = copy.deepcopy(data)
6.  
7. # 數據中的列要跟預測時用到的列一致
8. predict_cols = combos.columns[1:]
9.  
10. # 預測集也要添加intercept變量
11. combos['intercept'] = 1.0
12.  
13. # 進行預測，並將預測評分存入 predict 列中
14. combos['predict'] = result.predict(combos[predict_cols])
15.  
16. # 預測完成后，predict 的值是介於 [0, 1] 間的概率值
17. # 我們可以根據需要，提取預測結果
18. # 例如，假定 predict > 0.5，則表示會被錄取
19. # 在這邊我們檢驗一下上述選取結果的精確度
20. total = 0
21. hit = 0
22. for value in combos.values:
23. # 預測分數 predict, 是數據中的最后一列
24. predict = value[-1]
25. # 實際錄取結果
26. admit = int(value[0])
27.  
28. # 假定預測概率大於0.5則表示預測被錄取
29. if predict > 0.5:
30. total += 1
31. # 表示預測命中
32. if admit == 1:
33. hit += 1
34.  
35. # 輸出結果
36. print 'Total: %d, Hit: %d, Precision: %.2f' % (total, hit, 100.0*hit/total)
37. # Total: 49, Hit: 30, Precision: 61.22

Python

在這里，我是簡單的將原始數據再作為待預測的數據進行檢驗。通過上述步驟得到的是一個概率值，而不是一個直接的二分類結果（被錄取/不被錄取）。通常，我們可以設定一個閾值，若 predict 大於該閾值，則認為是被錄取了，反之，則表示不被錄取。

在上面的例子中，假定預測概率大於 0.5 則表示預測被錄取，一共預測有 49 個被錄取，其中有 30 個預測命中，精確度為 61.22%。

結果解釋

statesmodels提供了結果的摘要，如果你使用過R語言，你會發現結果的輸出與之相似。

1. # 查看數據的要點
2. print result.summary()

Python

1. Logit Regression Results
2. ==============================================================================
3. Dep. Variable: admit No. Observations: 400
4. Model: Logit Df Residuals: 394
5. Method: MLE Df Model: 5
6. Date: Sun, 03 Mar 2013 Pseudo R-squ.: 0.08292
7. Time: 12:34:59 Log-Likelihood: -229.26
8. converged: True LL-Null: -249.99
9. LLR p-value: 7.578e-08
10. ==============================================================================
11. coef std err z P>|z| [95.0% Conf. Int.]
12. ------------------------------------------------------------------------------
13. gre 0.0023 0.001 2.070 0.038 0.000 0.004
14. gpa 0.8040 0.332 2.423 0.015 0.154 1.454
15. prestige_2 -0.6754 0.316 -2.134 0.033 -1.296 -0.055
16. prestige_3 -1.3402 0.345 -3.881 0.000 -2.017 -0.663
17. prestige_4 -1.5515 0.418 -3.713 0.000 -2.370 -0.733
18. intercept -3.9900 1.140 -3.500 0.000 -6.224 -1.756
19. ==============================================================================

Python

你可以看到模型的系數，系數擬合的效果，以及總的擬合質量，以及一些統計度量。[待補充: 模型結果主要參數的含義]

當然你也可以只觀察結果的某部分，如置信區間(confidence interval)可以看出模型系數的健壯性。

1. # 查看每個系數的置信區間
2. print result.conf_int()
3. # 0 1
4. # gre 0.000120 0.004409
5. # gpa 0.153684 1.454391
6. # prestige_2 -1.295751 -0.055135
7. # prestige_3 -2.016992 -0.663416
8. # prestige_4 -2.370399 -0.732529
9. # intercept -6.224242 -1.755716

Python

在這個例子中，我們可以肯定被錄取的可能性與應試者畢業學校的聲望存在着逆相關的關系。

換句話說，高排名學校（prestige_1==True）的湘鄂生唄錄取的概率比低排名學校（prestige_4==True）要高。

相對危險度（odds ratio）

使用每個變量系數的指數來生成odds ratio，可知變量每單位的增加、減少對錄取幾率的影響。例如，如果學校的聲望為2，則我們可以期待被錄取的幾率減少大概50%。UCLA上有一個對odds ratio更為深入的解釋: 在邏輯回歸中如何解釋odds ratios?。

1. # 輸出 odds ratio
2. print np.exp(result.params)
3. # gre 1.002267
4. # gpa 2.234545
5. # prestige_2 0.508931
6. # prestige_3 0.261792
7. # prestige_4 0.211938
8. # intercept 0.018500

Python

我們也可以使用置信區間來計算系數的影響，來更好地估計一個變量影響錄取率的不確定性。

1. # odds ratios and 95% CI
2. params = result.params
3. conf = result.conf_int()
4. conf['OR'] = params
5. conf.columns = ['2.5%', '97.5%', 'OR']
6. print np.exp(conf)
7. # 2.5% 97.5% OR
8. # gre 1.000120 1.004418 1.002267
9. # gpa 1.166122 4.281877 2.234545
10. # prestige_2 0.273692 0.946358 0.508931
11. # prestige_3 0.133055 0.515089 0.261792
12. # prestige_4 0.093443 0.480692 0.211938
13. # intercept 0.001981 0.172783 0.018500

Python

更深入的挖掘

為了評估我們分類器的效果，我們將使用每個輸入值的邏輯組合（logical combination）來重新創建數據集，如此可以得知在不同的變量下預測錄取可能性的增加、減少。首先我們使用名為 cartesian 的輔助函數來生成組合值（來源於: 如何使用numpy構建兩個數組的組合）

我們使用 np.linspace 創建 “gre” 和 “gpa” 值的一個范圍，即從指定的最大、最小值來創建一個線性間隔的值的范圍。在本例子中，取已知的最大、最小值。

1. # 根據最大、最小值生成 GRE、GPA 均勻分布的10個值，而不是生成所有可能的值
2. gres = np.linspace(data['gre'].min(), data['gre'].max(), 10)
3. print gres
4. # array([ 220. , 284.44444444, 348.88888889, 413.33333333,
5. # 477.77777778, 542.22222222, 606.66666667, 671.11111111,
6. # 735.55555556, 800. ])
7. gpas = np.linspace(data['gpa'].min(), data['gpa'].max(), 10)
8. print gpas
9. # array([ 2.26 , 2.45333333, 2.64666667, 2.84 , 3.03333333,
10. # 3.22666667, 3.42 , 3.61333333, 3.80666667, 4. ])
11.  
12.  
13. # 枚舉所有的可能性
14. combos = pd.DataFrame(cartesian([gres, gpas, [1, 2, 3, 4], [1.]]))
15. # 重新創建啞變量
16. combos.columns = ['gre', 'gpa', 'prestige', 'intercept']
17. dummy_ranks = pd.get_dummies(combos['prestige'], prefix='prestige')
18. dummy_ranks.columns = ['prestige_1', 'prestige_2', 'prestige_3', 'prestige_4']
19.  
20. # 只保留用於預測的列
21. cols_to_keep = ['gre', 'gpa', 'prestige', 'intercept']
22. combos = combos[cols_to_keep].join(dummy_ranks.ix[:, 'prestige_2':])
23.  
24. # 使用枚舉的數據集來做預測
25. combos['admit_pred'] = result.predict(combos[train_cols])
26.  
27. print combos.head()
28. # gre gpa prestige intercept prestige_2 prestige_3 prestige_4 admit_pred
29. # 0 220 2.260000 1 1 0 0 0 0.157801
30. # 1 220 2.260000 2 1 1 0 0 0.087056
31. # 2 220 2.260000 3 1 0 1 0 0.046758
32. # 3 220 2.260000 4 1 0 0 1 0.038194
33. # 4 220 2.453333 1 1 0 0 0 0.179574

Python

現在我們已生成了預測結果，接着通過畫圖來呈現結果。我編寫了一個名為 isolate_and_plot 的輔助函數，可以比較給定的變量與不同的聲望等級、組合的平均可能性。為了分離聲望和其他變量，我使用了 pivot_table 來簡單地聚合數據。

1. def isolate_and_plot(variable):
2. # isolate gre and class rank
3. grouped = pd.pivot_table(combos, values=['admit_pred'], index=[variable, 'prestige'],
4. aggfunc=np.mean)
5.  
6. # in case you're curious as to what this looks like
7. # print grouped.head()
8. # admit_pred
9. # gre prestige
10. # 220.000000 1 0.282462
11. # 2 0.169987
12. # 3 0.096544
13. # 4 0.079859
14. # 284.444444 1 0.311718
15.  
16. # make a plot
17. colors = 'rbgyrbgy'
18. for col in combos.prestige.unique():
19. plt_data = grouped.ix[grouped.index.get_level_values(1)==col]
20. pl.plot(plt_data.index.get_level_values(0), plt_data['admit_pred'],
21. color=colors[int(col)])
22.  
23. pl.xlabel(variable)
24. pl.ylabel("P(admit=1)")
25. pl.legend(['1', '2', '3', '4'], loc='upper left', title='Prestige')
26. pl.title("Prob(admit=1) isolating " + variable + " and presitge")
27. pl.show()
28.  
29. isolate_and_plot('gre')
30. isolate_and_plot('gpa')

Python

結果圖顯示了 gre, gpa 和 prestige 如何影響錄取。可以看出，隨着 gre 和 gpa 的增加，錄取可能性如何逐漸增加，並且，不同的學校聲望對錄取可能性的增加程度相差很大。

結束語

邏輯回歸是用於分類的優秀算法，盡管有一些更加性感的，或是黑盒分類器算法，如SVM和隨機森林（RandomForest）在一些情況下性能更好，但深入了解你正在使用的模型是很有價值的。很多時候你可以使用隨機森林來篩選模型的特征，並基於篩選出的最佳的特征，使用邏輯回歸來重建模型。