R语言非线性回归入门

作者 Lionel Hertzog 　　　　　　　

在一簇散点中拟合一条回归线（即线性回归）是数据分析的基本方法之一。有时，线性模型能很好地拟合数据，但在某些（很多）情形下，变量间的关系未必是线性的。这时，一般有三类方法解决这个问题: (1) 通过变换数据使得其关系线性化, (2) 用多项式或者比较复杂的样条来拟合数据, (3) 用非线性函数来拟合数据

从标题你应该已经猜到非线性回归是本文的重点

什么是非线性回归

在非线性回归中，分析师通常采用一个确定的函数形式和相应的参数来拟合数据。最常用的参数估计方法是利用非线性最小二乘法（R中的nls函数）。该方法使用线性函数来逼近非线性函数，并且通过不断迭代这个过程来得到参数的最优解（本段来自维基百科）。非线性回归的良好性质之一是估计出的参数都有清晰的解释（如Michaelis-Menten模型的Vmax是指最大速率），而变换数据后得到的线性模型其参数往往难以解释。

非线性最小二乘拟合

首先，我们以Michaelis-Menten方程为例。

#生成一些仿真数据set.seed(20160227) x<-seq(0,50,1)y<-((runif(1,10,20)*x)/(runif(1,0,10)+x))+rnorm(51,0,1)#对于一些简单的模型，nls函数可以自动找到合适的参数初值m<-nls(y~a*x/(b+x))#计算模型的拟合优度cor(y,predict(m)) [1]0.9496598#将结果可视化plot(x,y) lines(x,predict(m),lty=2,col="red",lwd=3)123456789101112

输出的图片如下：

选择适宜的迭代初值

在非线性回归中，找到合适的迭代初值对于整个模型算法的收敛性而言至关重要。假如你设定的参数初值完全脱离了其潜在的取值范围，迭代算法可能不收敛或者返回一些没有意义的参数值。比如返回一个大小为1000的增长率，但其真值却是1.04。寻找合适初值的最好办法是“紧盯着”数据，绘制相应图表并结合你对方程的理解来确定参数的合适初值。

#生成仿真数据，并且此次对于参数没有先验信息y<-runif(1,5,15)*exp(-runif(1,0.01,0.05)*x)+rnorm(51,0,0.5)#可视化数据并选择一些参数初值plot(x,y)#通过这个散点图确定参数a,b的初值a_start<-8#参数a是x=0时y的取值b_start<-2*log(2)/a_start#b是衰减速率#拟合模型m<-nls(y~a*exp(-b*x),start=list(a=a_start,b=b_start))#计算拟合优度cor(y,predict(m)) [1]0.9811831#将结果可视化lines(x,predict(m),col="red",lty=2,lwd=3)1234567891011121314

输出的图片如下

使用自启动函数

不同的科学研究领域会对同一个模型设定不同的参数形式（即不同的方程），比如研究人口增长的逻辑斯蒂模型，在生态学中一般采用如下形式：

Nt=K∗N0∗ertK+N0∗(ert–1)

library(deSolve)#利用逻辑斯蒂模型生成人口增长的仿真数据，并用nls估计参数log_growth<-function(Time,State,Pars){ with(as.list(c(State,Pars)),{ dN<-R*N*(1-N/K)return(list(c(dN))) }) }#逻辑斯蒂增长的参数pars<-c(R=0.2,K=1000)#设定初值N_ini<-c(N=1)#常微分方程的时间阶段（下标t）times<-seq(0,50,by=1)#常微分方程out<-ode(N_ini,times,log_growth,pars)#添加一些随机波动N_obs<-out[,2]+rnorm(51,0,50)#个体数值不能小于1N_obs<-ifelse(N_obs<1,1,N_obs)#画图plot(times,N_obs)12345678910111213141516171819202122

这部分代码只是生成了带有随机误差的仿真数据，接下来的部分会展现估计参数初值的技巧。R语言中有一个估计逻辑斯蒂方程参数的内建函数（SSlogis），但它使用的是如下方程：

#寻找方程的参数SS<-getInitial(N_obs~SSlogis(times,alpha,xmid,scale),data=data.frame(N_obs=N_obs,times=times))12

我们可使用getInitial函数来对模型参数做一个基于数据的初步估计。然后把该函数的输出作为一个向量化参数传递给自启动函数（SSlogis），同时也将无引号的三个参数名赋值给逻辑斯蒂方程（译者注：即alpha，xmid，scale三个参数）。

然而，由于SSlogis的参数设定有些不同，我们需要对SSlogis的输出值做一些处理，使得其与逻辑斯蒂方程中的形式一致。

#改变参数形式K_start<-SS["alpha"] R_start<-1/SS["scale"] N0_start<-SS["alpha"]/(exp(SS["xmid"]/SS["scale"])+1)#构建模型的公式log_formula<-formula(N_obs~K*N0*exp(R*times)/(K+N0*(exp(R*times)-1)))#拟合模型m<-nls(log_formula,start=list(K=K_start,R=R_start,N0=N0_start))#估计参数summary(m) Formula:N_obs~K*N0*exp(R*times)/(K+N0*(exp(R*times)-1)) Parameters: EstimateStd.ErrortvaluePr(>|t|) K1.012e+033.446e+0129.366<2e-16*** R2.010e-011.504e-0213.360<2e-16*** N09.600e-014.582e-012.0950.0415* --- Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1  Residualstandarderror:49.01on48degreesoffreedom  Numberofiterationstoconvergence:1 Achievedconvergencetolerance:1.537e-06#计算拟合优度cor(N_obs,predict(m))  [1]0.9910316#结果可视化lines(times,predict(m),col="red",lty=2,lwd=3)12345678910111213141516171819202122232425262728293031

输出图形如下:

修改SSlogis输出的参数结构确实有些繁琐，不过值得一试。

在后续的文章中，除了非线性最小二乘法，我们将利用更为可靠和强大的极大似然估计法来拟合模型。它能使你构建你能想到的任何模型。

原文刊载于datascience+网站

链接:http://http://datascienceplus.com/first-steps-with-non-linear-regression-in-r/

PPV课转载的每一篇文章均来源于公开网络，仅供学习使用，不会用于任何商业用途，文章版权归原作者所有，如果侵犯到原作者的权益，请您与我们联系删除或者授权事宜，联系邮箱：149104196@qq.com。转载PPV课网站文章请注明原文章作者，否则产生的任何版权纠纷与PPV课无关。

PPV课-国内领先的大数据学习社区和职业培训平台

1、回复“数据分析师”查看数据分析师系列文章

2、回复“案例”查看大数据案例系列文章

3、回复“征信”查看相关征信的系列文章

4、回复“可视化”查看可视化专题系列文章

5、回复“SPPS”查看SPSS系列文章

6、回复“答案”查看hadoop面试题题目及答案

7、回复“爱情”查看大数据与爱情的故事

8、回复“笑话”查看大数据系列笑话

9、回复“大数据1、大数据2、大数据3、大数据4”查看大数据历史机遇连载

原文始发于微信公众号（PPV课数据科学社区）：R语言非线性回归入门