从泊松回归到伽马回归：用Python statsmodels库实战GLM处理非正态数据
当你的数据拒绝服从正态分布时，传统线性回归就像试图用螺丝刀敲钉子——不仅效率低下，还可能损坏工具。在真实业务场景中，我们常遇到计数数据（如网站点击量）、持续正数（如保险理赔金额）或比例数据（如转化率），这些数据往往呈现明显的右偏或离散特征。本文将带你用statsmodels库中的GLM模块，像专业数据科学家一样处理这些"叛逆"数据。
1. 理解GLM的核心武器库
广义线性模型(GLMs)是线性回归的瑞士军刀扩展版，通过三个关键组件解决非正态数据问题：
- 分布族(Family)：打破正态分布限制，支持泊松、伽马、负二项等分布
- 链接函数(Link Function)：建立线性预测与响应变量的非线性关系
- 方差函数(Variance Function)：描述均值与方差的关联方式
关键选择矩阵：
| 数据类型 | 典型分布族 | 常用链接函数 | 典型应用场景 |
|---|---|---|---|
| 计数数据 | Poisson/NegativeBinomial | log | 网站点击量分析 |
| 连续正数 | Gamma/InverseGaussian | log/inverse | 保险理赔建模 |
| 二元分类 | Binomial | logit/probit | 用户转化预测 |
| 比例数据 | Binomial | logit | 广告点击率分析 |
注意：选择链接函数时，需确保其能将线性预测值映射到响应变量的自然取值范围内。例如，对数链接确保伽马回归的输出保持正值。
2. 数据诊断：识别你的数据DNA
在构建模型前，我们需要像法医一样检验数据的分布特征：
import seaborn as sns import matplotlib.pyplot as plt from scipy import stats def diagnose_distribution(data, var): fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) # 分布直方图与密度曲线 sns.histplot(data[var], kde=True, ax=ax1) ax1.set_title(f'Distribution of {var}') # Q-Q图检验正态性 stats.probplot(data[var], plot=ax2) ax2.set_title(f'Q-Q Plot of {var}') plt.tight_layout() return fig # 示例：诊断保险理赔金额分布 claims = pd.read_csv('insurance_claims.csv') diagnose_distribution(claims, 'claim_amount')常见非正态数据特征及应对策略：
过度离散(Overdispersion)：方差明显大于均值时，泊松回归会低估标准误
- 解决方案：改用负二项回归
- 诊断代码：
print(f"离散系数：{data[var].var()/data[var].mean():.2f}")
零膨胀(Zero-inflation)：数据中零值比例异常高
- 解决方案：零膨胀泊松模型或 hurdle 模型
- 诊断代码：
print(f"零值占比：{(data[var]==0).mean()*100:.1f}%")
3. 泊松回归实战：点击量预测案例
当建模事件发生次数时（如每小时网站点击量），泊松回归是首选武器。以下完整流程展示如何用statsmodels实现：
import statsmodels.api as sm import statsmodels.formula.api as smf # 准备数据 clicks_data = pd.read_csv('website_clicks.csv') formula = 'clicks ~ time_of_day + page_type + user_segment' # 模型拟合 poisson_model = smf.glm( formula=formula, data=clicks_data, family=sm.families.Poisson(link=sm.families.links.log()) ).fit() # 过离散诊断 print(f"Pearson卡方统计量：{poisson_model.pearson_chi2/poisson_model.df_resid:.2f}") # 若存在过离散(>1.5)，改用负二项回归 if poisson_model.pearson_chi2/poisson_model.df_resid > 1.5: nb_model = smf.glm( formula=formula, data=clicks_data, family=sm.families.NegativeBinomial() ).fit() print(nb_model.summary()) else: print(poisson_model.summary())结果解释要点：
系数需按链接函数反向转换解释：
- 对数链接下：
exp(coefficient)表示倍数变化 - 示例：
time_of_day[T.Night] = 0.5表示夜间点击量是基准时段的exp(0.5)≈1.65倍
- 对数链接下：
模型评估关键指标：
- AIC/BIC：用于模型比较，值越小越好
- 残差分析：
model.resid_deviance检查模式结构 - 预测可视化：绘制实际值 vs 拟合值散点图
4. 伽马回归精解：处理右偏连续数据
保险理赔金额、服务器响应时间等持续正数常呈现右偏分布，伽马回归能有效处理这类数据：
# 伽马回归建模保险理赔金额 gamma_model = smf.glm( 'claim_amount ~ age + vehicle_type + claim_history', data=claims, family=sm.families.Gamma(link=sm.families.links.log()) ).fit() # 模型诊断图 def plot_gamma_diagnostics(model): fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) # 拟合值 vs 实际值 sns.scatterplot(x=model.fittedvalues, y=model.model.endog, ax=ax1) ax1.plot([0, max(model.fittedvalues)], [0, max(model.fittedvalues)], 'r--') ax1.set_xlabel('Predicted') ax1.set_ylabel('Actual') # 残差QQ图 stats.probplot(model.resid_pearson, plot=ax2) return fig plot_gamma_diagnostics(gamma_model)伽马回归调优技巧：
链接函数选择：
- log链接：默认选择，保证预测值为正
- inverse链接：当效应呈反比关系时更合适
形状参数估计：
# 估计伽马形状参数(alpha) alpha = len(gamma_model.params)/gamma_model.deviance print(f"Estimated shape parameter: {alpha:.2f}")处理极端值：
# 使用稳健标准误 gamma_model_robust = smf.glm( formula, data=claims, family=sm.families.Gamma(), var_weights=1/claims['claim_amount'] # 逆方差加权 ).fit()
5. 模型比较与生产部署
在实际业务中，我们常需要比较不同GLM配置：
# 定义候选模型 models = { 'Poisson': sm.families.Poisson(), 'NegativeBinomial': sm.families.NegativeBinomial(), 'Gamma': sm.families.Gamma(), 'Tweedie': sm.families.Tweedie(var_power=1.5) # 复合泊松-伽马 } # 自动化模型比较 results = [] for name, family in models.items(): try: model = smf.glm(formula, data, family=family).fit() results.append({ 'Model': name, 'AIC': model.aic, 'BIC': model.bic, 'Deviance': model.deviance, 'Params': len(model.params) }) except: continue pd.DataFrame(results).sort_values('AIC')生产部署检查清单：
性能优化：
# 使用稀疏矩阵处理高维分类变量 from patsy import dmatrices y, X = dmatrices(formula, data, return_type='sparse') sparse_model = sm.GLM(y, X.tocsc(), family=sm.families.Poisson()).fit()模型持久化：
import joblib joblib.dump(model, 'glm_model.pkl') # 加载时重新附加family类 loaded_model = joblib.load('glm_model.pkl') loaded_model.family = sm.families.Poisson() # 必须与保存时一致实时预测API示例：
from flask import Flask, request, jsonify app = Flask(__name__) model = joblib.load('glm_model.pkl') @app.route('/predict', methods=['POST']) def predict(): data = request.json X_new = pd.DataFrame([data]) pred = model.predict(X_new) return jsonify({'prediction': float(pred[0])})
在实际电商数据分析项目中，我发现伽马回归预测物流时间时，加入log(订单量)作为偏移量(offset)能显著提升模型效果——这相当于对"单位商品处理时间"建模，比直接预测总时间更具业务解释性。
