一、广义线性回归的概念

通过前面9讲的学习，我们了解了线性回归模型的结构，以及如何估计参数和假设检验。总结来说，线性回归结构的基本假定是，响应变量的期望值可以由自变量的线性预测子得到，即满足：

$EY=X\beta$

在现实应用当中，要满足这个假定其实很苛刻。一方面，某组自变量对应的响应变量取值的分布不一定是正态分布，可以是偏态的、拖尾的；另一方面，响应变量的取值可以是离散的。在这些情况下，线性回归模型就不再适用了。那应该如何处理呢？

既然是响应变量的问题，那么我们可以考虑，保留自变量的线性结构，但是和响应变量连接时套上一层函数进行转换，也就是满足：

$g(EY)=X\beta$

这里的 $g()$ 被称为联系函数(link function)，是已知且单调光滑的。联系函数，顾名思义，就是通过一个函数把响应变量和线性预测子联系起来，这样既处理了响应变量，也保留了线性成分。

显然，如果我们定义 $h=g^{-1}(t)$ ，那么就有 $EY=h(X\beta)$ 。不难猜测，如果 $g(t)=t$ ，那么就变成了一般的线性回归模型，通过这种方式我们就可以把广义线性回归和一般线性回归统一起来。

虽然这类模型的名称是广义线性回归，但是参数不是线性的，不过我们仍然称为广义线性回归模型，因为它可以和一般的线性回归模型联系起来。

解决了模型的结构，我们还需要确定响应变量的分布。通常我们假定响应变量服从指数分布族。如果你忘了第1讲学过的指数分布族，那么这里我再展示一下：

$f(y|\theta;\phi)=exp \left [ \frac{\theta y-b(\theta)}{\phi}+c(y;\phi) \right ]$

$\theta$ 称为自然参数(Natural Parameter)，是我们感兴趣的参数； $\phi$ 称为散布参数/尺度参数/讨厌参数(Dispersion Parameter)，在特定的分布下应当已知； $b(\theta)$ 和 $c(y;\phi)$ 在具体分布确定后已知表达式。

指数分布族有两条重要的性质：

$\begin{aligned} EY&=\dot{b(\theta)}\\ Var\ Y&=\phi \ddot{b(\theta)} \end{aligned}$

最后一个问题是，为什么要假定响应变量服从指数分布族？

这个问题其实无法回答得很完美，但是可以提供几个使用指数分布族的好处：

期望和方差的结构较为统一
估计方程具有封闭形式，求解高效
可选的连接函数非常自然且符合实际

二、广义线性回归的参数估计

假定模型之后，接下来的工作就是估计参数。由于响应变量服从指数分布族，所以不能直接使用最小二乘法。因此，我们考虑使用极大似然估计来完成。

为什么广义线性回归的参数估计不能使用最小二乘法？

线性回归中，最小二乘法的原理是最小化平方损失，也就是最小化 $||Y-X^T\beta||^2$ 。如果你使用极大似然估计，你会发现线性回归的对数似然函数是平方损失的一个线性函数，所以在求解优化时是等价的。而广义线性回归不满足这样的性质，且响应变量的方差要随均值变化而变化，所以只能使用极大似然估计。

2.1 极大似然估计与Fisher信息阵

根据极大似然估计的方法，首先我们需要构造对数似然函数：

$\begin{aligned} &l_i=\frac{\theta_iy_i-b(\theta_i)}{\phi_i}+c(y_i;\phi_i),i=1,\cdots,n\\ &l=\sum_{i=1}^nl_i=\sum_{i=1}^n \left[ \frac{\theta_iy_i-b(\theta_i)}{\phi_i}+c(y_i;\phi_i) \right] \end{aligned}$

为了求导后表示方便，我们定义得分统计量(score statistic)为：

$U_i=\frac{\partial l_i}{\partial\theta_i}=\frac{y_i-\dot{b(\theta_i)}}{\phi_i}=\frac{y_i-Ey_i}{\phi_i}$

不难看出，得分统计量的期望和方差满足：

$\begin{aligned} EU_i&=\frac{1}{\phi_i}E(y_i-Ey_i)=0\\ Var\ U_i&=EU_i^2=\frac{1}{\phi^2_i}E(y_i-Ey_i)^2=\frac{\ddot{b(\theta)}}{\phi_i} \end{aligned}$

我们称 $U_i$ 的方差为信息(information)，记作 $\mathfrak I_i=Var\ U_i$ 。注意到一个重要性质，即 $\mathfrak I=-E(U^{'})$ 。

既然我们的对数似然函数是要对 $\beta$ 求导，那么一定需要进行转换。由于对数似然函数无法直接对 $\beta$ 求导，那么我们可以根据链式法则展开：

$\frac{\partial l}{\partial\beta_j}=U_j=\sum_{i=1}^n\left[ \frac{\partial l_i}{\partial\beta_j} \right]=\sum_{i=1}^n\left[ \frac{\partial l_i}{\partial\theta_i}\cdot\frac{\partial \theta_i}{\partial \mu_i}\cdot\frac{\partial \mu_i}{\partial\beta_j} \right]$

其中 $\mu_i=Ey_i$ ， $g(\mu_i)=g(Ey_i)=x_i^T\beta$ ， $\mu_i=h(x^T_i\beta)$ 。

接下来我们依次计算每一项。首先，

$\frac{\partial l_i}{\partial\theta_i}=\frac{y_i-\dot{b(\theta_i)}}{\phi_i}=\frac{y_i-\mu_i}{\phi_i}$

其次：

$\frac{\partial \theta_i}{\partial\mu_i}=1/\left( \frac{\partial \mu_i}{\partial\theta_i} \right)=1/\left( \frac{\partial\dot{b(\theta_i)}}{\partial\theta_i} \right)=\frac{1}{\ddot{b(\theta_i)}}$

最后：

$\frac{\partial \mu_i}{\partial\beta_j}=\frac{\partial \mu_i}{\partial g(\mu_i)}\cdot\frac{\partial g(\mu_i)}{\partial\beta_j}=\frac{\partial \mu_i}{\partial g(\mu_i)}x_{ij}$

因此整合起来：

$U_j=\sum_{i=1}^n\left[ \frac{y_i-\mu_i}{Var\ y_i}x_{ij}\frac{\partial \mu_i}{\partial g(\mu_i)} \right]$

那么只需解方程 $U_j=0$ 即可求得回归系数的估计。

你一定注意到这里保留了 $\frac{\partial \mu_i}{\partial g(\mu_i)}$ 项，这是因为我们还不确定联系函数的表达式，一旦确定了联系函数就能求解了。

当然，这个似然方程不一定有解，或者解不唯一。理论上来说：

$\lim_{n\rightarrow\infty}P(有解)=1$
存在解 $\hat\beta_n$ 使得 $P(\lim_{n\rightarrow\infty}\hat\beta_n=\beta_0)=1$ ，其中 $\beta_0$ 是 $\beta$ 的真值

得到估计方程还不够。作为一种统计学的估计方法，我们还需考虑估计是否准确。对于广义线性回归来说，因为似然方程就是求解关于得分统计量的方程，因此我们可以使用得分统计量的方差来描述。这个方差又被称为信息，所有似然方程的联合信息就构成了Fisher信息矩阵，记作：

$\mathfrak I_{jk}=Cov(U_j,U_k)=E(U_jU_k)$

可以计算得到广义线性回归的Fisher信息阵为：

$\mathfrak I_{jk}=\sum_{i=1}^n\frac{x_{ij}x_{ik}}{Var(y_i)}\left( \frac{\partial \mu_i}{\partial g(\mu_i)} \right)^2\equiv F(\beta)$

由于在一定条件下有：

$F^{\frac{1}{2}}(\beta_0)(\hat\beta_n-\beta_0)\rightarrow N(0,I_p),p=\beta_0维数$

所以Fisher信息阵包含的方差越大，说明信息越多，那么参数估计也就越准确。

2.2 Fisher-Scoring算法

在2.1小节中，我们只搭建了估计方程，考察了估计准确性，但是还没有具体求解。事实上，对于这样的估计方程组来说，要求出解析解非常困难。不过，计算机的发展给我们带来了新的观点——迭代法。迭代思想一举成为几乎所有优化问题的解决办法。

如果你了解过优化问题的求解方法，你或许知道Newton-Raphson算法，这是非常经典的一种求解办法（这里不展开介绍）。这个算法更新很快，但是容易出现震荡或发散的问题。对于广义线性回归来说，统计学家开发了一种类似Newton-Raphson的算法，更新更平滑，收敛性更稳定——Fisher Scoring算法。

Fisher Scoring算法是一种基于得分函数和Fisher信息阵的方法，其基本思想是：

$b^{(m)}=b^{(m-1)}+\left[\mathfrak I^{(m-1)}\right]^{-1}U^{(m-1)}$

其中 $b^{(m)}$ 是指第 $m$ 次迭代时参数向量 $(\beta_1,\cdots,\beta_p)$ 的估计值。

为了求解方便，我们做一下简单变换：

$\mathfrak I^{(m-1)}b^{(m)}=\mathfrak I^{(m-1)}b^{(m-1)}+U^{(m-1)}$

信息阵可以改写为二次型的形式，即：

$\mathfrak I=X^TWX$

其中 $W$ 是一个 $n\times n$ 的对角阵，其对角元为：

$w_{ii}=\frac{1}{Var(y_i)}\left( \frac{\partial \mu_i}{\partial g(\mu_i)} \right)^2$

于是Fisher scoring算法的左侧可以改写为：

$X^TWXb^{(m)}$

同时，为了引入上面的 $W$ ，Fisher scoring算法的右侧可以改写为：

$\sum_{k=1}^p\sum_{i=1}^n\frac{x_{ij}x_{ik}}{Var(y_i)}\left( \frac{\partial \mu_i}{\partial g(\mu_i)} \right)^2b_{k}^{(m-1)}+ \sum_{i=1}^n\frac{(y_i-\mu_i)x_{ij}}{Var(y_i)}\left( \frac{\partial \mu_i}{\partial g(\mu_i)} \right)$

提取公因式 $x_{ij}w_{ii}$ ，我们再令：

$z_i=\sum_{k=1}^px_{ik}b_k^{(m-1)}+(y_i-\mu_i)\left( \frac{\partial g(\mu_i)}{\partial \mu_i} \right)$

于是右侧可以改写为：

$X^TWz$

最终，Fisher score算法可以改写为：

$X^TWXb^{(m)}=X^TWz$

你会惊奇地发现，从形式上看，这个方程类似于线性回归中的正规方程，仅仅只有两个区别：

带上了权重；
$z,W$ 都只与 $b$ 有关，可以参与迭代。

因此，上述这种改写后的方法又被称为迭代加权最小二乘法(Iterative weighted Least Square, IRWLS)，是Charnes等人在1976年提出来的方法，如今依然是广义线性模型参数估计当中最主流的方法。

如果进一步观察 $z_i$ ，我们还可以发现：

$z_i=g(\mu_i)+\dot{g(\mu_i)}(y_i-\mu_i)$

因此， $z_i$ 是 $g(y_i)$ 在 $\mu=Ey_i$ 处的一阶泰勒展开，于是我们把 $z_i$ 称为工作变量/调整响应变量，那么Fisher scoring算法也可以被看成调整响应变量转化版本的IRWLS，迭代步骤为：

每一步在当前 $\beta$ 下根据一阶泰勒展开计算新的 $z$ 以及 $w$ ；
以 $z$ 为新的响应变量，对 $X$ 做权重为 $w$ 的最小二乘法，更新 $\beta$ ；
重复上述步骤，直至收敛。

三、广义线性回归的假设检验

在得到参数估计后，就要对模型和参数进行统计推断和假设检验。

3.1 回归系数的检验统计量

我们先看回归系数如何检验。通过2.1小节的学习我们知道，广义线性回归的得分统计量为：

$U_j=\frac{\partial l}{\partial\beta_j}=\sum_{i=1}^n\left[ \frac{y_i-\mu_i}{Var(y_i)}x_{ij}(\frac{\partial\mu_i}{\partial g(\mu_i)}) \right]$

其中 $g(\mu_i)=g(Ey_i)=x_i^T\beta$ 。

我们也知道得分统计量的两个性质： $E(U_j)=0,j=1,\cdots,p$ ， $\mathfrak I_{jk}=E[U_jU_k]$ 。根据统计学的基本知识不难发现：

$\frac{U}{\sqrt{\mathfrak I}}\sim N(0,1)$

或者说：

$\frac{U^2}{\mathfrak I}\sim\chi^2_1$

如果写成向量形式，那么：

$\begin{aligned} &U\sim N_p(0,\mathfrak I)\\ &U^T\mathfrak I^{-1}U\sim\chi_p^2 \end{aligned}$

进一步，我们根据得分统计量，通过泰勒展开获取在极大似然估计 $\hat\beta=b$ 处的展开式：

$U(\beta)=U(b)+(\beta-b)U^{'}(b)$

对于 $U^{'}(b)$ ，我们使用期望值 $E(U^{'})=-\mathfrak I$ 来代替：

$U(\beta)=U(b)-\mathfrak I(b)(\beta-b)$

由于 $b$ 是极大似然估计，那么 $U(b)=\left.\frac{\partial l}{\partial \beta}\right|_{\beta=\hat\beta=b}=0$ ，所以化简得到：

$(b-\beta)=\mathfrak I^{-1}U$

假设 $\mathfrak I$ 非奇异且为常数，那么 $E(b-\beta)=0$ ，反映了 $b$ 就是 $\beta$ 的极大似然估计。又因为：

$E\left[ (b-\beta)(b-\beta)^T \right]=\mathfrak I^{-1}E(UU^T)\mathfrak I^{-1}=\mathfrak I^{-1}$

所以：

$(b-\beta)^T\mathfrak I(b)(b-\beta)\sim\chi_p^2$

左侧这个二次型被称为广义线性回归中的wald统计量。于是：

$b\sim N(\beta,\mathfrak I^{-1})$

你会发现，这个结论和线性回归的结论很像。事实上你可以自己尝试，令 $b=(X^TX)^{-1}X^Ty$ ，可以求解得到 $var(b)=\sigma^2(X^TX)^{-1}=\mathfrak I^{-1}$ ，和我们之前学过的结论完全一致。

3.2 模型整体显著性检验

在一般线性回归模型中，我们构造了F统计量来检验模型的显著性。而在广义线性回归中，由于响应变量服从指数分布族，F统计量便不再适用。但是，我们在构建F统计量时用到了一个思想，就是构造似然比。如果能够从似然比出发构建似然比检验(Likelihood Ratio Test, LRT)，就可以实现检验了。

事实上是可以的，我们来尝试导出。我们从对数似然函数出发，先用泰勒展开的前三项来表达 $\hat\beta=b$ 处的对数似然函数：

$l(\beta)=l(b)+(\beta-b)U(b)+\frac{1}{2}(\beta-b)^2U^{'}(b)$

对于 $U^{'}(b)$ ，我们使用期望值 $E(U^{'})=-\mathfrak I$ 来代替，写成矩阵形式：

$l(\beta)=l(b)+(\beta-b)^TU(b)-\frac{1}{2}(\beta-b)^T\mathfrak I(b)(\beta-b)$

3.1小节中我们提到 $U(b)=0$ ，所以：

$l(\beta)-l(b)=-\frac{1}{2}(\beta-b)^T\mathfrak I(b)(\beta-b)$

做一下变换：

$2[l(b)-l(\beta)]=(\beta-b)^T\mathfrak I(b)(\beta-b)\sim\chi_p^2$

也就是说对数似然函数的差的两倍服从卡方分布。我们假定似然比：

$\Lambda=\frac{\sup L(b_{H1};y)}{\sup L(b_{H0};y)}$

其中 $H0$ 代表一个更小的模型（原假设）， $H1$ 代表更大的模型（备择假设）。于是有：

$\begin{aligned} LRT&=2[l(b_{H1};y)-l(b_{H0};y)]\\ &=2[l(b_{H1};y)-l(\beta_{H1};y)]-2[l(b_{H0};y)-l(\beta_{H0};y)]+2[l(\beta_{H1};y)-l(\beta_{H0};y)] \end{aligned}$

其中第一项是备择假设下的统计量，服从 $\chi_p^2$ ， $p$ 是参数量；第二项是原假设下的统计量，服从 $\chi_m^2$ ， $m$ 是参数量；第三项是个未知的常数。因此，LRT应当渐近服从卡方分布：

$LRT=2\log\Lambda\sim\chi_{p-m,v}^2$

其中 $v$ 是一个非中心化参数。如果 $Y_i$ 服从正态分布，那么 $LRT$ 就精确服从卡方分布。

3.3 模型结构检验

在一般线性回归中，我们通过偏回归图和偏残差图检验了线性结构的存在。而在广义线性回归中，我们的主要检验目标是是否存在联系函数，即：

原假设 $H_0$ ：存在 $g$ 使得 $g(\mu_i)=x_i^T\beta$
备择假设 $H_1$ ：没有这样的 $g$ 使得 $g(\mu_i)=x_i^T\beta$

要想检验这一点，我们需要先介绍全模型的概念。全模型/饱和模型(saturated model)是指相对于数据集来说包含的参数量最多的模型。例如，如果 $N$ 个随机变量的样本 $Y_i$ 都服从正态分布 $N(\mu_i,\sigma^2)$ ，那么全模型就是参数量为 $N$ 的模型，因为每一个观测值 $y_i$ 都是 $\mu_i$ 的极大似然估计。如果拟合为线性回归模型，那么参数量就会缩减变为自变量的个数 $p$ ，因为 $\mu_i$ 的极大似然估计变为 $x_i^T\beta$ 。

也就是说，全模型是不考虑任何模型结构的模型，单纯从样本数值直接推断总体。

有了全模型的概念，那么模型结构检验也就很明显了：检验拟合模型和全模型有显著差异。那么很自然地想到用3.2小节学习的似然比检验方法。设全模型的回归系数是 $b_{max}$ ，拟合模型是 $b$ ，那么我们可以构造似然比：

$\Lambda=\frac{\sup L(b_{max};y)}{\sup L(b;y)}$

Nelder和Wedderburn在1972年定义了一个统计量Deviance，转换似然比使其服从卡方分布：

$D=2\log\Lambda=2\log\frac{L(b_{max};y)}{L(b;y)}=2[l(b_{max};y)-l(b;y)]$

这个检验被称为Deviance检验。不过你一定看出来了，这个Deviance统计量就是上面的似然比检验统计量，二者其实是一回事。

那么对于广义线性模回归来说，其全模型为 $g(\mu_i)=\psi_i$ ，即模型结构未知，那么似然比为：

$\Lambda=\frac{L(Y,\hat\psi)}{L(Y,\hat\beta)}$

因此有

$2\log\Lambda\sim\chi_r^2$

其中 $r$ 为自由度。当然，这种服从也是渐近的。

最后我们来看两个例子，学习如何计算Deviance。

例1：设一组独立样本 $y_1,\cdots,y_n$ ，且 $y_i\sim Possion(\lambda_i)$ ，即泊松分布。

对于广义线性回归来说，取 $\log\lambda_i=x^T_i\beta$ ，此时为我们后续要学习的泊松回归，那么对数似然为：

$l(Y,\beta)=\sum_{i=1}^ny_ix_i^T\beta-\sum_{i=1}^ne^{x_i^T\beta}-\sum_{i=1}^n\log y_i!$

并且在 $\beta=\hat\beta$ 处取得极大值。

对于全模型来说，对数似然为：

$l(Y,\psi)=\sum_{i=1}^ny_i\psi_i-\sum_{i=1}^ne^{\psi_i}-\sum_{i=1}^n\log y_i!$

由于在 $e^{\psi_i}=y_i$ 时取得极大值，于是当 $y_i\gt0$ 时，取 $\psi_i=\log y_i$ ；其他情况下 $e^{\psi_i}=0$ ，于是极大值为：

$\max l(Y,\psi)=\sum_{\{i|y_i\gt0\}}(y_i\log y_i-y_i)-\sum_{i=1}^n\log y_i!$

于是：

$D=2\sum_{i=1}^n\left[y_i\log\frac{y_i}{\hat\lambda_i}-(y_i-\hat\lambda_i)\right],\hat\lambda_i=e^{x_i^T\beta}$

例2：设一组独立样本 $y_1,\cdots,y_n$ ，且 $y_i\sim N(\mu_i,\sigma^2)$ ，即正态分布。

对于广义线性回归来说，以 $g(t)=t$ 为例，此时为一般线性回归模型，对数似然为：

$l=\sum_{i=1}^n\left[\frac{y_i(x_i^T\beta)-\frac{1}{2}(x_i^T\beta)^2}{\sigma^2}+c(y_i;\sigma^2)\right]$

并且在 $\beta=\hat\beta=(XX^T)^{-1}X^TY$ 处取得极大值。

对于全模型来说，每一个 $y_i$ 都是对 $\mu_i$ 的极大似然估计，因此似然函数为：

$l=\sum_{i=1}^n\left[\frac{y_i\psi_i-\frac{1}{2}\psi_i^2}{\sigma^2}+c(y_i;\sigma^2)\right]$

并且在 $\psi_i=y_i$ 处取得极大值。

于是：

$D=\frac{\sum_{i=1}^n(y_i-x_i^T\beta)^2}{\sigma^2}$

你会发现这个就是我们之前的结论：

$D=\frac{(n-r)\sum_{i=1}^n(y_i-x_i^T\beta)^2}{(n-r)\sigma^2}=\frac{(n-r)\hat\sigma^2}{\sigma^2}\sim\chi_{n-r}^2$