用R语言进行数据分析：最小二乘法和最大似然法

特定形式的非线性模型可以通过广义线性模型 (glm()) 拟合。但是许多时候，我们必须把非线性拟合的问题 作为一个非线性优化的问题解决。 R的非线性优化程序是 optim() 和 nlm()。 二者分别替换 S-Plus 的 ms() 和 nlminb()。我们通过搜寻 参数值使得缺乏度（lack-of-fit）最低，如 nlm() 就是通过循环调试各种参数值得到最优值。 和线性回归不同，程序不一定会 收敛到一个稳定值。 nlm() 需要设定参数搜索的起始值， 而参数估计是否收敛在很大程度上依赖于 起始值设置的质量。

最小二乘法

拟合非线性模型的一种办法就是使误差平方和（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 是误差的平方和（SSE）， out$estimate 是参数的最小二乘估计值。 为了得到参数估计过程中近似的标准误(SE)，我们可以：

     > sqrt(diag(2*out$minimum/(length(y) - 2) * solve(out$hessian)))

上述命令中的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.68370711   0.06412123
      residual sum-of-squares:  1195.449
     > 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
     
     Residual standard error: 10.93 on 10 degrees of freedom
     
     Correlation of Parameter Estimates:
           Vm
     K 0.7651

最大似然法

最大似然法（Maximum likelihood）也是一种非线性拟合办法。它甚至可以用于 误差非正态的数据。这种方法估计的参数 将会使得对数似然值最大或者负的对数似然值 最小。下面的例子来自 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 就是负对数似然值， out$estimate 就是最大似然拟合的参数值。 为了得到拟合过程近似的标准误，我们可以：

     > sqrt(diag(solve(out$hessian)))

一个 95% 信度期间可根据公式 +/- 1.96 SE 计算得到。