用R语言进行数据分析：常规和广义线性模型

线性模型

对于常规的多重模型（multiple model）拟合，最基本的函数是lm()。 下面是调用它的方式的一种改进版：

     > fitted.model <- lm(formula, data = data.frame)

例如

     > fm2 <- lm(y ~ x1 + x2, data = production)

将会拟合 y 对 x1 和 x2 的多重回归模型(和一个隐式的截距项)。

一个重要的(技术上可选)参数是 data = production。它指定任何构建这个模型的参数首先必须来自 数据框 production。 这里不需要考虑数据框 production 是否被绑定在搜索路径中。

广义线性模型

广义线性建模是线性模型在研究响应值的非正态分布以及非线性模型的简洁直接的线性转化 时的一种发展。 广义线性模型 是基于下面一系列 假设前提的：

• 有一个响应变量 y和一系列有趣的刺激变量（stimulus variable） x_1, x_2, …。 这些刺激变量决定响应变量的最终分布。
• 刺激变量仅仅通过一个线性函数影响响应值y 的分布。 该线性函数称为 线性预测器（linear predictor），常常写成

            eta = beta_1 x_1 + beta_2 x_2 + ... + beta_p x_p,

  因此 x_i 当且仅当 beta_i 等于0时对 y 的分布没有影响。

• y 分布的形式为

            f_Y(y; mu, phi)
              = exp((A/phi) * (y lambda(mu) - gamma(lambda(mu))) + tau(y, phi))

  其中 phi 是度量参数（scale parameter）(可能已知)，对所有观测 恒定；A 是一个先验的权重，假定知道但是 可能随观测不同有所不同；mu 是 y 的均值。 也就是说假定 y 的分布是由 均值和一个可能的度量参数决定的。

• 均值 mu 是线性预测器的平滑可逆函数（smooth invertible function）：

            mu = m(eta),    eta = m^{-1}(mu) = ell(mu)

  该可逆函数 ell() 称为 关联函数（link function）。

这些假定比较宽松，足以包括统计实践中大多数有用的统计模型， 但也足够严谨，使得可以发展计算和推论中一致的方法（ 至少可以近似一致）。 读者如果想了解这方面最新的进展，可以 参考 McCullagh & Nelder (1989) 或者 Dobson (1990)。

族

R 提供了一系列广义线性建模工具，从类型上来说包括 gaussian, 二项式, poisson, 反 gaussian 和 gamma 模型的响应变量分布以及 在响应变量分布没有明确给定时的逆似然（quasi-likelihood）模型。 在后者，方差函数（variance function） 可以认为是均值的函数，但是在另外一些情况下， 该函数可以由响应变量的分布得到。

每一种响应分布允许各种关联函数将均值和线性预测器关联起来。 这些自动关联函数如 下表所示：

这些用于模型构建过程中的响应分布，关联函数和各种 其他必要的信息统称为 广义线性模型的族（family）。

glm()函数

既然响应的分布仅仅通过单一的一个线性函数依赖于 刺激变量，那么用于线性模型的机制同样 可以用于指定一个广义模型的线性部分。 但是族必须以一种不同的方式指定。

R 用于广义线性回归的函数是 glm()， 它的使用形式为

     > fitted.model <- glm(formula, family=family.generator, data=data.frame)

和 lm() 相比，唯一的一个新特性就是描述族的参数family.generator。 它是产生函数和表达式列表的函数名字。这些函数 用于定义和控制模型的构建与计算过程。 尽管开始看起来有点复杂， 但它们非常容易 使用。

这些名字是标准的。程序给定的族生成器可以参见 Families 列表中 的“族名”。当选择一个关联函数时， 该关联函数名和族名可以同时在括弧里面作为 参数设定。在 拟（quasi） 家族里面，方差函数也是以这种方式设定。

一些例子可能会使这个过程更清楚。

gaussian族

命令

     > fm <- glm(y ~ x1 + x2, family = gaussian, data = sales)

和下面的命令结果一致。

     > fm <- lm(y ~ x1+x2, data=sales)

但是效率上，前者差一点。注意，gaussian 族没有相关参数， 因此它不提供关联函数的。 如一个问题需要用非标准关联函数的 gaussian 族， 那么只能采用我们后面讨论的 拟族。

二项式族

考虑 Silvey (1970) 提供的一个小的例子。

在 Kalythos 的 Aegean 岛上，男性居民常常患有 一种先天的眼科疾病，并且随着年龄的增长而变的愈显著。 现在搜集了各种年龄段岛上男性居民的样本，同时记录了盲眼的数目。 数据显示如下：

我们想知道的是这些数据是否吻合 logistic 和 probit 模型， 并且分别估计各个模型的 LD50，也就是一个男性居民盲眼的概率 为50%时候的年龄。

如果 y 和 n 是年龄为 x 时的盲眼数目和检测 样本数目，两种模型的形式都为 y ~ B(n, F(beta_0 + beta_1 x))， 其中在 probit 模型中， F(z) = Phi(z) 是标准的正态分布函数，而在 logit 模型 (默认)中， F(z) = e^z/(1+e^z)。 这两种模型中， LD50 = – beta_0/beta_1 ，即分布函数的参数为0时 所在的点。

第一步是把数据转换成数据框。

     > kalythos <- data.frame(x = c(20,35,45,55,70), n = rep(50,5),
                              y = c(6,17,26,37,44))

在 glm() 拟合二项式模型时，响应变量 有三种可能性：

• 如果响应变量是向量， 则假定操作二元（binary） 数据，因此要求是0/1向量。
• 如果响应变量是双列矩阵，则假定第一列 为试验成功的次数第二列 为试验失败的次数。
• 如果响应变量是因子，则第一水平作为失败 (0) 考虑而其他的作为`成功'(1) 考虑。

我们采用的是第二种惯例。我们在数据框中 增加了一个矩阵：

     > kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)

为了拟合这些模型，我们采用

     > fmp <- glm(Ymat ~ x, family = binomial(link=probit), data = kalythos)
     > fml <- glm(Ymat ~ x, family = binomial, data = kalythos)

既然 logit 的关联函数是默认的，因此我们可以在第二条命令中省略该参数。 为了查看拟合结果，我们使用

     > summary(fmp)
     > summary(fml)

两种模型都拟合的很好。为了计算 LD50，我们可以 利用一个简单的函数：

     > ld50 <- function(b) -b[1]/b[2]
     > ldp <- ld50(coef(fmp)); ldl <- ld50(coef(fml)); c(ldp, ldl)

从这些数据中得到的年龄分别是43.663年和 43.601年。

Poisson 模型

在 Poisson 族中，默认的关联函数是 log。在实际操作中， 这一族常常用于拟合计数资料的 Poisson 对数线性模型。 这些计数资料的实际分布往往符合二项式分布。 这是一个非常重要而又庞大的话题，我们不想在这里深入展开。 它构成了非-gaussian 广义模型内容 的很大一部分。

有时候，实践中产生的 Poisson 数据在对数或者平方根 转化后可当作正态数据处理。 作为后者的另一种选择是，一个 Poisson 广义线性模型可以通过下面的例子 拟合：

     > fmod <- glm(y ~ A + B + x, family = poisson(link=sqrt),
                   data = worm.counts)

拟似然模型

对于所有的族，响应变量的方差依赖于均值并且拥有 作为系数（multiplier）的尺度参数。 方差对均值的依赖方式是响应分布的一个特性； 例如对于poisson分布 Var(y) = mu。

对于拟似然估计和推断，我们不是设定精确的响应分布而是 设定关联函数和方差函数的形式。因为关联函数和方差函数都依赖于均值。 既然拟似然估计 和 gaussian 分布使用的技术非常相似， 因此这一族顺带提供了一种用非标准关联函数或者方差函数 拟合gaussian模型的 方法。

例如，考虑非线性回归的拟合 y = theta_1 z_1 / (z_2 – theta_2) + e 同样还可以写成 y = 1 / (beta_1 x_1 + beta_2 x_2) + e 其中 x_1 = z_2/z_1, x_2 = -1/x_1, beta_1 = 1/theta_1, and beta_2 = theta_2/theta_1。 假如有适合的数据框，我们可以如下 进行非线性拟合

     > nlfit <- glm(y ~ x1 + x2 - 1,
                    family = quasi(link=inverse, variance=constant),
                    data = biochem)