Python邏輯迴歸原理及實際案例應用
前言
上面我們介紹了線性迴歸, 嶺迴歸, Lasso迴歸, 今天我們來看看另外一種模型—"邏輯迴歸". 雖然它有"迴歸"一詞, 但解決的卻是分類問題
目錄
正文
在前面所介紹的線性迴歸, 嶺迴歸和Lasso迴歸這三種迴歸模型中, 其輸出變數均為連續型, 比如常見的線性迴歸模型為:
其寫成矩陣形式為:
現在這裡的輸出為連續型變數, 但是實際中會有"輸出為離散型變數"這樣的需求, 比如給定特徵預測是否離職(1表示離職, 0表示不離職). 顯然這時不能直接使用線性迴歸模型, 而邏輯迴歸就派上用場了.
引用百度百科定義
邏輯(logistic)迴歸, 又稱logistic迴歸分析，是一種廣義的線性迴歸分析模型，常用於資料探勘，疾病自動診斷，經濟預測等領域。
也就是說邏輯迴歸是從線性迴歸模型推廣而來的, 我們從假設函式開始說起.
1. 假設函式
現假設因變數取值0和1, 在自變數X的條件下因變數y=1的概率為p, 記作p=P(y=1|X), 那麼y=0的概率就為1-p, 把因變數取1和取0的概率比值p/(1-p)稱為優勢比, 對優勢比取自然對數, 則可以得到Sigmoid函式:
令Sigmoid(p)=z, 則有:
而Logistic迴歸模型則是建立在Sigmoid函式和自變數的線性迴歸模型之上(這可能就是為什麼帶有"迴歸"二字的原因吧), 那麼Logistic迴歸模型可以表示為:
上式也常常被稱為邏輯迴歸模型的假設函式, 其函式影象為:
通過影象可以看出
的取值範圍為
, h(x)的取值範圍為[0, 1], 對於二分類問題來說, h(x)>=0.5則y=1, h(x)<0.5則y=0, 而且通過影象得知: 當
時, h(x)>=0.5, 因此時y=1, 否則y=0.
模型的假設函式知道了, 接下來看看損失函式.
2. 損失函式
既然邏輯迴歸是建立線上性迴歸模型之上, 那麼我們先來回顧線性迴歸的損失函式:
如果將我們邏輯迴歸的假設函式代入上式損失函式, 繪製出來的影象則是非凸函式, 其很容易產生區域性最優解, 而非全域性最優解, 為了找到使得全域性最優解, 我們需要構造一個凸函式.
由於對數函式能夠簡化計算過程, 因此這裡也是通過對數函式來構建, 先來回歸下對數函式的影象(原圖來自百度百科):
通過上圖可以發現綠線部分與我們要構造的凸函式較接近. 當a=e時, 綠線部分可以分別表示為: -loge(x)和-loge(1-x). 現將x替換為h(x)並同時加入輸出變數y (取值為1或0), 則有:
當上式中的y=1時, 其結果為-logeh(x); 當y=0時, 其結果-loge[1-h(x)].
最後, 將上式代入我們的損失函式中, 則有:
當然, 也可以用統計學中的極大似然法構造出上式損失函式. 損失函式有了, 下一步則是求解損失函式最小的演算法了.
3. 演算法
常用的求解演算法有梯度下降法, 座標軸下降法, 擬牛頓法. 下面只介紹梯度下降法(其他方法還未涉及)
你也許會有疑問, 既然是線性迴歸模型推廣而來, 那麼為什麼沒有最小二乘法呢? 最小二乘法是用來求解最小誤差平方和的演算法, 而誤差的平方和正是我們上面提到的線性迴歸的損失函式, 通過其構造出來的邏輯迴歸的損失函式是非凸的不容易找到全域性最優解, 故不選用最小二乘法, 而通過極大似然法則可以構造出凸函式, 進而可以使用梯度下降法進行求解.
對於梯度下降法的理解在這節, 這裡直接給出其表示:
具體的求解過程:
因此, 我們的梯度下降法可以寫成(其中, x0=1):
上式也被稱為批量梯度下降法, 另外兩種: 隨機梯度下降法和小批量梯度下降法分別表示為:
2. 優缺點及優化問題
1. 優點
1) 模型簡單, 訓練速度快, 且對於輸出變數有很好的概率解釋
2) 可以適用連續型和離散型自變數.
3) 可以根據實際需求設定具體的閥值
2. 缺點
1) 只能處理二分類問題.
2) 適用較大樣本量, 這是由於極大似然估計在較小樣本量中表現較差.
3) 由於其是基於線性迴歸模型之上, 因此其同樣會出現多重共線性問題.
4) 很難處理資料不均衡問題
3. 優化
1) 可以在二分類上進行推廣, 將其推廣到多分類迴歸模型
2) 對於多重共線性問題, 也可以通過刪除自變數, 進行資料變換, 正則化, 逐步迴歸, 主成分分析等方法改善, 對於正則化邏輯迴歸同樣有兩種: L1和L2, 其分別表示為:
L1正則化
L2正則化
3. 實際案例應用
1. 資料來源及背景
資料來源: https://www.kaggle.com/jiangzuo/hr-comma-sep/version/1
該資料集包含14999個樣本以及10個特徵, 通過現有員工是否離職的資料, 建立模型預測有可能離職的員工.
2. 資料概覽
1) 檢視前2行和後2行資料
import pandas as pd df = pd.read_csv(r'D:\Data\HR_comma_sep.csv') pd.set_option('display.max_rows', 4) df
10個欄位分別是: 員工對公司滿意度, 最新考核評估, 專案數, 平均每月工作時長, 工作年限, 是否出現工作事故, 是否離職, 過去5年是否升職, 崗位, 薪資水平.
可以看到除過崗位以及薪資水平是字元型外, 其餘均是數值型.
2) 檢視資料型別等資訊
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 14999 entries, 0 to 14998 Data columns (total 10 columns): satisfaction_level 14999 non-null float64 last_evaluation 14999 non-null float64 number_project 14999 non-null int64 average_montly_hours 14999 non-null int64 time_spend_company 14999 non-null int64 Work_accident 14999 non-null int64 left 14999 non-null int64 promotion_last_5years 14999 non-null int64 sales 14999 non-null object salary 14999 non-null object dtypes: float64(2), int64(6), object(2) memory usage: 1.1+ MB
前兩個特徵為浮點型, 後兩個為字元型, 其餘為整型, 且均無缺失值.
3). 描述性統計
df.describe() df.describe(include=['O']).T
滿意度: 範圍 0.09~1, 中位數0.640, 均值0.613.
最新考核評估: 範圍 0.36~1, 中位數0.720, 均值0.716
專案數: 範圍 2~7個, 中位數4, 均值3.8
平均每月工作時長 範圍96~310小時, 中位數200, 均值201
工作年限: 範圍2~10年, 中位數3, 均值3.5.
工作中出現工作事故的佔14.46%.
已經離職的佔23.81%.
過去5年升職的佔2.13%.
員工崗位有10種, 其中最多的是銷售, 多達4140.
薪資水平共有3個等級, 最多的是低等, 多達7316.
3. 資料預處理
沒有缺失值, 因此不用處理缺失值. 對於記錄來說, 其沒有唯一標識的欄位, 因此會存在重複記錄, 這裡採取不處理.
1. 異常值
通過箱線圖檢視異常值.
import seaborn as sns fig, ax = plt.subplots(1,5, figsize=(12, 2)) sns.boxplot(x=df.columns[0], data=df, ax=ax[0]) sns.boxplot(x=df.columns[1], data=df, ax=ax[1]) sns.boxplot(x=df.columns[2], data=df, ax=ax[2]) sns.boxplot(x=df.columns[3], data=df, ax=ax[3]) sns.boxplot(x=df.columns[4], data=df, ax=ax[4])
除了工作年限外, 其他均無異常值. 該異常值也反映了該公司員工中以年輕人為主
4. 視覺化分析
1. 人力資源總體情況
from pyecharts import Pie attr = ["離職", "在職"] v1 =[df.left.value_counts()[1], df.left.value_counts()[0]] pie = Pie("該公司人力資源總體情況", title_pos='center') pie.add( "", attr, v1, radius=[35, 65], label_text_color=None, is_label_show=True, legend_orient="vertical", legend_pos="left", ) pie.render()
離職3571人,佔比23.81%; 在職11428人, 佔比76.19%
2. 對公司滿意度與是否離職的關係
from pyecharts import Boxplot #欄位重新命名 df.columns=['satisfaction', 'evaluation', 'project', 'hours', 'years_work','work_accident', 'left', 'promotion', 'department', 'salary'] #繪製箱線圖 boxplot = Boxplot("對公司滿意度與是否離職關係圖", title_pos='center') x_axis = ['在職', '離職'] y_axis = [df[df.left == 0].satisfaction.values, df[df.left == 1].satisfaction.values] boxplot.add("", x_axis, boxplot.prepare_data(y_axis)) boxplot.render()
就中位數而言, 離職人員對公司滿意度相對較低, 且離職人員對公司滿意度整體波動較大. 另外離職人員中沒有滿意度為1的評價.
3. 最新考核評估與是否離職的關係
boxplot = Boxplot("最新評估與是否離職關係圖", title_pos='center') x_axis = ['在職', '離職'] y_axis = [df[df.left == 0].evaluation.values, df[df.left == 1].evaluation.values] boxplot.add("", x_axis, boxplot.prepare_data(y_axis)) boxplot.render()
就中位數而言, 離職人員的最新考核評估相對較高, 但其波動也大.
4. 所參加專案與是否離職的關係
from pyecharts import Bar, Pie, Grid #按照專案數分組分別求離職人數和所有人數 project_left_1 = df[df.left == 1].groupby('project')['left'].count() project_all = df.groupby('project')['left'].count() #分別計算離職人數和在職人數所佔比例 project_left1_rate = project_left_1 / project_all project_left0_rate = 1 - project_left1_rate attr = project_left1_rate.index bar = Bar("所參加專案數與是否離職的關係圖", title_pos='10%') bar.add("離職", attr, project_left1_rate, is_stack=True) bar.add("在職", attr, project_left0_rate, is_stack=True, legend_pos="left", legend_orient="vertical") #繪製圓環圖 pie = Pie("各專案數所佔百分比", title_pos='center') pie.add('', project_all.index, project_all, radius=[35, 60], label_text_color=None, is_label_show=True, legend_orient="vertical", legend_pos="67%") grid = Grid(width=1200) grid.add(bar, grid_right="67%") grid.add(pie) grid.render()
通過下圖可以發現以下2點:
- 離職人員所佔比例隨著專案數的增多而增大, 2個專案數是特例
- 離職人員比例較高的專案數2, 6, 7在總專案數中所佔百分比相對較少. 專案數為2的這部分人可能是工作能力不被認可, 其離職人數也相對較高; 專案數為6, 7的這部分人則工作能力較強, 其可能在其他企業能有更好的發展, 自然離職比例也相對較高.
5. 平均每月工作時長和是否離職的關係
boxplot = Boxplot("平均每月工作時長與是否離職關係圖", title_pos='center') x_axis = ['在職', '離職'] y_axis = [df[df.left == 0].hours.values, df[df.left == 1].hours.values] boxplot.add("", x_axis, boxplot.prepare_data(y_axis)) boxplot.render()
通過下圖可以看到: 離職人員的平均每月工作時長相對較長, 每月按照22個工作日計算, 每日工作時數的中位數為10.18小時, 最大值為14.09小時.
6. 工作年限和是否離職的關係
from pyecharts import Bar, Pie, Grid #按照工作年限分別求離職人數和所有人數 years_left_0 = df[df.left == 0].groupby('years_work')['left'].count() years_all = df.groupby('years_work')['left'].count() #分別計算離職人數和在職人數所佔比例 years_left0_rate = years_left_0 / years_all years_left1_rate = 1 - years_left0_rate attr = years_all.index bar = Bar("工作年限與是否離職的關係圖", title_pos='10%') bar.add("離職", attr, years_left1_rate, is_stack=True) bar.add("在職", attr, years_left0_rate, is_stack=True, legend_pos="left" , legend_orient="vertical") #繪製圓環圖 pie = Pie("各工作年限所佔百分比", title_pos='center') pie.add('', years_all.index, years_all, radius=[35, 60], label_text_color=None, is_label_show=True, legend_orient="vertical", legend_pos="67%") grid = Grid(width=1200) grid.add(bar, grid_right="67%") grid.add(pie) grid.render()
通過下圖可以得出:
- 在各工作年限中, 離職人員較集中於3, 4, 5, 6年, 而6年以上則相對穩定
- 企業中3年人數所佔百分比最多, 其次是2年, 主要以年輕人為主
7. 是否發生工作事故與是否離職的關係
from pyecharts import Bar accident_left = pd.crosstab(df.work_accident, df.left) attr = accident_left.index bar = Bar("是否發生工作事故與是否離職的關係圖", title_pos='center') bar.add("離職", attr, accident_left[1], is_stack=True) bar.add("在職", attr, accident_left[0], is_stack=True, legend_pos="left" , legend_orient="vertical", is_label_show=True) bar.render()
可以看到少部分出現工作事故, 且其中有較少部分人離職.
8. 5年內是否升職與是否離職的關係
promotion_left = pd.crosstab(df.promotion, df.left) attr = promotion_left.index bar = Bar("5年內是否升職與是否離職的關係圖", title_pos='center') bar.add("離職", attr, promotion_left[1], is_stack=True) bar.add("在職", attr, promotion_left[0], is_stack=True, legend_pos="left" , legend_orient="vertical", is_label_show=True) bar.render()
5年內多數人沒有升職, 離職率就相對較高.
9. 崗位與是否離職的關係
#分別計算各崗位離職人員比例和各崗位佔總體百分比 department_left_0 = df[df.left == 0].groupby('department')['left'].count() department_all = df.groupby('department')['left'].count() department_left0_rate = department_left_0 / department_all department_left1_rate = 1 - department_left0_rate attr = department_all.index bar = Bar("崗位與離職比例的關係圖", title_top='40%') bar.add("離職", attr, department_left1_rate, is_stack=True) bar.add("在職", attr, department_left0_rate, is_stack=True, is_datazoom_show=True, xaxis_interval=0, xaxis_rotate=30, legend_top="45%", legend_pos="80%") #繪製圓環圖 pie = Pie("各個崗位所佔百分比", title_pos='left') pie.add('', department_all.index, department_all,center=[50, 23], radius=[18, 35], label_text_color=None, is_label_show=True, legend_orient="vertical", legend_pos="80%", legend_top="4%") grid = Grid(width=1200, height=700) grid.add(bar, grid_top="50%", grid_bottom="25%") grid.add(pie) grid.render()
通過下圖可以看出:
- 銷售崗位所佔百分比最多, 達到27.6%, 最少是管理層, 其所佔百分比是4.2%
- 令人意外的是hr崗位離職比例最大.
10. 薪資水平和是否離職的關係
from pyecharts import Bar #按照薪資水平分別求離職人數和所有人數 salary_left = pd.crosstab(df.salary, df.left).sort_values(0, ascending = False) attr = salary_left.index bar = Bar("薪資水平和是否離職的關係圖", title_pos='center') bar.add("離職", attr, salary_left[1], is_stack=True) bar.add("在職", attr, salary_left[0], is_stack=True, legend_pos="left" , legend_orient="vertical", is_label_show=True) bar.render()
薪資分為三個水平: 低等, 中等, 高等. 低等水平離職人數最多, 所佔比例也最大, 而高等則最少.
5. 特徵工程
1. 離散型資料處理
離散型資料可分為兩種: 一種是定序, 一種是定類.
1) 定序
薪資水平其含有順序意義, 因此將其字元型轉化為數值型
df['salary'] = df.salary.map({"low": 0, "medium": 1, "high": 2}) df.salary.unique()
array([0, 1, 2], dtype=int64)
2) 定類
崗位是定型別變數, 對其進行one-hot編碼, 這裡直接利用pandas的get_dummies方法.
df_one_hot = pd.get_dummies(df, prefix="dep") df_one_hot.shape
(14999, 19)
2. 連續型資料處理
邏輯迴歸模型能夠適應連續型變數, 因此可以不用進行離散化處理, 又由於多個特徵之間差異差異較大會造成梯度下降演算法收斂速度變慢, 故進行歸一化處理
#採用max-min歸一化方法 hours = df_one_hot['hours'] df_one_hot['hours'] = df_one_hot.hours.apply(lambda x: (x-hours.min()) / (hours.max()-hours.min()))
3. 相關係數
兩個變數均是連續型且具有線性關係, 則可以使用皮爾遜相關係數, 否則使用斯皮爾曼相關係數, 這裡採用斯皮爾曼相關係數
#計算相關係數 correlation = df_one_hot.corr(method = "spearman") plt.figure(figsize=(18, 10)) #繪製熱力圖 sns.heatmap(correlation, linewidths=0.2, vmax=1, vmin=-1, linecolor='w',fmt='.2f', annot=True,annot_kws={'size':10},square=True)
6. 邏輯迴歸模型
1. 劃分資料集
from sklearn.model_selection import train_test_split #劃分訓練集和測試集 X = df_one_hot.drop(['left'], axis=1) y = df_one_hot['left'] X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)
2. 訓練模型
from sklearn.linear_model import LogisticRegression LR = LogisticRegression() print(LR.fit(X_train, y_train)) print("訓練集準確率: ", LR.score(X_train, y_train)) print("測試集準確率: ", LR.score(X_test, y_test))
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
intercept_scaling=1, max_iter=100, multi_class='warn',
n_jobs=None, penalty='l2', random_state=None, solver='warn',
tol=0.0001, verbose=0, warm_start=False)
訓練集準確率: 0.7978998249854155
測試集準確率: 0.7966666666666666
參考官方文件說明, 引數C是正則化項引數的倒數, C的數值越小, 懲罰的力度越大. penalty可選L1, L2正則化項, 預設是L2正則化.
引數solver可選{‘newton-cg’, ‘lbfgs’, ‘liblinear’, ‘sag’, ‘saga’}這5個優化演算法:
newton-cg, lbfgs是擬牛頓法, liblinear是座標軸下降法, sag, saga是隨機梯度下降法, saga可以適用於L1和L2正則化項, 而sag只能用於L2正則化項.
#指定隨機梯度下降優化演算法 LR = LogisticRegression(solver='saga') print(LR.fit(X_train, y_train)) print("訓練集準確率: ", LR.score(X_train, y_train)) print("測試集準確率: ", LR.score(X_test, y_test))
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
intercept_scaling=1, max_iter=100, multi_class='warn',
n_jobs=None, penalty='l2', random_state=None, solver='saga',
tol=0.0001, verbose=0, warm_start=False)
訓練集準確率: 0.7980665055421285
測試集準確率: 0.7973333333333333
在選擇隨機梯度下降法後, 訓練集和測試集準確率均略有提升.
3. 調參
#用準確率進行10折交叉驗證選擇合適的引數C from sklearn.linear_model import LogisticRegressionCV Cs = 10**np.linspace(-10, 10, 400) lr_cv = LogisticRegressionCV(Cs=Cs, cv=10, penalty='l2', solver='saga', max_iter=10000, scoring='accuracy') lr_cv.fit(X_train, y_train) lr_cv.C_
array([25.52908068])
用該引數進行預測
LR = LogisticRegression(solver='saga', penalty='l2', C=25.52908068) print("訓練集準確率: ", LR.score(X_train, y_train)) print("測試集準確率: ", LR.score(X_test, y_test))
訓練集準確率: 0.7984832069339112 測試集準確率: 0.798
訓練集和測試集準確率均有所提升, 對於二分類問題, 準確率有時不是很好的評估方法, 這時需要用到混淆矩陣
4. 混淆矩陣
from sklearn import metrics X_train_pred = LR.predict(X_train) X_test_pred = LR.predict(X_test) print('訓練集混淆矩陣:') print(metrics.confusion_matrix(y_train, X_train_pred)) print('測試集混淆矩陣:') print(metrics.confusion_matrix(y_test, X_test_pred))
訓練集混淆矩陣: [[8494 647] [1771 1087]] 測試集混淆矩陣: [[2112 175] [ 431 282]]
from sklearn.metrics import classification_report print('訓練集:') print(classification_report(y_train, X_train_pred)) print('測試集:') print(classification_report(y_test, X_test_pred))
訓練集:
precision recall f1-score support
0 0.83 0.93 0.88 9141
1 0.63 0.38 0.47 2858
micro avg 0.80 0.80 0.80 11999
macro avg 0.73 0.65 0.67 11999
weighted avg 0.78 0.80 0.78 11999
測試集:
precision recall f1-score support
0 0.83 0.92 0.87 2287
1 0.62 0.40 0.48 713
micro avg 0.80 0.80 0.80 3000
macro avg 0.72 0.66 0.68 3000
weighted avg 0.78 0.80 0.78 3000
在訓練集有0.83的精準率和0.93的召回率, 在測試集上有0.83的精準率和0.92的召回率.
7. 樸素貝葉斯模型
樸素貝葉斯模型是基於特徵條件獨立假設和貝葉斯理論.
from sklearn.naive_bayes import GaussianNB from sklearn.model_selection import cross_val_score #構建高斯樸素貝葉斯模型 gnb = GaussianNB() gnb.fit(X_train, y_train) print("訓練集準確率: ", gnb.score(X_train, y_train)) print("測試集準確率: ", gnb.score(X_test, y_test)) X_train_pred =gnb.predict(X_train) X_test_pred = gnb.predict(X_test) print('訓練集混淆矩陣:') print(metrics.confusion_matrix(y_train, X_train_pred)) print('測試集混淆矩陣:') print(metrics.confusion_matrix(y_test, X_test_pred)) print('訓練集:') print(classification_report(y_train, X_train_pred)) print('測試集:') print(classification_report(y_test, X_test_pred))
訓練集準確率: 0.7440620051670973
測試集準確率: 0.741
訓練集混淆矩陣:
[[6791 2350]
[ 721 2137]]
測試集混淆矩陣:
[[1680 607]
[ 170 543]]
訓練集:
precision recall f1-score support
0 0.90 0.74 0.82 9141
1 0.48 0.75 0.58 2858
micro avg 0.74 0.74 0.74 11999
macro avg 0.69 0.75 0.70 11999
weighted avg 0.80 0.74 0.76 11999
測試集:
precision recall f1-score support
0 0.91 0.73 0.81 2287
1 0.47 0.76 0.58 713
micro avg 0.74 0.74 0.74 3000
macro avg 0.69 0.75 0.70 3000
weighted avg 0.80 0.74 0.76 3000
可以看到其準確率較邏輯迴歸低, 但是精準率高於邏輯迴歸.
8. ROC曲線
from sklearn import metrics from sklearn.metrics import roc_curve #將邏輯迴歸模型和高斯樸素貝葉斯模型預測出的概率均與實際值通過roc_curve比較返回假正率, 真正率, 閾值 lr_fpr, lr_tpr, lr_thresholds = roc_curve(y_test, LR.predict_proba(X_test)[:,1]) gnb_fpr, gnb_tpr, gnb_thresholds = roc_curve(y_test, gnb.predict_proba(X_test)[:,1]) #分別計算這兩個模型的auc的值, auc值就是roc曲線下的面積 lr_roc_auc = metrics.auc(lr_fpr, lr_tpr) gnb_roc_auc = metrics.auc(gnb_fpr, gnb_tpr) plt.figure(figsize=(8, 5)) plt.plot([0, 1], [0, 1],'--', color='r') plt.plot(lr_fpr, lr_tpr, label='LogisticRegression(area = %0.2f)' % lr_roc_auc) plt.plot(gnb_fpr, gnb_tpr, label='GaussianNB(area = %0.2f)' % gnb_roc_auc) plt.xlim([0.0, 1.0]) plt.ylim([0.0, 1.0]) plt.title('ROC') plt.xlabel('FPR') plt.ylabel('TPR') plt.legend() plt.show()
ROC曲線越靠近左上角說明分類效果越好, 與之對應的auc的值就越大. 對於該資料集來說, 高斯樸素貝葉斯模型略優於邏輯迴歸模型
總結
本文從邏輯迴歸模型的原理開始介紹, 並通過實際案例對邏輯迴歸模型進行應用, 但是結果還不是很好, 一方面是模型表現不是很好(當然, 也僅僅用到了邏輯迴歸和樸素貝葉斯) ; 另一方面是特徵工程沒有處理好(樸素貝葉斯模型對特徵條件獨立假設較敏感), 應當進行特徵選擇, 比如主成分分析.
參考資料:
網易雲課堂《吳恩達機器學習》
《Python資料分析與挖掘實戰》
宣告: 僅用於學習交流