机器学习——Logit模型

机器学习——Logit模型

1 OLS 缺陷

传统线性回归的响应变量为连续变量，使用最小二乘法或极大似然估计可以获得各变量的回归系数。当响应变量为二值(如企业是否采取投资策略，个人选择工作还是深造等)，使用普通最小二乘法进行估计存在严重的问题：给定如下解释二值响应变量 y i ∈ { 0 , 1 } y_i\in\{0,1\} yi​∈{0,1}的模型

$$y_i={\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_p{x}_{ip}+{\epsilon }_i={\mathbf{x}}_i^{\prime }\mathbf{\beta }+{\epsilon }_i$$

• 当$y_i=1$，$1-{\mathbf{x}}_i^{\prime }\mathbf{\beta }={\epsilon }_i$；
• 当$y_i=0$, $-{\mathbf{x}}_i^{\prime }\mathbf{\beta }={\epsilon }_i$

$cov({\epsilon }_i,{\mathbf{x}}_i)\ne 0$，$E({\epsilon }_i^2|{\mathbf{x}}_i)$与${\mathbf{x}}_i$有关，因此使用普通最小二乘法存在内生性和异方差问题，难以进行预测。

2 估计策略

通过对不同特征变量进行线性组合，如果存在某一阈值$c\in (0,1)$和单调递增的连接函数$G({\mathbf{x}}_i,\mathbf{\beta })$，当连接函数$G({\mathbf{x}}_i)>c$时，选择$y_i=1$，当$G({\mathbf{x}}_i,\mathbf{\beta })<c$，选择$y_i=0$。这里要求连接函数$G({\mathbf{x}}_i,\mathbf{\beta })$值域为[0,1]。可见连续随机变量的累积概率分布函数满足上述条件。累积概率分布函数一般采用正态累积分布函数$\Phi$和逻辑累积分布函数$F$。当使用$\Phi$时，则为Probit模型，使用$F$则为Logit模型。由于$\Phi$不存在显式，不利于后续计算，因此一般采用Logit模型。函数为$F$具有显性表达式
$$F=\frac{1}{1+e^{-z}},z\in R$$
该函数连续单调递增，关于点$(0,0.5)$对称。当$z\to -\infty$，$F\to 0$；当$z\to \infty$，$F\to 1$。令$z={\mathbf{x}}_i^{\prime }\mathbf{\beta }$，当特征变量线性组合使$F>0.5$时，$y=1$为正例；反之$F<0.5$时，$y=0$为反例；$F=0.5$，$y$位于$y=0,1$的边界。换言之特征变量组合能映射到$F\in (0,1)$中。因此$y$的条件期望可以写作
$$\{\begin{array}{l}\mathrm{P}(y=1\mid \mathbf{x})=F(\mathbf{x},\mathbf{\beta })=\frac{1}{1+e^{-z}}\\ \mathrm{P}(y=0\mid \mathbf{x})=1-F(\mathbf{x},\mathbf{\beta })=\frac{e^{-z}}{1+e^{-z}}\end{array}$$
即$y$服从两点分布。进一步将上述两点概率分布写成统一形式
$$P\left(y_i\mid {\mathbf{x}}_i,\mathbf{\beta }\right)={\left[F\left({\mathbf{x}}_i^{\prime }\mathbf{\beta }\right)\right]}^{y_i}{\left[1-F\left({\mathbf{x}}_i^{\prime }\mathbf{\beta }\right)\right]}^{1-y_i}$$
其似然函数为
$$L(\mathbf{\beta }\mid \mathbf{y},\mathbf{X})=\prod _{i=1}^n{\left[F\left({\mathbf{x}}_i^{\prime }\mathbf{\beta }\right)\right]}^{y_i}{\left[1-F\left({\mathbf{x}}_i^{\prime }\mathbf{\beta }\right)\right]}^{1-y_i}$$
对上式取对数，通过确定参数$\mathbf{\beta }$最大化似然函数
$$\mathbf{\beta }=arg\underset{\mathbf{\beta }}{\max }\ln L(\mathbf{\beta }\mid \mathbf{y},\mathbf{X})=\sum _{i=1}^n{y}_i\ln \left[F\left({\mathbf{x}}_i^{\prime }\mathbf{\beta }\right)\right]+\sum _{i=1}^n\left(1-y_i\right)\ln \left[1-F\left({\mathbf{x}}_i^{\prime }\mathbf{\beta }\right)\right]$$
由于似然函数是非线性函数，通过最优化方法(梯度下降法、牛顿法)数值求解。

3 模型解释

3.1 边际影响

参数$\mathbf{\beta }$并非特征变量的边际影响，因为
$$\frac{\partial \mathrm{P}(y=1\mid \mathbf{x})}{\partial x_k}=\frac{\partial F\left({\mathbf{x}}^{\prime }\mathbf{\beta }\right)}{\partial \left({\mathbf{x}}^{\prime }\mathbf{\beta }\right)}\cdot \frac{\partial \left({\mathbf{x}}^{\prime }\mathbf{\beta }\right)}{\partial x_k}=F\left({\mathbf{x}}^{\prime }\mathbf{\beta }\right)\cdot {\mathbf{\beta }}_k$$
可见特征变量$x_k$对$y=1$(选择正例)的边际影响并非常数，取决于给定的不同特征变量线性组合下的$F$值与第$k$个特征变量的系数${\beta }_k$的乘积。${\beta }_k$符号与$x_k$边际影响符号一致($F\left({\mathbf{x}}^{\prime }\mathbf{\beta }\right)>0$)。

3.2 相对风险

令$p=\mathrm{P}(y=1\mid \mathbf{x})=\frac{1}{1+e^{-z}}$表示$y=1$的概率，则$1-p=\frac{e^{-z}}{1+e^{-z}}$表示$y=0$的概率，则
$$\frac{p}{1-p}=\exp \left({\mathbf{x}}^{\prime }\mathbf{\beta }\right)$$
两边取对数
$$ln\frac{p}{1-p}={\mathbf{x}}^{\prime }\mathbf{\beta }$$
其中$odds=ln\frac{p}{1-p}$表示几率比或相对风险，表示正例和反例概率之比的对数。若选择$y=1$的概率$p$越大，则几率比越高。显然
$$\frac{\partial odds}{\partial x_k}={\beta }_k$$
在使用logit模型时，${\beta }_k$即为特征变量增加一单位，$y=1$的相对风险增加${\beta }_k$单位。当$x_k$是离散的(年龄)，则无法使用导数，$x_k$变化一单位后的相对风险与变化前的相对风险之比为
$$\frac{\frac{p^{\ast }}{1-p^{\ast }}}{\frac{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_k^{\beta }$$

4 拟合优度

Logit模型拟合效果用Pseudo$R^2$表示，计算公式为
$$Pseudo{R}^2\equiv \frac{\ln L_0-\ln L_1}{\ln L_0}$$
$ln{L}_1$ 为原模型的对数似然函数之最大值，$ln{L}_0$仅以常数项为变量的对数似然函数之最大值。

---

5 代码实现

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.]]

基于混淆矩阵即可计算出预测准确率、错分率、真阳率等指标。

---

-END-

陈强,《机器学习及Python应用》高等教育出版社, 2021年3月