Lecture 09 广义线性模型（GLM）

我们已经见过了广义线性模型的一些特例，其中，响应变量 $Y$ 服从伯努利分布。现在我们可以看一下更加广义的情况，响应变量 $Y$ 不再是正态的。本章将包含以下三部分：

• 介绍
• 估计
• 推断

1. 介绍

令 $Y_1,Y_2,\dots,Y_n$ 为随机变量 $Y$ 的独立观测，其关联的协变量向量为 $\mathbf x_1,\mathbf x_2,\dots,\mathbf x_n$。

如果随机变量 $Y$ 的这些独立观测服从正态分布，那么我们可以采用线性模型来建立随机变量 $Y$ 和其关联的协变量之间的联系。对于 $i=1,\dots,n$，我们有

其中，响应变量 $Y_i$ 的期望 $\mu_i$ 是关于协变量 $\mathbf x_i$ 的线性函数。

现在，我们试图将该模型从正态分布推广到 指数族（exponential famil y） 分布：我们假设

其中，$\eta_i=\mathbf x_i^{\mathrm{T}}\boldsymbol \beta$ 是 线性预测变量（linear predictor）。

函数 $g(.)$ 被称为 连接函数（link function）：它提供了线性预测变量和响应 $Y$ 的均值之间的连接。

如果 $g(\mu_i)$ 被设置为等于指数族中的 典范参数（canonical parameter），又称 自然参数（natural parameter），那么，$g(.)$ 被称为 自然连接（natural link） 或者 典范连接（canonical link）。

如果 $Y$ 的分布是 典范形式（canonical form），那么，它的概率密度函数（pdf）$f(y\mid \theta,\phi)$ 满足

其中，$\theta$ 被称为 典范参数（canonical parameter）或者 自然参数（natural parameter），代表了 位置（location）；而 $\phi$ 被称为 散布参数（dispersion parameter），代表了 尺度（scale）。我们可以通过指定函数 $a,b,c$ 来定义不同的指数族分布。通常，$a(\phi)=\phi \big/ w$，其中，$w$ 是已知的 权重（weight）。

可以推导出

另一方面，

注意，这里有个技巧是利用 $(\ln f)’=f’\big/f$ 的性质。

结合上面的推导，可以得到

• 例子：泊松分布
  如果 $Y\overset{\mathrm d}{=}\mathrm{Poi}(\lambda)$，那么 $\ln f(y\mid \lambda)=-\lambda + y\ln \lambda - \ln y!$，所以 $\theta=\ln \lambda$。

  因此，我们有 $\ln f(y\mid \lambda)=-e^{\theta}+y\theta-\ln y!$；所以，$b(\theta)=e^{\theta}$，并且 $a(\phi)=1$。

  由上面的结果可得

  \[\begin{align} \mu &= \mathbb{E}(Y)=b'(\theta)=e^{\theta}=\lambda \\\\ \sigma^2 &= \mathrm{var}(Y)=a(\phi)b''(\theta)=e^{\theta}=\lambda \end{align}\]

• 例子：二项分布
  如果 $Y\overset{\mathrm d}{=}\mathrm{Bin}(m,p)$，那么 $\ln f(y\mid m,p)=\text{const}+y\ln \frac{p}{1-p} + m\ln(1-p)$，所以 $\theta=\ln \frac{p}{1-p}$，$p=\frac{e^{\theta}}{1+e^{\theta}}$。

  因此，我们有 $b(\theta)=-m\ln(1-p)=m\ln(1+e^{\theta})$，并且 $a(\phi)=1$。

  由上面的结果可得

  \[\begin{align} \mu &= \mathbb{E}(Y)=b'(\theta)=\dfrac{me^{\theta}}{1+e^{\theta}}=mp \\\\ \sigma^2 &= \mathrm{var}(Y)=a(\phi)b''(\theta)=\dfrac{me^{\theta}}{(1+e^{\theta})^2}=mp(1-p) \end{align}\]

2. 估计

2.1 例子：具有对数连接的泊松回归模型

我们更详细地考虑一个相当标准的案例 —— 具有对数连接的泊松回归模型（Poisson regression model with log link），该模型可以视为 GLM 的一个模板，其连接函数等于自然参数：

例子：具有对数连接的泊松回归模型

• $Y_i\overset{\mathrm d}{=}\mathrm{Poi}(\lambda_i)$，其中 $i=1,2,\dots,n$；$Y_i$ 之间彼此相互独立。
• 均值参数 $\lambda$ 依赖于协变量 $\mathbf x_1,\mathbf x_2,\dots,\mathbf x_q$。
• 自然参数 $\theta_i=\ln \lambda_i$。自然连接 $g(\lambda_i)=\ln \lambda_i$。
• 这提供了一个对数线性模型 \(\eta_i=\ln \lambda_i=\mathbf x_i^{\mathrm{T}}\boldsymbol \beta=\sum_{\ell=1}^{q}x_{i\ell}\beta_{\ell}\)
• 其矩阵形式为：$\ln \boldsymbol{\lambda}=\boldsymbol{\eta}=X\boldsymbol{\beta}$，其中 $\boldsymbol{\eta}=(\eta_1,\dots,\eta_n)^{\mathrm{T}},\boldsymbol{\lambda}=(\lambda_1,\dots,\lambda_n)^{\mathrm{T}}$。

• 联合对数似然函数（Joint log-likelihood function）

  \[\begin{align} \ell (\boldsymbol \beta)=\ln f &= -\sum_{i=1}^{n}\lambda_i+\sum_{i=1}^{n}y_i\ln \lambda_i-\sum_{i=1}^{n}\ln y! \\ &= -\sum_{i=1}^{n}e^{\sum_{\ell=1}^{q}x_{i\ell}\beta_{\ell}}+\sum_{i=1}^{n}y_i\sum_{\ell=1}^{q}x_{i\ell}\beta_{\ell}+\mathrm{const} \end{align}\]

• 记分函数（score function）

  \[\mathbf u(\boldsymbol \beta)=\left(\dfrac{\partial \ln f}{\partial \beta_j}\right)=\left(-\sum_{i=1}^{n}x_{ij}e^{\sum_{\ell=1}^{q}x_{i\ell}\beta_{\ell}}+\sum_{i=1}^{n}x_{ij}y_i\right)=X^{\mathrm{T}}(\mathbf y-\boldsymbol \lambda)\]

• 海森函数（Hessian function）

  \[H(\boldsymbol \beta)=\left(\dfrac{\partial^2 \ln f}{\partial \beta_j \partial \beta_k}\right)=\left(-\sum_{i=1}^{n}x_{ij}x_{ik}e^{\sum_{\ell=1}^{q}x_{i\ell}\beta_{\ell}}\right)=-X^{\mathrm T}\Lambda X\]

  其中，$\Lambda=\Lambda(\boldsymbol \beta)=\mathrm{diag}(\lambda_1,\lambda_2,\dots,\lambda_n)$。

因此，

• 观测信息（Observed information）

  \[J(\boldsymbol \beta)=-H(\boldsymbol \beta)=X^{\mathrm T}\Lambda X\]

• 费雪信息（Fisher information）

  \[I(\boldsymbol \beta)=-\mathbb{E}[H(\boldsymbol \beta)]=X^{\mathrm T}\Lambda X\]

• 记分函数的方差（Variance of score function）

  \[\mathrm{var}(\mathbf u(\boldsymbol \beta))=I(\boldsymbol \beta)=X^{\mathrm T}\Lambda X \overset{\mathrm{denoted}}{=}V(\boldsymbol \beta)\]

2.2 记分法（Method of scoring）

$\boldsymbol \beta$ 的 极大似然估计 MLE $\hat {\boldsymbol \beta}$ 可以通过 牛顿-拉弗森（Newton-Raphson）或者 费雪记分法（Fisher scoring）得到。

其中，

这本质上是 加权最小二乘等式（weighted least squares equation）：

因此，记分法中的 迭代步（iterative step）相当于在 $\mathbf x_i$ 上拟合一个 $\hat z_i^{(k)}$ 的加权回归，权重为 $\hat \lambda_i^{(k)}$。

注意，在每次迭代中，$\hat \lambda_i^{(k)}$ 和 $\hat z_i^{(k)}$ 的值都需要更新，因为它们都依赖于 $\hat{\boldsymbol \beta}^{(k)}$。

在广义线性模型中，我们总是可以将记分法中的迭代步像这样表示为加权回归：该过程称为 迭代重新加权最小二乘法（iteratively re-weighted least squares，IRWLS）。

在 R 语言中，函数 glm() 中实现了 IRWLS 方法：

R 代码：

1
2
3
4
glm(formula, family = poisson, data, weights, subset, 
      na.action, start = NULL, etastart, mustart, offset, 
      control = list(...), model = TRUE, method = "glm.fit", 
      x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...)

例子：拟合模型 $Y_i\overset{\mathrm d}{=}\mathrm{Poi}(\lambda_i)$，其中 $\ln \lambda_i=\beta_0+\beta_1x_i$，数据如下：

1
2
x :  -1   -1   0   0   0   0   1   1   1
y :   2    3   6   7   8   9  10  12  15

R 代码：

1
2
3
4
x <- c(-1, -1, 0, 0, 0, 0, 1, 1, 1)
y <- c(2, 3, 6, 7, 8, 9, 10, 12, 15)
pr.1 <- glm(y ~ x, family = poisson)
summary(pr.1)

输出结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Call:
glm(formula = y ~ x, family = poisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8472  -0.2601  -0.2137   0.5214   0.8788  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.8893     0.1421  13.294  < 2e-16 ***
x             0.6698     0.1787   3.748 0.000178 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 18.4206  on 8  degrees of freedom
Residual deviance:  2.9387  on 7  degrees of freedom
AIC: 41.052

Number of Fisher Scoring iterations: 4

2.3 具有自然连接的 GLM 估计

具有自然连接函数的广义指数族分布和上面的泊松分布的推导过程略有不同。

在具有自然连接的广义情况下，$a(\phi)=\phi \big / w$，我们有

其中，

令 $V=X^{\mathrm T}W\Sigma X$，那么 MLE $\hat {\boldsymbol \beta}$ 可以通过记分法求解：

其中，

在上面最后两行等式中，记分法变成了 IRWLS：

其中，

• IRWLS 过程产生了收敛的 MLE $\hat {\boldsymbol \beta}$。

• $\hat {\boldsymbol \beta}$ 的方差矩阵可以被估计为：

  \[\widehat {\mathrm{var}}(\hat {\boldsymbol \beta})=\hat \phi \left(X^{\mathrm T}W\Sigma(\hat {\boldsymbol \beta}) X \right)^{-1}\]

3. 推断

3.1 渐近似然理论

为了进行推断，我们需要了解有关统计分布的一些信息。GLM 的基础是 渐近似然理论（asymptotic likelihood theory）：

其中，$q=\mathrm {dim}(\mathbf x)$。

对于 GLM，$V$ 中并不包含 $\mathbf y$，所以 $I(\boldsymbol \beta)=\phi^{-1}V$，我们有：

更进一步，最大似然估计量满足：

通常，$V$ 中包含 $\boldsymbol \beta$，所以我们需要利用 $V(\hat{\boldsymbol \beta})$ 来计算 标准误（se）：

标准误（se）为：$\mathrm{se}(\hat \beta_j)=\sqrt{\phi v^{jj}(\hat{\boldsymbol \beta})}$ $\;\;;\;\;$近似的 $95\%$ 置信区间（CI）为：$\hat \beta_j\pm 1.96\cdot \mathrm{se}(\hat \beta_j)$

3.2 检查模型的充分性

评估模型 $M$ 的 充分性（adequacy）：比较模型 $M$ 的似然和 饱和模型（full/saturated model）$F$ 的似然。其中，饱和模型 $F$ 是一个和模型 $M$ 具有相同分布和相同连接函数，但是参数数量为 $n$ 的 GLM 模型。

如果模型 $M$ 是一个好的模型，那么 $L(\boldsymbol \beta_M)$ 将非常接近 $L(\boldsymbol \beta_F)$：它能够解释数据中的大部分变化。

我们定义检验统计量：

这被称为 剩余偏差（residual deviance）或者简称 偏差（deviance），它类似于线性模型中的 残差平方和（residual sum of squares），但是这里不需要估计 $\sigma^2$，所以现在如果 $D$ 很大的话，原因就是模型拟合不够好。

检验基于以下结果，如果模型 $M$ 是正确的，那么：

否则，$D$ 将趋于更大（这表明模型 $M$ 拟合不够好）。

因此，如果 $D<\chi_{0.95}^{2}(n-q)$，那么，模型 $M$ 是可以接受的。

回忆前面 泊松分布的例子，现在，我们对 R 代码输出中剩余偏差的意义更清楚了：

1
Residual deviance:  2.9387  on 7  degrees of freedom

它表明 $D=2.9387$，并且 $n-q=7$（因为 $n=9,\;q=2$），因此，我们拟合的模型是相当可以接受的：$\chi_{0.95}^{2}(7)=14.07>D$。

3.3 模型充分性检验背后的理论

$D$ 的分布结果来自渐近似然结果：

注意，对于 MLE $\hat {\boldsymbol \beta}$，$\mathbf u(\hat {\boldsymbol \beta})=\mathbf 0$。

由上式可得：

关于模型充分性的标准检查，更进一步的做法是查看 残差（residuals），我们可以利用前面章节中介绍过的诊断工具。

我们可以将 $D$ 表示为 $D=\sum_{i=1}^{n}d_i^2$，其中 $d_i^2$ 是来自第 $i$ 个观测的作用。那么，观测 $i$ 的 偏常残差（deviance residual）定义为：

3.4 比较嵌套模型

关于 嵌套模型（nested models）比较的检验，处理方法与广义线性模型理论的检验基本相同，唯一区别是现在不需要估计 $\sigma^2$。

所以，检验可以直接基于类似于线性模型中残差平方和的变化（即，GLM 中的剩余偏差的变化）。

假设模型 $M_0$ 嵌套在模型 $M_1$ 中。假设模型 $M_1$（$H_1$ 假设下的模型）为真。然后，

如果模型 $M_0$ 也为真，那么：

如果模型 $M_0$ 不为真，那么 $\Delta D$ 将趋近于更大。

因此，如果 $\Delta D>\chi_{0.95}^{2}(q_1-q_0)$，那么，我们拒绝模型 $M_0$。

R 代码：

1
anova(pr.1, test = "Chi")

输出结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Analysis of Deviance Table

Model: poisson, link: log

Response: y

Terms added sequentially (first to last)


     Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
NULL                     8    18.4206             
x     1   15.482         7     2.9387 8.33e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

结果显示，$D_0=18.4206$，$D_1=2.9387$，所以

如果 $M_0$ 为真，那么它将会是 $\chi^{2}(1)$ 上的一个观测。

由于 $\chi_{0.95}^{2}(1)=3.841<\Delta D$，我们拒绝 $M_0$；并且得出结论 $\beta_1\neq 0$，就像我们已经从考虑 $\hat \beta_1$ 和 $\mathrm{se}(\hat \beta_1)$ 中得到的一样。

注意：比较 $M_0\subset M_1$ 和 $M_1$ 等价于检验线性假设 $H_0: \boldsymbol \beta_{M_1-M_0}=\mathbf 0$。对于这种比较，我们已经学习了三种检验方法：LR 检验、Wald 检验 和 Score 检验。而 $\Delta D=D_0-D_1$ 检验对应的是 LR 检验。

下节内容：逻辑回归模型分析分类数据

本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。 欢迎转载，并请注明来自：YEY 的博客 同时保持文章内容的完整和以上声明信息！

• Previous
  自然语言处理 08：NLP 深度学习：循环网络
• Next
  自然语言处理 09：词汇语义学