R官方入门教程(11)R中的统计模型

R是一门著名的可用于数据和统计分析的程序语言，本文翻译自R软件官方文档教程《An Introduction to R》，仅供学习和参考。

11 R中的统计模型

本节假定读者对统计方法有一定的了解，特别是回归分析和方差分析。后面我们还会假定读者对广义线性模型和非线性模型也有所了解。

拟合统计模型的基本要求在R中已经得到了充分的定义，从而使得我们可以轻松构建适用于广泛问题的通用工具。

R 提供了一系列紧密联系的统计模型拟合的工具，使得拟合工作变得简单。正如我们在预备与简介中提到的一样，基本的屏幕输出是简洁的，因此用户需要调用一些辅助函数来提取更多的细节信息。

11.1 定义统计模型中的公式

以下统计模型的模板是具有独立、均方差误差的线性回归模型:
$$y_i=\sum _{j=0}^p{\beta }_j{x}_{ij}+e_i,e_i\sim N(0,{\sigma }^2),i=1,\dots ,n$$
改写成矩阵形式:
$$y=X\beta +e$$
其中 $y$ 是响应向量，$X$ 是模型矩阵或者设计矩阵。$X$ 的列 $x_0,x_1,\dots ,x_p$ 是决定变量。通常，$x_0$ 列都是1，用来定义截距项。

例子

在给出正式的说明之前，先举几个例子来说明情况。

假设 y,x,x0,x1,x2,$\dots$ 是数值变量，X 是矩阵，A,B,C,$\dots$ 是因子。

下面的例子中，上面给出公式，下面给出该公式的统计模型的描述。

y~x

y~x+1

二者都反映了一个相同的 y 对 x 的简单线性模型。第一个公式包含了一个隐式的截距项，而第二个则是一个显式的截距项。

y~0+x

y~-1+x

y~x-1

y 对 x 过原点的简单线性模型(也就是说，没有截距项)。

log(y) ~ x1 + x2

y 的变换形式log(y) 对x1 和x2 进行的多重回归(有一个隐式的截距项)。

y ~ poly(x,2)

y ~ 1 + x + I(x^2)

y 对x 的 二 次 多 项 式 回 归 。 第 一 种是使用正交多项式，第二种则显式地注明各项的幂次。

y ~ X + poly(x,2)

y 利用模型矩阵 X 和二次多项式项x进行多重回归。

y~A

y 的单因素方差分析模型，类别由A 决定。

y ~ A + x

y 的单因素协方差分析模型，类别由A 决定，协方差项为x

y ~ A*B

y ~ A + B + A:B

y ~ B %in% A

y ~ A/B

y 对A 和B的非可加两因子方差分析模型。前两个公式表示相同的交叉分类设计，后两个公式表示相同的嵌套分类设计。抽象一点说，这四个公式指明同一个模型子空间。

y ~ (A + B + C)^2

y ~ A*B*C - A:B:C

三因子实验。该模型只包含一个主效应和两个因子的交互效应。这两个公式等价。

y ~ A * x

y ~ A/x

y ~ A/(1 + x) - 1

在A的各个Level独立拟合y 对x 的简单线性回归。三个公式的编码不一样。最后一个公式会对A 各个Level分别估计截距项和斜率项。

y ~ A*B + Error(C)

一个包含两个处理因子A 和B以及由因子C 决定的误差分层的实验。如在裂区实验设计中，所有区组（还包括子区组）都是由因子C决定的。

操作符 ~ 用来定义 R 的模型公式model formula。一个普通的线性模型公式可以表示为

response ~ op_1 term_1 op_2 term_2 op_3 term_3 ...

其中,

• response是一个作为响应变量的向量或者矩阵，或者是一个值为向量/矩阵的表达式
• op_i是一个操作符。它要么是+ 要么是-，分别代表在一个模型中加入或者去掉某一项(公式第一项的操作符是可选的）
• term_i可以
• • 是一个向量/矩阵表达式或者1
  • 因子
  • 是一个由因子、向量或矩阵通过公式操作符连接产生的公式表达式

基本上，公式中的项决定了模型矩阵中的列要么被加入要么被去除。1 表示截距项，并且默认就已加入模型矩阵，除非显式地去除这一选项。

R中的公式操作符与用于Glim 和Genstat程序中的Wilkinson& Rogers 标记符相似。但操作符. 在 R 里面变成了:,这是一个不可避免的改变,因为.在R 里面是可以出现在变量名中的。

这些公式操作符总结如下(参考Chambers & Hastie, 1992, p.29):

Y ~ M Y被建模为M。

M_1 + M_2 同时包括M_1和M_2项

M_1 - M_2 包括M1项但忽略M_2项

M_1 : M_2 M_1和M_2项的张量积。如果两项都是因子，那么产生“子类”因子。

M_1 %in% M_2 和M_1 : M_2 类似，但编码方式不一样。

M_1 * M_2 M_1 + M_2 + M_1:M_2
M_1 / M_2 M_1 + M_2 %in% M_1

M^n M 的所有各项以及所有到n阶为止的“交互作用”项

I(M) 隔离M。M内所有操作符当一般的算术运算符处理。并且该项出现在模型矩阵中。

请注意，在通常括起函数参数的圆括号内，所有操作符都具有其正常的算术含义。I() 是一个恒等函数，它使得常规的算术运算符可以用在模型公式中。

还要特别注意模型公式仅仅指定了模型矩阵的列项，但暗含了对参数项的指定。在某些情况下可能不是这样的，如非线性模型的参数指定。

11.1.1 对照

我们至少需要知道模型公式是如何指定模型矩阵的列的。

如果我们面临的是连续变量，这很容易，因为每个变量提供模型矩阵的一列(如果截距项包含在模型中，将额外增加一列1)。

那么如果我们面临的是一个k-levels的因子a呢?对于无序因子和有序因子来说，答案是不同的。对于无序因子，因子的第二级、…、第k级指标生成模型矩阵中的k-1列。((因此隐含的参数设置就是把其他水平和第一个水平的响应程度进行比较)。对于有序因子，k − 1列是在1*, . . . , k* 上的正交项，并且忽略常数项。

虽然目前答案已经很复杂了，但这并不是全部。首先，如果在包含因子项的模型中省略了截距项，这一项将会被编入指示所有因子水平的k 列中。其次整个行为可以通过options 设置参数contrasts 而改变。R 的默认设置为

options(contrasts = c("contr.treatment", "contr.poly"))

提这些内容的主要原因是R 和 S 对无序因子采用不同的默认值。 S采用Helmert对照。因此，当你需要比较你的结果和某本书上或论文上用SPLUS 代码的结果时，你必须设置

options(contrasts = c("contr.helmert", "contr.poly"))

这是一个经过认真考虑的改变。因为处理对照(R默认) 对于新手是比较容易理解的。

还没有结束，因为可以使用函数contrast和C为模型中的每个项设置要使用的对照方案。

虽然细节很复杂，但在正常情况下，R中的模型公式会生成统计学专家所期望的各种模型。例如，利用关联项而非主要效应的模型拟合常常会产生令人惊讶的结果，不过这些仅仅为统计学专家设计。

11.2 线性模型

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

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

例如

> df <- data.frame(x1=c(1,2,3,4,5),x2=c(2,3,4,2,2),y=c(1,3,2,4,2))
> fm2 <- lm(y ~ x1 + x2, data = df)
> fm2
Call:
lm(formula = z ~ x + y, data = df)

Coefficients:
(Intercept)            x            y  
    1.40645      0.30323      0.03226 

将拟合 y 对 x1 和 x2的多元回归模型(包含隐式截距项)。

重要的(但在技术上是可选的)参数data = df指定构造模型所需的任何变量都应该首先来自df数据框。无论在搜索路径上是否附加了数据框df，情况都是如此。

11.3 提取模型信息的泛型函数

lm() 的返回值是一个模型拟合结果对象；技术上就是属于类lm的一个结果列表。与拟合模型相关的信息可以用适合于lm对象的泛型函数来进行显示、提取、绘图等等。这包括

add1 coef effects kappa predict residuals
alias deviance family labels print step
anova drop1 formula plot proj summary

其中一些常用的泛型函数有

函数	说明
anova(object_1, object_2)	将子模型与外部模型进行比较，生成方差分析表。
coef(object)	提取回归系数(矩阵)。全称：coefficients(object).
deviance(object)	残差平方和，若有权重可加权。
formula(object)	提取模型公式信息。
plot(object)	产生四个图，显示残差，拟合值和一些诊断。
predict(object, newdata=data.frame)	提供的数据框必须有同原始变量一样标签的变量。结果是对应于data.frame中决定变量预测值的向量或矩阵。
print(object)	简要打印一个对象的内容。常常隐式使用。
residuals(object)	提取残差(矩阵)，有权重时可加权。缩写形式:resid(object)。
step(object)	通过增加或者减少模型中的项并且保留层次来选择合适的模型。在逐步搜索过程中，AIC (Akaike信息规范)值最大的模型将会被返回
summary(object)	打印回归分析结果的综合摘要
vcov(object)	返回拟合模型对象的主要参数的方差-协方差矩阵。

11.4 方差分析和模型比较

模型拟合函数aov(formula, data=data.frame) 和函数lm() 非常相似，上面列出的泛型函数对它同样适用。

需要注意的是aov() 还允许分析多误差层次的模型，如裂区实验设计，利用区组内信息进行的平衡不完全区组设计等。模型公式

response ~ mean.formula + Error(strata.formula)

指定了一个多层次实验设计，误差层由strata.formula 定义。最简单的情况是，strata.formula只是一个因子。它定义了一个双层次的实验，也就是研究在这些因子的水平内或者水平间的实验响应。

例如，对于所有决定变量因子，模型公式如下:

> fm <- aov(yield ~ v + n*p*k + Error(farms/blocks), data=farm.data)

这通常用于描述具有均值模型v + n*p*k和三个误差层次的实验，即“农场之间”，“农场内”，“区组之间”和“区组内”。

11.4.1 方差分析表

需要注意，方差分析表是针对一系列拟合模型的。

在模型序列的特定位置增加特定的项 会使残差平方和降低。因此仅仅在正交实验中，模型中增加项的次序是没有影响的。

在多层实验设计中，程序首先把响应值依次投射到各个误差层次上，并且用均值模型去拟合各个投射。细节内容可以参考Chambers &Hastie (1992)。

与默认的完整方差分析表相比，更灵活的替代方法是直接使用anova()函数比较两个或多个模型。

> anova(fitted.model.1, fitted.model.2, ...)

结果将是一个方差分析表以显示依次加入的拟合模型的差异。需要比较的拟合模型常常是等级序列（hierarchical sequence）。这个和默认的设置实际上没有差别，只是使它更容易理解和控制。

11.5 更新拟合模型

函数update() 是一个非常便利的函数。它允许拟合一个比原先模型只增加或减少一个项的模型。它的调用形式为

> new.model <- update(old.model, new.formula)

在new.formula 中，公式中包含的.， 表示old.model模型中的对应部分。例如

> df <- data.frame(x1 = c(1, 2, 3, 4, 3, 2, 9, 10), x2 = c(2, 4, 4, 5, 3, 4, 5, 6), x3 = c(1, 4, 5, 6, 1, 3, 4, 8),
                 x4 = c(4, 5, 6, 9, 4, 5, 6, 7), x5 = c(1, 6, 7, 8, 2, 3, 5, 6), x6 = c(6, 4, 8, 9, 4, 4, 3, 1),
                 y = c(8, 9, 10, 11, 12, 21, 11, 2))
> fm05 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = df)
> fm05
Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4           x5  
    -8.5068      -2.0530      10.3460      -3.7048       0.6450      -0.8073  

> fm6 <- update(fm05, . ~ . + x6)
> fm6
Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4           x5  
    -24.947       -1.774       19.058       -2.719       -4.643       -3.094  
         x6  
      3.413  

> smf6 <- update(fm6, sqrt(.) ~ .)
> smf6
Call:
lm(formula = sqrt(y) ~ x1 + x2 + x3 + x4 + x5 + x6, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4           x5  
    -2.0647      -0.2389       2.7631      -0.5625      -0.6620      -0.3400  
         x6  
     0.5106  

这将分别拟合从数据框df中得到的五个变量的多重回归模型，拟合额外增加一个变量的六个回归量的模型，和进一步对响应值进行平方根变换后的模型拟合。

注意参数data= 在最开始调用模型拟合函数的时候指定。这个信息将会通过拟合模型对象传递给函数update() 及其相关函数。

符号. 也可以用在别的情况下，不过含义会稍微不同。如

> fmfull <- lm(y ~ . , data = df)
> fmfull

Call:
lm(formula = y ~ ., data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4           x5  
    -24.947       -1.774       19.058       -2.719       -4.643       -3.094  
         x6  
      3.413  

它将会拟合y 对数据框df中所有其他变量回归的模型.

其他研究逐步回归的函数有add1(), drop1() 和step()。 从字面上就可以看出这些函数的意义，更多细节可以参考帮助文档。

11.6 广义线性模型

广义线性模型是线性模型的一种拓展，它通过一种简洁而又直接的方式使得线性模型既适合非正态分布的响应值又可以进行线性变换。广义线性模型基于下面的一系列假设前提：

• 有一个待研究的响应变量 $y$ 和一系列刺激变量$x_1,x_2,...$。这些刺激变量决定响应变量的最终分布。
• 刺激变量仅仅通过一个线性函数影响响应变量 $y$ 的分布。该线性函数称为线性预测器，常常写成

$$\eta ={\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_p{x}_p$$

因此 $x_i$ 当且仅当 ${\beta }_i$ = 0 时对 $y$ 的分布没有影响。

• $y$ 分布的形式为

f Y ( y ; μ , φ ) = exp ⁡ [ A φ { y λ ( μ ) − γ ( λ ( μ ) ) } + τ ( y , φ ) ] f_{Y}(y;\mu,\varphi)=\exp\left[\frac{A}{\varphi}\left\{y\lambda(\mu)-\gamma\left(\lambda(\mu)\right)\right\}+\tau(y,\varphi)\right] fY​(y;μ,φ)=exp[φA​{yλ(μ)−γ(λ(μ))}+τ(y,φ)]

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

• 均值 $\mu$ 是线性预测器的平滑可逆函数：

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

其中的反函数 $\ell (\ast )$ 被称为关联函数。

这些假定比较宽松，足以囊括统计实践中大多数有用的统计模型，同时也足够严谨，使得可以开发出与参数估计和统计推论中一致的方法（至少可以近似一致）。

11.6.1 族

R 提供了一系列广义线性建模工具，从类型上来说包括高斯，二项式，泊松，逆高斯和伽马模型的响应变量分布以及响应变量分布无须明确给定的拟似然模型。在后者，方差函数可以由均值的函数指定，但在其它情况下，该函数可以由响应变量的分布得到。

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

族名	关联函数
binomial	logit, probit, log, cloglog
gaussian	identity, log, inverse
Gamma	identity, inverse, log
inverse.aussian	1/mu^2, identity, inverse, log
poisson	identity, log, sqrt
quasi	logit, probit, cloglog, identity,inverse, log, 1/mu^2, sqrt

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

11.6.2 glm()函数

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

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

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

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

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

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

gaussian族

> fm <- glm(y ~ x1 + x2, family = gaussian, data = df)
> fm
Call:  glm(formula = y ~ x1 + x2, family = gaussian, data = df)

Coefficients:
(Intercept)           x1           x2  
      9.610       -1.183        1.435  

Degrees of Freedom: 7 Total (i.e. Null);  5 Residual
Null Deviance:	    194 
Residual Deviance: 140.9 	AIC: 53.65

和

> fm <- lm(y ~ x1+x2, data=df)
> fm
Call:
lm(formula = y ~ x1 + x2, data = df)

Coefficients:
(Intercept)           x1           x2  
      9.610       -1.183        1.435 

结果一致.

但是效率上，前者差一点。注意，高斯族没有自动提供关联函数设定的选项，因此不允许设置参数。如一个问题需要用非标准关联函数的高斯族，那么只能采用我们后面讨论的拟族。

二项式族

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

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

Age: 20 35 45 55 70
No.: tested: 50 50 50 50 50
No.: blind: 6 17 26 37 44

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

如果 $y$ 和 $n$ 分别是年龄为 $x$ 时的盲眼数目和检测样本数目，两种模型的形式都为
$$y\sim \mathrm{B}(n,F({\beta }_0+{\beta }_{1}x))$$
其中在probit模型中，$F(z)=\Phi (z)$ 是标准的正态分布函数，而在logit 模型(默认)中，$F(z)=\frac{e^z}{(1+e^z)}$ 。这两种模型中，
$$\mathrm{LD}50=-\frac{{\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() 拟合二项式模型时，响应变量有三种可能性：

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

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

> kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)
> kalythos
   x  n  y Ymat.1 Ymat.2
1 20 50  6      6     44
2 35 50 17     17     33
3 45 50 26     26     24
4 55 50 37     37     13
5 70 50 44     44      6

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

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

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

> summary(fmp)

Call:
glm(formula = Ymat ~ x, family = binomial(link = probit), data = kalythos)

Deviance Residuals: 
       1         2         3         4         5  
-0.15582   0.02545  -0.08009   0.51246  -0.40097  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.102270   0.276287  -7.609 2.76e-14 ***
x            0.048147   0.005885   8.181 2.82e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 82.14455  on 4  degrees of freedom
Residual deviance:  0.45473  on 3  degrees of freedom
AIC: 24.27

Number of Fisher Scoring iterations: 4


> summary(fml)
 
 
Call:
glm(formula = Ymat ~ x, family = binomial, data = kalythos)

Deviance Residuals: 
      1        2        3        4        5  
-0.1797   0.1157  -0.1182   0.3791  -0.3372  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.53778    0.50232  -7.043 1.88e-12 ***
x            0.08114    0.01082   7.498 6.47e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 82.14455  on 4  degrees of freedom
Residual deviance:  0.31707  on 3  degrees of freedom
AIC: 24.132

Number of Fisher Scoring iterations: 4

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

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

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

Poisson 模型

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

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

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

拟似然模型

对于所有的族，响应变量的方差依赖于均值并且拥有作为乘数的尺度参数。方差对均值的依赖方式是响应分布的一个特性；例如对于poisson分布 $\mathrm{Var}[y]=\mu$。

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

例如，考虑非线性回归的拟合
$$y=\frac{{\theta }_1{z}_1}{z_2-{\theta }_2}+e$$
同样还可以写成
$$y=\frac{1}{{\beta }_1{x}_1+{\beta }_2{x}_2}+e$$
其中 $x_1=\frac{z_2}{z_1},x_2=-\frac{1}{x_1},{\beta }_1=\frac{1}{{\theta }_1},{\beta }_2=\frac{{\theta }_2}{{\theta }_1}$ 。假如有适合的数据框，我们可以如下进行非线性拟合

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

11.7 非线性最小二乘法和最大似然法模型

拟合非线性模型的一种办法就是使误差平方和(SSE)或残差平方和最小。如果观测到的误差极似正态分布，这种方法是非常有效的。

下面的例子来自Bates & Watts (1988)，51页。数据是：

> x <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.56, 0.56,
1.10, 1.10)
> y <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200)

被拟合的模型是：

> fn <- function(p) sum((y - (p[1] * x)/(p[2] + x))^2)

为了进行拟合，我们需要估计参数初始值。一种寻找合理初始值的办法是把数据图形化，然后估计一些参数值，并且利用这些值初步添加模型曲线。

> plot(x, y)
> xfit <- seq(.02, 1.1, .05)
> yfit <- 200 * xfit/(0.1 + xfit)
> lines(spline(xfit, yfit))

当然，我们可以做的更好，但是初始值200 和0.1 应该足够了。现在做拟合：

> out <- nlm(fn, p = c(200, 0.1), hessian = TRUE)
> out
$minimum
[1] 1195.449

$estimate
[1] 212.68384222   0.06412146

$gradient
[1] -0.000153498  0.093420975

$hessian
            [,1]        [,2]
[1,]    11.94725   -7661.319
[2,] -7661.31875 8039421.153

$code
[1] 3

$iterations
[1] 26

模型拟合之后，out$minimum就是误差的平方和(SSE)，out$estimate 就是参数的最小二乘估计值。为了得到参数估计过程中近似的标准误(SE)，我们可以：

> sqrt(diag(2*out$minimum/(length(y) - 2) * solve(out$hessian)))
[1] 7.173465199 0.008744815

上述命令中的2表示参数的个数。一个95% 的置信区间可以通过±1.96 SE 计算得到。我们可以把这个最小二乘拟合曲线画在一个新的图上：

> plot(x, y)
> xfit <- seq(.02, 1.1, .05)
> yfit <- 212.68384222 * xfit/(0.06412146 + xfit)
> lines(spline(xfit, yfit))

标准包stats 提供了许多用最小二乘法拟合非线性模型的扩展工具。我们刚刚拟合过的模型是Michaelis-Menten模型，因此可以利用下面的代码得到类似的结果。

> df <- data.frame(x=x, y=y)
> fit <- nls(y ~ SSmicmen(x, Vm, K), df)
> fit
Nonlinear regression model
  model: y ~ SSmicmen(x, Vm, K)
   data: df
       Vm         K 
212.68371   0.06412 
 residual sum-of-squares: 1195

Number of iterations to convergence: 0 
Achieved convergence tolerance: 1.937e-06

> summary(fit)
Formula: y ~ SSmicmen(x, Vm, K)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
Vm 2.127e+02  6.947e+00  30.615 3.24e-11 ***
K  6.412e-02  8.281e-03   7.743 1.57e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.93 on 10 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 1.937e-06

11.7.2 最大似然法

最大似然法也是一种非线性拟合方法。它甚至可以用在误差非正态的数据中。这种方法估计的参数将会使得对数似然值最大或者负的对数似然值最小。下面的例子来自Dobson (1990), pp. : 108–111。这个例子对剂量-响应数据拟合logistic模型（当然也可以用glm() 拟合）。数据是：

> x <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113,
1.8369, 1.8610, 1.8839)
> y <- c( 6, 13, 18, 28, 52, 53, 61, 60)
> n <- c(59, 60, 62, 56, 63, 59, 62, 60)

要使负对数似然值最小，则编写函数如下：

> fn <- function(p)
sum( - (y*(p[1]+p[2]*x) - n*log(1+exp(p[1]+p[2]*x))
+ log(choose(n, y)) ))

我们选择一个适当的初始值，开始拟合：

> out <- nlm(fn, p = c(-50,20), hessian = TRUE)
> out
$minimum
[1] 18.71513

$estimate
[1] -60.71727  34.27021

$gradient
[1] 1.456959e-08 2.643526e-08

$hessian
          [,1]     [,2]
[1,]  58.48405 103.9787
[2,] 103.97873 184.9662

$code
[1] 1

$iterations
[1] 21

拟合后，out$minimum 就是负对数似然值，out$estimate 就是最大似然拟合的参数值。为了得到拟合过程近似的标准误差，我们可以：

> sqrt(diag(solve(out$hessian)))
[1] 5.554567 3.123365

参数估计的95% 信度期间可由估计值 ± 1.96 SE 计算得到。

11.8 一些非标准模型

在本章的最后，我们简单介绍一下R 里面某些用于特殊回归和数据分析问题的工具。

• 混合模型。用户捐献包nlme里面提供了函数lme() 和nlme()。这些函数可以用于混合效应模型的线性和非线性回归。也就是说在线性和非线性回归中，一些系数受随机因素的影响。这些函数需要充分利用公式来描述模型。
• 局部近似回归。函数loess() 利用局部加权回归进行一个非参数回归。这种回归对显示一组凌乱数据的趋势和描述大数据集的整体情况非常有用。函数loess 和投影跟踪回归的代码都放在标准包stats 中。
• 稳健回归。有多个函数可以用于拟合回归模型，同时尽量不受数据中极端值的影响。在推荐包MASS 中的函数lqs() 为高稳健性的拟合提供了最新的算法。另外，稳健性较低但统计学上高效的方法同样可以在包MASS中得到，如函数rlm。
• 累加模型。这种技术期望可以通过决定变量的平滑累加函数构建回归函数。一般来说，每个决定变量都有一个平滑累加函数。第三方包acepack 里面的函数avas() 和ace() 以及包mda 里面的函数bruto() 和mars() 为这种技术提供了一些例子。这种技术的一个扩展是第三方包gam 和mgcv 里面实现的广义累加模型。
• 树型模型。除了利用外在的全局线性模型预测和解释数据，还可以利用树型模型递归地在决定性变量的判断点上将数据的分叉分开。这样做会把数据最终分成几个不同组，使得组内尽可能相似而组间尽可能差异。这样常常会得到一些其他数据分析方法不能产生的的信息。模型可以用一般的线性模型形式指定。该模型拟合函数是tree()， 而且许多泛型函数，如plot() 和text() 都可以很好的用于树型模型拟合结果的图形显示。R 里面的树型模型函数可以通过第三方包rpart 和tree 得到。