广义线性模型（GLIM，Generalized Linear Model）

在本篇文章中，我们将按照以下流程来介绍广义线性模型。

• 指数族和最大熵模型
• 基于指数族的广义线性模型——广义线性模型的构成、极大似然估计和求解算法
• 广义线性模型的偏差和分析

• 广义线性模型的特征选择

1、指数族和最大熵模型

1.1 指数族的形式

指数族是概率统计中最重要的一类分布族。具有以下的一般形式：

$$f(y)=exp[\{\frac {\theta y-b(\theta)}{a(\phi)}\}+c(y,\phi)]$$

这里指数族标准型相较于我们常见的指数族的形式有所不同。之前我们接触到的指数族的标准形式如下：

$$f(y|\theta)=h(y)c(\theta )exp\{\ \sum_{i=1}^k\omega_i(\theta)t_i(y)\}$$

这里我们将$\theta$和$\phi$作为两类不同的参数加以区分。$\theta$反映指数族随机变量的集中趋势,$\phi$反映指数族随机变量的离散趋势。

我们称$\theta$为指数族的典型参数($Canonical\ Parameter$)，$\phi$为指数族的离散参数($Dispersion \ Parameter$)，在给定指数族分布时，函数$a(.),b(.),c(.,.)$是已知函数。

1.2 指数族的性质

设$l(\theta)$为对数似然函数，指数分布族有以下的性质：

• $E(l'(\theta))=0; E(Y)=b'(\theta)$
• $-E(l''(\theta))=E(l'(\theta)^2); Var(Y)=a(\phi)b''(\theta)$

证明：
(1) $$l(\theta)={\frac {\theta y-b(\theta)}{a(\phi)}}+c(y,\phi)$$

$$E(l'(\theta))=E(\frac{\partial f(y;\theta)}{\partial \theta} /f(y;\theta))=\int_{y\in D} \frac{\partial}{\partial \theta} f(y;\theta)dy$$

由莱布尼茨法则，这里的积分和求导可以互换：

$$E(l'(\theta))=\frac{\partial }{\partial \theta} \int_{y\in D} f(y;\theta)dx=0$$

故有：

$$E( l'(\theta)) =E[\frac{(y-b'(\theta))}{a(\phi)}]=0$$

$$E(y)=b'(\theta)$$

(2) $$\frac{\partial}{\partial \theta} E(l'(\theta))=\frac{\partial}{\partial \theta} \int l'(\theta) f_y(y;\theta)dy$$
由莱布尼茨定理交换积分求导次序：

$$\int \frac{\partial}{\partial \theta} l'(\theta) f_y(y;\theta)dy=\int l''(\theta)f_y(y;\theta)dy+\int l'(\theta)\frac{\partial}{\partial \theta}f_y(y;\theta)dy$$

$$=\int l''(\theta)f_y(y;\theta)dy+\int l'(\theta)[\frac{\partial f_y(y;\theta)}{\partial \theta}/f_y(y;\theta)]f_y(y;\theta)dy$$

$$=\int l''(\theta)f_y(y;\theta)dy+\int l'(\theta)^2f_y(y;\theta)dy$$

$$=E(l''(\theta))+E(l'(\theta)^2)=0$$

又有：

$$-E(l''(\theta))=-E[\frac{\partial}{\partial \theta}\frac{(y-b'(\theta))}{a(\phi)}]=E(\frac{b''(\theta)}{a(\phi)})=\frac{b''(\theta)}{a(\phi)}$$

$$=E(l'(\theta)^2)=E(\frac{(y-b'(\theta))^2}{a(\phi)^2})=\frac{Var(y)}{a(\phi)^2}$$

故有：$Var(Y)=a(\phi)b''(\theta)$

下面给出后文构建广义线性模型时候的重要定义。由指数族性质知，函数$b'(.)$将前文提到的典型参数$\theta$和Y的期望$E(Y)=\mu$联系在一起，当$b'(.)$存在逆函数的时候，令$g(.)=(b')^{-1}(.)$有：

$$b'(\theta)=\mu \Leftrightarrow (b')^{-1}(\mu)=\theta\Leftrightarrow g(\mu)=\theta$$

称$g(\mu)$为典型链接函数（$Canonical\ Link\ function$），名如其意，是将随机变量的期望和典型参数链接起来的函数。此外，称$b''(\theta)$为 方差函数,记作$V(\theta)$，

1.3 常见指数族的标准型&典型链接函数

下面是常见的指数族分布随机变量的标准型。

一维正态分布$N(\mu,\sigma^2)$：

$$f(y;\mu,\sigma^2)=exp\{\frac{y\mu-\frac{\mu^2}{2}}{\sigma^2}-\frac{y^2}{2\sigma^2}-log(\sqrt{2\pi}\sigma)\}$$

二项分布$b(n,p)$：

$$P(Y=y)=exp\{ylog(\frac{p}{1-p})+mlog(1-p)+log\tbinom{m}{y}\}$$

泊松分布$P(\lambda)$：

$$P(Y=y)=exp\{ylog(\lambda)-\lambda-log(y!)\}$$

Gamma分布$Gamma(\alpha,\frac{\alpha}{\mu})$：

$$f(y;\alpha,\mu)=\frac{1}{\Gamma(\alpha)}(\frac{\mu y }{\alpha})^\alpha \frac{1}{y} exp\{-\frac{\mu y }{\alpha}\}$$

$$f(y;\alpha,\mu)=exp\{\alpha(-y/\mu)+(\alpha-1)log(y)+\alpha log(\alpha)-log(\Gamma(\alpha))\}$$

	$\theta$	$\phi$	$a(\phi)$	$b(\theta)$	$c(y,\theta)$	$Canonical\ Link$	$Range$
Normal	$\mu$	$\sigma^2$	$\sigma^2$	$\frac{\theta^2}{2}$	$-\frac{y^2}{2\sigma^2}-log(\sqrt{2\pi}\sigma)$	$g(\mu)=\mu$	$(-\infty,\infty)$
Binomial	$log(\frac{p}{1-p})$	1	1	$mlog(1-p)=-mlog(1+exp(\theta))$	$log\tbinom{m}{y}$	$g(\mu)=log(\mu)$	$(-\infty,\infty)$
Possion	$log(\lambda)$	1	1	$\lambda=exp(\theta)$	$log(y!)$	$g(\mu)=log(\frac{\mu}{1-\mu})$	$(-\infty,\infty)$
Gamma	$-\frac{1}{\mu}$	$\frac{1}{\alpha}$	$\frac{1}{\alpha}$	$log(u)=-log(-\theta)$	$(\alpha-1)log(y)+\alpha log(y)-log(\Gamma(\alpha))$	$g(\mu)=-\frac{1}{\mu}$	$(-\infty,0)$

1.4 常见指数族分布的关系

下面我们不假证明的说明指数族分布之间的一些联系。

• $Bernoulli\to Binomial$:概率为$p$的独立重复的伯努利试验

我们在独立重复实验的背景下做一个推广，我们假定一个独立同分布的实验场景：将单位时间或者空间n等分，每一份上随机事件X都独立等概率发生，且假定单位时间内随机事件X的平均发生次数为$\lambda$。后文我们称之为单位时间内推广的独立重复试验场景（很不严谨）。

• $Binomial\to Possion$:在单位时间内推广的独立重复试验场景下，单位时间内事件发生的次数服从参数为$\lambda$的$Possion$分布。时二项分布在n趋于无穷且$np=\lambda$的情形下的推广。
• $Possion \to Gamma(\alpha,\lambda=\alpha/\mu)$:依旧是在单位时间内推广的独立重复试验场景下，我们假定单位时间内事件X平均发生次数为$\lambda$，事件X第$\alpha$次发生的时间服从$Gamma(\alpha,\lambda=\alpha/\mu)$。当$\alpha=1$时，即事件第一次发生的时间服从指数分布$Exp(\lambda)$即$Gamma(1,\lambda)$。
• 二项分布、泊松分布的正态近似：极限性质下有中心极限性定理中最早的棣莫弗拉普拉斯定理，即二项分布的正态近似。一个更直观的例子是高尔顿板，即一维左右等概率的随机游走。
  
  下面是流传已久的机器学习常用分布的关系图：
  

1.5 指数族与最大熵模型

指数族分布被证实是满足最大熵的一族分布，在自然界中最常见也是我们研究中最偏好的一类概率模型。指数族分布广泛存在于自然界和统计应用当中是有原因的。由以上的关系可以直观体会到这一众指数族分布之间存在的关系，其中很核心的一点是“等概率”，从熵的视角去看，“等概率=最大熵”。

• 最大熵原理
  最大熵原理是概率学习模型的一个准则。最大熵原理认为，学习概率模型时，在模型假设空间里，熵最大的模型是最好的模型。
  熵是由概率定义的一个随机事件的混乱程度的度量，其在样本空间内的元素等可能发生时达到最大值。

$$H(X)=-\int_{X\in D}f_x(x)log(f_x(x))dx$$

模型若要熵最大，在已知条件之外应当是等可能地无知的。由此求最大熵模型的问题转化为了以模型熵值为目标函数，以已知条件为约束的优化问题。此处不加证明给出符合最大熵原理的概率模型经证明服从指数族分布。（详细证明见李航《统计学习方法》第六章）。下面是几个给定约束的最大熵分布的例子：

2、基于指数族的广义线性模型

首先回顾一下经典线性模型的假定：

经典线性回归模型的5个假定

$$Y=X^T\beta+\epsilon$$

• 假设1 线性假设$E(Y|X)=X\beta$
• 假设2 X严格外生，$E(\epsilon|X)=0$
• 假设3 $E(XX')$非奇异且最小特征值$\lambda_{min}\to +\infty(n\to +\infty)$
• 假设4 随机误差项条件同方差且条件不相关
• 假设5 随机误差项$\epsilon \sim N(0,\sigma^2I)$

当我们模型的响应变量$Y$是离散型的时候，贸然套用经典的线性回归模型是完全错误的。因此我们基于指数族（不限于离散型）设计了广义的线性模型。对于设计的广义线性模型应当满足以下基本要求：

• 经典的正态假定下线性模型应当是广义回归的特例
• 回归函数应当是关于$X\beta$的单调函数，以保证我们的系数有更好的解释性

2.1 广义线性模型构成

广义线性模型由三部分构成：响应变量分布、线性预测量和链接函数。

• 响应变量分布
  给定观测变量$X$的条件下，响应变量$Y$的概率分布属于指数族。 $$f(y|X)=exp([\theta(X)y-b(\theta(X))]/a(\phi)+c(y,\phi))$$ 这里我们令典型参数$\theta=\theta(X)$，表明典型参数$\theta$you观测变量决定
• 链接函数和线性预测量
  假定$\mu(X)=E(Y|X)$，定义链接函数$g(.)$ （区分这里的链接函数和前文的典型链接函数），$g(\mu(X))=X^T\beta,\mu(X)=g^{-1}(X^T\beta)$。给定链接函数$g(.)$,若链接函数满足$g(.)=(b')^{-1}$时，即链接函数是将条件期望和典型参数链接起来的典型链接函数，可以将典型参数$\theta(X)$写为: $$\theta (X) =(b')^{-1} (\mu (X))= (b')^{-1} (g^{-1} (X^T \beta))=(g\circ b')^{-1}( X^T \beta)=h( X^T\beta)$$ 此处$h( X^T\beta)$在$g(.)$是典型链接函数时，$h( X^T\beta)=X^T\beta=\theta(X)$.

2.2 广义线性模型的极大似然估计

以往我们处理经典线性模型的参数估计时，更多看到的是最小二乘法的视角下的参数估计，但是极大似然估计法是更加广泛的参数估计方法。在经典线性模型中，极大似然估计和最小二乘法是等价的。

经典线性模型MLE求解参数估计(考虑异方差存在的情形，$\epsilon\sim N(0,W\sigma^2)$)

$$l(\beta)=-\frac{n}{2}log(2\pi)-\frac{n}{2}log(\det{W})+\sum_{i=1}^n \frac{(Y_i-X_i'\beta)^2}{W_{i,i}\ \sigma^2}$$

$$\propto(Y-X\beta)'W^{-1}(Y-X\beta) \propto -2\beta'X'W^{-1}Y+\beta'X'W^{-1}X\beta$$

$$l'(\beta)\propto -2W^{-1}X'Y+2X'W^{-1}X\beta=0$$

解得：

$$\hat{\beta}=(X'W^{-1}X)^{-1}X'W^{-1}Y$$

此处略过对二阶条件得验证，等价于最小二乘法或者最大投影法的角度下的参数估计值

2.2.1 广义线性模型似然函数

• 似然函数
  单个样本的似然函数$f(Y_i|X,\theta,\phi)=exp\{(\theta_i y_i-b(\theta_i))/a_i(\phi)+c(y_i,\phi)\}$(注意广义线性模型的典型参数$\theta$ 由 $X_i$ 决定，故有下标)，我们将$\theta (X) =(g\circ b')^{-1}( X' \beta)=h( X'\beta)$带入，写出对数似然函数如下：

$$l_n(\beta,\phi)=\sum_{i=1}^n[Y_i h(X_i'\beta)-b(h(X'\beta))]/a(\phi)+\sum_{i=1}{n}c(Y_i,\phi)$$

最大似然估计即求解:

$$(\hat\beta_{mle} ,\hat\phi_{mle})=\mathop{argmax}\limits_{\beta,\phi} l_n(\beta,\phi)$$

• 关于离散参数$\phi$
  由前文中关于对数似然导数的期望： $$E( \frac{\partial}{\partial\theta}l(\theta,\phi)) =E[\frac{(y-b'(\theta))}{a(\phi)}]$$ - $$E( \frac{\partial^2}{\partial\theta\ \partial \phi}l(\theta,\phi)) =a'(\phi)E[-\frac{(y-b'(\theta))}{a(\phi)^2}]$$
  由于$E(y)=b'(\theta)$，故此处的混合偏导为0，由样本数据中，倘若对$\theta$有渐进无偏的估计，通过混合求导为0，可以说明由样本计算得到的$\hat\beta_{mle} ,\hat\phi_{mle}$是渐进不相关的。对典型参数$\beta$和离散参数$\phi$采取相对独立的估计方法是可行的。实际应用中参数$\phi$通常采用矩估计而非最大似然估计。
  此外为了设计可能存在的异方差问题，我们假定： $$a_i(\phi)=\frac{\phi}{\omega_i}$$ 故有关于$\phi$的矩估计推导如下：

$$Var(Y)=a(\phi)b''(\theta)=\frac{\phi b''(\theta)}{\omega_i}; S^2_Y=\frac{1}{n-d}\sum_{i=1}^n (Y_i-\mu_i)^2$$

$$\hat\phi=\frac{1}{n-d}\sum_{i=1}^n\frac{\omega_i (Y_i-\mu_i)^2}{b''(\theta_i)}$$

其中$d$为参数$\beta$的维度。此外这里暗含了经典线性模型中的随机扰动项条件不相关假设，以保证样本的方差矩阵是对角阵，才可以类比我们在经典线性回归中的异方差情形。
由此我们得到用于求解最大似然估计的对数似然函数的形式：

$$l_n(\beta,\phi)=\sum_{i=1}^n \frac{\omega_i [Y_i h(X_i'\beta)-b(h(X'\beta))]}{\phi}+\sum_{i=1}{n}\ c(Y_i,\phi)$$

2.2.2 迭代加权最小二乘法求解MLE

求解关于参数$\beta$的最大似然估计，等价于求解以下目标函数的最大值(即此时离散参数$\phi=1$的情形)：

$$l_n(\beta,\phi)=\sum_{i=1}^n \omega_i [Y_i h(X_i'\beta)-b(h(X'\beta))]$$

首先了解一下目标函数的性质：

$$l''(\beta)=-\sum_{i=1}^n \omega_i b''(\theta)X_iX_i'$$

目标函数关于参数$\beta$的二阶导为半负定的矩阵，即Hessian矩阵为半负定的。而此处参数$\beta$取值范围为$(-\infty,+\infty)$，由凸优化理论知，对于定义在凸集（$R为凸集$）上的Hessian矩阵为半负定的情形下，目标函数的极大值唯一，且为全局的最大值。因此此处求解最大似然估计是一个凸优化问题。

此处求解采用牛顿法迭代求解。求解思路如下:

$$\beta^{(t+1)}=\beta^{(t)}+[l''(\beta^{(t)})]^{-1}l'(\beta^{(t)})$$

为了将问题转化为矩阵运算的形式，我们定义矩阵$W_{n\times n}=diag(\omega_i[b''(\theta_i)g'(\mu_i)]^{-1})$,$\ G_{n\times n}=diag(g'(\mu))$。用矩阵形式表示$l'(\theta)，l''(\theta)$的流程如下：

为了方便后面直接使用我们给出$\theta_i(X_i)$的导数形式

$$\frac{\partial}{\partial \beta}\theta_i(X_i)=h(X_i'\beta)'=(g\ \circ\ b')'(X_i'\beta)=\frac{X_i}{b''(\theta_i(X))g'(\mu(X_i))}$$

$$l'(\beta)=\sum_{i=1}^n \frac{\omega_i[Y_i-b'(h(X_i'\beta))]X_i}{b''(\theta_i(X))g'(\mu(X_i))}$$

由前文$E(Y)=b'(\theta)=\mu$，此处$\mu(X_i)=b'(h(X_i'\beta))$，令$Y=[Y_1,Y_2...Y_n]',\mu(X)=[\mu_1(X_1),\mu_2(X_2)...\mu_n(X_n)]'$,$X_{n\times d}=[X_1,X_2...X_n]'$,$d$为随机变$X_i$的维数，也是$\beta$维数。

$$l'(\beta)=\sum_{i=1}^n \frac{\omega_i[Y_i-b'(h(X_i'\beta))]X_i}{b''(\theta_i(X))g'(\mu(X_i))}=X'WG(Y-\mu(X))$$

由指数函数的性质知，$E(l'(\theta)^2)=-E(l''(\theta))$，矩阵形式下有 $$E(l'(\theta)l'(\theta)^T)=-E(l''(\theta))$$ 故由链式法则:

$$E([l'(\theta(\beta))][l'(\theta(\beta)]')=E(\theta'(\beta)\ l'(\theta)l'(\theta)^T\ \theta'(\beta)))=E(\theta'(\beta)\ l''(\theta)\ \theta'(\beta)))=-E(l''(\beta))$$

故:

$$-E(l''(\beta))=E(l'(\beta)l'(\beta)^T)$$

$$=E(X'WG(Y-\mu(X))[X'WG(Y-\mu(X))]')$$

$$=X'WG\ E[(Y-\mu(X))^2]\ GWX$$

$$=X'WG\ diag(\frac{\phi b''(\theta_i)}{\omega_i}) GWX$$

由最大似然目标函数的等价形式下，离散参数$\phi=1$，且参照矩阵$W,G$的定义，可得：

$$-E(l''(\beta))=X'WG\ diag(\frac{b''(\theta_i)}{\omega_i}) GWX=X'WW^{-1}WX=X'WX$$

综上我们得到了$l'(\beta)$的形式和$E(l''(\beta))$的形式，值得注意的是我们这里没有采用$l''(\beta)$的具体形式，而是选取了$l''(\beta)$的期望这一渐进的估计形式，在样本量足够大的时候，有大数律保证这种替代是有效的。因此牛顿法迭代式可以写为：

$$\beta^{(t+1)}=\beta^{(t)}+[X'WX]^{-1}X'WG(Y-\mu)$$

$$=[X'WX]^{-1}X'W\{G(Y-\mu)+X\beta^{(t)}\}$$

2.2.3 迭代加权最小二乘法

• Step1 初始化 $\beta^{(0)},t=0$
• Step2 对于给定的$t$，用$\beta^{(t)}$替代$\beta$来计算此时的$G,\mu$，计算 $$z=G(Y-\mu)+X\beta^{(t)}$$
• Step3 对于给定的$t$，用$\beta^{(t)}$替代$\beta$来计算此时的$W$,给更新参数

$$\beta^{(t+1)}=[X'WX]^{-1}X'Wz$$

• Step4 重复Step2，Step3直至\beta^{(t)}收敛,记此时计算得到的$W$为$\hat W$,参数收敛后有

$$\hat \beta=(X'\hat WX)^{-1}X'\hat Wz$$

注意到 $$cov(z|X)=cov(G(Y-\mu)|X)=GCov(Y)G'$$
前文已证$Cov(Y)=diag(\frac{\phi b''(\theta_i)}{\omega_i})$，故有

$$cov(z|X)=\phi W^{-1}$$

将关于离散参数$\phi$的估计值$\hat \phi$一并带入，由$\hat \beta$和$z$的关系，此时的参数的分布我们也可以得到：

$$\hat cov(\hat \beta|X)=\hat \phi (X'\hat W X)^{-1}X'\ \hat W \hat W^{-1} \hat W \ X(X'\hat W X)^{-1}$$

$$=\hat \phi (X'\hat W X)^{-1}X'\hat WX(X'\hat W X)^{-1}$$

在实际操作中，对于每一时刻的$E(l''(\beta))$求逆会带来$O(d^2)$的运算复杂度，因此实际计算中更常用的是梯度上升算法来计算。

2.3 广义线性模型的偏差和残差

2.3.1 广义线性模型的偏差以及偏差分析

以上我们已经给出了如何在无参数约束的条件，即$\beta\in R^d$下来估计参数$\beta$，实际应用中通常这些参数取值存在着一些限制，导致参数$\beta$的估计值无法收敛到无约束时的全局最优，出现欠拟合问题。为了诊断由参数取值的约束带来的偏差，我们在广义线性模型中引入了偏差的概念。
此处为了方便表示，我们记无约束典型参数估计值为$\tilde\theta=\theta(\tilde \beta)$，有约束下的典型参数估计记作$\hat \theta=\theta(\hat \beta)$。下面定义偏差（Deviance） 为两类参数下对数似然之差。

$$D(Y;\hat \mu)=2\sum_{i=1}^n \omega_i \{Y_i(\tilde \theta-\hat \theta)-b(\tilde \theta)+b(\hat \theta)\}$$

这在经典的正态线性模型下，若$\theta(\mu)=\mu$，$b(\theta)=\frac{\theta^2}{2}$，有：

$$D(Y;\hat\mu)=\sum_{i=1}^n (Y_i-\hat Y_i)^2$$

即加权残差平方和。
换个角度看我们这里定义的偏差的概念，其实是在不同的参数空间下的似然比检验的对数结果。

$$LRT=\frac{\mathop{sup}\limits_{\theta\in \Theta_1 }L(Y;\theta)}{\mathop{sup}\limits_{\theta\in \Theta }L(Y;\theta)}$$

$$log(LRT)=2\phi(l(Y;\tilde \theta)-l(Y;\hat \theta))$$

$$=2\sum_{i=1}^n \omega_i \{Y_i(\tilde \theta-\hat \theta)-b(\tilde \theta)+b(\hat \theta)\}=D(Y;\hat \mu)$$

当比较两个不同大小（参数个数不同）的模型时，即$\hat \theta_1 \in \Theta_1$，$\hat \theta_2 \in \Theta_q$，两个参数空间的维度不同，在似然比检验的角度下：

$$D(Y;\hat \mu,\theta_1)-D(Y;\hat \mu,\theta_2)=2\phi\{\mathop{sup}\limits_{\theta\in \Theta_1 } l(\theta)-\mathop{sup}\limits_{\theta\in \Theta_2 } l(\theta)\}=2\phi log(LRT)$$

由似然比检验的渐进性质可知，$\phi D(Y;\hat \mu)\to \chi^2_{df}$，即服从参数维度为$df=dim(\Theta_1)-dim(\Theta_0)$的$\chi^2$分布.由此我们可以检验两个不同大小的模型的拟合能力。

2.3.2 广义线性模型的残差定义

2.4 广义线性模型的例子

2.4.1 logistic回归

2.4.2 Possion回归

3、广义线性模型的稀疏解

3.1 高置信区间的稀疏解

3.2 采用惩罚似然的特征选择