详解广义线性模型的来龙去脉

广义线性模型[generalize linear model]是线性模型的扩展，通过联结函数建立响应变量的数学期望值与线性组合的预测变量之间的关系。其特点是不强行改变数据的自然度量，数据可以具有非线性和非恒定方差结构。是线性模型在研究响应值的非正态分布以及非线性模型简洁直接的线性转化时的一种发展。

1.指数分布族

指数族分布是指概率密度函数（或是离散型随机变量的概率分布）具有特定形式的一类分布的集合，其概率密度函数的形式如下：
$f(y_i,{\theta }_i,\varphi )=\exp \left\{\frac{y_i{\theta }_i-b({\theta }_i)}{a(\varphi )}+h(y_i,\varphi )\right\}$
常见的指数分布族有高斯分布，伯努利分布，多项分布，泊松分布，$\Gamma$分布等。分别以高斯分布和伯努利分布为例来说明这个问题。
高斯分布：
$f(y_i,{\mu }_i,\sigma )=\exp \left\{-\frac{(y_i-{\mu }_i{)}^2}{2{\sigma }^2}-\ln (\sqrt{2\pi }\sigma )\right\}=\exp \left\{\frac{y_i{\mu }_i-\frac{{\mu }_i^2}{2}}{{\sigma }^2}-\frac{y_i^2}{2\sigma }-\ln (\sqrt{2\pi }\sigma )\right\}$
可以看出，${\theta }_i={\mu }_i,b({\theta }_i)=\frac{{\mu }_i^2}{2},a(\varphi )={\sigma }^2,h(y_i,\varphi )=-\frac{y_i^2}{2\sigma }-\ln (\sqrt{2\pi }\sigma )$
伯努利分布：
$p(y_i,{\pi }_i)={\pi }_i^{y_i}(1-{\pi }_i{)}^{1-y_i}=\exp \left\{y_i\ln \frac{{\pi }_i}{1-{\pi }_i}+\ln (1-{\pi }_i)\right\}$
显然，${\theta }_i=\ln \frac{{\pi }_i}{1-{\pi }_i},b({\theta }_i)=-\ln (1-{\pi }_i),a(\varphi )=1,h(y_i,\varphi )=0$
指数族分布的概率密度函数在$\mathbb{R}$上的积分应当为1，根据这个性质，进一步地能推出更多普适结论，比如，对密度函数关于$\theta$求导：
$0=\frac{\partial }{\partial {\theta }_i}{\int }_{\mathbb{R}}f(y_i,{\theta }_i,\varphi )d{y}_i={\int }_{\mathbb{R}}\frac{\partial }{\partial {\theta }_i}f(y_i,{\theta }_i,\varphi )d{y}_i=\frac{1}{a(\varphi )}{\int }_{\mathbb{R}}\left(y_i-\frac{db({\theta }_i)}{d{\theta }_i}\right)f(y_i,{\theta }_i,\varphi )d{y}_i=\frac{1}{a(\varphi )}\left({\mu }_i-\frac{db({\theta }_i)}{d{\theta }_i}\right)$
得到指数分布族的期望${\mu }_i=\frac{db({\theta }_i)}{d{\theta }_i}$;在此基础上再对${\theta }_i$求一次导数：
$0=\frac{{\partial }^2}{\partial {\theta }_i^2}{\int }_{\mathbb{R}}f(y_i,{\theta }_i,\varphi )d{y}_i=\frac{1}{[a(\varphi ){]}^2}\left(E(y_i^2)-{\mu }_i^2-a(\varphi )\frac{d^{2}b({\theta }_i)}{d{\theta }_i^2}\right)$
即$E(y_i^2)={\mu }_i^2+a(\varphi )\frac{d^{2}b({\theta }_i)}{d{\theta }_i^2}$,再由方差公式知 $\text{var}(y_i)=a(\varphi )\frac{d^{2}b({\theta }_i)}{d{\theta }_i^2}=a(\varphi )\frac{d{\mu }_i}{d{\theta }_i}$
整理一下上面所做的工作：通过推导，我们得到了指数分布族期望与方差的计算公式
${\mu }_i=\frac{db({\theta }_i)}{d{\theta }_i}$
$\text{var}(y_i)=a(\varphi )\frac{d^{2}b({\theta }_i)}{d{\theta }_i^2}=a(\varphi )\frac{d{\mu }_i}{d{\theta }_i}$
这两个公式作为指数分布族的重要性质将会在下文用到，至于指数分布族的其他一些性质，比如充分性和完备性，由于超出本文范畴，此处不再介绍。

2.广义线性模型

普通线性模型具有的模型形式为：$y_i={x\text{x}x}_i^T\beta \text{β}\beta +{\epsilon }_i$;建立普通线性模型是基于常数方差，独立性和高斯假设前提；如果误差项不满足常数方差和独立性假设，那么我们已经有一些以广义最小二乘法（GLS）为代表的解决方案；但是，如果数据不满足常数方差及高斯分布的假设呢？比如，如果响应变量 y i ∈ { 0 , 1 } y_i\in\left\{0,1\right\} yi​∈{0,1},这是一个只含二值的问题，我们建立模型：$y_i=E(y_i|{x\text{x}x}_i,\beta \text{β}\beta )+{\epsilon }_i$则误差项显然是离散型随机变量，而且它的方差随期望的变化而变化，当然，我们仍然可以用线性模型来拟合数据，然后用GLS学习参数取值，不过由于真实数据不再具有线性结构，所以由GLS强行用线性函数拟合得到的估计量不再是BLUE估计量，另一方面，这样得到的模型鲁棒性不够的，容易受新样本的影响，而且预测时，由线性模型计算出的预测值可能远远大于1或者小于0.
考虑到线性模型在解决非线性结构数据时存在的这些问题，可以利用一个连接函数$g(\cdot )$,然后建立模型$y_i=g({x\text{x}x}_i^T\beta \text{β}\beta )+{\epsilon }_i$,比如对于上文提到的二值问题，可以取$g(x)=\text{sigmoid}(x)$.指数分布族为这类广义线性模型（GLM）中连接函数的寻找提供了一个一般的思路。
为了说明这个问题，先以最简单的普通线性模型为例。当响应变量服从高斯分布时，我们假设响应变量的期望是回归变量的线性函数，即${\mu }_i={x\text{x}x}_i^T\beta \text{β}\beta$，将高斯分布的密度函数化成指数族的形式后可以看到${\theta }_i={\mu }_i$,即${\theta }_i={x\text{x}x}_i^T\beta \text{β}\beta$，将这一假设拓展到广义线性模型中也是可行的，因为在文章的第一部分已经推导出：${\mu }_i=\frac{db({\theta }_i)}{d{\theta }_i}$，这说明${\mu }_i$的取值由${\theta }_i$唯一确定而与尺度化参数$\varphi$无关，于是，可以取连接函数$g({\theta }_i)=\frac{db({\theta }_i)}{d{\theta }_i}$,根据线性假设有${\mu }_i=g({x\text{x}x}_i^T\beta \text{β}\beta )$,所以我们可以建立模型$y_i=g({x\text{x}x}_i^T\beta \text{β}\beta )+{\epsilon }_i$.
以上文中提到的二值问题为例，我们来建立它的广义线性模型。对于二值问题最常用的假设就是认为它服从伯努利分布，根据第一部分对伯努利分布指数族形式的推导有${\theta }_i=\ln \frac{{\pi }_i}{1-{\pi }_i}$，由于伯努利分布的期望等于取1的概率，所以${\theta }_i=\ln \frac{{\mu }_i}{1-{\mu }_i}$，再由${\theta }_i={x\text{x}x}_i^T\beta \text{β}\beta$，得${\mu }_i=\frac{1}{1+\exp (-{x\text{x}x}_i^T\beta \text{β}\beta )}$,所以合适模型可以为$y_i=\frac{1}{1+\exp (-{x\text{x}x}_i^T\beta \text{β}\beta )}+{\epsilon }_i$,这就得到了经典的logstic回归模型。
相信读者通过以上过程已经清楚了应该如何建立一个广义线性模型，模型建立之后，接下来的问题就是参数的学习了，广义线性模型用到负对数似然损失函数学习参数的最大似然估计量（MLE）,使用牛顿—拉普森算法进行数值求解，可以推导出GLM参数估计的迭代加权最小二乘（IRLS）算法。

3.迭代加权最小二乘（IRLS）算法

对数似然函数
$\ln L(\beta \text{β}\beta )=\frac{1}{a(\varphi )}\sum _i\left(y_i{\theta }_i-b({\theta }_i)\right)+\sum _{i}h(y_i,\varphi )$
令对数似然函数关于$\beta \text{β}\beta$的导数等于0
$\frac{\partial \ln L(\beta \text{β}\beta )}{\partial \beta \text{β}\beta }=\sum _i\frac{\partial \ln L(\beta \text{β}\beta )}{\partial {\theta }_i}\frac{\partial {\theta }_i}{\partial \beta \text{β}\beta }=\frac{1}{a(\varphi )}\sum _i\left(y_i-{\mu }_i\right){x\text{x}x}_i=0\text{0}0$
可以将方程写成矩阵的形式
${X\text{X}X}^T(y\text{y}y-\mu \text{μ}\mu )=0\text{0}0$
这个方程称为得分方程,如果是线性模型，即$\mu \text{μ}\mu =X\beta \text{Xβ}X\beta$，这时得分方程等价于正规方程，于是很容易对这个方程求解。但是对一般的广义线性模型，得分方程关于参数$\beta \text{β}\beta$是一个非线性方程组，因此只能用近似逼近的方法求解。
为了求解得分方程，我们希望用一个线性表达式来代替它的非线性部分，可以取${\eta }_i$满足：
$y_i-{\mu }_i=\frac{d{\mu }_i}{d{\theta }_i}({\eta }_i-{\theta }_i)$
根据第一部分的推导
$\frac{d{\mu }_i}{d{\theta }_i}=\frac{\text{var}(y_i)}{a(\varphi )}$
于是得分方程可化为
$\sum _i\frac{\text{var}(y_i)}{a(\varphi )}\left({\eta }_i-{\theta }_i\right){x\text{x}x}_i=0\text{0}0$
令$V\text{V}V=\frac{1}{a(\varphi )}diag(\text{var}(y_1),\cdots ,\text{var}(y_n))$,则
${X\text{X}X}^{T}V\text{V}V(\eta \text{η}\eta -X\beta \text{Xβ}X\beta )=0\text{0}0$
很明显$\widehat{\beta \text{β}\beta }=({X\text{X}X}^{T}VX\text{VX}VX{)}^{-1}{X\text{X}X}^{T}V\text{V}V\eta \text{η}\eta$. all right…看上去似乎已经估计出参数取值了，不过好像又有点不太对劲。仔细检查发现我们并不知道$V\text{V}V$和$\eta \text{η}\eta$，注意到${\eta }_i={\theta }_i+\frac{d{\theta }_i}{d{\mu }_i}(y_i-{\mu }_i)$，要求$V\text{V}V$和$\eta \text{η}\eta$，就必须得知道参数$\beta \text{β}\beta$的取值，这听起来令人沮丧，不过好在有迭代法这一神器，不管收不收敛，试一试总不会有问题。事实证明，只要初值选取合适，在本问题上使用迭代法是可行的。说了这么多废话，现将算法描述如下：
输入：样本$({x\text{x}x}_1,y_1),\cdots ,({x\text{x}x}_n,y_n)$
输出：参数估计$\widehat{\beta \text{β}\beta }$
1.初始化：$\widehat{\beta \text{β}\beta }\leftarrow {\beta \text{β}\beta }_0$；
2.for i in 1:n，计算：
${\hat{\theta }}_i\leftarrow {x_i\text{xi}{x}_i}^T\widehat{\beta \text{β}\beta }$,
${\hat{\mu }}_i\leftarrow g({\hat{\theta }}_i)$,
$V_{ii}\leftarrow \frac{{\text{var}}_{y_i|{\hat{\theta }}_i}(y_i)}{a(\varphi )}$;
3.计算：$\eta \text{η}\eta \leftarrow \widehat{\theta \text{θ}\theta }+{V\text{V}V}^{-1}(y\text{y}y-\widehat{\mu \text{μ}\mu })$；
4.更新：$\widehat{\beta \text{β}\beta }\leftarrow ({X\text{X}X}^{T}VX\text{VX}VX{)}^{-1}{X\text{X}X}^{T}V\text{V}V\eta \text{η}\eta$；
5.重复2到5步，直至算法收敛.

4.IRLS为什么有效？

不少读者读到这里可能会产生一个疑惑，那就是在推导IRLS算法过程中一个关键步骤$y_i-{\mu }_i=\frac{d{\mu }_i}{d{\theta }_i}({\eta }_i-{\theta }_i)$似乎显得有些天马行空，${\eta }_i-{\theta }_i$是关于$\beta \text{β}\beta$的线性函数通过变换后可以“解出”$\beta \text{β}\beta$，这点挺好理解，但是为什么前面的系数是$\frac{d{\mu }_i}{d{\theta }_i}$而不是其他，要知道这个系数直接关系到权重矩阵$V\text{V}V$,取不同的系数参数的估计结果就会有不同的值。直观上看，这个表达式形式上类似于连接函数在${\theta }_i$处的一阶泰勒展开，但是$g({\eta }_i)$一般是不等于$y_i$的，否则的话我们就可以认为这一步是利用泰勒公式取的一个线性近似。实际上IRLS可以看做由牛顿-拉普森算法推导出来的。
下面我们利用牛顿-拉普森算法来求解得分方程，具体的过程就不写了（网上一搜一堆介绍牛顿算法的，也可以参考我之前的博文），直接写出算法的更新步骤：
$\widehat{\beta \text{β}\beta }\leftarrow \widehat{\beta \text{β}\beta }+({X\text{X}X}^{T}VX\text{VX}VX{)}^{-1}{X\text{X}X}^T(y\text{y}y-\widehat{\mu \text{μ}\mu })$
在IRLS中将$\eta \text{η}\eta \leftarrow \widehat{\theta \text{θ}\theta }+{V\text{V}V}^{-1}(y\text{y}y-\widehat{\mu \text{μ}\mu })$代入到更新步骤$\widehat{\beta \text{β}\beta }\leftarrow ({X\text{X}X}^{T}VX\text{VX}VX{)}^{-1}{X\text{X}X}^{T}V\text{V}V\eta \text{η}\eta$中，利用$\widehat{\theta \text{θ}\theta }=X\text{X}X\widehat{\beta \text{β}\beta }$同样可以得到
$\widehat{\beta \text{β}\beta }\leftarrow \widehat{\beta \text{β}\beta }+({X\text{X}X}^{T}VX\text{VX}VX{)}^{-1}{X\text{X}X}^T(y\text{y}y-\widehat{\mu \text{μ}\mu })$
这就说明了IRLS算法和Newton-Raphson算法是等价的。

5.IRLS算法与GLS算法

记$W\text{W}W=diag(\text{var}({\eta }_1),\cdots ,\text{var}({\eta }_n))$,根据${\eta }_i$的定义
$\text{var}({\eta }_i)=\text{var}\left[\frac{a(\varphi )}{\text{var}(y_i)}(y_i-{\mu }_i)\right]=\frac{a(\varphi {)}^2}{\text{var}(y_i)}$
故$V\text{V}V=a(\varphi ){W\text{W}W}^{-1}$
${\widehat{\beta \text{β}\beta }}_{IRLS}=({X\text{X}X}^T{W\text{W}W}^{-1}X\text{X}X{)}^{-1}{X\text{X}X}^T{W\text{W}W}^{-1}\eta \text{η}\eta$
如果直接对广义线性模型的线性部分用广义最小二乘法学习参数结果会是怎样呢？
记${\theta }_i^{\ast }=g^{-1}(y_i)$,则$({x\text{x}x}_1,{\theta }_1^{\ast }),\cdots ,({x\text{x}x}_n,{\theta }_n^{\ast })$可以看成线性模型$\theta ={x\text{x}x}^T\beta \text{β}\beta +\epsilon$的$n$个观测，利用GLS算法估计$\beta \text{β}\beta$,记$U\text{U}U=diag(\text{var}({\theta }_1^{\ast }),\cdots ,\text{var}({\theta }_n^{\ast }))$,则：
${\widehat{\beta \text{β}\beta }}_{GLS}=({X\text{X}X}^T{U\text{U}U}^{-1}X\text{X}X{)}^{-1}{X\text{X}X}^T{U\text{U}U}^{-1}{\theta }^{\ast }\text{θ∗}{\theta }^{\ast }$
将$g^{-1}(y_i)$在${\mu }_i$处一阶泰勒展开有
${\theta }_i^{\ast }\approx {\theta }_i+\frac{\partial {\theta }_i}{\partial {\mu }_i}(y_i-{\mu }_i)={\eta }_i$
故${\widehat{\beta \text{β}\beta }}_{IRLS}\approx {\widehat{\beta \text{β}\beta }}_{GLS}$，即IRLS与GLS是渐进等价的。

6.后记

GLM的好处在于建立模型时结合了响应变量的原始分布信息，其响应变量的分布必须是指数族分布，它是当响应变量不满足正态分布和方差不变假设时，区别于GLS的另一种变量变换方法，当然，这种方法也不是万能的，面对具体问题时，GLM只是提供了一种建模思路，而模型是否适用还得做进一步检验，很多情况下，GLM对数据拟合效果不佳，这时候应该尝试使用其他模型。