星期六, 一月 07, 2012

用nls函数进行非线性回归

在许多实际问题中,回归模型中响应变量和预测变量之间的关系可能是复杂的非线性函数。有时候能通过变量变换的方法可以将其变为线性模型,有时则不能。在后一种情况下,就需要采取专门的非线性回归方法来建立模型。

非线性回归是在对变量的非线性关系有一定认识前提下,对非线性函数的参数进行最优化的过程,最优化后的参数会使得模型的RSS(残差平方和)达到最小。在R语言中最为常用的非线性回归建模函数是nls,下面以car包中的USPop数据集为例来讲解其用法。数据中population表示人口数,year表示年份。如果将二者绘制散点图可以发现它们之间的非线性关系。在建立非线性回归模型时需要事先确定两件事,一个是非线性函数形式,另一个是参数初始值。


一、模型拟合
对于人口模型可以采用Logistic增长函数形式,它考虑了初期的指数增长以及总资源的限制。其函数形式如下。

首先载入car包以便读取数据,然后使用nls函数进行建模,其中theta1、theta2、theta3表示三个待估计参数,start设置了参数初始值,设定trace为真以显示迭代过程。nls函数默认采用Gauss-Newton方法寻找极值,迭代过程中第一列为RSS值,后面三列是各参数估计值。然后用summary返回回归结果。
library(car)
pop.mod1 <- nls(population ~ theta1/(1+exp(-(theta2+theta3*year))),start=list(theta1 = 400, theta2 = -49, theta3 = 0.025), data=USPop, trace=T)
summary(pop.mod)
在上面的回归过程中我们直接指定参数初始值,另一种方法是采用搜索策略,首先确定参数取值范围,然后利用nls2包的暴力方法来得到最优参数。但这种方法相当费时。

还有一种更为简便的方法就是采用内置自启动模型(self-starting Models),此时我们只需要指定函数形式,而不需要指定参数初始值。本例的logistic函数所对应的selfstarting函数名为SSlogis
pop.mod2 <- nls(population ~ SSlogis(year,phi1,phi2,phi3),data=USPop)
二、判断拟合效果
非线性回归模型建立后需要判断拟合效果,因为有时候参数最优化过程会捕捉到局部极值点而非全局极值点。最直观的方法是在原始数据点上绘制拟合曲线。
library(ggplot2)
p <- ggplot(USPop,aes(year, population))
p+geom_point(size=3)+geom_line(aes(year,fitted(pop.mod1)),col='red')

若比较多个模型的拟合效果可使用AIC函数,取最小值为佳。

三、残差诊断
和线性回归类似,非线性回归假设误差是正态、独立和同方差性。为了检测这些假设是否成立我们用拟合模型的残差来代替误差进行判断。
plot(fitted(pop.mod1) , resid(pop.mod1),type='b')
同方差假设采用残差绝对值和拟合值的散点图判断,或是使用bartlett.test函数检测。正态性检测除了画QQ图外还可用shapiro.test函数。独立性检验则可以绘制滞后残差图或是使用acf函数。

其它可以使用的泛型函数还包括了anova、coef、confint、deviance、fitted、plot、predict、vcov。如果预测变量中包括了分类数据,这种情况下可利用nlme包中的nlsList函数进行非线性回归。nlstools包中也有许多对非线性回归有用的函数。

参考资料:Nonlinear Regression with R (use R)

3 条评论:

  1. 您好!请教个小问题,nls()拟合好,应该采取什么方法,将拟合后的theta保存到一个变量中?

    回复删除
    回复
    1. 估计系数已经存好了,你可以使用某个对象来保存,例如
      result = summary(modle)
      result$coef里头就是所有的系数。

      删除
  2. 您好!!!非常感谢您的文章。但我想问getInitial函数是不是可以用来得到估计的参数啊?具体应该怎么用呢?

    如果您肯抽时间回答,感激不尽!

    回复删除