R语言中级篇-PPT课件.ppt

Section II 基基础统计础统计与数学模型与数学模型数据的数据的产生生 dat1-rnorm(20)#一一组数的平均数的平均值dat1dat20.05时候，我候，我们可以可以认为这组数据是正数据是正态的的，但是适用，但是适用条件是所条件是所测的数据个数在的数据个数在 3 and 5000 之之间格式格式 test.nor-dat$biomass shapiro.test(test.nor)格式格式binom.test(x,n,p=0.5,alternative=c(“two.sided”,“less”,“greater”),conf.level=0.95)#数据是否是二数据是否是二项分布的分布的binom.test(c(682,243),p=3/4)binom.test(682,682+243,p=3/4)#The same 1 基基础统计求一求一组数据的数据的95%置信区置信区间 xbar-mean(x)#一一组数的平均数的平均值sigma-sd(x)#标准差准差sem-sigma/sqrt(length(x)#平均平均标准准误 low-xbar+sem*qnorm(0.025)#置信区置信区间下限下限 high-xbar+sem*qnorm(0.975)#置信区置信区间上限上限1 基基础统计随机抽随机抽样（random sampling）sample(x,size,replace=FALSE)a1 T H T H T H H H H T1 基基础统计概率概率统计（Probability）a-1:40asample(a,6)40个数字中，抽取想要的个数字中，抽取想要的6个数字，由于是非放回抽个数字，由于是非放回抽样，每个抽，每个抽取数字的概率分取数字的概率分别为:1/40,1/39,1/38,1/37,1/36,1/351基基础统计整体的概率：Pro-1/prod(40:35)如果是仅仅选取5个数字，那么就有不同组合，总体的概率为prod(5:1)/prod(40:35)或者是：1/choose(40,5)数据分布数据分布-密度密度图（Densities）x-seq(-4,4,0.1)plot(x,dnorm(x),type=l)1基基础统计数据分布数据分布-密度密度图（Densities）x-0:50plot(x,dbinom(x,size=50,prob=.33),type=h)#prob解释#probability of success on each trial 1基基础统计基基础统计练习 x-rnorm(20)#产生一生一组随机随机树Exercise:#编辑一个能一个能够返回所有基本返回所有基本统计量的函数量的函数var(x)mean(x)sd(x)1 基基础统计R语言中求言中求导#R语语言求言求导导y=D(expression(sin(x)+x3),x)x=2;eval(y);#就可以就可以计计算出算出f(2);值值y=deriv(expression(sin(x)+x3),x,function.arg=T)y(2)1 基基础统计Apply函数的函数的应用用 函数格式：apply(X,MARGIN,FUN,.)#MARGIN可以选择，1为按行，2为按列dat-matrix(1:6,2,3)#创建一个23的矩阵apply(dat,1,mean)apply(dat,2,mean)apply(dat,1,sd)Ex:a-tapply(dat$biomass,dat$disturb,mean)1基基础统计2、一般、一般线性模型（性模型（General Linear models）2.1 一般一般线性模型（性模型（General Linear models）一般一般线线性回性回归归模型模型 general linear models#Note1:适用的假适用的假设前提是前提是1)变量的独立性，非独立的量的独立性，非独立的变量会有很多意量会有很多意想不到的想不到的结果出果出现；2）正）正态的残差分布；的残差分布；3）方差）方差齐性；性；4）小的采）小的采样误差差#Note2；线性模型性模型 lm:自自变量是量是连续变量（即量（即x是是1,2,3,4,4.5这种）；种）；ANOVA:变量是因子型的（比如加温、森林干量是因子型的（比如加温、森林干扰类型等）；型等）；ANCOVA:既既有因子有因子变量又有量又有连续变量量2.1 一般一般线性模型（性模型（General Linear models）线线性回性回归归模型模型 lm()lm(formula,data,subset,weights,na.action,method=qr,model=TRUE,x=FALSE,y=FALSE,qr=TRUE,singular.ok=TRUE,contrasts=NULL,offset,.)Ex1:dat-read.csv(data_eg1.csv,head=T)head(dat)lm.model|t|)(Intercept)1.834e+02 6.954e-01 263.78 2e-16*num.sp 6.356e-02 5.011e-03 12.68 2e-16*-Signif.codes:0*0.001*0.01*0.05.0.1 1 Residual standard error:10.19 on 598 degrees of freedomMultiple R-squared:0.212,Adjusted R-squared:0.2107 F-statistic:160.9 on 1 and 598 DF,p-value:F)num.sp 1 16690 16689.5 160.87 2.2e-16*Residuals 598 62039 103.7 -Signif.codes:0*0.001*0.01*0.05.0.1 2.1 一般一般线性模型（性模型（General Linear models）对于要于要执行很多个模型，而且你行很多个模型，而且你还期望能期望能够自自动提取提取这些模型的参数，特些模型的参数，特别是在大量的随机模是在大量的随机模拟中，那么中，那么这些参数是以些参数是以list的形式存的形式存储在在summary(lm.model)中的，可以中的，可以调用和存用和存储的，如：的，如：res.summary-summary(lm.model)mode(res.summary)a-res.summary1 a a|t|)(Intercept)183.42426737 0.695375275 263.77738 0.000000e+00num.sp 0.06356016 0.005011243 12.68351 8.146262e-33 r2-res.summary8#调取R2#或者是 ls(res.summary)r2-res.summary$r.squared r22.1 一般一般线性模型（性模型（General Linear models）一个模型的一个模型的结结果里面存果里面存储储了不止了不止这这些些显显示的参数，示的参数，还还有很多的有很多的结结果可以果可以调调用：用：其他的一些其他的一些专门专门的看的函数，比如的看的函数，比如 coef,effects,residuals,fitted模型模型预测预测 predict()模型模型预测函数函数predict,主要主要传递的参数名字和模型一致，比如我的参数名字和模型一致，比如我们有模型有模型 lm.model和一列新的物种数和一列新的物种数new.sp，要，要预测基于基于这个物种数的生物量个物种数的生物量 new.sp-c(90,45,26,79,357)new.sp new.bio-predict(lm.model,data.frame(num.sp=new.sp)new.bio 2.1 一般一般线性模型（性模型（General Linear models）一个模型的构建完一个模型的构建完毕毕后，我后，我们们用什么方法来用什么方法来评评价模型的好坏呢？价模型的好坏呢？有很多的有很多的评评判判标标准，其中比如高的准，其中比如高的拟拟合合R平方，相平方，相对较少的因子，低的少的因子，低的AIC值（相（相对于同一的于同一的应变量量 Y值）fitvalue-fitted(lm.model)res-resid(lm.model)#获取残差。（残差即：模型的取残差。（残差即：模型的预测值与与实际值之之间的差的差值，即余下模型没，即余下模型没有解有解释的部分）的部分）plot(fitvalue,res)#可看下残差的和可看下残差的和预测值关系，来关系，来评判其模型判其模型拟合效果，合效果，较好的残差分布好的残差分布 abline(h=0,col=red)#abline(h=0)，好的模型，残差的散点，好的模型，残差的散点图是符合正是符合正态分布的分布的#标标准化的残差准化的残差图图 fit2-rstandard(lm.model)plot(fit2,res)或者可以看自或者可以看自变量和量和x的残差的残差图2.1 一般一般线性模型（性模型（General Linear models）#模型的残差正模型的残差正态检测 shapiro.test(lm.model$resid)qqnorm(lm.model$resid)#如果一个模型不是很好，那么我如果一个模型不是很好，那么我们就可能会用就可能会用glm模型来模型来处理理 AIC(lm.model)#但是注意，但是注意，AIC是只是只针对一一组相同的相同的y值才有比才有比较意意义。AIC：Akaike Information Criterion(1973)，后面会，后面会详细讲 1 4491.881#如果模型不如果模型不够好，那么我好，那么我们就准就准备对其其进行行转换，如，如log,sqrt,2.1 一般一般线性模型（性模型（General Linear models）较好的模型残差分布2.1 一般一般线性模型（性模型（General Linear models）残差残差标标准残差准残差2.2 方差分析方差分析Analysis of Variance（ANOVA）方差分析方差分析Analysis of Variance（ANOVA）单单因素方差分析因素方差分析dat$Increase_T-as.factor(dat$Increase_T)#注意，我们做的是水平处理的话，那么应该隐去相关的数值信息，做因子转换dis.aovF)Increase_T 2 507.07 253.54 22.903 2.728e-06*Residuals 24 265.68 11.07 -Signif.codes:0*0.001*0.01*0.05.0.1 1 2.2 方差分析方差分析Analysis of Variance（ANOVA）假假设前提是前提是这几种加温几种加温处理理对群落的生物量作用无明群落的生物量作用无明显差异，而差异，而结果果显示示这种假种假设的概率低于的概率低于p0.001,说明我明我们不可接受不可接受这个假个假设，即，即结果有明果有明显差异。差异。2.2 方差分析方差分析Analysis of Variance（ANOVA）箱箱线图线图作作图图-顺便便讲下作下作图出出图语句形式tiff()par()plot()dev.off 2.2 方差分析方差分析Analysis of Variance（ANOVA）箱箱线图线图作作图图-顺便便讲下作下作图出出图1)确定确定输出的出的图的文件名，大小等信息的文件名，大小等信息tiff(file=fig.tif,width=72,height=72,units=mm,res=1200,compression=lzw,pointsize=4)2）确定画）确定画图到到边界的距离界的距离par(mar=c(5,6,2,2)#分别表示边界的下、左、上、右的距离，防止输出的字不完整，可以测试下，如果不用这个语句的结果是什么3）作）作图的主体函数的主体函数boxplot(dat$biomassdat$Increase_T,xlab=expression(paste(Increased temperature (,degree,C),ylab=expression(paste(Biomass g/,m2),lwd=0.5,col=grey,cex.axis=2,cex.lab=2,cex=2)dev.off()#输出出结束束语 2.2 方差分析方差分析Analysis of Variance（ANOVA）均均值值的多重比的多重比较较 我们刚才得到的结果表示的是在目前的3个水平上具有差异，但并不是指他们两两具有差异，我们有时候想要的是他们两两间是否有差异#Note1；其原理用的是每组数据进行t检验那样，但是当因素水平角多的时候，检验同时进行，容易造成误差，就需要对显著性进行调整，即P值进行调整#Note2:在R中有自动调整P值的函数，在pairwise.t.test()内可以设置pairwise.t.test(x,g,p.adjust.method=p.adjust.methods,pool.sd=!paired,paired=FALSE,alternative=c(two.sided,less,greater),.)#x为输入的数值变量，g为分组变量p.adjust.methods用缺省设置值holm“-Holm(1979)2.2 方差分析方差分析Analysis of Variance（ANOVA）多重比多重比较较tem.t-pairwise.t.test(dat$biomass,dat$Increase_T)tem.t结果显示：Pairwise comparisons using t tests with pooled SD data:dat$biomass and dat$Increase_T 2 4 4 0.38 -6 5.4e-06 3.4e-05P value adjustment method:holm 2.2 方差分析方差分析Analysis of Variance（ANOVA）因为我们的假设前提是数据正态和方差齐性有时候我们可以做个检验shapiro.test()#对每一组数据是否符合正态分布bartlett.test(exp.gen,organ.nam)2.2 方差分析方差分析Analysis of Variance（ANOVA）练习EX1：目前有某种基因在各个器官中的表达水平，求各个器官间该基因的表达水平是否有差异organ.nam-factor(rep(c(root,stem,leaf,flower),c(4,6,8,6)exp.gen-c(0.5,0.46,0.23,0.29,0.2,0.11,0.14,0.14,0.13,0.19,0.19,0.19,0.34,0.43,0.57,0.75,0.78,0.71,0.88,0.39,0.45,0.28,0.53,0.88)pairwise.t.test(exp.gen,organ.nam)2.2 方差分析方差分析Analysis of Variance（ANOVA）双因素方差分析双因素方差分析 比如全球比如全球变变化控制化控制实验实验研究中的加温，添加研究中的加温，添加N，加水分析，往往需要多个因素考，加水分析，往往需要多个因素考虑，并且存在交互作用，并且存在交互作用1）不考虑交互作用dat$Increase_T-as.factor(dat$Increase_T)#因子因子转换dat$Increase_T-as.factor(dat$add_N)dis.aov-aov(biomassIncrease_T+add_N,data=dat)#不考不考虑交互作用交互作用summary(dis.aov)加温有显著影响加N没有显著影响 2.2 方差分析方差分析Analysis of Variance（ANOVA）2）考）考虑交互作用交互作用加温有显著影响加N有显著影响dis.aov.inter-aov(biomassIncrease_T+add_N+Increase_T:add_N,data=dat)#考虑交互作用summary(dis.aov.inter)并且交互存在交互作用Note,由于交互作用的存在，影响到了N添加的作用解释，相比前面的结果 2.2 方差分析方差分析Analysis of Variance（ANOVA）Ancova分析分析既有非既有非连续变量又有量又有连续变量量*号Eg:见原始数据eg1.csv中，我们收集了全国600个样点的生物量及气候因子及干扰方式，而干扰行为在空间上可能与气候又存在一定的关系执行语句dat-read.csv(data_eg1.csv,head=T)dat$disturb-as.factor(dat$disturb)dis.aov-lm(biomasstemperature*disturb,data=dat)summary(dis.aov)#结果显示：2.2 方差分析方差分析Analysis of Variance（ANOVA）Ancova分析分析 2.2 方差分析方差分析Analysis of Variance（ANOVA）Ancova分析分析解读：干扰不显著，温度和干扰的交互作用对biomass有影响，降水非常显著 2.2 方差分析方差分析Analysis of Variance（ANOVA）#我我们可以看到干可以看到干扰影响不明影响不明显，那么我，那么我们期望看看去掉期望看看去掉disturb的因子，的因子，模型会如何反模型会如何反应，是否可以去掉？，是否可以去掉？dis.aov1-lm(biomasstemperature*disturb,data=dat)dis.aov2-update(dis.aov1,.-disturb-disturb:temperature)#update()的作用是去掉一个因子和dis.aov2-lm(biomasstemperature,data=dat)等效AIC(dis.aov1,dis.aov2)2.3 多元多元线性回性回归、模型、模型筛选、方差分解、方差分解模型模型筛选筛选AIC准准则 step(yx1+x2+x3+x4)AIC准准则是在因子个数与模型解是在因子个数与模型解释率之率之间找一个平衡的参找一个平衡的参数，来确定和数，来确定和选择一个最一个最优模型。去掉具有和其他因子共模型。去掉具有和其他因子共线性性部分解部分解释的因子。的因子。2.3 多元多元线性回性回归、模型、模型筛选、方差分解、方差分解多元多元线线性回性回归归 max.model-lm(num.sp plot_coverage+temperature+precipitation+elevation,data=dat)step(max.model)summary(step(max.model)模型模型筛选筛选函数函数step()step(object,scope,scale=0,direction=c(both,backward,forward),trace=1,keep=NULL,steps=1000,k=2,.)缺省值当k=2是基于AIC选择；steps,最大接受筛选的步骤，所以对于因子庞大的回归模型，这个就要小心，其因子组合数为2n 2.3 多元多元线性回性回归、模型、模型筛选、方差分解、方差分解 2.3 多元多元线性回性回归、模型、模型筛选、方差分解、方差分解#对于一些回于一些回归系数系数检验结果不果不显著的一些因子，我著的一些因子，我们可能考可能考虑进行行删除操作，在除操作，在R中，提供了两个函数中，提供了两个函数drop1()和和add1()#利用利用drop1函数看各个因子舍取的指函数看各个因子舍取的指标标，对对模型因子的修正可以用上面模型因子的修正可以用上面提到的提到的update()函数，其格式函数，其格式为为update(model1,.+another)max.model-lm(num.sp plot_coverage+temperature+precipitation+elevation,data=dat)drop.mod-drop1(max.model,test=“Chisq”)#基于卡方基于卡方检验的的结果果结结果果显显示：示：2.3 多元多元线性回性回归、模型、模型筛选、方差分解、方差分解虽然我们看起来，这个因子去掉AIC不会变化很多，但是基于drop1两个模型检验的结果，我们可以发现，plot_coverage还不能去掉 2.3 多元多元线性回性回归、模型、模型筛选、方差分解、方差分解方差分解方差分解mod$part$indfract 3 广广义线性模型性模型 Generalized linear model 3.1 广广义线性模型性模型 Generalized linear model#glm函数的格式函数的格式glm(formula,family=gaussian,data,weights,subset,na.action,start=NULL,etastart,mustart,offset,control=list(.),model=TRUE,method=glm.fit,x=FALSE,y=TRUE,contrasts=NULL,.)相关定相关定义广广义线性模型由性模型由Nelder和和Wedderburn在在1972年的年的时候提出。候提出。其主体思想是在不改其主体思想是在不改变应变量量y和自和自变量量x之之间的的线性关系，通性关系，通过一个一个连接函接函数数linkfunction 将自将自变量量线性性连接到接到应变量量y的关系中。的关系中。另外，在很多研究另外，在很多研究中中还还是要看模型的解是要看模型的解释释能力作能力作为为一个重要的指一个重要的指标标，而，而选择选择哪种分布族函数哪种分布族函数比如：比如：y=ax+b 变换为 y=log(ax+b)关键的参数 广广义线义线性模型，是在生性模型，是在生态态学学应应用中及其用中及其重要的模型，也是重点重要的模型，也是重点讲讲解的解的模型模型。广广义线义线性模型有很多好性模型有很多好处处，适用于，适用于连续连续和非和非连续连续数据，并且数据，并且应变应变量量y可可以是任何形式的分布，比如常以是任何形式的分布，比如常见的有高斯分布、二的有高斯分布、二项式分布、泊松分式分布、泊松分布等等。主要可以解决数据特布等等。主要可以解决数据特别离散（非正离散（非正态）数据的模型构建。）数据的模型构建。但是注意一点：但是注意一点：GLM模型的模型的获得的系数和通得的系数和通过summary()获获得的得的结结果果并不是原始的响并不是原始的响应变应变量，所量，所获获得的参数，而是得的参数，而是转换转换后的模型所后的模型所获获得的得的参数，因此，需要通参数，因此，需要通过过你所用的你所用的连连接函数的反函数去接函数的反函数去获获得得对应对应的原始的原始Y值的参数，或者使用的参数，或者使用predict()里面的参数里面的参数设设置成置成为为type=“respones”3.1 广广义线性模型性模型 Generalized linear model广广义线义线性模型（性模型（Generalized linear models）1)分布族与连接函数分布族（family）连接函数binomiallogit,probit,cogloggaussianidentityGamamaIdentity,inverse,logInverse.gaussian1/mu2poissonIdentity,log,sqrtquasilogit,probit,cloglog,identity,inverse,log,1/mu2,sqrt广广义线义线性模型（性模型（Generalized linear models）2)主要分布族的适用的条件分布族链接函数反函数适用条件高斯Identity1正态分布泊松logexp计数，很多的0或者变量是整数二项式logit1(1+1/exp(x)比例或者有无（0,1）数据分布inverse1/x连续变量不全一定用缺省设置的函数来作为连接函数，比如对0,1的二元判断数据就可以使用其提供的family=binomial(link=cloglog)广广义线义线性模型（性模型（Generalized linear models）1)正态分布（高斯分布gaussian）EX1：我们看物种的盖度、年均温、年降水联合是否可以较好解释物种多样性mod.glm-glm(num.sp plot_coverage+temperature+precipitation,family=gaussian(link=identity),data=dat)mod.lm-lm(num.sp plot_coverage+temperature+precipitation,data=dat)summary(mod.glm)#仍然可以用summary和anova函数这个来看结果，但是主要里面的参数并非是原始的参数，而需要转换，对于identity无所谓summary(mod.lm)广广义线义线性模型（性模型（Generalized linear models）正态分布结果解读其他和线性模型都类似，但是就这里有区别，这就是我们所生态最想要的，模型解释能力，但没直接提供广广义线义线性模型（性模型（Generalized linear models）正态分布结果解读由于GLM提供了零模型的总体方差和模型解释的剩余方差，那么我们就可以计算模型的，后面的分布类型一致语句：mod.res-summary(mod.glm)mod.power-100*(mod.res$null.deviance-mod.res$deviance)/mod.res$null.deviancemod.power广广义线义线性模型（性模型（Generalized linear models）1)泊松分布（Poisson）泊松分布的假设是平均值和方差相等log转换poisson分布的使用条件，y值是计数变量，即整数。但很多计数变量是很难转变为正态的 进行泊松分布可以考虑过度离散现象。其检测方法同样是残差除以其自由度。若过度离散存在，则要将family参数设置为准泊松分布(quasipoisson)例子，样地的盖度是否能够通过GLM模型构建出能够很好解释生物多样性的模型，注意假设生物量是个计数变量（整数）。广广义线义线性模型（性模型（Generalized linear models）广广义线义线性模型（性模型（Generalized linear models）#执行语句 mod.pos-glm(num.sp plot_coverage,family=poisson(link=log),data=dat)summary(mod.pos)plot(dat$plot_coverage,dat$num.sp)#进行画图，讲解下画图 x-seq(0,500,len=500)#先生成一系列密集的x,y值 pre-predict(mod.pos,data.frame(plot_coverage=x)#x里面的参数 p-exp(pre)#一定要进行转换 lines(x,p,col=red,lwd=2)#这个例子看起来还不是很完美，我们可以看一个随机的例子广广义线义线性模型（性模型（Generalized linear models）训练一个随机例子，并操作卡方模型检验#分步执行以下语句 x-rnorm(200)y-rpois(200,exp(1+x)#生成泊松分布的随机数 mod.pos2-glm(yx,family=poisson(link=log)plot(x,y)x.p-seq(-2,3,len=200)pre-predict(mod.pos2,data.frame(x=x.p)#x里面的参数 p-exp(pre)pre2-predict(mod.pos2,data.frame(x=x.p),type=response)#或者添加response head(p)head(pre2)lines(x.p,p)#模型的卡方检验 mod.pos33），我们一般认为两个模型（统一应变量），Delta AIC小于3那么，这两个模型的解释能力等价广广义线义线性模型（性模型（Generalized linear models）二项式分布（binomial）二项式分布中Logistic回归时最重要也是最典型的回归模型，特别是对应的响应变量是二分类的，比如，这是著名的伯努利试验，就是什么比如抛硬币的实验。对应生态学中，有些物种的出现与否，也与可能一些环境因子具有很重要的关系，这可能就需要我们利用逻辑斯蒂回归来解决这种问题广广义线义线性模型（性模型（Generalized linear models）模型参数获得后，通过反函数获取模型形式方程为：y=exp(34.678368-0.026798*x)/(1+exp(34.678368-0.026798*x)广广义线义线性模型（性模型（Generalized linear models）Logistic回归作图#自己分步练习下，我休息会x-seq(1000,1500,len=200)xpre-predict(mod.bin1,data.frame(elevation=x),type=response)preplot(dat$elevation,dat$species1)lines(x,pre,lwd=2,col=blue)广广义线义线性模型（性模型（Generalized linear models）Logistic回归，出现与不出现的概率，基于一定调查样本Eg2,在海拔梯度上观测一种高山分布的鸟类出现情况-逻辑斯蒂回归(用中文大家看起来方便点)dat-read.csv(eg6.csv,head=T)head(dat)#看下数据Y-cbind(dat$某鸟的出现,dat$鸟类观测次数-dat$某鸟的出现)#获取Y/N分类变量X-dat$观测海拔mod.bin2-glm(YX,family=binomial(link=logit)#作图x-seq(0,3500,len=500)pre1-predict(mod.bin2,data.frame(X=x)#主要大X和模型中的要一致pre2-exp(pre1)/(1+exp(pre1)广广义线义线性模型（性模型（Generalized linear models）plot(X,dat$出现比例,xlab=Altitude(m),ylab=Probability of occurrence,ylim=c(0,1),lwd=3,col=blue,cex.axis=1.5,cex.lab=1.5,cex=3)#注意前后循序lines(x,pre2,col=red,lty=2,lwd=2)#添加拟合线自己练习：#这种鸟在该区域海拔2750m出现的概率是？非非线线性回性回归归模型（模型（non-linear regression models）多项式回归函数形式lm(yx+I(x2)+I(x3)例子，由于受到中位效应或者面积效应的影响，往往物种在中海拔地区比较丰富。dat-read.csv(data_eg1.csv,head=T)y-dat$num.spx-dat$elevationplot(x,y,xlab=Altitude(m),ylab=Species richness)#看下散点图的分布lm.obj-lm(yx+I(x2)summary(lm.obj)#结果显示回归方程为：y=-2299+4.049x-1658x2非非线线性回性回归归模型（模型（non-linear regression models）作图x.pre-seq(min(x),max(x),len=1000)#进行分段y.pre-predict(lm.obj,data.frame(x=x.pre)#lm.obj函数形式，得到预测值,模型都可以预测plot(yx,lwd=0.2,cex=1.5,xlab=Altitude(m),ylab=Species richness)lines(x.pre,y.pre,col=red,lwd=3)也可以加根对称线y=-b/(2a)xa-sum.obj$coefficients3b-sum.obj$coefficients2sym.lin-(-b)/(2*a)abline(v=sym.lin,col=blue,lwd=3,lty=2)非非线线性回性回归归模型（模型（non-linear regression models）(Toms et al.Ecology,2003)非非线线性回性回归归模型（模型（non-linear regression models）分段回归分段回归有时候是获得一些生态边界、阈值、转折点的很好方法 (Rincon et al.Ecology,2003)非非线线性回性回归归模型（模型（non-linear regression models）非非线线性回性回归归模型（模型（non-linear regression models）函数形式piecewise.linear(x,y,middle=1,CI=FALSE,bootstrap.samples=1000,sig.level=0.05)(0,1),0代表断点只能发生在x的中点值范围内，1代表可以任何地方是否进行bootstrap操作获取不确定区间，对原始的点数据，可以自己操作下CI=TRUE非非线线性回性回归归模型（模型（non-linear regression models）非非线线性回性回归归模型（模型（non-linear regression models）#但是我们有时候想获得断点等相关参数，并且精细编辑相关的线段#获取想要的参数 cor.test(x,y)r1-x.yx.y,1=model1,cor.test(r1,1,r1,2)r2 model1,cor.test(r2,1,r2,2)#也可以用nls实现，后面讲非非线线性回性回归归模型（模型（non-linear regression models）#其他函数类似#有时候我们不知道大概用多少个多项式进行回归合适，那么我们利用局部多项式回归函数可以看看大致的格局#判断曲线样式可以用lowess函数plot(x,y)lines(lowess(yx),lty=1,lwd=2,col=red)lowess(x,y=NULL,f=2/3,iter=3,delta=0.01*diff(range(xy$xo)数据平滑参数非非线线性回性回归归模型（模型（non-linear regression models）非非线线性回性回归归模型（模型（non-linear regression models）非线性回归的一般函数nls#非线性回归模型-nls应用#化学计量dat-read.csv(eg7.csv,head=T)plot(dat$biomass,dat$element)#大致看下散点图，那么比较适合指数模型，但是我们要获取一个初始值#模型公式形式为 a*exp(b*x)#看形式b为负数,a即为初始值，当x=0的时候的值#构建的模型为：mod.nls-nls(elementa*exp(b*biomass),data=dat,start=list(a=12,b=0.1)summary(mod.nls)#我们可以获取参数，模型的结果为 y=13.182765*exp(-0.081681*x)非非线线性回性回归归模型（模型（nonlinear regression models）#作图X-seq(0,35,len=100)pre-predict(mod.nls,data.frame(biomass=X)line