基于泊松回归的计数数据建模：链接函数、极大似然估计与优化过程
1. 题目背景
在实际问题中，我们经常遇到因变量是计数（Count）类型的数据，例如：
- 一个小时内到达某个网站的访客数。
- 某条道路上一天内发生的交通事故数。
- 一本书中某章节的单词数。
这类数据的特点是其取值是非负整数。传统的线性回归（假设误差服从高斯分布）不适用于此，因为其预测值可能是负数，且对计数数据的离散性和方差特性建模不佳。泊松回归是专门为建模计数型响应变量而设计的回归模型。
2. 模型描述
泊松回归是广义线性模型（GLM）的一个特例。其核心包含三个部分：
A. 随机成分（响应分布）
假设响应变量 \(Y\) 服从泊松分布：
\[Y_i | \mathbf{x}_i \sim \text{Poisson}(\lambda_i) \]
其中，\(\lambda_i > 0\) 是给定特征 \(\mathbf{x}_i\) 时 \(Y_i\) 的期望（均值）。泊松分布的概率质量函数为：
\[P(Y_i = y_i) = \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!} \]
泊松分布的一个重要性质是期望等于方差：\(E[Y_i] = \text{Var}(Y_i) = \lambda_i\)。
B. 系统成分（线性预测器）
与线性回归类似，我们构造一个关于特征 \(\mathbf{x}_i = (1, x_{i1}, x_{i2}, \dots, x_{ip})^T\) 的线性组合：
\[\eta_i = \mathbf{x}_i^T \boldsymbol{\beta} = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip} \]
其中，\(\boldsymbol{\beta} = (\beta_0, \beta_1, \dots, \beta_p)^T\) 是待估计的系数向量。
C. 链接函数（Link Function）
为了确保预测的均值 \(\lambda_i\) 是非负的，我们不能直接将 \(\lambda_i\) 建模为线性预测器（因为 \(\eta_i\) 可以是任意实数）。因此，我们使用一个链接函数 \(g(\cdot)\)，将均值 \(\lambda_i\) 映射到线性预测器 \(\eta_i\)。
泊松回归使用对数链接函数：
\[\eta_i = g(\lambda_i) = \log(\lambda_i) \]
反之：
\[\lambda_i = e^{\eta_i} = \exp(\mathbf{x}_i^T \boldsymbol{\beta}) \]
这个指数变换确保了 \(\lambda_i > 0\)。
综合模型：将三部分结合，泊松回归模型表示为：
\[Y_i | \mathbf{x}_i \sim \text{Poisson}(\lambda_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta})) \]
这意味着，特征 \(x_j\) 增加一个单位，响应变量 \(Y\) 的期望值 \(\lambda\) 会变为原来的 \(e^{\beta_j}\) 倍（发生率比）。
3. 参数估计：极大似然估计（MLE）
给定 \(n\) 个独立同分布的观测数据 \(\{ (\mathbf{x}_i, y_i) \}_{i=1}^n\)，我们需要估计参数 \(\boldsymbol{\beta}\)。
A. 构造似然函数
由于观测独立，联合概率（似然函数）是所有观测概率的乘积：
\[L(\boldsymbol{\beta}) = \prod_{i=1}^{n} P(Y_i = y_i) = \prod_{i=1}^{n} \frac{e^{-\lambda_i} \lambda_i^{y_i}}{y_i!} \]
其中 \(\lambda_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta})\)。
B. 构造对数似然函数
为简化计算，通常取对数：
\[\ell(\boldsymbol{\beta}) = \log L(\boldsymbol{\beta}) = \sum_{i=1}^{n} [-\lambda_i + y_i \log(\lambda_i) - \log(y_i!)] \]
代入 \(\lambda_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta})\)：
\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} [y_i (\mathbf{x}_i^T \boldsymbol{\beta}) - \exp(\mathbf{x}_i^T \boldsymbol{\beta}) - \log(y_i!)] \]
由于 \(\log(y_i!)\) 与 \(\boldsymbol{\beta}\) 无关，在优化中可以忽略。因此，我们需要最大化的目标函数是：
\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} [y_i (\mathbf{x}_i^T \boldsymbol{\beta}) - \exp(\mathbf{x}_i^T \boldsymbol{\beta})] \]
C. 求解极大似然估计
我们通过求解对数似然函数对参数 \(\boldsymbol{\beta}\) 的导数（梯度）为零的方程来估计参数。首先计算梯度：
\[\frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} = \sum_{i=1}^{n} [y_i \mathbf{x}_i - \exp(\mathbf{x}_i^T \boldsymbol{\beta}) \mathbf{x}_i] = \sum_{i=1}^{n} (y_i - \lambda_i) \mathbf{x}_i \]
令梯度为零，得到得分方程（Score Equations）：
\[\sum_{i=1}^{n} (y_i - \exp(\mathbf{x}_i^T \boldsymbol{\beta})) \mathbf{x}_i = 0 \]
这是一个关于 \(\boldsymbol{\beta}\) 的非线性方程组，没有解析解。因此，我们必须使用数值优化方法。
4. 优化求解：迭代重加权最小二乘法（IRLS）
对于广义线性模型，最常用的优化算法是牛顿-拉弗森法（Newton-Raphson Method）的一种特殊形式——迭代重加权最小二乘法。
A. 推导过程
牛顿-拉弗森法的迭代更新公式为：
\[\boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - [H^{(t)}]^{-1} \cdot \mathbf{g}^{(t)} \]
其中，\(\mathbf{g}^{(t)}\) 是梯度向量，\(H^{(t)}\) 是海森矩阵（Hessian Matrix，即二阶导数矩阵）。
-
计算梯度（已得）：
\(\mathbf{g} = \sum_i (y_i - \lambda_i) \mathbf{x}_i = \mathbf{X}^T (\mathbf{y} - \boldsymbol{\lambda})\)，其中 \(\mathbf{X}\) 是设计矩阵，\(\boldsymbol{\lambda} = (\lambda_1, \dots, \lambda_n)^T\)。 -
计算海森矩阵：
\(H = \frac{\partial^2 \ell}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^T} = \sum_i \frac{\partial}{\partial \boldsymbol{\beta}}[(y_i - \lambda_i) \mathbf{x}_i]\)。
由于 \(\frac{\partial \lambda_i}{\partial \boldsymbol{\beta}} = \lambda_i \mathbf{x}_i\)，可得：
\[ H = -\sum_i \lambda_i \mathbf{x}_i \mathbf{x}_i^T = -\mathbf{X}^T \mathbf{W} \mathbf{X} \]
这里 $\mathbf{W}$ 是一个 $n \times n$ 的对角权重矩阵，其第 $i$ 个对角元素为 $w_i = \lambda_i = \exp(\mathbf{x}_i^T \boldsymbol{\beta})$。
- 代入牛顿法更新公式：
\[ \begin{aligned} \boldsymbol{\beta}^{(t+1)} &= \boldsymbol{\beta}^{(t)} - [-\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X}]^{-1} \cdot [\mathbf{X}^T (\mathbf{y} - \boldsymbol{\lambda}^{(t)})] \\ &= \boldsymbol{\beta}^{(t)} + (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T (\mathbf{y} - \boldsymbol{\lambda}^{(t)}) \end{aligned} \]
B. 重新表述为迭代加权最小二乘（IRLS）
上述公式可以巧妙变形。定义一个“工作响应变量” \(z_i\)：
\[z_i = \mathbf{x}_i^T \boldsymbol{\beta}^{(t)} + \frac{y_i - \lambda_i^{(t)}}{\lambda_i^{(t)}} = \eta_i^{(t)} + \frac{y_i - \lambda_i^{(t)}}{\lambda_i^{(t)}} \]
令 \(\mathbf{z} = (z_1, \dots, z_n)^T\)。则可以验证，上述牛顿更新公式等价于求解以下加权最小二乘问题：
\[\boldsymbol{\beta}^{(t+1)} = \arg\min_{\boldsymbol{\beta}} (\mathbf{z} - \mathbf{X}\boldsymbol{\beta})^T \mathbf{W}^{(t)} (\mathbf{z} - \mathbf{X}\boldsymbol{\beta}) \]
其解析解为：
\[\boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z} \]
IRLS算法步骤：
- 初始化：设定初始参数 \(\boldsymbol{\beta}^{(0)}\)（例如，全零或用简单线性回归的结果）。
- 迭代直至收敛（第 \(t\) 次迭代）：
a. 使用当前参数 \(\boldsymbol{\beta}^{(t)}\) 计算：
* 预测均值：\(\lambda_i^{(t)} = \exp(\mathbf{x}_i^T \boldsymbol{\beta}^{(t)})\)
* 权重矩阵：\(\mathbf{W}^{(t)} = \text{diag}(\lambda_1^{(t)}, \dots, \lambda_n^{(t)})\)
* 工作响应：\(z_i^{(t)} = \mathbf{x}_i^T \boldsymbol{\beta}^{(t)} + (y_i - \lambda_i^{(t)}) / \lambda_i^{(t)}\)
b. 求解加权最小二乘问题，更新参数：
\(\boldsymbol{\beta}^{(t+1)} = (\mathbf{X}^T \mathbf{W}^{(t)} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}^{(t)} \mathbf{z}^{(t)}\) - 终止条件：当参数的变化 \(\|\boldsymbol{\beta}^{(t+1)} - \boldsymbol{\beta}^{(t)}\|\) 小于预设阈值时，停止迭代。
5. 模型评估与过离散
- 模型评估：可以使用偏差（Deviance）或皮尔逊卡方统计量来评估模型拟合优度，并与空模型（仅截距）进行比较。
- 过离散（Overdispersion）：这是泊松回归一个关键假设检验。泊松分布要求均值等于方差。如果数据中观测到的方差明显大于均值，则存在“过离散”，可能导致标准误低估。解决方案包括：
- 使用准泊松回归，它引入一个离散参数来缩放方差。
- 使用负二项回归，它为计数数据提供了更大的灵活性，允许方差大于均值。
总结
泊松回归通过对数链接函数，将计数数据的期望值映射到特征的线性组合上，并通过极大似然估计参数。其优化通过高效的IRLS算法实现。理解此过程的关键在于掌握广义线性模型的三个组成部分，以及如何从泊松分布的假设出发，推导出用于数值优化的得分方程和迭代算法。在实际应用中，务必检查过离散假设。