如何在Statsmodels逻辑回归中获取实际预测边际而非dy/dx？

问题场景 #

你使用smf.logit拟合了如下逻辑回归模型：

model_5a_1 = smf.logit('''fluid ~ C(examq3_n, Treatment(reference = 2.0)) + C(pmhq3_n) + C(fluidq3_n) + C(mapq3_n, Treatment(reference = 3.0)) + 
                  C(examq6_n, Treatment(reference = 2.0)) + C(pmhq6_n) + C(fluidq6_n) + C(mapq6_n, Treatment(reference = 3.0)) +
                  + C(case, Treatment(reference = 2))''',
                 data = case1_2_vars).fit()
print(model_5a_1.summary())

尝试用get_margeff()得到的是dy/dx形式的边际效应，但你需要的是预测边际（即给定自变量取值时的平均预测概率），类似Stata的margin命令；使用get_prediction()时又出现PatsyError: categorical data cannot be >1-dimensional错误。

解决方案 #

要实现类似Stata margin的预测边际计算，核心是先构造符合要求的参考数据集，再用模型计算预测概率，具体步骤如下：

1. 构造参考数据集 #

创建一个数据集，将目标分类变量遍历所有水平，其他变量固定在均值（或你需要的参考值，比如中位数）：

import pandas as pd

# 获取模型用到的自变量列表
var_names = model_5a_1.model.exog_names
# 生成基础数据：非分类变量取均值，分类变量先保留一个基准水平
base_data = case1_2_vars[var_names].mean().to_frame().T

# 定义所有分类变量
categorical_vars = ['examq3_n', 'pmhq3_n', 'fluidq3_n', 'mapq3_n', 'examq6_n', 'pmhq6_n', 'fluidq6_n', 'mapq6_n', 'case']

# 生成包含所有分类水平组合的边际数据集
margins_data = pd.DataFrame()
for var in categorical_vars:
    # 获取当前变量的所有唯一取值
    levels = case1_2_vars[var].unique()
    # 复制基础数据并替换当前变量的取值
    temp_data = pd.concat([base_data]*len(levels), ignore_index=True)
    temp_data[var] = levels
    margins_data = pd.concat([margins_data, temp_data], ignore_index=True)

# 去重避免重复组合
margins_data = margins_data.drop_duplicates().reset_index(drop=True)

2. 计算预测边际概率 #

用构造好的数据集调用get_prediction()，即可得到带置信区间的预测边际：

# 计算预测结果及置信区间
predictions = model_5a_1.get_prediction(margins_data)
pred_summary = predictions.summary_frame()

# 合并变量取值与预测结果
final_result = pd.concat([margins_data, pred_summary[['mean', 'mean_ci_lower', 'mean_ci_upper']]], axis=1)

# 示例：按case变量分组展示预测边际
print(final_result.groupby('case')[['mean', 'mean_ci_lower', 'mean_ci_upper']].first())

3. 简化版：针对单个变量计算边际 #

如果只需要计算某一个分类变量的预测边际，可以用自定义函数简化：

def get_single_var_margins(model, var_name, data):
    levels = data[var_name].unique()
    base_data = data.drop(var_name, axis=1).mean().to_frame().T
    result_list = []
    
    for level in levels:
        temp_data = base_data.copy()
        temp_data[var_name] = level
        pred = model.get_prediction(temp_data).summary_frame()
        result_list.append({
            var_name: level,
            'pred_prob': pred['mean'].iloc[0],
            'lower_ci': pred['mean_ci_lower'].iloc[0],
            'upper_ci': pred['mean_ci_upper'].iloc[0]
        })
    return pd.DataFrame(result_list)

# 计算case变量的预测边际
case_margins = get_single_var_margins(model_5a_1, 'case', case1_2_vars)
print(case_margins)

错误原因说明 #

之前调用get_prediction()报错，是因为直接传入的原数据或格式错误的数据集，包含多维分类变量结构，Patsy无法正确解析；而构造的参考数据集是单条基准记录复制替换分类水平，格式符合模型的输入要求。

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