2014-05-03 (updated: 2015-05-22) 3097
本文基于yhat上
Logistic Regression in Python，作了中文翻譯，并相應(yīng)補(bǔ)充了一些內(nèi)容。本文并不研究邏輯回歸具體算法實(shí)現(xiàn)，而是使用了一些算法庫，旨在幫助需要用Python來做邏輯回歸的訓(xùn)練和預(yù)測的讀者快速上手。
邏輯回歸是一項(xiàng)可用于預(yù)測二分類結(jié)果(binary outcome)的統(tǒng)計(jì)技術(shù)，廣泛應(yīng)用于金融、醫(yī)學(xué)、犯罪學(xué)和其他社會(huì)科學(xué)中。邏輯回歸使用簡單且非常有效，你可以在許多機(jī)器學(xué)習(xí)、應(yīng)用統(tǒng)計(jì)的書中的前幾章中找到個(gè)關(guān)于邏輯回歸的介紹。邏輯回歸在許多統(tǒng)計(jì)課程中都會(huì)用到。
我們不難找到使用R語言的高質(zhì)量的邏輯回歸實(shí)例，如UCLA的教程
R Data Analysis Examples: Logit Regression就是一個(gè)很好的資源。Python是機(jī)器學(xué)習(xí)領(lǐng)域最流行的語言之一，并且已有許多Python的資源涵蓋了
支持向量積和
文本分類等話題，但少有關(guān)于邏輯回歸的資料。
本文介紹了如何使用Python來完成邏輯回歸。
簡介
示例代碼中使用了一些算法包，請(qǐng)確保在運(yùn)行這些代碼前，你的電腦已經(jīng)安裝了如下包:
numpy: Python的語言擴(kuò)展，定義了數(shù)字的數(shù)組和矩陣
pandas: 直接處理和操作數(shù)據(jù)的主要package
statsmodels: 統(tǒng)計(jì)和計(jì)量經(jīng)濟(jì)學(xué)的package，包含了用于參數(shù)評(píng)估和統(tǒng)計(jì)測試的實(shí)用工具
pylab: 用于生成統(tǒng)計(jì)圖
可參考
Setting Up Scientific Python來安裝所需要的環(huán)境。
邏輯回歸的實(shí)例
在此使用與
Logit Regression in R相同的數(shù)據(jù)集來研究Python中的邏輯回歸，目的是要辨別不同的因素對(duì)研究生錄取的影響。
數(shù)據(jù)集中的前三列可作為預(yù)測變量(predictor variables)：
gpa
gre分?jǐn)?shù)
rank表示本科生母校的聲望
第四列admit則是二分類目標(biāo)變量(binary target variable)，它表明考生最終是否被錄用。
加載數(shù)據(jù)
使用 pandas.read_csv加載數(shù)據(jù)，這樣我們就有了可用于探索數(shù)據(jù)的DataFrame。
import pandas as pd
import statsmodels.api as sm
import pylab as pl
import numpy as np
# 加載數(shù)據(jù)
# 備用地址: http://cdn.powerxing.com/files/lr-binary.csv
df = pd.read_csv("http://www.ats.ucla.edu/stat/data/binary.csv")
# 瀏覽數(shù)據(jù)集
print df.head()
# admit gre gpa rank
# 0 0 380 3.61 3
# 1 1 660 3.67 3
# 2 1 800 4.00 1
# 3 1 640 3.19 4
# 4 0 520 2.93 4
# 重命名'rank'列，因?yàn)閐ataframe中有個(gè)方法名也為'rank'
df.columns = ["admit", "gre", "gpa", "prestige"]
print df.columns
# array([admit, gre, gpa, prestige], dtype=object)
注意到有一列屬性名為rank，但因?yàn)閞ank也是pandas dataframe中一個(gè)方法的名字，因此需要將該列重命名為”prestige”.
統(tǒng)計(jì)摘要(Summary Statistics) 以及 查看數(shù)據(jù)
現(xiàn)在我們就將需要的數(shù)據(jù)正確載入到Python中了，現(xiàn)在來看下數(shù)據(jù)。我們可以使用pandas的函數(shù)describe來給出數(shù)據(jù)的摘要–describe與R語言中的summay類似。這里也有一個(gè)用于計(jì)算標(biāo)準(zhǔn)差的函數(shù)std，但在describe中已包括了計(jì)算標(biāo)準(zhǔn)差。
我特別喜歡pandas的pivot_table/crosstab聚合功能。crosstab可方便的實(shí)現(xiàn)多維頻率表(frequency tables)(有點(diǎn)像R語言中的table)。你可以用它來查看不同數(shù)據(jù)所占的比例。
# summarize the data
print df.describe()
# admit gre gpa prestige
# count 400.000000 400.000000 400.000000 400.00000
# mean 0.317500 587.700000 3.389900 2.48500
# std 0.466087 115.516536 0.380567 0.94446
# min 0.000000 220.000000 2.260000 1.00000
# 25% 0.000000 520.000000 3.130000 2.00000
# 50% 0.000000 580.000000 3.395000 2.00000
# 75% 1.000000 660.000000 3.670000 3.00000
# max 1.000000 800.000000 4.000000 4.00000
# 查看每一列的標(biāo)準(zhǔn)差
print df.std()
# admit 0.466087
# gre 115.516536
# gpa 0.380567
# prestige 0.944460
# 頻率表，表示prestige與admin的值相應(yīng)的數(shù)量關(guān)系
print pd.crosstab(df['admit'], df['prestige'], rownames=['admit'])
# prestige 1 2 3 4
# admit
# 0 28 97 93 55
# 1 33 54 28 12
# plot all of the columns
df.hist()
pl.show()
運(yùn)行代碼后，繪制的柱狀統(tǒng)計(jì)圖如下所示：
使用pylab繪制的柱狀統(tǒng)計(jì)圖
虛擬變量(dummy variables)
虛擬變量，也叫啞變量，可用來表示分類變量、非數(shù)量因素可能產(chǎn)生的影響。在計(jì)量經(jīng)濟(jì)學(xué)模型，需要經(jīng)?？紤]屬性因素的影響。例如，職業(yè)、文化程度、季節(jié)等屬性因素往往很難直接度量它們的大小。只能給出它們的“Yes—D=1”或”No—D=0”，或者它們的程度或等級(jí)。為了反映屬性因素和提高模型的精度，必須將屬性因素“量化”。通過構(gòu)造0-1型的人工變量來量化屬性因素。
pandas提供了一系列分類變量的控制。我們可以用get_dummies來將”prestige”一列虛擬化。
get_dummies為每個(gè)指定的列創(chuàng)建了新的帶二分類預(yù)測變量的DataFrame，在本例中，prestige有四個(gè)級(jí)別：1，2，3以及4（1代表最有聲望），prestige作為分類變量更加合適。當(dāng)調(diào)用get_dummies時(shí)，會(huì)產(chǎn)生四列的dataframe，每一列表示四個(gè)級(jí)別中的一個(gè)。
# 將prestige設(shè)為虛擬變量
dummy_ranks = pd.get_dummies(df['prestige'], prefix='prestige')
print dummy_ranks.head()
# prestige_1 prestige_2 prestige_3 prestige_4
# 0 0 0 1 0
# 1 0 0 1 0
# 2 1 0 0 0
# 3 0 0 0 1
# 4 0 0 0 1
# 為邏輯回歸創(chuàng)建所需的data frame
# 除admit、gre、gpa外，加入了上面常見的虛擬變量（注意，引入的虛擬變量列數(shù)應(yīng)為虛擬變量總列數(shù)減1，減去的1列作為基準(zhǔn)）
cols_to_keep = ['admit', 'gre', 'gpa']
data = df[cols_to_keep].join(dummy_ranks.ix[:, 'prestige_2':])
print data.head()
# admit gre gpa prestige_2 prestige_3 prestige_4
# 0 0 380 3.61 0 1 0
# 1 1 660 3.67 0 1 0
# 2 1 800 4.00 0 0 0
# 3 1 640 3.19 0 0 1
# 4 0 520 2.93 0 0 1
# 需要自行添加邏輯回歸所需的intercept變量
data['intercept'] = 1.0
這樣，數(shù)據(jù)原本的prestige屬性就被prestige_x代替了，例如原本的數(shù)值為2，則prestige_2為1，prestige_1、prestige_3、prestige_4都為0。
將新的虛擬變量加入到了原始的數(shù)據(jù)集中后，就不再需要原來的prestige列了。在此要強(qiáng)調(diào)一點(diǎn)，生成m個(gè)虛擬變量后，只要引入m-1個(gè)虛擬變量到數(shù)據(jù)集中，未引入的一個(gè)是作為基準(zhǔn)對(duì)比的。
最后，還需加上常數(shù)intercept，statemodels實(shí)現(xiàn)的邏輯回歸需要顯式指定。
執(zhí)行邏輯回歸
實(shí)際上完成邏輯回歸是相當(dāng)簡單的，首先指定要預(yù)測變量的列，接著指定模型用于做預(yù)測的列，剩下的就由算法包去完成了。
本例中要預(yù)測的是admin列，使用到gre、gpa和虛擬變量prestige_2、prestige_3、prestige_4。prestige_1作為基準(zhǔn)，所以排除掉，以防止
多元共線性(multicollinearity)和引入分類變量的所有虛擬變量值所導(dǎo)致的陷阱(
dummy variable trap)。
# 指定作為訓(xùn)練變量的列，不含目標(biāo)列`admit`
train_cols = data.columns[1:]
# Index([gre, gpa, prestige_2, prestige_3, prestige_4], dtype=object)
logit = sm.Logit(data['admit'], data[train_cols])
# 擬合模型
result = logit.fit()
在這里是使用了statesmodels的
Logit函數(shù)，更多的模型細(xì)節(jié)可以查閱statesmodels的
文檔使用訓(xùn)練模型預(yù)測數(shù)據(jù)
（本小節(jié)是博主補(bǔ)充的）通過上述步驟，我們就得到了訓(xùn)練后的模型?；谶@個(gè)模型，我們就可以用來預(yù)測數(shù)據(jù)，代碼如下：
# 構(gòu)建預(yù)測集
# 與訓(xùn)練集相似，一般也是通過 pd.read_csv() 讀入
# 在這邊為方便，我們將訓(xùn)練集拷貝一份作為預(yù)測集（不包括 admin 列）
import copy
combos = copy.deepcopy(data)
# 數(shù)據(jù)中的列要跟預(yù)測時(shí)用到的列一致
predict_cols = combos.columns[1:]
# 預(yù)測集也要添加intercept變量
combos['intercept'] = 1.0
# 進(jìn)行預(yù)測，并將預(yù)測評(píng)分存入 predict 列中
combos['predict'] = result.predict(combos[predict_cols])
# 預(yù)測完成后，predict 的值是介于 [0, 1] 間的概率值
# 我們可以根據(jù)需要，提取預(yù)測結(jié)果
# 例如，假定 predict > 0.5，則表示會(huì)被錄取
# 在這邊我們檢驗(yàn)一下上述選取結(jié)果的精確度
total = 0
hit = 0
for value in combos.values:
# 預(yù)測分?jǐn)?shù) predict, 是數(shù)據(jù)中的最后一列
predict = value[-1]
# 實(shí)際錄取結(jié)果
admit = int(value[0])
# 假定預(yù)測概率大于0.5則表示預(yù)測被錄取
if predict > 0.5:
total += 1
# 表示預(yù)測命中
if admit == 1:
hit += 1
# 輸出結(jié)果
print 'Total: %d, Hit: %d, Precision: %.2f' % (total, hit, 100.0*hit/total)
# Total: 49, Hit: 30, Precision: 61.22
在這里，我是簡單的將原始數(shù)據(jù)再作為待預(yù)測的數(shù)據(jù)進(jìn)行檢驗(yàn)。通過上述步驟得到的是一個(gè)概率值，而不是一個(gè)直接的二分類結(jié)果（被錄取/不被錄?。?。通常，我們可以設(shè)定一個(gè)閾值，若 predict 大于該閾值，則認(rèn)為是被錄取了，反之，則表示不被錄取。
在上面的例子中，假定預(yù)測概率大于 0.5 則表示預(yù)測被錄取，一共預(yù)測有 49 個(gè)被錄取，其中有 30 個(gè)預(yù)測命中，精確度為 61.22%。
結(jié)果解釋
statesmodels提供了結(jié)果的摘要，如果你使用過R語言，你會(huì)發(fā)現(xiàn)結(jié)果的輸出與之相似。
# 查看數(shù)據(jù)的要點(diǎn)
print result.summary()
Logit Regression Results
==============================================================================
Dep. Variable: admit No. Observations: 400
Model: Logit Df Residuals: 394
Method: MLE Df Model: 5
Date: Sun, 03 Mar 2013 Pseudo R-squ.: 0.08292
Time: 12:34:59 Log-Likelihood: -229.26
converged: True LL-Null: -249.99
LLR p-value: 7.578e-08
==============================================================================
coef std err z P>|z| [95.0% Conf. Int.]
------------------------------------------------------------------------------
gre 0.0023 0.001 2.070 0.038 0.000 0.004
gpa 0.8040 0.332 2.423 0.015 0.154 1.454
prestige_2 -0.6754 0.316 -2.134 0.033 -1.296 -0.055
prestige_3 -1.3402 0.345 -3.881 0.000 -2.017 -0.663
prestige_4 -1.5515 0.418 -3.713 0.000 -2.370 -0.733
intercept -3.9900 1.140 -3.500 0.000 -6.224 -1.756
==============================================================================
你可以看到模型的系數(shù)，系數(shù)擬合的效果，以及總的擬合質(zhì)量，以及一些統(tǒng)計(jì)度量。[待補(bǔ)充: 模型結(jié)果主要參數(shù)的含義]
當(dāng)然你也可以只觀察結(jié)果的某部分，如置信區(qū)間(confidence interval)可以看出模型系數(shù)的健壯性。
# 查看每個(gè)系數(shù)的置信區(qū)間
print result.conf_int()
# 0 1
# gre 0.000120 0.004409
# gpa 0.153684 1.454391
# prestige_2 -1.295751 -0.055135
# prestige_3 -2.016992 -0.663416
# prestige_4 -2.370399 -0.732529
# intercept -6.224242 -1.755716
在這個(gè)例子中，我們可以肯定被錄取的可能性與應(yīng)試者畢業(yè)學(xué)校的聲望存在著逆相關(guān)的關(guān)系。
換句話說，高排名學(xué)校（prestige_1==True）的湘鄂生唄錄取的概率比低排名學(xué)校（prestige_4==True）要高。
相對(duì)危險(xiǎn)度（odds ratio）
odds ratio
OR值，是相對(duì)危險(xiǎn)度，又稱比值比、優(yōu)勢比。
使用每個(gè)變量系數(shù)的指數(shù)來生成odds ratio，可知變量每單位的增加、減少對(duì)錄取幾率的影響。例如，如果學(xué)校的聲望為2，則我們可以期待被錄取的幾率減少大概50%。UCLA上有一個(gè)對(duì)odds ratio更為深入的解釋:
在邏輯回歸中如何解釋odds ratios?。
# 輸出 odds ratio
print np.exp(result.params)
# gre 1.002267
# gpa 2.234545
# prestige_2 0.508931
# prestige_3 0.261792
# prestige_4 0.211938
# intercept 0.018500
我們也可以使用置信區(qū)間來計(jì)算系數(shù)的影響，來更好地估計(jì)一個(gè)變量影響錄取率的不確定性。
# odds ratios and 95% CI
params = result.params
conf = result.conf_int()
conf['OR'] = params
conf.columns = ['2.5%', '97.5%', 'OR']
print np.exp(conf)
# 2.5% 97.5% OR
# gre 1.000120 1.004418 1.002267
# gpa 1.166122 4.281877 2.234545
# prestige_2 0.273692 0.946358 0.508931
# prestige_3 0.133055 0.515089 0.261792
# prestige_4 0.093443 0.480692 0.211938
# intercept 0.001981 0.172783 0.018500
更深入的挖掘
為了評(píng)估我們分類器的效果，我們將使用每個(gè)輸入值的邏輯組合（logical combination）來重新創(chuàng)建數(shù)據(jù)集，如此可以得知在不同的變量下預(yù)測錄取可能性的增加、減少。首先我們使用名為 cartesian 的輔助函數(shù)來生成組合值（來源于:
如何使用numpy構(gòu)建兩個(gè)數(shù)組的組合）
我們使用 np.linspace 創(chuàng)建 “gre” 和 “gpa” 值的一個(gè)范圍，即從指定的最大、最小值來創(chuàng)建一個(gè)線性間隔的值的范圍。在本例子中，取已知的最大、最小值。
# 根據(jù)最大、最小值生成 GRE、GPA 均勻分布的10個(gè)值，而不是生成所有可能的值
gres = np.linspace(data['gre'].min(), data['gre'].max(), 10)
print gres
# array([ 220. , 284.44444444, 348.88888889, 413.33333333,
# 477.77777778, 542.22222222, 606.66666667, 671.11111111,
# 735.55555556, 800. ])
gpas = np.linspace(data['gpa'].min(), data['gpa'].max(), 10)
print gpas
# array([ 2.26 , 2.45333333, 2.64666667, 2.84 , 3.03333333,
# 3.22666667, 3.42 , 3.61333333, 3.80666667, 4. ])
# 枚舉所有的可能性
combos = pd.DataFrame(cartesian([gres, gpas, [1, 2, 3, 4], [1.]]))
# 重新創(chuàng)建啞變量
combos.columns = ['gre', 'gpa', 'prestige', 'intercept']
dummy_ranks = pd.get_dummies(combos['prestige'], prefix='prestige')
dummy_ranks.columns = ['prestige_1', 'prestige_2', 'prestige_3', 'prestige_4']
# 只保留用于預(yù)測的列
cols_to_keep = ['gre', 'gpa', 'prestige', 'intercept']
combos = combos[cols_to_keep].join(dummy_ranks.ix[:, 'prestige_2':])
# 使用枚舉的數(shù)據(jù)集來做預(yù)測
combos['admit_pred'] = result.predict(combos[train_cols])
print combos.head()
# gre gpa prestige intercept prestige_2 prestige_3 prestige_4 admit_pred
# 0 220 2.260000 1 1 0 0 0 0.157801
# 1 220 2.260000 2 1 1 0 0 0.087056
# 2 220 2.260000 3 1 0 1 0 0.046758
# 3 220 2.260000 4 1 0 0 1 0.038194
# 4 220 2.453333 1 1 0 0 0 0.179574
現(xiàn)在我們已生成了預(yù)測結(jié)果，接著通過畫圖來呈現(xiàn)結(jié)果。我編寫了一個(gè)名為 isolate_and_plot 的輔助函數(shù)，可以比較給定的變量與不同的聲望等級(jí)、組合的平均可能性。為了分離聲望和其他變量，我使用了 pivot_table 來簡單地聚合數(shù)據(jù)。
def isolate_and_plot(variable):
# isolate gre and class rank
grouped = pd.pivot_table(combos, values=['admit_pred'], index=[variable, 'prestige'],
aggfunc=np.mean)
# in case you're curious as to what this looks like
# print grouped.head()
# admit_pred
# gre prestige
# 220.000000 1 0.282462
# 2 0.169987
# 3 0.096544
# 4 0.079859
# 284.444444 1 0.311718
# make a plot
colors = 'rbgyrbgy'
for col in combos.prestige.unique():
plt_data = grouped.ix[grouped.index.get_level_values(1)==col]
pl.plot(plt_data.index.get_level_values(0), plt_data['admit_pred'],
color=colors[int(col)])
pl.xlabel(variable)
pl.ylabel("P(admit=1)")
pl.legend(['1', '2', '3', '4'], loc='upper left', title='Prestige')
pl.title("Prob(admit=1) isolating " + variable + " and presitge")
pl.show()
isolate_and_plot('gre')
isolate_and_plot('gpa')
結(jié)果圖顯示了 gre, gpa 和 prestige 如何影響錄取?？梢钥闯觯S著 gre 和 gpa 的增加，錄取可能性如何逐漸增加，并且，不同的學(xué)校聲望對(duì)錄取可能性的增加程度相差很大。
結(jié)束語
邏輯回歸是用于分類的優(yōu)秀算法，盡管有一些更加性感的，或是黑盒分類器算法，如SVM和隨機(jī)森林（RandomForest）在一些情況下性能更好，但深入了解你正在使用的模型是很有價(jià)值的。很多時(shí)候你可以使用隨機(jī)森林來篩選模型的特征，并基于篩選出的最佳的特征，使用邏輯回歸來重建模型。
其他資源
UCLA的R教程scikit-learn文檔純Python實(shí)現(xiàn)邏輯回歸邏輯回歸在線交互式例子http://www.powerxing.com/logistic-regression-in-python/筆記Python,
數(shù)據(jù)挖掘,
邏輯回歸相關(guān)文章
Python網(wǎng)絡(luò)爬蟲入門教程 Python使用pyQuery解析HTML內(nèi)容 買彩票機(jī)選和守號(hào)哪個(gè)中獎(jiǎng)概率高?
