GLM原理与代码实例：从理论到实践的全面解析

2025年11月26日 互联网

GLM原理与代码实例讲解

一、广义线性模型（GLM）的数学基础

1.1 传统线性模型的局限性

经典线性回归模型（OLS）假设因变量服从正态分布，且预测值与误差项独立。但在实际场景中，如二分类问题（是否购买）、计数数据（日访问量）等非连续变量，OLS的假设无法满足。例如预测用户点击率时，因变量取值范围为[0,1]，直接使用线性回归会导致预测值超出有效范围。

1.2 GLM的三大核心组件

GLM通过随机分量、系统分量和连接函数三要素突破传统限制：

  • 随机分量：定义因变量的概率分布，支持二项分布（逻辑回归）、泊松分布（计数数据）、伽马分布（正偏态连续数据）等
  • 系统分量：保持线性预测形式 $X^T\beta$，其中$X$为特征矩阵，$\beta$为参数向量
  • 连接函数：建立系统分量与随机分量期望的映射关系，常用形式包括：
    • Logit连接（逻辑回归）：$g(\mu)=\ln(\frac{\mu}{1-\mu})$
    • Probit连接（分位数回归）：$g(\mu)=\Phi^{-1}(\mu)$
    • 对数连接（泊松回归）：$g(\mu)=\ln(\mu)$

1.3 最大似然估计的优化过程

GLM采用迭代加权最小二乘法（IRLS）进行参数估计。以逻辑回归为例，似然函数为：
<br>L(β)=∏<em>i=1n[πiyi(1−πi)1−yi]<br></em><br>L(\beta)=\prod<em>{i=1}^n [\pi_i^{y_i}(1-\pi_i)^{1-y_i}]<br></em>
其中$\pi_i=g^{-1}(X_i^T\beta)$。通过牛顿-拉夫森迭代更新参数：
<br>β(k+1)=β(k)+(XTWX)−1XT(y−μ)<br><br>\beta^{(k+1)}=\beta^{(k)}+(X^TWX)^{-1}X^T(y-\mu)<br>
其中$W$为对角权重矩阵，元素$w{ii}=\pi_i(1-\pi_i)$。

二、Python实现全流程解析

2.1 环境准备与数据加载

  1. import numpy as np
  2. import pandas as pd
  3. import statsmodels.api as sm
  4. from sklearn.datasets import make_classification
  5. from sklearn.model_selection import train_test_split
  7. # 生成模拟二分类数据
  8. X, y = make_classification(n_samples=1000, n_features=5, 
  9.                           n_classes=2, random_state=42)
  10. X = pd.DataFrame(X, columns=[f'x{i}' for i in range(5)])
  11. X['intercept'] = 1  # 添加截距项
  13. # 划分训练集/测试集
  14. X_train, X_test, y_train, y_test = train_test_split(
  15.     X, y, test_size=0.3, random_state=42)

2.2 逻辑回归模型构建

  1. # 使用statsmodels构建GLM
  2. logit_model = sm.GLM(
  3.     endog=y_train,          # 因变量
  4.     exog=X_train[['x0', 'x1', 'x2', 'x3', 'x4', 'intercept']],  # 自变量
  5.     family=sm.families.Binomial()  # 指定二项分布
  6. ).fit()
  8. # 输出模型摘要
  9. print(logit_model.summary())

输出结果包含：

  • 参数估计值（coef）及其显著性（P>|z|）
  • 对数似然值（-2LL）
  • AIC/BIC模型选择指标
  • 伪R方（McFadden’s R²）

2.3 模型诊断与结果解释

2.3.1 参数显著性检验

  1. # 获取参数估计与置信区间
  2. params = logit_model.params
  3. conf_int = logit_model.conf_int()
  4. result_df = pd.DataFrame({
  5.     'Coefficient': params,
  6.     '95% Lower': conf_int[:,0],
  7.     '95% Upper': conf_int[:,1]
  8. })
  9. print(result_df)

示例输出：

  1.           Coefficient  95% Lower  95% Upper
  2. x0            0.8234     0.6123     1.0345
  3. x1           -0.4567    -0.7210    -0.1924
  4. intercept     0.1234    -0.0456     0.2924

解释：x0每增加1个单位，对数几率比增加0.8234（Z=6.78, P<0.001）

2.3.2 预测与阈值优化

  1. # 生成预测概率
  2. y_pred_prob = logit_model.predict(X_test[['x0','x1','x2','x3','x4','intercept']])
  4. # 寻找最优分类阈值（通过ROC曲线）
  5. from sklearn.metrics import roc_curve, auc
  6. fpr, tpr, thresholds = roc_curve(y_test, y_pred_prob)
  7. optimal_idx = np.argmax(tpr - fpr)
  8. optimal_threshold = thresholds[optimal_idx]
  10. # 应用最优阈值进行分类
  11. y_pred = (y_pred_prob > optimal_threshold).astype(int)

三、工程化实践建议

3.1 特征工程关键点

  • 分箱处理：对连续变量进行分箱（如年龄分为5组）可提升模型可解释性
  • 交互项引入：通过X['x0_x1'] = X['x0'] * X['x1']添加特征交互
  • 标准化：对正态分布特征进行Z-score标准化（$\frac{x-\mu}{\sigma}$）

3.2 模型调优策略

  1. 正则化实现：

    1. from sklearn.linear_model import LogisticRegression
    2. # L2正则化示例
    3. lr_model = LogisticRegression(
    4.  penalty='l2', 
    5.  C=0.1,  # 1/λ，值越小正则化越强
    6.  solver='lbfgs'
    7. ).fit(X_train.drop('intercept', axis=1), y_train)

  2. 类别不平衡处理：

    1. # 使用statsmodels的权重参数
    2. sample_weights = np.where(y_train==1, 0.7, 0.3)  # 正样本权重更高
    3. weighted_model = sm.GLM(
    4.  y_train, X_train[['x0','x1','x2','x3','x4','intercept']],
    5.  family=sm.families.Binomial(),
    6.  freq_weights=sample_weights  # 频率权重
    7. ).fit()

3.3 生产环境部署要点

  • 模型序列化：使用pickle或joblib保存训练好的模型 
    1. import joblib
    2. joblib.dump(logit_model, 'glm_model.pkl')
    3. loaded_model = joblib.load('glm_model.pkl')
  • API封装：通过FastAPI创建预测服务
    ```python
    from fastapi import FastAPI
    import pandas as pd

app = FastAPI()

@app.post(“/predict”)
async def predict(features: dict):
df = pd.DataFrame([features])
df[‘intercept’] = 1
prob = loaded_model.predict(df[[‘x0’,’x1’,’x2’,’x3’,’x4’,’intercept’]])[0]
return {“probability”: float(prob)}

  2. ## 四、常见问题解决方案
  3. ### 4.1 收敛失败处理
  4. 当遇到`ConvergenceWarning`时，可尝试：
  5. 1. 增加最大迭代次数：`max_iter=1000`
  6. 2. 调整优化算法：`solver='newton'`或`'bfgs'`
  7. 3. 特征缩放：对大范围变量进行对数变换
  9. ### 4.2 完全分离问题
  10. 当某个特征能完美预测类别时（如所有正样本的x1>5），会导致参数估计无限大。解决方案：
  11. - 使用Firth的惩罚似然法
  12. - 合并相关特征或移除问题特征
  13. - 增加正则化强度
  15. ### 4.3 多重共线性诊断
  16. ```python
  17. from statsmodels.stats.outliers_influence import variance_inflation_factor
  19. # 计算VIF
  20. vif_data = pd.DataFrame()
  21. vif_data["feature"] = X_train.columns[:-1]  # 排除intercept
  22. vif_data["VIF"] = [variance_inflation_factor(
  23.     X_train.values, i) for i in range(len(X_train.columns[:-1]))]
  25. print(vif_data[vif_data["VIF"] > 5])  # 输出VIF>5的特征

五、扩展应用场景

5.1 泊松回归实现

  1. # 生成计数数据
  2. np.random.seed(42)
  3. counts = np.random.poisson(lam=np.exp(0.5*X['x0'] + 0.3*X['x1']), size=1000)
  5. # 构建泊松GLM
  6. poisson_model = sm.GLM(
  7.     counts, X[['x0','x1','intercept']],
  8.     family=sm.families.Poisson()
  9. ).fit()
  11. print(poisson_model.summary())

5.2 Gamma回归用于偏态连续数据

  1. # 生成右偏连续数据
  2. gamma_data = np.random.gamma(shape=2, scale=np.exp(0.2*X['x0']), size=1000)
  4. # 构建Gamma GLM（对数连接）
  5. gamma_model = sm.GLM(
  6.     gamma_data, X[['x0','intercept']],
  7.     family=sm.families.Gamma(sm.families.links.log())
  8. ).fit()
  10. print(gamma_model.summary())

六、总结与最佳实践

  1. 模型选择流程：

    • 确认因变量类型（连续/二分类/计数）
    • 检查分布假设（正态性检验、Q-Q图）
    • 选择对应的分布族和连接函数

  2. 解释性优先场景：

    • 使用statsmodels获取详细统计量
    • 生成部分依赖图（PDP）解释特征影响

  3. 预测性能优先场景：

    • 使用scikit-learn的LogisticRegression
    • 结合网格搜索进行超参数调优

  4. 持续监控指标：

    • 预测准确率/AUC（分类）
    • 残差分布（回归）
    • 参数稳定性（时间序列数据）

通过系统掌握GLM的原理与实现细节，数据分析师能够更精准地构建符合业务需求的统计模型，同时工程师可实现高效可靠的预测系统部署。建议读者结合实际业务数据，按照本文提供的代码框架进行实践验证，逐步构建对GLM的深度理解。

最新文章