广义线性模型（GLMs）及算法介绍

  一般我们了解的线性模型是针对连续性变量，并且服从正态分布的，但是在实际应用上显得非常的局限。因为我们我看到的数据很多都是离散的，而且不是服从正态分布的。针对这种情况，对传统线性模型进行推广，行成了现在的广义线性模型。广义线性模型使得变量从正态分布拓展到指数分布族，从连续型变量拓展到离散型变量，这就使得在现实中有着很好的运用，下面开始介绍广义线性模型。

广义线性模型（GLM）定义

由以下三部分组成：
1 随机部分
  随机样本$Y_1,Y_2,...,Y_n$服从的分布来自指数分布族，即$Y_i$的分布形式为$$f\left(y_i;{\theta }_i,\varphi \right)=exp\left\{\frac{y_i{\theta }_i-b\left({\theta }_i\right)}{a\left(\varphi \right)}+c\left(y_i,\varphi \right)\right\}$$
其中参数${\theta }_i$叫做正则参数，并且随着指数$i（i=1,2,...,n）$而变化，但是扰乱因子$\varphi$ 是个常数。
2 系统部分
  对于第$i$个观测$Y_i$，我们有一个称为系统部分的线性预测值，即所研究变量的线性组合，即$${\eta }_i=x_i^T\beta =\sum _{j=1}^p{x}_{ij}{\beta }_j,i=1,2,...n$$
3 连接函数
  有一个单调可微函数$g\left(\right)$称为连接函数，它将随机部分的期望和系统部分连接起来，通过下面的等式$$g\left({\mu }_i\right)={\eta }_i=x_i^T\beta ,i=1,2,...n,$$其中${\mu }_i=E\left(Y_i\right)$是$Y_i$的期望。
  矩阵表示：
$$\eta ={\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\\ ⋮\\ {\eta }_n\end{array}\right]}_{n\times 1},\mu ={\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ {\mu }_n\end{array}\right]}_{n\times 1},X={\left[\begin{array}{c}{x}_1^{{}^{\prime }}\\ x_2^{{}^{\prime }}\\ ⋮\\ x_n^{{}^{\prime }}\end{array}\right]}_{n\times p}$$
那么连接函数可以用矩阵形式表示$$g\left(\mu \right)=\eta =X\beta$$

连接函数介绍

1、正如$GLMs$的定义所指出的那样，连接函数是单调可微的，用于连接随机部分的期望和系统部分的线性预测值
2、选择与分布相关的自然参数作为连接函数，在这种情况下，它被称为点则连接。具体而言，如果连接函数$g()$采用与自然参数相同的函数形式，那么它被称为点则连接函数。
3、点则连接的优点是它可以带来非常好的统计特性，并且使用起来很方便。例如，对于最常用的分布，我们有以下点则连接函数。

Normal	$\mu =\eta$	(identity-link)
Poisson	$log\mu =\eta$	(log-link)
Bernoulli	$log\frac{\pi }{1-\pi }=\eta$	(logistic-link)
Binomial	$log\frac{\pi }{1-\pi }=\eta$	(logistic-link)

4、然而，点则连接并不是连接函数的唯一选择。其他可能的$GLMs$连接函数包括
（a）二项分布的Probit连接：$\eta ={\Phi }^{-1}\left(\pi \right)$;$0<\pi <1$,其中$\Phi ()$叫做累计分布函数（不是点则连接呦）
（b）补充的二项分布的log-log连接 η = l o g { − l o g ( 1 − π ) } , 0 < π < 1 \eta =log\left \{ -log\left ( 1-\pi \right ) \right \},0<\pi<1 η=log{−log(1−π)},0<π<1
（c）属于任何幂族分布的连接$$\eta =\{\begin{array}{cc}{\mu }^{\lambda },if\lambda \ne 0& \\ log\mu ,if\lambda =0& \end{array}$$

最大似然估计（MLE）的一般原则

  假设我们对未知参数$\theta$的对数似然函数，比如说$l\left(\theta \right)$我们想找出$\theta$的最大似然估计（MLE） $\hat{\theta }$，即 θ ^ ≡ a r g m a x θ ⊂ Ω { l ( θ ) } \hat{\theta }\equiv arg \underset{\theta \subset \Omega }{max}\left \{ l\left ( \theta \right ) \right \} θ^≡argθ⊂Ωmax​{l(θ)}
这是估计方程的解$$\frac{\partial l\left(\theta \right)}{\partial \theta }=0$$
1、在这种情况下，例如，对于正态分布参数$\theta$的最大似然估计 $\hat{\theta }$可以有一个明确的数学表达式$（\mu =\frac{1}{n}{\sum }_{i=1}^{n}ln{x}_i）$
2、一般来说，最大似然估计$\theta$没有没有明确的数学解。相反，需要某些数值优化方法。
3、统计学中最常用的两种优化方法是Newton-Raphson算法和Fisher得分算法，他们都涉及计算$l\left(\theta \right)$对$\theta$的2次偏导数。

Newton-Raphson算法

  该算法计算最大似然估计 $\hat{\theta }$，通过下面的迭代：
$${\theta }^m={\theta }^{m-1}+{\left[-l^{{}^{\prime \prime }}\left({\theta }^{\left(m-1\right)}\right)\right]}^{-1}\left[l^{{}^{\prime }}\left({\theta }^{\left(m-1\right)}\right)\right]（1）$$
其中$m=1,2,...$这里，$$l^{{}^{\prime }}\left({\theta }^{\left(m-1\right)}\right)=\frac{\partial l\left(\theta \right)}{\partial \theta }{|}_{\theta ={\theta }^{\left(m-1\right)}}$$$$l^{{}^{\prime \prime }}\left({\theta }^{\left(m-1\right)}\right)=\frac{{\partial }^{2}l\left(\theta \right)}{\partial \theta \partial {\theta }^T}{|}_{\theta ={\theta }^{\left(m-1\right)}}$$是$p\times 1$和$p\times p$的向量和矩阵$（p是\theta 的维数）$
注1 $l^{{}^{\prime }}\left(\theta \right)被称为\theta 的得分函数。-l^{{}^{\prime \prime }}\left(\theta \right)被称为\theta 的观测信息矩阵$
注2 算法（1）需要初始值，例如${\theta }^0$，以开始迭代过程。初始值的选择通常需要经验。
注3 算法（1）迭代直到收敛。例如，当迭代结果满足$$\frac{\Vert {\theta }^{\left(m\right)}-{\theta }^{\left(m-1\right)}\Vert }{\Vert {\theta }^{\left(m-1\right)}\Vert }\le 10^{-5}$$则迭代停止。${\theta }^{\left(m\right)}$可以认为是最大似然估计 $\hat{\theta }$。

Fisher得分算法

  Fisher得分算法与Newton-Raphson算法相同，只是（1）式中的观测矩阵$-l^{{}^{\prime \prime }}\left(\theta \right)$被Fisher信息矩阵所代替$$I\left(\theta \right)=E\left[-l^{{}^{\prime \prime }}\left(\theta \right)\right]=-\int l^{{}^{\prime \prime }}\left(\theta |Y\right)f_Y\left(Y|\theta \right)dY$$
注释 Fisher得分算法和Newton-Raphson算法一般收敛于同一解。对于前者，在某些情况下，在信息矩阵的解析式上可能比后者更简单。例如Fisher信息矩阵可能是对角阵或者块对角阵，二观测信息矩阵可能不是。其次作为副产物，这两种算法都提供了极大似然估计的协方差矩阵。

广义线性模型（GLMs）中的最大似然估计（MLE）

  首先，GLMs中的对数似然函数具有这样的形式$$l=\sum _{i=1}^{n}logf\left(y_i;{\theta }_i,\varphi \right)=\sum _{i=1}^n\frac{\left(y_i{\theta }_i-b\left({\theta }_i\right)\right)}{a_i\left(\varphi \right)}+\sum _{i=1}^{n}c\left(y_i,\varphi \right)$$其中，$\beta ={\left({\beta }_1,{\beta }_2,...,{\beta }_p\right)}^T$，${\beta }_j$的得分函数为$$U_j=\frac{\partial l}{\partial {\beta }_j}=\sum _{i=1}^n\frac{\left(y_i-b^{{}^{\prime }}\left({\theta }_i\right)\right)}{a_i\left(\varphi \right)}\frac{\partial {\theta }_i}{\partial {\beta }_j}=\sum _{i=1}^n\frac{\left(y_i-\mu \right)}{a_i\left(\varphi \right)}\frac{\partial {\theta }_i}{\partial {\beta }_j}（2）$$
其中，${\mu }_i=E\left(Y_i\right)=b^{{}^{\prime }}\left({\theta }_i\right)，Var\left(Y_i\right)=b^{{}^{\prime \prime }}\left({\theta }_i\right)a\left(\varphi \right)$，我们使用链式法则进行差异化$$\frac{\partial {\theta }_i}{\partial {\beta }_j}=\frac{\partial {\theta }_i}{\partial {\mu }_i}\frac{\partial {\mu }_i}{\partial {\beta }_j}$$因为$$\frac{\partial {\theta }_i}{\partial {\mu }_i}=\frac{1}{\frac{\partial {\mu }_i}{\partial {\theta }_i}}=\frac{1}{b^{{}^{\prime \prime }}\left({\theta }_i\right)}=\frac{a_i\left(\varphi \right)}{b^{{}^{\prime \prime }}\left({\theta }_i\right)a_i\left(\varphi \right)}=\frac{a_i\left(\varphi \right)}{Var\left(Y_i\right)}$$
并且$$\frac{\partial {\mu }_i}{\partial {\beta }_j}=\frac{\partial {\mu }_i}{\partial {\eta }_i}\frac{\partial {\eta }_i}{\partial {\beta }_j}=\frac{\partial {\mu }_i}{\partial {\eta }_i}{x}_{ij}$$其中$x_{ij}$是$x_i的第j个分量$，我们知道$$\frac{\partial {\theta }_i}{\partial {\beta }_j}=\frac{a_i\left(\varphi \right)}{Var\left(Y_i\right)}\frac{\partial {\mu }_i}{\partial {\eta }_j}{x}_{ij}$$因此（2）式就化为了$$U_j=\sum _{i=1}^n\left[\frac{\left(y_i-{\mu }_i\right)}{Var\left(Y_i\right)}{x}_{ij}\left(\frac{\partial {\mu }_i}{\partial {\eta }_i}\right)\right]=\sum _{i=1}^n\frac{\left(y_i-{\mu }_i\right)}{g^{{}^{\prime }}\left({\mu }_i\right)V_i}{x}_{ij}（3）$$其中$V_i=Var\left(Y_i\right)$，并且$$\frac{\partial {\mu }_i}{\partial {\eta }_i}=\frac{1}{\frac{\partial {\eta }_i}{\partial {\mu }_i}}=\frac{1}{g^{{}^{\prime }}\left({\mu }_i\right)}$$由于${\eta }_i=g\left({\mu }_i\right)$，因此$\beta$的得分向量是$$U\equiv U\left(\beta \right)=\sum _{i=1}^n\frac{\left(y_i-{\mu }_i\right)}{g^{{}^{\prime }}\left({\mu }_i\right)V_i}{x}_i（4）$$另一方面，（3）式对${\beta }_j$求偏导得$$\frac{{\partial }^{2}l}{\partial {\beta }_j\partial {\beta }_k}=\frac{\partial U_j}{\partial {\beta }_k}=\sum _{i=1}^n\left(-\frac{\partial {\mu }_i}{\partial {\beta }_k}\right)\frac{1}{g^{{}^{\prime }}\left({\mu }_i\right)V_i}{x}_{ij}+\sum _{i=1}^n\left(y_i-{\mu }_i\right)\frac{\partial \left[\frac{1}{g^{{}^{\prime }}\left({\mu }_i\right)V_i}\right]}{\partial {\beta }_k}{x}_{ij}（5）$$由于$E\left(Y_i-{\mu }_i\right)=0$，所以（5）式的第二项在进行期望时就消失了。即Fisher信息阵的矩阵形式就变成了$$I\left(\beta \right)=E\left(\frac{{\partial }^{2}l}{\partial \beta \partial {\beta }^T}\right)=\sum _{i=1}^n\frac{1}{g^{{}^{\prime }}{\left({\mu }_i\right)}^2{V}_i}{x}_{ij}{x}_{ik}$$因此当我们表示$W_i=\frac{1}{g^{{}^{\prime }}{\left({\mu }_i\right)}^2{V}_i}$，并且$$W=diag\left(W_1,W_2,...,W_n\right)=\left(\begin{array}{cccc}{W}_1& 0& \cdots & 0\\ 0& W_2& \cdots & 0\\ ⋮& 0& \ddots & 0\\ 0& \cdots & 0& W_n\end{array}\right)$$则Fisher信息阵就可以表示为$$I\left(\beta \right)=X^{T}WX$$令$D=diag\left(g^{{}^{\prime }}\left({\mu }_1\right),g^{{}^{\prime }}\left({\mu }_2\right),...,g^{{}^{\prime }}\left({\mu }_n\right)\right)$，这样（4）式就可以写成$$U=U\left(\beta \right)=X^{T}WD\left(y-\mu \right)$$

计算最大似然估计（MLE）参数$\beta$的算法

  假设我们有一个估计${\beta }^{\left(m-1\right)}$，基于这个估计我们计算$${\mu }^{\left(m-1\right)}=\mu \left({\beta }^{\left(m-1\right)}\right)，W^{\left(m-1\right)}=W\left({\beta }^{\left(m-1\right)}\right)$$
并且有$$D^{\left(m-1\right)}=D\left({\beta }^{\left(m-1\right)}\right)$$  那么Fisher得分算法就会显示$\beta$的下一次迭代$${\beta }^{\left(m\right)}={\beta }^{\left(m-1\right)}+{\left[I\left({\beta }^{\left(m-1\right)}\right)\right]}^{-1}\left[U\left({\beta }^{\left(m-1\right)}\right)\right]={\beta }^{\left(m-1\right)}+{\left[X^T{W}^{\left(M-1\right)}X\right]}^{-1}\left[X^T{W}^{\left(M-1\right)}{D}^{\left(M-1\right)}\left(y-{\mu }^{\left(m-1\right)}\right)\right]$$可以写成$${\beta }^{\left(m\right)}={\left[X^{T}WX\right]}^{-1}{X}^T{W}^{\left(m-1\right)}\left[X{\beta }^{\left(m-1\right)}+D^{\left(m-1\right)}\left(y-{\mu }^{\left(m-1\right)}\right)\right]$$令$$Z^{\left(m-1\right)}=X{\beta }^{\left(m-1\right)}+D^{\left(m-1\right)}\left(y-{\mu }^{\left(m-1\right)}\right)$$然后它又可以被写成$${\beta }^{\left(m\right)}={\left(X^{\left(T\right)}{W}^{\left(m-1\right)}X\right)}^{-1}{X}^T{W}^{m-1}{Z}^{\left(m-1\right)}（6）$$
注释 （6）式意味着，给定参数$\beta$的解，我们需要计算“工作权重矩阵”$W$和“工作响应向量”$Z$，然后利用广义加权最小二乘法得到$\beta$的更新解。

广义线性模型实例解析

  下表中的ARPI事物数据在协变量X的不同处观察到Y,并且数据是服从Poisson分布的。我们利用GLM来解决这个问题。

$Y_i$	2	3	6	7	8	9	10	12	15
$x_i$	-1	-1	0	0	0	0	1	1	1

数据即探索Y和X之间的关系。设$Y_i$为变量$y$的第$i$个数，表示$E\left(Y_i\right)={\mu }_i$。我们通过建立关系$$g\left({\mu }_i\right)=x_i^{{}^{\prime }}\beta$$  对于这个Poisson数据集，点则连接是对数连接函数。
$$log{\mu }_i={\beta }_0+{\beta }_1{x}_i=\left(1,x_i\right)\left(\begin{array}{c}{\beta }_0\\ {\beta }_1\end{array}\right)=x_i^T\beta$$接下来我们要来求$W和Z$的表达式。
  我们已知的条件有$g^{{}^{\prime }}\left({\mu }_i\right)=\frac{1}{{\mu }_i}$，对于Poisson分布显然有$V_i=E\left(Y_i\right)={\mu }_i$，所以$$W_i={\left[\left(g^{{}^{\prime }}{\left({\mu }_i\right)}^2\right)V_i\right]}^{-1}=exp\left\{x_i^T\beta \right\}$$并且$$Z_i=x_i^T\beta +g^{{}^{\prime }}\left({\mu }_i\right)\left(y_i-{\mu }_i\right)=x_i^T\beta +\frac{\left(y_i-{\mu }_i\right)}{{\mu }_i}$$
  我们选择$\beta 的初始值{\beta }_0=2，{\beta }_1=1$。结合Fisher迭代算法，代入数据。这个过程一直持续到收敛。结果如下表

m	0	1	2	3	4
${\beta }_0^m$	2	1.9150	1.8902	1.8892	1.8892
${\beta }_1^m$	1	0.7235	1.8902	1.8892	1.8892

因此$\beta 的MLE是{\beta }_0=1.8892，{\beta }_1=0.6697$

R语言代码

y <- c(2,3,6,7,8,9,10,12,15); 
x <- c(-1,-1,0,0,0,0,1,1,1)
X <- cbind(rep(1,9),x); beta_0 <- c(2,1)
for (i in 1:100){
beta <- beta_0
eta <- X %*% beta
mu <- exp(eta)
W <- diag(as.vector(mu))
Z <- X %*% beta + ((y-mu)*mu^(-1))
XWX <- t(X) %*% W %*% X
XWZ <- t(X) %*% W %*% Z
Cov <- solve(XWX)
beta_0 <- Cov %*% XWZ}
testdata<-data.frame(y,x)
summary(glm(y~x,family=poisson,data=testdata))