9.1 广义加性模型

回归模型在很多数据分析中扮演着重要的角色，它不仅给出了预测和分类规则，并为理解不同输入变量的重要程度提供了数据分析工具。

传统的线性模型尽管非常简洁，但通常会在这些分析中失效：在现实生活中，影响往往不是线性的。之前的章节介绍的方法通过使用预先定义的基函数来达到非线性关系。本节介绍可用于识别和刻画非线性回归关系的更自动和灵活的统计方法。这些方法被称为 广义加性模型（generalized additive model）。

在回归问题中，广义加性模型的形式为：

同以往一样，令 $X_1$、$X_2$、……、$X_p$ 代表自变量，$Y$ 代表输出变量；$f_j$ 代表一些未指定的平滑（非参数）函数。假如使用基函数展开来为这些函数建模（第五章），则得出的模型就可以用简单的最小二乘来拟合。我们在这里采用不同的方法：先用散点图平滑器（例如三次平滑样条或核平滑器）拟合每个函数，然后通过一个算法同时对所有 $p$ 个函数进行估计（第 9.1.1 节）。

回顾第 4.4 节中解决二分类问题的对数几率（逻辑）回归模型。其中将二元输出变量的均值 $\mu(X)=\operatorname{Pr}(Y=1|X)$ 与自变量的关系表达为线性回归模型和对数几率（logit）的 联系函数（link function）：

加性对数几率回归模型 将每个线性项替换为更一般性的函数形式：

其中每个 $f_j$ 同样是未指定的平滑函数。非参数形式的函数 $f_j$ 使模型更灵活，同时也保持了“加性”，从而我们可以跟和之前基本相同的方式来理解模型。加性对数几率回归模型是广义加性模型的一个例子。一般地，一个输出变量 $Y$ 的条件均值 $\mu(X)$ 可以与自变量的加性函数之间通过 联系函数（link function） $g$ 来建立联系：

以下为一些经典的联系函数的例子：

• $g(\mu)=\mu$，恒等联系函数，用于高斯输出变量的线性和加性模型。
• $g(\mu)=\operatorname{logit}(\mu)$，如上文；或 概率单位（probit） 函数 $g(\mu) = \text{probit}(\mu)$，用于二项分布概率模型。概率单位函数是高斯累积分布函数的反函数：$\operatorname{probit}(\mu)=\Phi^{-1}(\mu)$。
• $g(\mu)=\log(\mu)$，用于泊松计数样本的对数线性或对数加性模型。

这三个都属于指数类型的抽样概率模型，指数类型中另外还包括了伽玛（gamma）分布和负二项分布。这些概率模型被用于广义线性模型中一些常用的模型，也都可以用同样的方式扩展到广义加性模型。

函数 $f_j$ 的估计方法比较灵活，所使用的算法的基础构成是散点图平滑器。估计出的函数 $\hat{f}_j$ 可用于展现出 $X_j$ 中可能存在的非线性影响。但并不是所有的函数 $f_j$ 都需要是非线性的。广义加性模型可以方便地将线性以及其他参数形式与非线性项混合在一起，这也是当输入变量中存在定性类型（因子）变量时仍可进入模型的一个必要条件。非线性关系也不局限于输入变量对应的函数 $f_j$ 上；比如可以有两个或更多变量的（交叉）非线性成分，或根据因子变量 $X_k$ 取不同类别时 $X_j$ 采用不同的函数曲线。因此以下描述的都是非线性关系：

• $g(\mu)=X^T\beta+\alpha_k+f(Z)$
  半参数（semiparametric）模型，其中的 $X$ 为在模型中呈线性的自变量向量，$\alpha_k$ 为分类输入变量 $V$ 取值第 k 类别时对结果的影响，而自变量 $Z$ 的影响则是非参数的形式。
• $g(\mu)=f(X)+g_k(Z)$
  其中的 $k$ 同样代表了分类输入变量 $V$ 的取值类别，因此产生了一个 $V$ 和 $Z$ 影响的交叉项 $g(V,Z)=g_k(Z)$。
• $g(\mu)=f(X)+g(Z,W)$
  其中的 $g$ 为两个特征变量的非参数函数。

加性模型在很多种场景中可以替换线性模型，例如可对一个时间序列进行加性分解：

其中的 $S_t$ 为周期成分，$T_t$ 为趋势项，$\varepsilon_t$ 为误差项。

9.1.1 拟合加性模型

本节介绍一个拟合加性模型及其推广形式的模块化算法。算法的基础模块是灵活地拟合非线性影响的散点平滑器。具体来说，这里用的是第五章中介绍的三次平滑样条作为散点平滑器。

加性模型的表达式为：

其中的误差项 $\varepsilon$ 均值为零。给定观测样本 $(x_i,y_i)$，可以为这个问题制定一个类似于第 5.4 节中的惩罚残差平方和（式 5.9）的准则：

其中的 $\lambda_j\geq0$ 为调节参数。可证明式 9.7 的最小值解是一个加性三次样条模型；每个函数 $f_j$ 都是对应的 $X_j$ 的三次样条，其结点为 $x_{ij}$ 的唯一值，$i=1,\dots,N$。然而在不对模型做进一步假设时，这个解不唯一。常数 $\alpha$ 就无法唯一求解，因为可以对每个函数 $f_j$ 加减任意常数，并相应地调整 $\alpha$。标准惯例是假定 $\sum_{i=1}^N f(x_{ij})=0\,\forall j$，即函数在数据样本上的平均值为零。在这个惯例下可见 $\hat{\alpha}=\operatorname{ave}(y_i)$。如果在这个约束外，再要求输入变量的矩阵（第 ij 个元素为 $x_{ij}$）为满列秩矩阵，那么式 9.7 是一个严格凸性的准则函数，有唯一的最小值解。如果这个矩阵是一个奇异（singular）矩阵，则无法唯一确定函数成分 $f_j$ 中的线性部分（但仍然可确定其非线性部分）（Buja et al., 1989）。

此外，这个问题的求解存在一个简单迭代计算过程。令 $\hat{\alpha}=\operatorname{ave}(y_i)$，并固定不变。然后将 $\{y_i-\hat{\alpha}-\sum_{k\ne j}\hat{f}_k(x_{ik})\}_1^N$ 作为目标变量，将第 j 个变量 $\{x_{ij}\}_{i=1}^N$ 作为自变量，使用三次平滑样条 $\mathbf{S}_j$ 获得一个新的函数估计 $\hat{f}_j$。轮流对每个自变量进行这个操作，使用其他函数当前的估计 $\hat{f}_k$ 来计算 $y_i-\hat{\alpha}-\sum_{k\ne j}\hat{f}_k(x_{ik})$。持续这个过程，直到函数估计 $\hat{f}_j$ 达到稳定状态。算法 9.1 给出了这个过程的细节，它也称为 回修（backfitting），并且它的拟合可类比于线性模型中的多次回归。¹

算法 9.1：加性模型的回修（backfitting）算法

1. 初始化：$\hat{\alpha}=\frac{1}{N}\sum_{i=1}^Ny_i$，$\hat{f}_j\equiv 0$，$\forall i,j$。
2. 循环迭代: $j=1,2,\dots,p,\dots,1,2,\dots,p,\dots$, $$\begin{align} & \hat{f}_j \leftarrow \mathcal{S}_j \left[ \{y_i - \hat{\alpha} - \sum_{k \ne j}\hat{f}_k(x_{ik})\}_1^N \right] \\ & \hat{f}_j \leftarrow \hat{f}_j - \frac{1}{N}\sum_{i=1}^N \hat{f}_j(x_{ij}) \end{align}$$ 直到函数 $\hat{f}_j$ 的变动小于某个预设的阈值。

理论上来说算法 9.1 第（2）步中对 $\hat{f}_j$ 的第二次更新是没必要的，因为对零均值输出变量拟合的平滑样条，其均值也为零（练习 9.1）。在实践中，运算中的近似有可能造成结果的偏移，所以建议进行相应调整。

这个算法也可以完全相同的方式应用于其他拟合方法，只需要指定合适的平滑运算 $\mathcal{S}_j$：

• 其他单变量回归平滑器，比如局部多项回归和核方法；
• 可以产生多项式拟合、分段常数拟合、参数样条拟合、序列和傅立叶拟合的线性回归运算；
• 其他更复杂的运算，比如捕捉二次或更高次交叉作用的曲面平滑器，或捕捉周期性作用的周期平滑器。

如果只是在训练集上考虑平滑器运算 $\mathcal{S}_j$，则它可表达为一个 $N\times N$ 的算子矩阵（见第 5.4.1 节）。那么与第五和第六章讨论的平滑器的自由度相类比，第 j 项的自由度可（近似地）计算为 $\text{df}_j=\operatorname{trace}[\mathbf{S}_j]-1$。

对很多类型的线性平滑器 $\mathbf{S}_j$，回修算法在求解一个线性方程组时与高斯-赛德尔（Gauss-Seidel）算法是等价的。练习 9.2 给出了相关细节。

对于对数几率回归模型和其他部分广义加性模型而言，惩罚对数似然度是一个合理的准则。在对其最大化中，同时使用了回修过程和似然函数最大化方法。可将常规的最大化广义线性模型的对数似然函数的牛顿-拉弗森（Newton-Raphson）方法重改为一个 迭代重加权最小二乘算法（iteratively reweighted least squares，IRLS）。这个过程是重复地进行当前输出变量对协变量（covariate）的加权线性回归；每个回归会产生新的参数估计值，这反过来也会得出新的当前输出变量和权重，并迭代这个过程（见第 4.4.1 节）。在广义加性模型中，就是简单地将其中的加权线性回归替换成加权回修算法。下面更详细地介绍这个算法在对数几率回归中的应用，并且在 Hastie and Tiebshirani (1990) 的第六章有更一般性地介绍。

9.1.2 示例：加性对数几率回归

二分类数据的对数几率模型可能是医学研究中使用最广泛的模型。这个模型中的输出变量 $Y$ 可被编码为 0 或 1，用 1 表示一个事件（如死亡或疾病复发）的发生，而用 0 表示没有事件发生。我们要在给定的预后（prognostic）因子的取值上建立一个事件发生概率 $\operatorname{Pr}(Y=1|X)$ 的模型。通常的目标是为了理解这些预后因子的作用，而不是对新个体进行分类。对数几率回归也适用于估计类别概率从而筛查风险的应用中。除医学应用外，模型也在信用风险识别里普遍使用。

广义加性对数几率模型的表达式为：

函数 $f_1$，$f_2$，……，$f_p$ 在牛顿-拉弗森（Newton–Raphson）过程中通过回修（backfitting）算法估计出，详见算法 9.2。

算法 9.2：加性对数几率回归模型的局部分数（scoring）算法

1. 计算初始值：$\hat{\alpha}=\log[\bar{y}/(1−\bar{y})]$,其中的 $\bar{y}=\operatorname{ave}(y_i)$，即样本中类别一的比例，并设置函数 $\hat{f}_j\equiv 0\forall j$。
2. 定义 $\hat{\eta}_i=\hat{\alpha}+\sum_j\hat{f}_j(x_{ij})$ 以及 $\hat{\rho}_i=1/[1+\exp(−\hat{\eta}_ i)]$。进行迭代：
   1. 构建当前目标变量： $$z_i = \hat{\eta} + \frac{(y_i - \hat{p}_i)}{\hat{p}_i(1-\hat{p}_i)}$$
   2. 构建权重 $w_i = \hat{p}_i (1-\hat{p}_i)$。
   3. 通过一个加权回修算法，对目标变量 $z_i$ 和权重 $w_i$ 拟合一个加性模型。结果得出新的估计 $\hat{\alpha}$，$\hat{f}_j$，$\forall j$。
3. 持续步骤 2，直到函数的变动小于某个预设的阈值。

算法 9.2 的步骤（2.2）中的加性模型拟合需要一个加权的散点图平滑器。大多的平滑过程都适用于带权重的观测样本（练习 5.12）；Hastie and Tibshirani (1990) 的第三章提供了更多细节。

利用第 4.4 节中简述的多项对数几率（multilogit）表达式，可将加性对数几率回归模型进一步推广到两个以上类别的问题中。尽管其表达式是对式 9.8 的直接推广，但拟合这个模型的算法要更复杂。更多细节见 Yee and Wild (1996)；另外软件 VGAM 可见下面的网站：

http://www.stat.auckland.ac.nz/∼yee.

例：预测垃圾邮件

我们在第一章中介绍的垃圾邮件数据上应用广义加性模型。数据为 4601 个电子邮件中的信息，来自于一个筛查垃圾（spam or junk）邮件的研究。数据可从 ftp.ics.uci.edu 公开获取，数据由加利福尼亚州（California）帕洛阿尔托市（Palo Alto）惠普（Hewlett-Packard）实验室的 George Forman 捐赠。

输出变量是二分类，取值为“正常邮件”或“垃圾邮件”，并且有 57 个自变量：

• 48 个数值自变量——某个单词在一封邮件所有单词中所占的百分比。例如 business、address、internet、free和george。这是考虑到这些单词可能可以区分不同的用户。
• 6 个数值自变量——某个符号在一封邮件所有符号中所占的百分比。这几个符号为 ch;、ch(、ch[、ch!、ch$ 和 ch#。
• 连续无间断的大写字母字符串的平均长度：CAPAVE。
• 连续无间断的大写字母字符串的平均长度：CAPMAX。
• 连续无间断的大写字母字符串的长度之和：CAPTOT。

将“垃圾邮件”编码为 1，“正常邮件”编码为 0。随机选出 1536 个样本作为测试集，留在训练集中的观测样本有 3065 个。在拟合广义加性模型中，使用了三次平滑样条，每个自变量有四个名义的（nominal）自由度。也就是说对每个自变量 $X_j$，选择使得 $\operatorname{trace}[\mathbf{S}_j(\lambda_j)]-1=4$ 的平滑样条参数 $\lambda_j$，其中的 $\mathbf{S}_j(\lambda)$ 为观测值 $x_{ij}$，$i=1,\dots,N$ 构建的平滑样条算子矩阵。这是一个在这样复杂模型中设定平滑程度的一个简便的方法。

大部分自变量的分布都有很重的长尾（long-tailed）。在拟合广义加性模型前，先对每个变量进行了对数转换（实际的转换是 $\log(x+0.1)$），但图 9.1 中的曲线仍是对原始变量的函数。

表 9.1：用垃圾邮件训练数据拟合的加性对数几率回归模型，在测试集上的混淆（confusion）矩阵。总体测试误差率为 5.5%。

表 9.1 给出了测试误差率；总体的误差率为 5.3%。与之相比，线性对数几率回归的测试误差率为 7.6%。表 9.1 给出了加性模型中最显著的一些自变量。

表 9.2：用垃圾邮件训练数据拟合的加性模型中的显著自变量。系数、标准误差和 Z 评分代表的是 $\hat{f}_j$ 的线性部分。非线性 P 值来自对 $\hat{f}_j$ 的非线性检验。

为了易于理解，表 9.2 中每个变量的贡献被分解为线性成分和剩余的非线性成分。上半区为与垃圾邮件有正相关性的自变量，而下半区为有负相关性的自变量。线性的成分是来自拟合的曲线对自变量的加权最小二乘线性拟合，而非线性的成分就是其残差。系数、标准误差和 Z 分数概括了估计函数的线性成分的信息；Z 分数是系数与其标准误差的比值，若它超出某个标准正太分布的分位数则认为对应变量是显著的。“非线性 P 值”一列是对估计函数的非线性检验结果。然而需要注意的是，每个自变量的影响中已经排除了其他自变量整体的影响，而不仅仅是它们影响中线性的部分。表中展示的为至少有一个（线性或非线性）检验中在（双边） $p=0.01$ 水平上呈显著的自变量。

图 9.1 展示了表 9.2 中出现的显著自变量的函数估计。大部分的非线性影响看起来是由于在零点处的不连续性。例如，随着 george 的频率从零点增加，垃圾邮件的概率显著地下降，但之后却几乎不再变动。据此可尝试将每个频率自变量替换为零出现的指示变量，然后实施线性对数几率模型。这样得到的测试误差率为 7.4%；同时加入频率的线性作用后，测试误差率降为 6.6%。这样看来，加性模型中的非线性成分可提供额外的预测能力。

相比之下将一个正常邮件分类为垃圾邮件是后果更严重的错误，因为这样会使一封正常的邮件被过滤掉而无法发送给用户。可以通过改变损失函数来改变不同分类误差的轻重程度。如果将真实类别为 0 的样本预测为 1 的损失设为 $L_{01}$，将真实类别为 1 的样本预测为 0 的损失设为 $L_{10}$，那么根据估计的贝叶斯规则，当类别 1 的概率大于 $L_{01}/(L_{01}+L_{10})$ 时分类为类别 1。例如，如果令 $L_{01}=10$ 和 $L_{10}=1$，则（实际的）类别 0 和 类别 1 的误差率变为 0.8% 和 8.7%。

更进一步，为了促使模型对类别 0 这部分数据的拟合结果更好，可以对类别 0 的观测值使用权重 $L_{01}$并对类别 1 的观测值使用权重 $L_{10}$。再用同上的估计贝叶斯规则进行预测。这样得到的（实际的）类别 0 和 类别 1 的误差率分别为 1.2% 和 8.0%。后续会在树结构（tree-based）模型的场景中会进一步讨论不相等损失的话题。

在拟合了加性模型后，需要核查是否加入某些交叉项可能会显著地改善拟合。可通过手动的方式进行，即加入某些或全部显著输入变量的乘积；或者可通过 MARS 过程自动地完成（第 9.4 节）。

这个例子以自动化的方式应用加性模型。而作为一个数据分析工具，加性模型常常会以更交互的方式使用，即逐步增加或去除某些项来查看它们的作用。另外通过调整 $\text{df}_j$ 从而调节平滑程度，可以使拟合结果在线性模型（$\text{df}_j=1$）与部分线性模型（让某些项有更灵活的拟合）之间无缝地切换。更多细节见 Hastie and Tibshirani (1990)。

9.1.3 总结

加性模型是对线性模型的一个有效的推广，使模型更加灵活却仍基本上保持了其可解释性。线性模型中熟悉的建模和推断的工具仍然适用于加性模型，可参考表 9.2 的例子。拟合这些模型的回修过程简便而且模块化，即可为每个输入变量选择相应合适的拟合方法。因此加性模型在统计学中被广泛地使用。

然而加性模型在大型的数据挖掘应用中可能存在缺陷。回修算法会拟合所有的自变量，当存在很多自变量时，这个算法不太理想甚至不再可行。BRUTO 过程（Hastie and Tibshirani, 1990，第九章）在回修算法中加入了变量选择，但其设计并没有考虑大型的数据挖掘问题。近期也有一些用套索（lasso）类型的惩罚项来估计稀疏加性模型的研究，例如 Lin and Zhang (2006) 的 COSSO 过程以及 Ravikumar et al. (2008) 提出的 SpAM 方法。对大型的问题，一个诸如提升方法（第十章）的前向分段方法更有效，并且同时可在模型中加入交叉项。

本节练习

练习 9.1

证明 $x_i$ 对 $y_i$ 的一个平滑样条拟合保持了（preserve）拟合中的线性部分。或者说，如果 $y_i=\hat{y}_i+r_i$，其中的 $\hat{y}_i$ 表示线性回归的拟合，并且 $\mathbf{S}$ 是平滑矩阵，那么 $\mathbf{S}\mathbf{y}=\hat{\mathbf{y}}+\mathbf{S}\mathbf{r}$。

证明这个结论对局部线性回归（第 6.1.1 节）也成立。

因此可说明算法 9.1 中步骤（2）第二行的调节过程不是必需的。

练习 9.2

令 $\mathbf{A}$ 代表一个已知的 $k\times k$ 矩阵，$\mathbf{b}$ 代表一个已知的 $k$ 维向量，并且 $\mathbf{z}$ 代表一个未知的 $k$ 维向量。求解方程组 $\mathbf{A}\mathbf{z}=\mathbf{b}$ 的高斯-赛德尔（Gauss-Seidel）算法是一个迭代方法，将所有其他 $z_j$ 的值固定在它们的当前步骤取值上，然后在方程 $j$ 中求解元素 $z_j$。重复这个过程 $j=1,2,\dots,k,1,2,\dots,k\dots$，直到结果收敛。（Golub and Van Loan, 1983）。

1. 考虑一个有 $p$ 项和 $N$ 个样本的加性模型，它的第 $j$ 项会使用一个线性平滑器 $\mathbf{S}_j$ 来拟合。以下为一个方程组：

   $$\begin{pmatrix} \mathbf{I} & \mathbf{S}_1 & \mathbf{S}_1 & \cdots & \mathbf{S}_1 \\ \mathbf{S}_2 & \mathbf{I} & \mathbf{S}_2 & \cdots & \mathbf{S}_2 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{S}_p & \mathbf{S}_p & \mathbf{S}_p & \cdots & \mathbf{I} \end{pmatrix} \begin{pmatrix} \mathbf{f}_1 \\ \mathbf{f}_2 \\ \vdots \\ \mathbf{f}_p \end{pmatrix}= \begin{pmatrix} \mathbf{S}_1\mathbf{y} \\ \mathbf{S}_2\mathbf{y} \\ \vdots \\ \mathbf{S}_p\mathbf{y} \end{pmatrix} \tag{9.33}$$

   其中的每个 $\mathbf{f}_j$ 都是函数 $j$ 在数据样本点上取值的一个 $N$ 维向量，而 $\mathbf{y}$ 是 $N$ 维的数据变量取值向量。证明回修算法就是求解这个方程组的分块高斯-赛德尔算法。

2. 令 $\mathbf{S}_1$ 和 $\mathbf{S}_2$ 为特征值在 $[0,1)$ 的对称平滑算子（矩阵）。考虑一个在输出变量 $\mathbf{y}$ 和平滑器 $\mathbf{S}_1$、$\mathbf{S}_2$ 上的回修算法。证明在任意的初始值条件下，算法都会收敛，并给出最后一次迭代的表达式。

练习 9.3

回修方程组

Consider a backfitting procedure with orthog- onal projections, and let D be the overall regression matrix whose columns span V = Lcol (S1 ) ⊕ Lcol (S2 ) ⊕ · · · ⊕ Lcol (Sp ), where Lcol (S) denotes the column space of a matrix S. Show that the estimating equations

are equivalent to the least squares normal equations DT Dβ = DT y where β is the vector of coefficients.

练习 9.4

Suppose the same smoother S is used to estimate both terms in a two-term additive model (i.e., both variables are identical). Assume that S is symmetric with eigenvalues in [0, 1). Show that the backfitting residual converges to (I + S)−1 (I − S)y, and that the residual sum of squares con- verges upward. Can the residual sum of squares converge upward in less structured situations? How does this fit compare to the fit with a single term fit by S? [Hint: Use the eigen-decomposition of S to help with this comparison.]