广义线性模型（Generalized Linear Model）之一：指数分布族与glm函数

---

通常情况下，我们是使用回归分析和方差分析去探究线性模型，它们可以通过一系列连续型和/或类别型预测变量来预测正态分布的响应变量。但在许多情况下，假设因变量为正态分布（甚至连续型变量）并不合理，例如下面这几种情况：

1）结果变量可能是类别型的

• 二值变量（比如：是/否、通过/未通过、活着/死亡）
• 多分类变量（比如：差/良好/优秀）都显然不是正态分布

2）结果变量可能是计数型的

• 比如，一周交通事故的数目，每日酒水消耗的数量
• 这类变量都是非负的有限值，而且它们的均值和方差通常都是相关的（正态分布变量间不是如此，而是相互独立）

广义线性模型扩展了线性模型的框架，它包含了非正态因变量的分析。本文旨在简要概述广义线性模型，并梳理其基本理论和R软件实现方法。

一、广义线性模型

许多广泛应用的、流行的数据分析方法其实都归属于广义线性模型框架。所以，在第一部分，笔者将简短回顾这些方法背后的理论。

（一）经典线性模型

现假设我们想要对响应变量 Y 和 p 个预测变量 X1…Xp 间的关系进行建模。在标准线性模型中，可以假设Y呈正态分布，关系的形式为：
$${\mu }_Y={\beta }_0+{\sum }_{j=1}^p{\beta }_j{X}_j$$

该等式表明响应变量的条件均值是预测变量的线性组合。参数 βj 指一单位 Xj 的变化造成的 Y 预期的变化，β0 指当所有预测变量都为 0 时 Y 的预期值。

对于该等式，可以通俗地理解为：给定一系列X变量的值，赋予X变量合适的权重，然后将它们加起来，便可预测Y观测值分布的均值。

值得注意的是，你并没有对预测变量Xj做任何分布的假设。与Y不同，它们不需要呈正态分布。实际上，它们常为类别型变量（比如方差分析设计）。另外，对预测变量使用非线性函数也是允许的，比如你常会使用预测变量 $X_{}^2$ 或者 $X_1\times X_2$ ，只要等式的参数（β0,β1,…,βp）为线性即可。

总的来说，经典的线性模型 (LM) 可以分解为三个要素：

• LM1：随机要素，即Y服从正态分布，${\mu }_Y=E\left(y_i\right)$
• LM2：系统要素，$\eta ={\beta }_j{X}_j$
• LM3：连接要素，$\eta =\mu$

（二）广义线性模型

在上述经典线性回归模型中，通常假设因变量服从正态分布，方差为常数，解释变量通过线性相加关系直接影响因变量本身。

而广义线性模型（Generalized Linear Model，GLM）则假设因变量来自于指数型分布族，其方差随着均值而变化，解释变量通过线性相加关系对因变量的期望值的某种变换产生影响。

广义线性模型扩展了经典的线性模型，因此它适用于更广范围的数据分析问题。广义线性模型拟合的形式为：

$$g\left({\mu }_Y\right)={\beta }_0+{\sum }_{j=1}^p{\beta }_j{X}_j$$

其中 $g\left({\mu }_Y\right)$ 是条件均值的函数（称为连接函数）。另外，可放松Y为正态分布的假设，改为Y服从指数分布族中的一种分布即可。设定好连接函数和概率分布后，便可以通过最大似然估计的多次迭代推导出各参数值。

广义线性模型 (GLM) 可以分解为三个要素：

• GLM1：随机要素，因变量Y（或误差项）服从指数型分布族中的一个分布。指数型分布族包括许多常见分布，如正态分布、泊松分布、逆高斯分布、二项分布、伽玛分布等
• GLM2：系统要素，即解释变量的线性组合，表示为 $\eta ={\beta }_j{X}_j$ ，与经典线性回归模型没有区别，仍保持线性结构
• GLM3：连接要素，$g\left({\mu }_Y\right)=\eta$ 或 $E\left(y_i\right)={\mu }_Y=g_{}^{-1}\left(\eta \right)$。其中 g 严格单调且可导，称为联系函数，它建立了随机要素与系统要素之间的关系。

可见，在广义线性模型中，对因变量的预测值并不直接等于解释变量的线性组合，而是该线性组合的一个函数变换。

二、指数分布族

为便于建立各个变量与参数之间的关系（求导、求期望等），广义线性模型要求随机变量 Y 的分布来自指数分布族.。

指数分布族(Exponential Family) 是一系列分布的统称，包含连续和离散的相关分布。例如，正态分布(Gaussian)、泊松分布（Poisson）、二项分布(Bernoulli)、指数分布(exponential)、Gamma分布、多项式分布(multivariate)等。

（一）定义

单参数指数族是一组概率分布，其概率密度函数（或概率质量函数，对于离散分布的情况）可以表示为:

$$f_X(x\mid \theta )=h(x)\exp [\eta (\theta )\cdot T(x)-A(\theta )]$$

通常给出的另一种等效形式是:

$$f_X(x\mid \theta )=h(x)g(\theta )\exp [\eta (\theta )\cdot T(x)]f_X(x\mid \theta )=h(x)g(\theta )\exp [\eta (\theta )\cdot T(x)]$$

或者也等价为：

$$f_X(x\mid \theta )=\exp [\eta (\theta )\cdot T(x)-A(\theta )+B(x)]f_X(x\mid \theta )=\exp [\eta (\theta )\cdot T(x)-A(\theta )+B(x)]$$

总的来说，凡是概率（密度）函数能写成如下形式的分布都属于指数族分布：
$$p(x\mid \eta )=h(x)\exp [{\eta }_{}^T\cdot T(x)-A(\eta )]$$

其中 T(x), h(x), η, A(η) 是已知函数：

$\eta$ 为自然参数(natural parameter)，可以是向量形式
$T(x)$ 为充分统计量(sufficient statistic)
$A(\eta )$ 为累计函数(cumulant function)，作用是确保概率和为1
$h(x)$ 为underlying measure，必须是非负的

（二）典型分布转化

1. The Bernoulli distribution

Bernoulli分布的转化：

相关参数为：

2. The Poisson distribution

泊松分布的标准形式为：

相关参数为：

3. The Gaussian distribution

正态分布：

相关参数为：

4. The multinomial distribution

多项式分布：

rewrite the multinomial distribution as follows:

三、R语言：glm()函数

R中可通过glm()函数（还可用其他专门的函数）拟合广义线性模型。它的形式与lm()类似，只是多了一些参数。函数的基本形式为：

glm(formula, family, data, weights, subset, ...)

唯一新的参数是family，用来指定方差和链接函数。family有六种选择：

以下列出了概率分布（family）和相应默认的连接函数（function）：

• Binomial(link=“logit”)
• Gaussian(link=“identity”)·
• Gamma(link=“inverse”)·
• Inverse.gaussian (link=“1/mu^2”)·
• Poisson(link=“log”)·
• Quasi(link=“identity”, variance=“constant”)·
• Quasibinomial (link=“logit”)·
• Quasipoisson (link=“log”)

以Dobson (1990) 第 93 页的随机对照试验为例：

# NOT RUN {
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
anova(glm.D93)
# }
# NOT RUN {
summary(glm.D93)
# }
# NOT RUN {
## Computing AIC [in many ways]:
(A0 <- AIC(glm.D93))
(ll <- logLik(glm.D93))
A1 <- -2*c(ll) + 2*attr(ll, "df")
A2 <- glm.D93$family$aic(counts, mu=fitted(glm.D93), wt=1) +
  2 * length(coef(glm.D93))
stopifnot(exprs = {
  all.equal(A0, A1)
  all.equal(A1, A2)
  all.equal(A1, glm.D93$aic)
})


# }
# NOT RUN {
## an example with offsets from Venables & Ripley (2002, p.189)
utils::data(anorexia, package = "MASS")

anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
                family = gaussian, data = anorexia)
summary(anorex.1)
# }
# NOT RUN {
# A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
clotting <- data.frame(
  u = c(5,10,15,20,30,40,60,80,100),
  lot1 = c(118,58,42,35,27,25,21,19,18),
  lot2 = c(69,35,26,21,18,16,13,12,12))
summary(glm(lot1 ~ log(u), data = clotting, family = Gamma))
summary(glm(lot2 ~ log(u), data = clotting, family = Gamma))
## Aliased ("S"ingular) -> 1 NA coefficient
(fS <- glm(lot2 ~ log(u) + log(u^2), data = clotting, family = Gamma))
tools::assertError(update(fS, singular.ok=FALSE), verbose=interactive())
## -> .. "singular fit encountered"

# }
# NOT RUN {
## for an example of the use of a terms object as a formula
demo(glm.vr)
# }

Quality of Fit of a GLM in R

一些衡量标准可以说明 GLM 对于给定概率分布的拟合程度：

1）Deviance

偏差是衡量 GLM 拟合质量或优度的指标。偏差值过大就表示不合适。
偏差度量通常有两种类型：

• Null deviance： shows how well the response variable predicted by a model that includes only the intercept.
• Residual deviance ： shows the same but when the independent variables are included.

2）Fisher scoring

Fisher 评分使用 Newton-Raphson 算法来找到模型的最佳拟合。算法首先根据估计值拟合模型，然后它检查是否可以通过使用不同的估计来改进模型。这个过程一直持续到算法发现新值不会导致任何改进。

3）Hosmer-Lemeshow test

Hosmer-Lemeshow 检验用来检查模型与观测数据之间的相似性。如果模型与观测数据之间的差异不显着，我们可以说模型拟合得很好。

4）Akaike Information criterion

AIC 有助于避免不相关的预测变量不必要地增加模型的复杂性。选择具有最小 AIC 的模型会产生最佳拟合。