机器学习——Logit模型(python)
1 OLS 缺陷
传统线性回归的响应变量为连续变量,使用最小二乘法或极大似然估计可以获得各变量的回归系数。当响应变量为二值(如企业是否采取投资策略,个人选择工作还是深造等),使用普通最小二乘法进行估计存在严重的问题:给定如下解释二值响应变量 y i ∈ { 0 , 1 } y_i\in\{0,1\} yi∈{0,1}的模型
y i = β 1 x i 1 + β 2 x i 2 + ⋯ + β p x i p + ε i = x i ′ β + ε i y_i=\beta_1 x_{i 1}+\beta_2 x_{i 2}+\cdots+\beta_p x_{i p}+\varepsilon_i=\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}+\varepsilon_i yi=β1xi1+β2xi2+⋯+βpxip+εi=xi′β+εi
- 当 y i = 1 y_i=1 yi=1, 1 − x i ′ β = ε i 1 -\boldsymbol{x}_i^{\prime} \boldsymbol{\beta} = \varepsilon_i 1−xi′β=εi;
- 当 y i = 0 y_i=0 yi=0, − x i ′ β = ε i -\boldsymbol{x}_i^{\prime} \boldsymbol{\beta} = \varepsilon_i −xi′β=εi
c o v ( ε i , x i ) ≠ 0 cov(\varepsilon_i,\boldsymbol{x}_i)\ne0 cov(εi,xi)=0, E ( ε i 2 ∣ x i ) E(\varepsilon_i^2|\boldsymbol{x}_i) E(εi2∣xi)与 x i \boldsymbol{x}_i xi有关,因此使用普通最小二乘法存在内生性和异方差问题,难以进行预测。
2 估计策略
通过对不同特征变量进行线性组合,如果存在某一阈值
c
∈
(
0
,
1
)
c\in(0,1)
c∈(0,1)和单调递增的连接函数
G
(
x
i
,
β
)
G(\boldsymbol{x}_i,\boldsymbol{\beta})
G(xi,β),当连接函数
G
(
x
i
)
>
c
G(\boldsymbol{x}_i)>c
G(xi)>c时,选择
y
i
=
1
y_i = 1
yi=1,当
G
(
x
i
,
β
)
<
c
G(\boldsymbol{x}_i,\boldsymbol{\beta}) <c
G(xi,β)<c,选择
y
i
=
0
y_i = 0
yi=0。这里要求连接函数
G
(
x
i
,
β
)
G(\boldsymbol{x}_i,\boldsymbol{\beta})
G(xi,β)值域为[0,1]。可见连续随机变量的累积概率分布函数满足上述条件。累积概率分布函数一般采用正态累积分布函数
Φ
\Phi
Φ和逻辑累积分布函数
F
F
F。当使用
Φ
\Phi
Φ时,则为Probit模型,使用
F
F
F则为Logit模型。由于
Φ
\Phi
Φ不存在显式,不利于后续计算,因此一般采用Logit模型。函数为
F
F
F具有显性表达式
F
=
1
1
+
e
−
z
,
z
∈
R
F =\dfrac{1}{1+e^{-z}},z\in R
F=1+e−z1,z∈R
该函数连续单调递增,关于点
(
0
,
0.5
)
(0,0.5)
(0,0.5)对称。当
z
→
−
∞
z\to-\infty
z→−∞,
F
→
0
F\to 0
F→0;当
z
→
∞
z\to\infty
z→∞,
F
→
1
F\to 1
F→1。令
z
=
x
i
′
β
z = \boldsymbol{x}_i^{\prime} \boldsymbol{\beta}
z=xi′β,当特征变量线性组合使
F
>
0.5
F>0.5
F>0.5时,
y
=
1
y=1
y=1为正例;反之
F
<
0.5
F<0.5
F<0.5时,
y
=
0
y=0
y=0为反例;
F
=
0.5
F=0.5
F=0.5,
y
y
y位于
y
=
0
,
1
y=0,1
y=0,1的边界。换言之特征变量组合能映射到
F
∈
(
0
,
1
)
F\in(0,1)
F∈(0,1)中。因此
y
y
y的条件期望可以写作
{
P
(
y
=
1
∣
x
)
=
F
(
x
,
β
)
=
1
1
+
e
−
z
P
(
y
=
0
∣
x
)
=
1
−
F
(
x
,
β
)
=
e
−
z
1
+
e
−
z
\left\{\begin{array}{l} \mathrm{P}(y=1 \mid \boldsymbol{x})=F(\boldsymbol{x}, \boldsymbol{\beta})=\dfrac{1}{1+e^{-z}} \\ \mathrm{P}(y=0 \mid \boldsymbol{x})=1-F(\boldsymbol{x}, \boldsymbol{\beta})=\dfrac{e^{-z}}{1+e^{-z}} \end{array}\right.
⎩
⎨
⎧P(y=1∣x)=F(x,β)=1+e−z1P(y=0∣x)=1−F(x,β)=1+e−ze−z
即
y
y
y服从两点分布。进一步将上述两点概率分布写成统一形式
P
(
y
i
∣
x
i
,
β
)
=
[
F
(
x
i
′
β
)
]
y
i
[
1
−
F
(
x
i
′
β
)
]
1
−
y
i
P\left(y_i \mid \boldsymbol{x}_i,\boldsymbol{\beta}\right)=\left[F\left(\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}\right)\right]^{y_i}\left[1-F\left(\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}\right)\right]^{1-y_i}
P(yi∣xi,β)=[F(xi′β)]yi[1−F(xi′β)]1−yi
其似然函数为
L
(
β
∣
y
,
X
)
=
∏
i
=
1
n
[
F
(
x
i
′
β
)
]
y
i
[
1
−
F
(
x
i
′
β
)
]
1
−
y
i
L(\boldsymbol{\beta} \mid \mathbf{y}, \mathbf{X})=\prod_{i=1}^n\left[F\left(\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}\right)\right]^{y_i}\left[1-F\left(\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}\right)\right]^{1-y_i}
L(β∣y,X)=i=1∏n[F(xi′β)]yi[1−F(xi′β)]1−yi
对上式取对数,通过确定参数
β
\boldsymbol{\beta}
β最大化似然函数
β
=
a
r
g
max
β
ln
L
(
β
∣
y
,
X
)
=
∑
i
=
1
n
y
i
ln
[
F
(
x
i
′
β
)
]
+
∑
i
=
1
n
(
1
−
y
i
)
ln
[
1
−
F
(
x
i
′
β
)
]
{\boldsymbol{\beta}} = arg\max _{\boldsymbol{\beta}} \ln L(\boldsymbol{\beta} \mid \boldsymbol{y}, \mathbf{X})=\sum_{i=1}^n y_i \ln \left[F\left(\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}\right)\right]+\sum_{i=1}^n\left(1-y_i\right) \ln \left[1-F\left(\boldsymbol{x}_i^{\prime} \boldsymbol{\beta}\right)\right]
β=argβmaxlnL(β∣y,X)=i=1∑nyiln[F(xi′β)]+i=1∑n(1−yi)ln[1−F(xi′β)]
由于似然函数是非线性函数,通过最优化方法(梯度下降法、牛顿法)数值求解。
3 模型解释
3.1 边际影响
参数
β
\boldsymbol{\beta}
β并非特征变量的边际影响,因为
∂
P
(
y
=
1
∣
x
)
∂
x
k
=
∂
F
(
x
′
β
)
∂
(
x
′
β
)
⋅
∂
(
x
′
β
)
∂
x
k
=
F
(
x
′
β
)
⋅
β
k
\frac{\partial \mathrm{P}(y=1 \mid \boldsymbol{x})}{\partial x_k}=\frac{\partial F\left(\boldsymbol{x}^{\prime} \boldsymbol{\beta}\right)}{\partial\left(\boldsymbol{x}^{\prime} \boldsymbol{\beta}\right)} \cdot \frac{\partial\left(\boldsymbol{x}^{\prime} \boldsymbol{\beta}\right)}{\partial x_k}=F\left(\boldsymbol{x}^{\prime} \boldsymbol{\beta}\right) \cdot \boldsymbol{\beta}_k
∂xk∂P(y=1∣x)=∂(x′β)∂F(x′β)⋅∂xk∂(x′β)=F(x′β)⋅βk
可见特征变量
x
k
x_k
xk对
y
=
1
y=1
y=1(选择正例)的边际影响并非常数,取决于给定的不同特征变量线性组合下的
F
F
F值与第
k
k
k个特征变量的系数
β
k
\beta_k
βk的乘积。
β
k
\beta_k
βk符号与
x
k
x_k
xk边际影响符号一致(
F
(
x
′
β
)
>
0
F\left(\boldsymbol{x}^{\prime} \boldsymbol{\beta}\right)>0
F(x′β)>0)。
3.2 相对风险
令
p
=
P
(
y
=
1
∣
x
)
=
1
1
+
e
−
z
p = \mathrm{P}(y=1 \mid \boldsymbol{x})=\dfrac{1}{1+e^{-z}}
p=P(y=1∣x)=1+e−z1表示
y
=
1
y=1
y=1的概率,则
1
−
p
=
e
−
z
1
+
e
−
z
1-p = \dfrac{e^{-z}}{1+e^{-z}}
1−p=1+e−ze−z表示
y
=
0
y=0
y=0的概率,则
p
1
−
p
=
exp
(
x
′
β
)
\frac{p}{1-p}=\exp \left(\boldsymbol{x}^{\prime} \boldsymbol{\beta}\right)
1−pp=exp(x′β)
两边取对数
l
n
p
1
−
p
=
x
′
β
ln\frac{p}{1-p}=\boldsymbol{x}^{\prime} \boldsymbol{\beta}
ln1−pp=x′β
其中
o
d
d
s
=
l
n
p
1
−
p
odds= ln\dfrac{p}{1-p}
odds=ln1−pp表示几率比或相对风险,表示正例和反例概率之比的对数。若选择
y
=
1
y=1
y=1的概率
p
p
p越大,则几率比越高。显然
∂
o
d
d
s
∂
x
k
=
β
k
\dfrac{\partial odds}{\partial x_k} = \beta_k
∂xk∂odds=βk
在使用logit模型时,
β
k
\beta_k
βk即为特征变量增加一单位,
y
=
1
y=1
y=1的相对风险增加
β
k
\beta_k
βk单位。当
x
k
x_k
xk是离散的(年龄),则无法使用导数,
x
k
x_k
xk变化一单位后的相对风险与变化前的相对风险之比为
p
∗
1
−
p
∗
p
1
−
p
=
exp
[
β
1
x
1
+
⋯
+
β
k
(
x
k
+
1
)
+
⋯
+
β
p
x
p
]
exp
(
β
1
x
1
+
⋯
+
β
k
x
k
+
⋯
+
β
p
x
p
)
=
e
k
β
\dfrac{\dfrac{p^*}{1-p^*}}{\dfrac{p}{1-p}}=\frac{\exp \left[\beta_1 x_1+\cdots+\beta_k\left(x_k+1\right)+\cdots+\beta_p x_p\right]}{\exp \left(\beta_1 x_1+\cdots+\beta_k x_k+\cdots+\beta_p x_p\right)}=e^\beta_k
1−pp1−p∗p∗=exp(β1x1+⋯+βkxk+⋯+βpxp)exp[β1x1+⋯+βk(xk+1)+⋯+βpxp]=ekβ
4 拟合优度
Logit模型拟合效果用Pseudo
R
2
R^2
R2表示,计算公式为
P
s
e
u
d
o
R
2
≡
ln
L
0
−
ln
L
1
ln
L
0
Pseudo R^2 \equiv \frac{\ln L_0-\ln L_1}{\ln L_0}
PseudoR2≡lnL0lnL0−lnL1
l
n
L
1
ln L_1
lnL1 为原模型的对数似然函数之最大值,
l
n
L
0
ln L_0
lnL0仅以常数项为变量的对数似然函数之最大值。
5 Python模拟
import numpy as np
from scipy.stats import norm
import pandas as pd
import statsmodels.api as sm
# DGP
np.random.seed(123456)
n = 30000
x = norm.rvs(loc=0, scale=1, size=n).reshape(10000, 3)
x = pd.DataFrame(x)
x.columns = ['x1', 'x2', 'x3']
# 误差项
u = pd.DataFrame(norm.rvs(loc=0, scale=0.5, size=10000).reshape(10000, 1))
y = 2 * x['x1'] - 3* x['x2'] + 4 * x['x3'] +u.iloc[1:10000,0]
y = (y > 0)
# logit估计
results = sm.Logit(y, x).fit()
print(results.summary())
# Optimization terminated successfully.
# Current function value: 0.488599
# Iterations 6
# Logit Regression Results
# ==============================================================================
# Dep. Variable: y No. Observations: 10000
# Model: Logit Df Residuals: 9997
# Method: MLE Df Model: 2
# Date: Wed, 05 Oct 2022 Pseudo R-squ.: 0.2951
# Time: 23:00:29 Log-Likelihood: -4886.0
# converged: True LL-Null: -6931.5
# Covariance Type: nonrobust LLR p-value: 0.000
# ==============================================================================
# coef std err z P>|z| [0.025 0.975]
# ------------------------------------------------------------------------------
# x1 0.6953 0.027 25.644 0.000 0.642 0.748
# x2 -0.9989 0.029 -34.333 0.000 -1.056 -0.942
# x3 1.3071 0.032 40.867 0.000 1.244 1.370
# ==============================================================================
# 边际系数
margin = results.get_margeff()
print(margin.summary())
# Logit Marginal Effects
# =====================================
# Dep. Variable: y
# Method: dydx
# At: overall
# ==============================================================================
# dy/dx std err z P>|z| [0.025 0.975]
# ------------------------------------------------------------------------------
# x1 0.1123 0.004 29.165 0.000 0.105 0.120
# x2 -0.1613 0.004 -44.566 0.000 -0.168 -0.154
# x3 0.2111 0.003 62.726 0.000 0.205 0.218
# ==============================================================================
# 相对风险、几率比
odd = np.exp(results.params)
print('odds:\n',odd)
# odds:
# x1 2.004249
# x2 0.368295
# x3 3.695308
# dtype: float64
# 混淆矩阵
table = results.pred_table()
print(f'混淆矩阵\n{table}\n')
# [[3778. 1217.]
# [1217. 3788.]]
基于混淆矩阵即可计算出预测准确率、错分率、真阳率等指标。
陈强,《机器学习及Python应用》高等教育出版社, 2021年3月