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

特定形式的非线性模型可以通过广义线性模型 (glm()) 拟合。但是许多时候,我们必须把非线性拟合的问题 作为一个非线性优化的问题解决。 R的非线性优化程序是 optim()nlm()二者分别替换 S-Plusms()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 计算得到。

本文采用「CC BY-SA 4.0 CN」协议转载自互联网、仅供学习交流,内容版权归原作者所有,如涉作品、版权和其他问题请给「我们」留言处理。

返回顶部