4.4 对数几率回归

对数几率回归（logistic regression） ¹模型的思路为用 $x$ 的线性函数为 $K$ 个类型的后验概率建立模型，同时保证产生的概率和为 1 并且取值在 $[0, 1]$ 区间上。模型的形式为：

\begin{matrix} (4.17) & \begin{aligned} \log \frac{\Pr (G = 1 | X = x)}{\Pr (G = K | X = x)} & = β_{10} + β_{1}^{T} x \\ \log \frac{\Pr (G = 2 | X = x)}{\Pr (G = K | X = x)} & = β_{20} + β_{2}^{T} x \\ ⋮ \\ \log \frac{\Pr (G = K - 1 | X = x)}{\Pr (G = K | X = x)} & = β_{(K - 1) 0} + β_{K - 1}^{T} x \end{aligned} \end{matrix}

模型由 $K - 1$ 个对数几率（log-odds 或 logit）转换，减去一个是为了满足其和为 1 的约束条件。上式中用最后一个类型作为几率比的分母，但选择任意一个类型作为分母对最后的结果没有影响。经过简单的转换后：

\begin{aligned} \Pr (G = k | X = x) & = \frac{\exp (β_{k 0} + β_{k}^{T} x)}{1 + \sum_{l = 1}^{K - 1} \exp (β_{l 0} + β_{l}^{T} x)}, k = 1, \dots, K - 1 \\ (4.18) & \Pr (G = K | X = x) & = \frac{1}{1 + \sum_{l = 1}^{K - 1} \exp (β_{l 0} + β_{l}^{T} x)} \end{aligned}

显然这些概率的和为 1。为了凸显这个模型依赖于整个参数集合 $θ = β_{10}, β_{1}^{T}, \dots, β_{(K - 1) 0}, β_{K - 1}^{T}$ ，通常将概率写为 $\Pr (G = k | X = x) = p_{k} (x; θ)$ 。

当 $K = 2$ 时，模型非常简洁，只有一个线性函数。在生物统计中二元输出变量（两个类型）的场景很多，这个模型使用非常广泛。例如，患者是否生存、是否患有心脏病、或某个症状是否存在等等。

4.4.1 对数几率回归模型的拟合

对数几率回归模型通常对条件于 $X$ 的 $G$ 的似然度使用最大似然估计拟合。 $\Pr (G | X)$ 可完整地确定条件概率分布，故使用多项式分布是合理的选择。 $N$ 个样本的对数似然度为：

\begin{matrix} (4.19) & ℓ (θ) = \sum_{i = 1}^{N} \log p_{g_{i}} (x_{1}; θ) \end{matrix}

其中 $p_{k} (x_{i}; θ) = \Pr (G = k | X = x_{i}; θ)$ 。

下面详细介绍两个类型的情况，其算法简化了很多。简便起见，将两个类型 $g_{i}$ 编码为取值为 0 或 1 的输出变量 $y_{i}$ ，当 $g_{i} = 1$ 时 $y_{i} = 1$ ，当 $g_{i} = 2$ 时 $y_{i} = 0$ 。记 $p_{1} (x; θ) = p (x; θ)$ ， $p_{2} (x; θ) = 1 - p (x; θ)$ 。对数似然度可写为：

\begin{aligned} ℓ (β) & = \sum_{i = 1}^{N} {y_{i} \log p (x_{i}; β) + (1 - y_{i}) \log (1 - p (x_{i}; β))} \\ (4.20) & = \sum_{i = 1}^{N} {y_{i} β^{T} x_{i} - \log (1 + e^{β^{T} x_{i}})} \end{aligned}

其中 $β = {β_{10}, β_{1}}$ ，并假设输入变量向量 $x_{i}$ 包含了对应截距项的全为 1 的列。

通过设置一阶导数为 0 来求对数似然度的最大化，得到 评分（score） 等式

\begin{matrix} (4.21) & \frac{\partial ℓ (β)}{\partial β} = \sum_{i = 1}^{N} x_{i} (y_{i} - p (x_{i}; β)) = 0 \end{matrix}

为 $p + 1$ 个 $β$ 的非线性方程式。注意输入变量 $x_{i}$ 中的第一个元素为常数 1，所以第一个评分等式为 $\sum_{i = 1}^{N} y_{i} = \sum_{i = 1}^{N} p (x_{i}; β)$ ，即样本中类型一的个数的期望值与实际值相等，因此类型二同样成立。

运用牛顿-拉弗森（Newton–Raphson）求解评分等式 4.21，需要用到海森（Hessian）矩阵，即二阶导数矩阵：

\begin{matrix} (4.22) & \frac{\partial^{2} ℓ (β)}{\partial β \partial β^{T}} = - \sum_{i = 1}^{N} x_{i} x_{i}^{T} p (x_{i}; β) (1 - p (x_{i}; β)) \end{matrix}

若从 $β^{old}$ 出发，一次牛顿算法更新为：

\begin{matrix} (4.23) & β^{new} = β^{old} - {(\frac{\partial^{2} ℓ (β)}{\partial β \partial β^{T}})}^{- 1} \frac{\partial ℓ (β)}{\partial β} \end{matrix}

其中的一阶和二阶导数在 $β^{old}$ 处取值。

用矩阵的形式标记评分和海森可以更简便地进行推导。让 $y$ 表示 $y_{i}$ 输出变量的向量， $X$ 为输入变量 $x_{i}$ 组成的 $N \times (p + 1)$ 的矩阵， $p$ 为拟合的概率值向量，其第 i 个元素为 $p (x_{i}; β^{old})$ ， $W$ 为 $N \times N$ 的加权对角矩阵，其第 i 个对角线元素为 $p (x_{i}; β^{old}) (1 - p (x_{i}; β^{old}))$ 。则导数可写为：

\begin{aligned} (4.24) & \frac{\partial ℓ (β)}{\partial β} & = X^{T} (y - p) \\ (4.25) & \frac{\partial^{2} ℓ (β)}{\partial β \partial β^{T}} & = - X^{T} W X \end{aligned}

牛顿算法更新可写为：

\begin{aligned} β^{new} & = β^{old} + (X^{T} W X)^{- 1} X^{T} (y - p) \\ = (X^{T} W X)^{- 1} X^{T} W (X β^{old} + W^{- 1} (y - p)) \\ (4.26) & = (X^{T} W X)^{- 1} X^{T} W z \end{aligned}

从第二和第三行可见牛顿算法的步骤可以被视为一个加权最小二乘，其使用的输出变量为：

\begin{matrix} (4.27) & z = X β^{old} + W^{- 1} (y - p) \end{matrix}

其有时称之为 调整输出变量（adjusted response）。在每一次更新中， $p$ 、 $W$ 、和 $z$ 均不同，因此在每一步中都要进行重新计算。由于在每一步迭代都是解一个加权最小二乘问题，这个算法也被称为 迭代重加权最小二乘（iteratively reweighted least squares，IRLS）：

\begin{matrix} (4.28) & β^{new} \leftarrow \underset{β}{\arg min} (z - X β)^{T} W (z - X β) \end{matrix}

一般 $β = 0$ 是一个不错的迭代过程初始点，但理论上并不能保证最终的收敛。由于对数似然度函数为凹函数，这个算法通常会收敛，但也有过调（overshooting）的可能性。更新中很少会出现对数似然度下降的情况，若有可通过将步长减半来解决。

多类型问题（ $K \geq 3$ ）的牛顿算法也可被表达成迭代重加权最小二乘方法，其输出变量个数为 $K - 1$ ，而且每个观测的加权矩阵为非对角矩阵。这样的加权矩阵使简化算法不再适用，此时直接对综合的参数向量 $θ$ 进行数值计算反而更方便（练习 4.4）。或者，第 3.8.6 节介绍的坐标下降方法可快速地最大化对数似然度。R 程序包 glmnet （Friedman et al., 2010）可以快速地拟合较大 $N$ 和 $p$ 的对数几率回归问题。其为拟合正则化的模型而设计，但通过选项也可拟合没有正则化的模型。

对数几率回归模型大多作为一个数据分析和推断的工具，其目的在于理解输入变量对输出变量的影响方式和程度。很多时候会通过对数几率回归的结果来选择输入变量以及可能有交叉项，最终建立一个简洁的模型。

4.4.2 例：南非心脏病数据

下面通过一个二元输出变量数据的分析，来演示对数几率回归模型在传统统计学中的应用。图 4.12 中展示的是在南非西开普省的三个农村地区进行的（Rousseauw et al., 1983）“冠心病风险因子研究（Coronary Risk Factor Study, CORIS）”基础研究的部分数据。这项研究的目的是探究缺血性心脏病在高发地区中的各种风险因子的严重程度。数据中包括 15 到 64 岁的白人男性，输出变量为在调查时是否存在心肌梗塞（MI）（区域整体的患病率为 5.1%）。图中的数据集中包括了 160 个患病案例和 302 个未患病对照组。对数据更详细的说明可参考 Hastie and Tibshirani (1987)。

**图 4.12**：南非心脏病数据的散点图矩阵。每个图为两个风险因子的散点图，红色点为患病样本，蓝色为未患病样本。心脏病家族史的输入变量（famhist）是一个二元变量（是与否）。

表 4.2 为用最大似然估计法拟合的对数几率回归模型的结果。表中包含了模型中每个系数的 Z 分数（系数值与其标注误差的比值）；Z 分数不显著代表其对应的输入变量可以从模型中排除。每个 Z 分数对应的是在该系数为 0 的原假设下的检验统计量（沃德检验）。在 5% 的置信水平上显著的条件大约为 Z 评分的绝对值大于 2。

	系数值	标准误差	Z 分数
(Intercept)	−4.130	0.964	−4.285
sbp	0.006	0.006	1.023
tobacco	0.080	0.026	3.034
ldl	0.185	0.057	3.219
famhist	0.939	0.225	4.178
obesity	-0.035	0.029	−1.187
alcohol	0.001	0.004	0.136
age	0.043	0.010	4.184

表 4.2：南非心脏病数据的对数几率回归拟合结果。

结果中有些出乎意料的发现，需要谨慎地审视。收缩压（高压，sbp）不是显著变量！肥胖度（obesity）也不显著，并且其系数方向为负。这种违背常理的结果是由系数之间的相关性造成的。在单独使用时，两个变量 sbp 和 obesity 均为显著的，且符号为正。但在存在其他与之相关的变量时，这个结论不一定继续成立（甚至会产生符号相反的系数）。

根据上面的分析结果，可对模型做一些筛选，即挑选一部分可以较好的对是否患病（chd）做出解释的输入变量。一种方法是排除最不显著的变量，对模型重新拟合，再进行排除；直到模型中所有的变量均为显著变量。这个方法的结果为表 4.3：

	系数值	标准误差	Z 分数
(Intercept)	−4.204	0.498	−8.45
tobacco	0.081	0.026	3.16
ldl	0.168	0.054	3.09
famhist	0.924	0.223	4.14
age	0.044	0.010	4.52

表 4.3：南非心脏病数据的逐步对数几率回归拟合结果。

一个更优但更费时的策略是轮流去掉每一个变量后重新拟合模型，通过分析偏差来绝对排除哪一个变量。拟合的残差偏差为 $- 2$ 乘以拟合后的对数似然度，两个模型的差别为两者的残差偏差的差（类比于残差平方和）。这个策略最终会给出和表 4.3 中一样的结果。

以吸烟量（tobacco）为例，如何解释一个系数值 0.081 和标准误差 0.026 的影响？吸烟量衡量的是以公斤为单位的一生吸烟量，未患病组的中位数为 1 公斤，患病组的中位数为 4.1 公斤。因此，每增加 1 公斤的一生吸烟量，可增加冠心病患病几率 $\exp (0.081) = 1.084$ ，或 8.4%。考虑进标准误差后，可得到大约 95% 的置信区间为 $\exp (0.081 \pm 2 \times 0.026) = (1.03, 1.14)$ 。

第五章会继续以这个数据集为例，届时可看到其中某些变量的影响是非线性的，在合适的模型中不用被排除模型。

4.4.3 二次项近似和推断

最大似然的参数估计 $\hat{β}$ 满足一个自洽的关系：它是加权最小二乘估计的参数，其输出变量为：

\begin{matrix} (4.29) & z_{i} = x_{i}^{T} \hat{β} + \frac{(y_{i} - {\hat{p}}_{i})}{{\hat{p}}_{i} (1 - {\hat{p}}_{i})} \end{matrix}

其权重为 $w_{i} = {\hat{p}}_{i} (1 - {\hat{p}}_{i})$ ，输出变量和权重均依赖于 $\hat{β}$ 本身²。与最小二乘的关联不仅为其提供了简便的算法，同时也使其有以下性质：

加权残差平方和即为常见的皮尔森卡方检验统计量： {« math »} $\begin{matrix} (4.30) & \sum_{i = 1}^{N} \frac{(y_{i} - {\hat{p}}_{i})^{2}}{{\hat{p}}_{i} (1 - {\hat{p}}_{i})} \end{matrix}$ {« /math »} 也是对偏差的二次项近似。
似然度的渐进理论可证明，如果模型的选择正确，则 $\hat{β}$ 为一致估计量（即随着样本量增大，会收敛到真实的参数 $β$ ）。
利用中央极限定理可证明 $\hat{β}$ 的分布收敛于 $N (β, (X^{T} W X)^{- 1})$ 。这些渐进性质均可从加权最小二乘拟合推导出来，其过程类似于正态分布下的推断理论。
由于拟合中需要迭代计算，对数几率回归模型的计算可能耗时较多。一些常见的加速技巧包括用 拉奥评分检验（Rao score test） 判断是否包含某个变量，用 沃德检验（Wald test） 判断是否排除某个变量。两者均基于当前模型的最大似然估计，不涉及迭代拟合。两者的结果均为从加权最小二乘法中添加或去除变量，而维持原有权重值。这样的计算过程不用反复地计算整个加权最小二乘拟合，因此可快速完成。

软件中的实现会利用上述的性质。例如，R 中的广义线性模型（其以对残差为二项分布的模型支持覆盖了对数几率回归）即利用了这些性质。拟合的结果为 GLM 类型（generalized linear model）的对象，其继承了线性模型对象的性质，所以可应用于线性模型的工具函数仍然适用。

4.4.4 $L_{1}$ 正则化对数几率回归

套索回归（第 3.4.2 节）中的 $L_{1}$ 惩罚项可起到变量选择和参数收缩的作用，可用在任意线性回归模型中。在对数几率回归中，加上正则化的式 4.20 为：

\begin{matrix} (4.31) & max_{β_{0}, β} {\sum_{i = 1}^{N} [y_{i} (β_{0} + β^{T} x_{i}) - \log (1 + e^{β_{0} + β^{T} x_{i}})] - λ \sum_{j = 1}^{p} | β_{j} |} \end{matrix}

与套索回归类似，通常截距项不被加入到惩罚项中，并且事先对输入变量进行标准化处理。式 4.31 中的最大化准则为凹函数，可通过非线性规划方法求解（例如 Koh et al., 2007）。或者，利用类似于第 4.4.1 节中牛顿算法的二次项近似，用迭代的加权套索算法求解式 4.31。有趣的是，系数非零的输入变量的评分等式（对比式 4.24）为：

\begin{matrix} (4.32) & x_{j}^{T} (y - p) = λ \cdot sign (β_{j}) \end{matrix}

其为第 3.4.4 节式 3.58 的推广；可理解为，所有选入模型的变量与残差的广义相关性相同。

此时系数的曲线不再是分段线性而是分段曲线，类似于套索回归中的最小角回归的路径算法变得更困难。不过仍可利用二次项近似进行迭代计算。

**图 4.13**：南非心脏病数据的 $L_1$ 正则化对数几率回归系数，曲线图的横轴为系数的 $L_1$ 范式。输入变量均已标准化为方差 1。曲线图为逐点进行计算后描绘的。 — **图 4.13**：南非心脏病数据的 $L_{1}$ 正则化对数几率回归系数，曲线图的横轴为系数的 $L_{1}$ 范式。输入变量均已标准化为方差 1。曲线图为逐点进行计算后描绘的。

图 4.13 为建立在第 4.4.2 节中的南非心脏病数据上的 $L_{1}$ 正则化系数路径。计算使用的是 R 程序包 glmpath （Park and Hastie, 2007），其中运用了凸优化的 预测-校正（predictor-corrector） 方法来确定当被选入模型的变量集合发生变动时（图中竖线处） $λ$ 的取值。这个例子中，曲线看起来近乎于分段线性；在其他例子中，曲线的非线性会更明显。

坐标下降方法（第 3.8.6 节）可以快速地计算很多 $λ$ 的取值网格点上的系数值。R 程序包 glmnet （Friedman et al., 2010）可快速地计算出较大 $N$ 或较大 $p$ 的对数几率回归问题的系数路径。他们的算法利用了自变量矩阵 $X$ 的稀疏性，因而可适用于更大维度的问题。第 18.4 节会介绍更深入细节，并讨论 $L_{1}$ 正则化的多项分布模型。

4.4.5 对数几率回归与线性判别分析

在第 4.3 节中，可推导出类型 k 和类型 K 之间的对数后验几率比为 $x$ 的线性函数（式 4.9）：

\begin{aligned} \log \frac{\Pr (G = k | X = x)}{\Pr (G = K | X = x)} & = \log \frac{π_{k}}{π_{K}} - \frac{1}{2} (μ_{k} + μ_{K})^{T} Σ^{- 1} (μ_{k} - μ_{K}) \\ + x^{T} Σ^{- 1} (μ_{k} - μ_{K}) \\ (4.33) & = α_{k 0} + α_{k}^{T} x \end{aligned}

这个线性关系来自于类型条件密度函数的高斯分布假设和具有相同协方差矩阵的假设。而从定义上（式 4.17），对数几率回归模型中的对数几率为线性函数：

\begin{matrix} (4.34) & \log \frac{\Pr (G = k | X = x)}{\Pr (G = K | X = x)} = β_{k 0} + β_{k}^{T} x \end{matrix}

两个模型看起来形式是一样的。尽管如此，两者系数的估计方法不同。对数几率回归模型更广泛，即它所依赖的假设更少。 $X$ 与 $G$ 的联合分布可写为：

\begin{matrix} (4.35) & \Pr (X, G = k) = \Pr (X) \Pr (G = k | X) \end{matrix}

其中 $\Pr (X)$ 为输入变量 $X$ 的边际密度函数。线性判别分析和对数几率回归中，右手边的第二项均为线性的对数几率形式：

\begin{matrix} (4.36) & \Pr (G = k | X = x) = \frac{e^{β_{k 0} + β_{k}^{T} x}}{1 + \sum_{ℓ = 1}^{K - 1} e^{β_{ℓ 0} + β_{ℓ}^{T} x}} \end{matrix}

其中任意地选择了最后一个类型 K 作为参考类型³。

对数几率回归将 $X$ 的边际密度函数处理为任意的密度函数 $\Pr (X)$ ，而只通过最大化条件似然度——由概率 $\Pr (G = k | X)$ 构成的多项分布似然度——来拟合 $\Pr (G | X)$ 中的参数。尽管完全忽视了 $\Pr (X)$ 部分，但也可以理解为是在用完全的非参数无约束的方式来估计这个边际密度函数，即根据经验分布函数给每个样本基于 $1 / N$ 的概率。

而在线性判别分析中，通过基于联合分布的对数似然度的最大化里拟合参数：

\begin{matrix} (4.37) & \Pr (X, G = k) = ϕ (X; μ_{k}, Σ) π_{k} \end{matrix}

其中 $ϕ$ 为高斯密度函数。根据“标准正态理论（standard normal theory）”，可推导出同第 4.3 节中的估计量 ${\hat{μ}}_{k}$ 、 $\hat{Σ}$ 、和 ${\hat{π}}_{k}$ 。由于在对数几率等式 4.33 中，线性模型的参数为高斯分布参数的函数，带入高斯分布参数的估计量后，即可得到线性参数的最大似然估计值。然而与条件分布的最大似然结果不同，这里的“似然”中包含了边际密度函数 $\Pr (X)$ 的影响。边际密度函数为一个混合密度：

\begin{matrix} (4.38) & \Pr (X) = \sum_{k = 1}^{K} π_{k} ϕ (X; μ_{k}, Σ) \end{matrix}

可见其中包含了模型的参数。

似然度中额外的因子或约束会对拟合产生什么影响？附加的模型假设包含了关于参数的更多信息，因此加以利用后可更有效地估计参数（更低的方差）。若真实的 $f_{k} (x)$ 即为高斯分布，则在似然度中忽略边际部分，会在渐进误差率中丢失高至 30% 的信息（Efron, 1975）。原结论复述为：基于条件似然度的估计在样本量多 30% 时能达到同样的拟合程度。

远离判别边界的样本点，在对数几率回归中会赋予低权重，但在线性判别分析中仍被用于估计协方差矩阵。但这也会使得线性判别分析对离群值不够稳定。

从混合分布的式 4.38 中可见，没有被标记类型的观测点仍包含了参数的信息。实践中样本的类型标记通常比较珍贵，而获得无标记的样本点相对容易很多。通过附加比较强模型假设，例如上述，可以同时利用有标记和无标记的数据中关于参数的信息。

边际似然度可被理解为某种正则条件（regularizer），从某种程度上要求边际分布包含了类型的概率分布的一些信息。例如，若在可以被一个超平面完全区分开的两个类型的数据中，对数几率回归的最大似然估计无法定义（有无穷多个解，见练习 4.5）。而在同样的数据上，线性判别分析的系数是可计算的，其边际似然度避免了信息维度的降低（degeneracy）。

在实践中不存在正确的假设，而且通常输入变量 $X$ 中存在分类的变量。一般会认为对数几率回归比线性判别分析更安全、稳健，依赖于更少的假设。以作者的经验来说，甚至在诸如存在分类自变量的并不适合线性判别分析的情景下，两个模型通常会得到非常的近似的结果。

本节练习

练习 4.4

Consider the multilogit model with K classes (4.17). Let β be the (p + 1)(K − 1)-vector consisting of all the coefficients. Define a suitably enlarged version of the input vector x to accommodate this vectorized coefficient matrix. Derive the Newton-Raphson algorithm for maximizing the multinomial log-likelihood, and describe how you would implement this algorithm.

练习 4.5

Consider a two-class logistic regression problem with x ∈ IR. Characterize the maximum-likelihood estimates of the slope and intercept parameter if the sample xi for the two classes are separated by a point x 0 ∈ IR. Generalize this result to (a) x ∈ IR p (see Figure 4.16), and (b) more than two classes.

印象中似乎“逻辑回归”的叫法更多，但与西瓜书的观点一致，这个音译的方式没有逻辑。 ↩︎
即为“自洽（self consistency）”的含义，中文的译法参考了“诺维科夫自洽性原则”。 ↩︎
即构建线性的对数几率模型形式时，作为几率的分母。 ↩︎