R包mgcv的广义加模型函数gam

R包mgcv的广义加模型函数gam（）library(mgcv)set.seed(0)##simulatesomedata...[＃模拟一些数据...]dat-gamSim(1,n=400,dist=normal,scale=2)b-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)summary(b)plot(b,pages=1,residuals=TRUE)##showpartialresiduals[＃显示部分残差]plot(b,pages=1,seWithMean=TRUE)##`withintercept'CIs[＃拦截“证明书]##samefitintwoparts.....[＃同样适合两部分......]G-gam(y~s(x0)+s(x1)+s(x2)+s(x3),fit=FALSE,data=dat)b-gam(G=G)print(b)##changethesmoothnessselectionmethodtoREML[＃REML法改变平滑的选择方法]b0-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,method=REML)plot(b0,pages=1,scheme=1)##Wouldasmoothinteractionofx0andx1bebetter?[＃X0和X1的流畅交互更好？]##Usetensorproductsmoothofx0andx1,basis[＃使用张量积X0和X1顺利，基础]##dimension49(see?tefordetails,also?t2).[＃尺寸49（见详情TE，也T2）。]bt-gam(y~te(x0,x1,k=7)+s(x2)+s(x3),data=dat,method=REML)plot(bt,pages=1)plot(bt,pages=1,scheme=2)##alternativevisualization[＃替代的可视化]AIC(b0,bt)##interactionworsethanadditive[＃交互比添加剂差]##Ifitisbelievedthatx0andx1arenaturallyon[＃如果它被认为X0和X1自然]##thesamescale,andshouldbetreatedisotropically[＃相同的规模，应被视为各向同性]##thencouldtry...[＃然后可以尝试...]bs-gam(y~s(x0,x1,k=50)+s(x2)+s(x3),data=dat,method=REML)plot(bs,pages=1)AIC(b0,bt,bs)##additivestillbetter.[＃添加剂更好。]##Nowdoautomatictermsselectionaswell[＃现在做选择自动条款，以及]b1-gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat,method=REML,select=TRUE)plot(b1,pages=1)##setthesmoothingparameterforthefirstterm,estimaterest...[＃设置的第一任期内，估计截断的平滑参数...]bp-gam(y~s(x0)+s(x1)+s(x2)+s(x3),sp=c(0.01,-1,-1,-1),data=dat)plot(bp,pages=1,scheme=1)##alternatively...[＃或者...]bp-gam(y~s(x0,sp=.01)+s(x1)+s(x2)+s(x3),data=dat)#setlowerboundsonsmoothingparameters....[平滑参数设定的下限....]bp-gam(y~s(x0)+s(x1)+s(x2)+s(x3),min.sp=c(0.001,0.01,0,10),data=dat)print(b);print(bp)#samewithREML[同与REML法]bp-gam(y~s(x0)+s(x1)+s(x2)+s(x3),min.sp=c(0.1,0.1,0,10),data=dat,method=REML)print(b0);print(bp)##nowaGAMwith3dfregressionsplineterm&2penalizedterms[＃现在与的3DF回归样条任期和2的GAM处罚条款]b0-gam(y~s(x0,k=4,fx=TRUE,bs=tp)+s(x1,k=12)+s(x2,k=15),data=dat)plot(b0,pages=1)##nowsimulatepoissondata...[＃现在模拟泊松资料的...]dat-gamSim(1,n=4000,dist=poisson,scale=.1)##usecrbasistosavetime,with4000data...[＃使用“CR”的基础上，与4000数据保存时间，...]b2-gam(y~s(x0,bs=cr)+s(x1,bs=cr)+s(x2,bs=cr)+s(x3,bs=cr),family=poisson,data=dat,method=REML)plot(b2,pages=1)##dropx3,butinitializesp'sfrompreviousfit,to[＃下降X3，但初始化SP，从以前的契合，]##savemoretime...[＃节省更多的时间......]b2a-gam(y~s(x0,bs=cr)+s(x1,bs=cr)+s(x2,bs=cr),family=poisson,data=dat,method=REML,in.out=list(sp=b2$sp[1:3],scale=1))par(mfrow=c(2,2))plot(b2a)par(mfrow=c(1,1))##similarexampleusingperformanceiteration[＃类似的例子，使用性能迭代]dat-gamSim(1,n=400,dist=poisson,scale=.25)b3-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,data=dat,optimizer=perf)plot(b3,pages=1)##repeatusingGACVasinWood2008...[＃重复使用GACV作为2008年在木材...]b4-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,data=dat,method=GACV.Cp,scale=-1)plot(b4,pages=1)##repeatusingREMLasinWood2011...[＃重复使用REML法作为2011年在木材...]b5-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,data=dat,method=REML)plot(b5,pages=1)##abinaryexample(seelaterforlargedatasetversion)...[＃二进制的例子（见稍后为大型数据集版本）...]dat-gamSim(1,n=400,dist=binary,scale=.33)lr.fit-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=binomial,data=dat,method=REML)##plotmodelcomponentswithtruthoverlaidinred[覆盖在红色的图与真理的模型组件]op-par(mfrow=c(2,2))fn-c(f0,f1,f2,f3);xn-c(x0,x1,x2,x3)for(kin1:4){plot(lr.fit,residuals=TRUE,select=k)ff-dat[[fn[k]]];xx-dat[[xn[k]]]ind-sort.int(xx,index.return=TRUE)$ixlines(xx[ind],(ff-mean(ff))[ind]*.33,col=2)}par(op)anova(lr.fit)lr.fit1-gam(y~s(x0)+s(x1)+s(x2),family=binomial,data=dat,method=REML)lr.fit2-gam(y~s(x1)+s(x2),family=binomial,data=dat,method=REML)AIC(lr.fit,lr.fit1,lr.fit2)##AGammaexample,bymodify`gamSim'output...[＃伽玛例如，通过修改gamSim的输出...]dat-gamSim(1,n=400,dist=normal,scale=1)dat$f<-dat$f/4##truelinearpredictor[＃真正的线性预测]Ey<-exp(dat$f);scale<-.5##meanandGLMscaleparameter[＃的意思，GLM尺度参数]##Notethat`shape'and`scale'in`rgamma'arealmost[＃注意形和规模rgamma“几乎]##oppositeterminologytothatusedwithGLM/GAM...[＃对面的GLM/自由亚齐运动所使用的术语来...]dat$y-rgamma(Ey*0,shape=1/scale,scale=Ey*scale)bg-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=Gamma(link=log),data=dat,method=REML)plot(bg,pages=1,scheme=1)##ForinverseGaussian,see?rig[＃逆高斯，看到了什么？钻机]##nowa2Dsmoothingexample...[＃现在的二维平滑的例子...]eg-gamSim(2,n=500,scale=.1)attach(eg)op-par(mfrow=c(2,2),mar=c(4,4,1,1))contour(truth$x,truth$z,truth$f)##contourtruth[＃轮廓真相]b4<-gam(y~s(x,z),data=data)##fitmodel[＃适合模型]fit1-matrix(predict.gam(b4,pr,se=FALSE),40,40)contour(truth$x,truth$z,fit1)##contourfit[＃轮廓适合]persp(truth$x,truth$z,truth$f)##persptruth[＃persp真相]vis.gam(b4)##perspfit[＃persp适合]detach(eg)par(op)##################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]##largishdatasetexamplewithuserdefinedknots[＃稍大例如与用户定义的节集]##################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]par(mfrow=c(2,2))eg-gamSim(2,n=10000,scale=.5)attach(eg)ind-sample(1:10000,1000,replace=FALSE)b5-gam(y~s(x,z,k=50),data=data,knots=list(x=data$x[