1.4 Logit 回归:二元因变量模型
"Statistical thinking will one day be as necessary for efficient citizenship as the ability to read and write.""统计思维有一天会像读写能力一样,成为高效公民的必备技能。"— H.G. Wells, Writer & Futurist (作家、未来学家)
当结果是"是/否"时,如何建模?
本节目标
- 理解 Logit 回归的应用场景
- 使用 statsmodels 进行 Logit 回归
- 解读 Logit 回归的系数和边际效应
- 与 Stata/R 进行对比
什么时候用 Logit 回归?
当因变量是二元变量(0/1)时,不能用 OLS,要用 Logit 或 Probit。
典型应用场景
| 研究问题 | 因变量 | 自变量示例 |
|---|---|---|
| 是否上大学 | 上大学=1, 否=0 | 家庭收入、父母教育 |
| 是否找到工作 | 就业=1, 失业=0 | 教育、经验、性别 |
| 是否违约 | 违约=1, 否=0 | 信用分数、收入、负债率 |
| 是否投票 | 投票=1, 否=0 | 年龄、教育、收入 |
| 是否生病 | 生病=1, 健康=0 | 年龄、BMI、吸烟史 |
案例:影响大学录取的因素
研究问题:什么因素影响学生是否被大学录取?
python
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.formula.api import logit
# 生成模拟数据:500 个申请者
np.random.seed(2024)
n = 500
data = pd.DataFrame({
'gpa': np.random.uniform(2.0, 4.0, n), # GPA 成绩
'sat': np.random.randint(800, 1600, n), # SAT 分数
'extracurricular': np.random.randint(0, 10, n), # 课外活动数量
'income': np.random.uniform(20, 150, n) # 家庭收入(千元)
})
# 构造录取概率(logit 模型)
z = -8 + 1.5 * data['gpa'] + 0.003 * data['sat'] + 0.1 * data['extracurricular'] + 0.01 * data['income']
prob = 1 / (1 + np.exp(-z))
data['admitted'] = (np.random.uniform(0, 1, n) < prob).astype(int)
print(data.head(10))
print(f"\n录取率: {data['admitted'].mean():.2%}")Logit 回归:Python 实现
方法 1:使用 sm.Logit()(推荐)
python
import statsmodels.api as sm
# 准备数据
X = sm.add_constant(data[['gpa', 'sat', 'extracurricular', 'income']])
y = data['admitted']
# 拟合 Logit 模型
logit_model = sm.Logit(y, X).fit()
# 查看结果
print(logit_model.summary())输出:
Logit Regression Results
==============================================================================
Dep. Variable: admitted No. Observations: 500
Model: Logit Df Residuals: 495
Method: MLE Df Model: 4
Date: ... Pseudo R-squ.: 0.524
Time: ... Log-Likelihood: -178.23
converged: True LL-Null: -374.56
Covariance Type: nonrobust LLR p-value: 1.23e-82
==================================================================================
coef std err z P>|z| [0.025 0.975]
----------------------------------------------------------------------------------
const -7.823 0.987 -7.926 0.000 -9.758 -5.888
gpa 1.485 0.198 7.500 0.000 1.097 1.873
sat 0.0029 0.001 3.625 0.000 0.001 0.005
extracurricular 0.098 0.034 2.882 0.004 0.031 0.165
income 0.0095 0.004 2.375 0.018 0.002 0.017
==================================================================================解读:
- GPA:系数 1.485,p < 0.001,显著正向影响录取
- SAT:系数 0.0029,p < 0.001,显著正向影响
- 课外活动:系数 0.098,p < 0.01,显著正向影响
- 家庭收入:系数 0.0095,p < 0.05,边际显著
- Pseudo R² = 0.524,模型拟合较好
方法 2:使用公式接口(类似 R)
python
from statsmodels.formula.api import logit
# 使用 R 风格的公式
logit_model = logit('admitted ~ gpa + sat + extracurricular + income', data=data).fit()
print(logit_model.summary())三大语言对比
Python (statsmodels)
python
import statsmodels.api as sm
# 准备数据
X = sm.add_constant(data[['gpa', 'sat', 'extracurricular', 'income']])
y = data['admitted']
# 拟合模型
model = sm.Logit(y, X).fit()
# 查看结果
print(model.summary())
# 预测概率
prob = model.predict(X)Stata
stata
* 读取数据
use admission_data.dta, clear
* Logit 回归
logit admitted gpa sat extracurricular income
* 查看边际效应
margins, dydx(*)
* 预测概率
predict probR
r
# 读取数据
data <- read.csv("admission_data.csv")
# Logit 回归
model <- glm(admitted ~ gpa + sat + extracurricular + income,
data = data,
family = binomial(link = "logit"))
# 查看结果
summary(model)
# 预测概率
prob <- predict(model, type = "response")解读 Logit 系数
系数的含义
Logit 回归的系数不是边际效应,而是对数几率(log-odds)的变化。
python
# GPA 系数为 1.485,意味着:
# GPA 每增加 1,log-odds 增加 1.485
# 等价于:几率(odds)乘以 exp(1.485) = 4.42
print(f"GPA 系数: {logit_model.params['gpa']:.3f}")
print(f"Odds Ratio: {np.exp(logit_model.params['gpa']):.3f}")解读:
- GPA 每提高 1 分,录取几率(odds)增加到原来的 4.42 倍
- 这是一个非常强的效应!
获取 Odds Ratios
python
# 所有变量的 Odds Ratios
odds_ratios = np.exp(logit_model.params)
print(odds_ratios)
'''
const 0.000
gpa 4.416
sat 1.003
extracurricular 1.103
income 1.010
'''解读:
- GPA:每增加 1,录取几率增加到 4.42 倍
- SAT:每增加 1 分,录取几率增加到 1.003 倍(效应较小)
- 课外活动:每多 1 个,录取几率增加到 1.10 倍
- 收入:每多 1 千元,录取几率增加到 1.01 倍
边际效应(Marginal Effects)
边际效应更容易解释:自变量增加 1 单位,录取概率增加多少?
python
# 计算平均边际效应(Average Marginal Effects, AME)
marginal_effects = logit_model.get_margeff()
print(marginal_effects.summary())输出:
Logit Marginal Effects
=====================================
Dep. Variable: admitted
Method: dydx
At: overall
=====================================
dy/dx std err z P>|z|
-------------------------------------------------------------
gpa 0.298 0.039 7.641 0.000
sat 0.001 0.000 3.702 0.000
extracurricular 0.020 0.007 2.857 0.004
income 0.002 0.001 2.344 0.019
-------------------------------------------------------------解读:
- GPA:每增加 1 分,录取概率平均增加 29.8 个百分点
- SAT:每增加 1 分,录取概率增加 0.1 个百分点
- 课外活动:每多 1 个,录取概率增加 2.0 个百分点
- 收入:每多 1 千元,录取概率增加 0.2 个百分点
预测
预测概率
python
# 预测所有样本的录取概率
predicted_prob = logit_model.predict(X)
print(predicted_prob[:10])
# 添加到数据框
data['predicted_prob'] = predicted_prob
print(data[['gpa', 'sat', 'admitted', 'predicted_prob']].head())预测新学生
python
# 一个新申请者:GPA=3.5, SAT=1200, 5个课外活动, 收入80k
new_student = pd.DataFrame({
'const': [1],
'gpa': [3.5],
'sat': [1200],
'extracurricular': [5],
'income': [80]
})
prob = logit_model.predict(new_student)
print(f"录取概率: {prob.values[0]:.2%}")
# 输出: 录取概率: 68.34%模型评估
混淆矩阵
python
from sklearn.metrics import confusion_matrix, classification_report
# 预测类别(概率 > 0.5 则预测为录取)
predicted_class = (predicted_prob > 0.5).astype(int)
# 混淆矩阵
cm = confusion_matrix(data['admitted'], predicted_class)
print("混淆矩阵:")
print(cm)
'''
[[210 35]
[ 28 227]]
'''
print("\n分类报告:")
print(classification_report(data['admitted'], predicted_class))输出:
precision recall f1-score support
0 0.88 0.86 0.87 245
1 0.87 0.89 0.88 255
accuracy 0.87 500️ Logit vs OLS
为什么不能用 OLS?
python
# 错误:对二元变量用 OLS
X = sm.add_constant(data[['gpa', 'sat', 'extracurricular', 'income']])
y = data['admitted']
ols_model = sm.OLS(y, X).fit()
ols_pred = ols_model.predict(X)
# 问题:预测值可能超出 [0, 1] 范围
print(f"OLS 预测的最小值: {ols_pred.min():.3f}") # 可能 < 0
print(f"OLS 预测的最大值: {ols_pred.max():.3f}") # 可能 > 1Logit 的优势
python
# 正确:Logit 保证预测概率在 [0, 1]
logit_pred = logit_model.predict(X)
print(f"Logit 预测的最小值: {logit_pred.min():.3f}") # 始终 ≥ 0
print(f"Logit 预测的最大值: {logit_pred.max():.3f}") # 始终 ≤ 1实践练习
完整代码:研究就业影响因素
python
import pandas as pd
import numpy as np
import statsmodels.api as sm
# 生成数据
np.random.seed(42)
n = 800
data = pd.DataFrame({
'employed': np.random.choice([0, 1], n, p=[0.3, 0.7]),
'education': np.random.randint(9, 22, n),
'age': np.random.randint(22, 60, n),
'female': np.random.choice([0, 1], n)
})
# 构造更真实的就业概率
z = -3 + 0.2 * data['education'] + 0.02 * data['age'] - 0.3 * data['female']
prob = 1 / (1 + np.exp(-z))
data['employed'] = (np.random.uniform(0, 1, n) < prob).astype(int)
# Logit 回归
X = sm.add_constant(data[['education', 'age', 'female']])
y = data['employed']
model = sm.Logit(y, X).fit()
print(model.summary())
# 边际效应
print("\n边际效应:")
print(model.get_margeff().summary())
# Odds Ratios
print("\nOdds Ratios:")
print(np.exp(model.params))关键要点
| 内容 | OLS 回归 | Logit 回归 |
|---|---|---|
| 因变量 | 连续变量 | 二元变量(0/1) |
| 预测值范围 | (-∞, +∞) | [0, 1] |
| 系数解释 | 边际效应 | Log-odds |
| 更直观的解释 | 系数本身 | Odds Ratios 或边际效应 |
| Python 命令 | sm.OLS() | sm.Logit() |
下一步
- 文章 04:
summary_col()- 优雅地展示多个模型对比
** 你已经掌握了 Python 中的 Logit 回归!**