pyGAM相关知识

pyGAM相关知识

0.随笔序

第一次写博客，尝试一下用csdn做笔记。
最近做毕设要学习可解释的模型相关知识。正好学到了pyGAM包，所以写一下试试。

1.GAMs

背景

学习过机器学习的同学应该对线性回归与逻辑回归都比较熟悉了。GAM和他们两个有很多的相似之处。

1.1 线性回归

我们都知道线性回归的基本形式：
$${\hat{y}}^{(i)}={\beta }_0+{\beta }_1{x}_1^{(i)}+\dots +{\beta }_p{x}_p^{(i)}$$
其中${\hat{y}}^{(i)}$为模型对样本$i$的预测值，${\beta }_0$为模型截距，${\beta }_p$为每个特征$x_p^{(i)}$的系数，我们也称其为特征的权重。
这很好理解，我们选定一个特征$x_p^{(i)}$，将其他的特征值固定，只将该特征增大一个单位，那么模型的预测值${\hat{y}}^{(i)}$就会变化其相应的系数值${\beta }_p$。
换句话说，${\beta }_p$的大小可以代表该特征对模型预测值的影响程度的大小。

1.2 逻辑回归

让我们再看看逻辑回归。
线性回归可以很好的解决很多回归问题，但是它不能解决分类问题。
拿二分类为例，我们可以把两类分别标记为0，1，并使用线性回归。从技术上讲，这是可行的，你确实能够得到一个模型，但这种方法存在着一个问题：
线性回归输出的不是概率，而是一个连续的预测值。 你不能简单地就因为它的值大于0.5就把它分类为1，小于0.5就分类为0。
而且线性回归不能扩展到多分类问题中，你不可避免地会给其他类标记上2，3等等数值，这些标签原本之间不存在着大小上的联系，但是当我们将其标记上了数值后，我们会将其强加上一个数值大小的相关性，这对我们的模型是有一定影响的。
由此，我们就要引入逻辑函数来构建逻辑回归。逻辑函数：
$$\text{logistic}(\eta )=\frac{1}{1+exp(-\eta )}$$
此式可以将线性方程的输出收缩到0与1之间。
我们将线性回归的方程带入$\eta$。如下式：
$$P(y^{(i)}=1)=\frac{1}{1+exp(-({\beta }_0+{\beta }_1{x}_1^{(i)}+\dots +{\beta }_p{x}_p^{(i)}))}$$
等式右侧为我们的逻辑函数，等式左侧为样本数据$i$被分类为1的概率。
为了便于观察，我们将公式进行变换。
$$log\left(\frac{P(y=1)}{1-P(y=1)}\right)=log\left(\frac{P(y=1)}{P(y=0)}\right)={\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_p{x}_p$$
我们想办法将线性项单独提取出来，得到上式。我们可以发现线性项之和就是二分类概率比的对数。我们将log()函数中的项称为“odds(发生比)”(事件的概率除以没有事件的概率)，并将其封装在对数中，称为log(odds)。

该公式表明，logistic回归模型是对数概率的线性模型。稍微调整一下方程，您就可以知道当某个局部特征$x_p$改变1个单位时，预测会发生什么变化。为此，我们可以首先将exp()函数应用到等式两边:
$$\frac{P(y=1)}{1-P(y=1)}=odds=exp\left({\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_p{x}_p\right)$$
然后比较其中一个特征值增加1时的情况。但我们不是看它们的区别，而是看这两种预测的比例:
$$\frac{odd{s}_{x_j+1}}{odds}=\frac{exp\left({\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_j(x_j+1)+\dots +{\beta }_p{x}_p\right)}{exp\left({\beta }_0+{\beta }_1{x}_1+\dots +{\beta }_j{x}_j+\dots +{\beta }_p{x}_p\right)}$$
简化一下右式：
$$\frac{odd{s}_{x_j+1}}{odds}=exp\left({\beta }_j(x_j+1)-{\beta }_j{x}_j\right)=exp\left({\beta }_j\right)$$
最后，我们得到了像exp()这样简单的特征权重。一个特征变化一个单位,比值比(odds)改变$\exp ({\beta }_j)$倍(乘法)。

1.3 GAM

线性回归模型最大的优点也是最大的缺点是预测被建模为特征的加权和。
特征效应是加性的，意味着特征间没有交互作用；关系是线性的，意味着一个单位特征的增加可以直接转化为预测结果的增加/减少。
线性模型允许我们将一个特征和预期结果之间的关系压缩成一个单一的数字，即估计的权重。但是简单的加权和对于现实世界中的许多预测问题来说限制太大。

接下来我们将简单讲一下关于经典线性回归模型的三个问题和解决方法。

问题: 给定特征的目标输出y不服从高斯分布。

例如: 假设我想要预测某一天我骑自行车的时间。作为特征，我有日期，天气等等。如果我用一个线性模型，它是可以预测到负分钟的，因为它假设了一个高斯分布，它不会在0分钟停止。如果我想用线性模型预测概率，我可以得到负的或大于1的概率。

解决方案: 广义线性模型(GLMs)。 Generalized Linear Models (GLMs).

问题: 特征之间相互影响。

例如： 一般来说，小雨对我骑车的欲望有轻微的负面影响。但是在夏天，高峰时间，我欢迎下雨，因为那时所有的好天气骑自行车的人都呆在家里，我可以拥有自己的自行车道!这是一种时间和天气之间的相互作用，不能用纯粹的相加模型来描述。

解决方案: 手动添加相关性。 Adding interactions manually.

问题: 特征和y之间的真实关系不是线性的。

例如: 在0到25摄氏度之间，温度对我骑自行车欲望的影响可能是线性的，这意味着从0到1摄氏度的升高与从20到21摄氏度的升高产生的骑车欲望是一样的。但是在较高的温度下，我骑车的动力减弱了，甚至下降了——我不喜欢在太热的时候骑车。

解决方案: 广义加性模型(GAMs);特征转换。Generalized Additive Models (GAMs); transformation of
features.

这里我们主要针对第三个问题：特征和y之间的真实关系不是线性的。
这种情况在真实世界中其实是很常见的。

世界不是线性的。 线性模型中的线性意味着无论一个实例在一个特定的特征中具有什么值，将该值增加一个单位对预测结果的影响总是相同的。假设气温在10摄氏度的情况下升高1度，对出租自行车数量的影响与气温已经40度相同时升高1度相同，这合理吗?直觉上，人们认为把温度从10摄氏度提高到11摄氏度会对自行车租赁产生积极的影响，而把温度从40摄氏度提高到41摄氏度则会产生消极的影响，正如你在书中看到的许多例子一样。温度特征对租用自行车的数量有一个线性的、积极的影响，但在某一时刻它会变平，甚至在高温下有一个消极的影响。线性模型并不关心，它会尽职地找到最好的线性平面(通过最小化欧几里得距离)。

我们可以通过广义可加性模型(GAMs)来引入非线性关系。

为什么不“简单”地允许(广义)线性模型学习非线性关系?这就是GAMs出现的背后动机。GAMs放宽了线性模型的特征关系必须是一个简单加权和的限制，取而代之的是假设结果可以由每个特征的任意函数的和来建模。数学上，GAM的关系是这样的:
$$g(E_Y(y|x))={\beta }_0+f_1(x_1)+f_2(x_2)+\dots +f_p(x_p)$$
GAM的核心仍然是特征效果的总和，但是您可以允许某些特征和输出之间拥有非线性关系。线性效应也包含在框架中，因为对于要被线性处理的特征们，你可以通过仅仅采取${\beta }_p{x}_p$的形式限制他们的函数，即，使$f_p(x_p)$=${\beta }_p{x}_p$。

问题是我们如何构建非线性函数。答案是使用“样条函数”。样条函数是可以用来近似任意函数的函数。

例如，为了使用样条函数来建模之前自行车问题的温度特征，我们就从数据中删除了温度特征，并将其替换为4列，每一列表示一个样条函数。通常在实际工作中你会用到更多的样条函数，我只是为了说明而减少了数量。这些新样条特征的每个实例的值$f_p^{(i)}(x_p)$取决于实例的温度值$x_p$。连同所有的线性效应，GAM也会估计这些样条的权值。GAMs还引入了一个惩罚项来保持权重接近于零。这有效地降低了样条的灵活性，减少了过度拟合。然后，通过交叉验证调整通常用于控制曲线灵活性的平滑度参数。忽略惩罚项，用样条进行非线性建模是一项奇特的特征工程。

在我们仅使用温度特征的GAM来预测自行车数量的例子中，模型特征矩阵如下所示:

(Intercept)	s(temp).1	s(temp).2	s(temp).3	s(temp).4
1	0.93	-0.14	0.21	-0.83
1	0.83	-0.27	0.27	-0.72
1	1.32	0.71	-0.39	-1.63
1	1.32	0.70	-0.38	-1.61
1	1.29	0.58	-0.26	-1.47
1	1.32	0.68	-0.36	-1.59

每一行表示数据中的一个单独实例(一天)。每个样条列包含在特定温度值下的样条函数值。下图显示了这些样条函数的样子:

为了平滑地模拟温度效应，我们使用了4个样条函数。每个温度值映射到(这里)4个样条值。如果一个实例的温度为30°C，则第一个样条特征值为-1，第二个为0.7，第三个为-0.8，第四个为1.7。

GAM为每个温度样条特征赋予权重:

	weight
(Intercept)	4504.35
s(temp).1	-989.34
s(temp).2	740.08
s(temp).3	2309.84
s(temp).4	558.27

实际的曲线是由样条函数的和与估计的权值加权得到的，它是这样的:

温度的GAM特征效应，用于预测租用自行车的数量(温度作为唯一的特征)。

平滑效果的解释需要对拟合曲线进行视觉检查。样条曲线通常以均值预测为中心，因此曲线上的点是与均值预测的差值。例如，在0摄氏度时，自行车的预测数量比平均预测低3000辆。

2.pyGAM导读

2.1 简介

广义加性模型(GAM)是以下形式的平滑半参数模型
$$g(\mathbb{E}[y|X])={\beta }_0+f_1(X_1)+f_2(X_2,X3)+\dots +f_M(X_N)$$
其中 X.T = [X_1, X_2, ..., X_N] 是自变量, y 是因变量, g() 是连接我们预测变量和因变量预期值的连接函数。
特征函数f_i()有惩罚过的B样条函数构造，该样条函数允许我们自动地建模非线性关系，不需要再手动为每个变量做血多不同的转换。

GAMs通过在保持可加性的同时，允许特征的非线性函数，来扩展线性模型。由于模型是可加的，我们很容易就可以解释每个特征X_i在其他所有预测常量保持不变的情况下单独对y的影响。
总结起来就是一个非常灵活的模型，可以很容易地合并先验知识和控制过度拟合。

2.2 GAMs模型实体

$$y\sim ExponentialFamily(\mu |X)$$
y服从的是指数族的分布。其中
$$g(\mu |X)={\beta }_0+f_1(X_1)+f_2(X_2,X3)+\dots +f_M(X_N)$$
由此我们可以找到GAM模型的三要素

• distribution —来自指数族的分布
• link function—连接函数$g(\cdot )$
• functional form—可加性的模型结构${\beta }_0+f_1(X_1)+f_2(X_2,X3)+\dots +f_M(X_N)$

分布

用GAM(distribution='...')来选定：

• 'normal'
• 'binomial'
• 'poisson'
• 'gamma'
• 'inv_gauss'

连接函数

用GAM(link='...')来选定
连接函数将分布均值用于线性预测。到目前为止，以下内容可用：

• 'identity'
• 'logit'
• 'inverse'
• 'log'
• 'inverse-squared'

函数形式

用GAM(terms=...)或者GAM(...)来选定：
在pyGAM中，我们用以下几个术语来选定函数形式：

• l()线性形式
• s()样条函数
• f()因素项（factor terms？）
• te()张量积
• intercept
  利用以上这些知识，我们可以快速的构建一个模型：

from pygam import GAM, s, te

GAM(s(0, n_splines=200) + te(3,1) + s(2), distribution='poisson', link='log')

结果为

GAM(callbacks=[‘deviance’, ‘diffs’], distribution=‘poisson’, fit_intercept=True, link=‘log’, max_iter=100, terms=s(0) + te(3, 1)+ s(2), tol=0.0001, verbose=False)

这里我们选择的为：

• 针对特征0的样条函数，有200个基本函数(s(0, n_splines=200))
• 特征1和特征3上的张量样条相互作用(te(3,1))
• 针对特征2的样条函数(s(2))

注意：
GAM(..., intercept=True)，这样模型会默认包含一个截距项。

实际

在pyGAM中，你可以依赖这三个元素来构建你自己的模型，或者你也可以从以下几个公共模型中选择：

• LinearGAM 恒等连接和正态分布（identity link and normal distribution）
• LogisticGAM 对数连接和二项分布（logit link and binomial distribution）
• PoissonGAM 对数链接和泊松分布（log link and Poisson distribution）
• GammaGAM 对数连接和伽马分布（log link and gamma distribution）
• InvGauss 对数连接和反高斯分布（log link and inv_gauss distribution）

公共模型的好处就在于，出了减少样板代码量以外，还具有一些额外的功能。

术语和交互影响

pyGAM可以通过使用te()来使用张量积

from pygam import PoissonGAM, s, te
from pygam.datasets import chicago

X, y = chicago(return_X_y=True)

gam = PoissonGAM(s(0, n_splines=200) + te(3, 1) + s(2)).fit(X, y)

画一个3D图

import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

plt.ion()
plt.rcParams['figure.figsize'] = (12, 8)

XX = gam.generate_X_grid(term=1, meshgrid=True)
Z = gam.partial_dependence(term=1, X=XX, meshgrid=True)

ax = plt.axes(projection='3d')
ax.plot_surface(XX[0], XX[1], Z, cmap='viridis')

对于一些比较简单的交互影响，有时加一条计量语句到命令里就很有用了。

from pygam import LinearGAM, s
from pygam.datasets import toy_interaction

X, y = toy_interaction(return_X_y=True)

gam = LinearGAM(s(0, by=1)).fit(X, y)
gam.summary()

回归

由于毕设主要是搞得面向分类数据。这部分以后再补充。

分类

对于二分类问题来说，我们可以用一个logistic GAM，其模型为：
$$log\left(\frac{P(y=1|X)}{P(y=0|X)}\right)={\beta }_0+f_1(X_1)+f_2(X_2,X3)+\cdots +f_M(X_N)$$
代码如下：

from pygam import LogisticGAM, s, f
from pygam.datasets import default

X, y = default(return_X_y=True)

gam = LogisticGAM(f(0) + s(1) + s(2)).gridsearch(X, y)

fig, axs = plt.subplots(1, 3)
titles = ['student', 'balance', 'income']

for i, ax in enumerate(axs):
    XX = gam.generate_X_grid(term=i)
    pdep, confi = gam.partial_dependence(term=i, width=.95)

    ax.plot(XX[:, i], pdep)
    ax.plot(XX[:, i], confi, c='r', ls='--')
    ax.set_title(titles[i]);

检查准确度：

gam.accuracy(X, y)

0.9739

本文参考文献：
[1]Christoph Molnar.“Interpretable Machine Learning”.https://christophm.github.io/interpretable-ml-book/index.html
[2]pyGAM的GitHub文档：https://pygam.readthedocs.io/en/latest/notebooks/tour_of_pygam.html

本文严禁转载，只用交流学习，不作商用。