扫码加入数据分析学习群

R语言怎么使用mgcv包做面板数据的非参数广义可加模型的估计？

 library(mgcv)
 set.seed(0)
 n<-400
 sig2<-4
 x0 <- runif(n, 0, 1)
 x1 <- runif(n, 0, 1)
 x2 <- runif(n, 0, 1)
 x3 <- runif(n, 0, 1)
 pi <- asin(1) * 2
 f <- 2 * sin(pi * x0)
 f <- f + exp(2 * x1) - 3.75887
 f <- f+0.2*x2^11*(10*(1-x2))^6+10*(10*x2)^3*(1-x2)^10-1.396
 e <-rnorm(n, 0, sqrt(abs(sig2)))
 y <- f + e
 b<-gam(y~s(x0)+s(x1)+s(x2)+s(x3))

###
 summary(b)
 plot(b,pages=1)
 # an extra ridge penalty (useful with convergence problems) ....
 bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),H=diag(0.5,41))
 print(b);print(bp);rm(bp)
 # set the smoothing parameter for the first term, estimate rest ...
 bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),sp=c(0.01,-1,-1,-1))
 plot(bp,pages=1);rm(bp)
 # set lower bounds on smoothing parameters ....
 bp<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),min.sp=c(0.001,0.01,0,10))
 print(b);print(bp);rm(bp)
 

 
 

###### now a GAM with 3df regression spline term & 2 penalized terms
 b0<-gam(y~s(x0,k=4,fx=TRUE,bs="tp")+s(x1,k=12)+s(x2,15))
 plot(b0,pages=1)
 # now fit a 2-d term to x0,x1
 b1<-gam(y~s(x0,x1)+s(x2)+s(x3))
 par(mfrow=c(2,2))
 plot(b1)
 par(mfrow=c(1,1))
 # now simulate poisson data
 g<-exp(f/5)
 y<-rpois(rep(1,n),g)
 b2<-gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson)
 plot(b2,pages=1)
 # and a pretty 2-d smoothing example....
 test1<-function(x,z,sx=0.3,sz=0.4)  
 { (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
   0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
 }
 n<-500
 old.par<-par(mfrow=c(2,2))
 x<-runif(n);z<-runif(n);
 xs<-seq(0,1,length=30);zs<-seq(0,1,length=30)
 pr<-data.frame(x=rep(xs,30),z=rep(zs,rep(30,30)))
 truth<-matrix(test1(pr
 x,prx,prz),30,30)
 contour(xs,zs,truth)
 y<-test1(x,z)+rnorm(n)*0.1
 b4<-gam(y~s(x,z))
 fit1<-matrix(predict.gam(b4,pr,se=FALSE),30,30)
 contour(xs,zs,fit1)
 persp(xs,zs,truth)
 vis.gam(b4)
 par(old.par)
 # very large dataset example using knots
 n<-10000
 x<-runif(n);z<-runif(n);
 y<-test1(x,z)+rnorm(n)
 ind<-sample(1:n,1000,replace=FALSE)
 b5<-gam(y~s(x,z,k=50),knots=list(x=x[ind],z=z[ind]))
 vis.gam(b5)
 # and a pure "knot based" spline of the same data
 b6<-gam(y~s(x,z,k=100),knots=list(x= rep((1:10-0.5)/10,10),
         z=rep((1:10-0.5)/10,rep(10,10))))
 vis.gam(b6,color="heat")