如何在Python中计算Logit模型代表性值边际效应（仿Stata margins）

在Python的statsmodels中没有直接对应Stata margins命令的现成函数，但可以通过手动构造目标数据集、结合Logit模型的预测与边际效应公式来实现。以下是具体步骤和代码示例：

1. 拟合基础Logit模型 #

首先导入依赖库并拟合Logit模型，这里用模拟数据演示，你可以替换成自己的真实数据集：

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import Logit

# 构造模拟数据集
np.random.seed(42)
n = 500
age = np.random.randint(20, 60, n)
income = np.random.normal(50000, 10000, n)
education = np.random.randint(12, 20, n)
# 生成二元因变量
y = np.random.binomial(
    1, 
    sm.tools.tools.predict(
        sm.Logit(endog=None, exog=sm.add_constant(np.column_stack((age, income, education)))),
        params=[-10, 0.05, 0.0001, 0.2]
    ), 
    n
)
data = pd.DataFrame({'y': y, 'age': age, 'income': income, 'education': education})

# 拟合Logit模型
X = sm.add_constant(data[['age', 'income', 'education']])
model = Logit(data['y'], X)
results = model.fit()
print(results.summary())

2. 构造目标计算数据集 #

模拟Stata at(age=(30(1)35))的逻辑：固定其他变量为代表性值（这里用均值，可替换为中位数/众数），生成age从30到35的序列：

# 获取其他变量的均值（可替换为median()取中位数）
other_vars_mean = data[['income', 'education']].mean().to_dict()
# 生成age序列：30到35（步长1）
age_range = np.arange(30, 36)
# 构造margins计算用数据集
margins_data = pd.DataFrame({
    'age': age_range,
    'income': [other_vars_mean['income']] * len(age_range),
    'education': [other_vars_mean['education']] * len(age_range)
})
# 添加模型所需的常数项
margins_data = sm.add_constant(margins_data)

3. 计算边际效应 #

有两种常用方法，分别对应解析解和数值近似解：

方法一：解析法（精确导数） #

Logit模型中，连续自变量的边际效应公式为：

边际效应 = 预测概率p * (1 - p) * 自变量系数

用该公式直接计算：

# 计算每个age对应的预测概率
margins_data['pred_prob'] = results.predict(margins_data)
# 获取age变量的系数
beta_age = results.params['age']
# 计算解析法边际效应
margins_data['mfx_analytic'] = margins_data['pred_prob'] * (1 - margins_data['pred_prob']) * beta_age

方法二：数值法（离散变化近似） #

模拟Stata中对连续变量的离散边际效应计算（即age增加1单位时，预测概率的变化量）：

# 构造age+1的数据集
margins_data_plus1 = margins_data.copy()
margins_data_plus1['age'] += 1
# 计算age+1时的预测概率
margins_data_plus1['pred_prob_plus1'] = results.predict(margins_data_plus1)
# 数值边际效应 = 概率差值
margins_data['mfx_numeric'] = margins_data_plus1['pred_prob_plus1'] - margins_data['pred_prob']

4. 查看结果 #

输出最终的age、预测概率和两种方法的边际效应：

print(margins_data[['age', 'pred_prob', 'mfx_analytic', 'mfx_numeric']])

说明 #

• 两种方法的结果非常接近，解析法是精确的导数，数值法是离散单位变化的近似值；
• 如果需要固定其他变量为非均值的代表性值（如中位数），只需将mean()替换为median()或其他统计量；
• 若自变量是离散变量，只需调整数值法的步长（如分类变量取不同类别）即可。

内容的提问来源于stack exchange，提问作者dataLearner