在格子中绘制 gam 模型

Question

我正在使用 mgcv 包中的 gam 拟合模型。我将结果存储在model 中，到目前为止，我一直在使用plot(model) 查看平滑组件。我最近开始使用lattice 并喜欢它的输出。所以我想知道，是否可以使用lattice 绘制这些图表？

这是我的数据集：https://gist.github.com/plxsas/fcef4a228c18c772b4f3

m2<- gam(TotalInd ~ s(dayinyear, by=as.numeric(Site=="1"), bs="cr")
  +s(dayinyear, by=as.numeric(Site=="2"), bs="cr") + s(dayinyear, 
   by=as.numeric(Site=="3"), bs="cr"), random=list(Replicate=~ 1), data=data)

我怎样才能在lattice 包中绘制这个模型，其中三个面板代表我的三个站点更流畅，拜托？

您可能还注意到我使用了dayinyear 而不是正确的月份格式（数据中的第一列）。这是因为广义加性模型不处理分类变量。但是，我想用月份的名称表示图表中的时间（如第一列），有人知道lattice 图中的前进方向吗？

Answer 1

这是使用一些虚假数据的一般方法。您需要对此进行调整以确保名称符合您的喜好，

library(reshape)
library(mgcv)
library(lattice)

X1<-rnorm(100)   # Make some fake data
X2<-rnorm(100)
X3<-rnorm(100)
Y<-rnorm(100)

Mod<-gam(Y~s(X1,bs="cr")+s(X2,bs="cr")+s(X3, bs="cr")) # make a model

Z<- predict(Mod,type="terms", se.fit=T)  #Z is the predicted value 
                               #for each smooth term, se.fit give you SE

Z2<-melt(Z$fit)                     #Z was in wide form, Z2 is long form
Z2$XX<-c(X1,X2,X3)            #add the original values for he predictors 
Z2$SE<-melt(Z$se.fit)$value  #add SE
Z2$UP<-Z2$value+2*Z2$SE      #+2 SE
Z2$Low<-Z2$value-2*Z2$SE     # - 2 SE
Z2<-Z2[order(Z2$XX),]

xyplot(value~XX|X2,data=Z2,type="l",col="black",as.table=T,
     prepanel=function (x,y,...)list(ylim=c(min(Z2$Low),max(Z2$UP))),
     panel=function(x,y,groups,subscripts,...){
       panel.xyplot(x,y,...)
       panel.lines(Z2$UP[subscripts]~Z2$XX[subscripts],lty=2, col="red")
       panel.lines(Z2$Low[subscripts]~Z2$XX[subscripts],lty=2, col="red")
     }
 ) 

value 是每个预测变量的预测值，X2是分组变量所在的位置（指示哪些数据属于每个预测变量）。如果您经常为我们工作，您应该重命名这些内容以使其更清晰。 order 部分只是避免了意大利面条情节

您可以使用 scales 参数中的 x 轴的 at 和 labels 参数来控制 x 轴的标记方式。详情见?xyplot

更新 - 这是一个适用于此数据的版本

m2<- gam(TotalInd ~ s(dayinyear, by=as.numeric(Site=="1"), bs="cr")
     +s(dayinyear, by=as.numeric(Site=="2"), bs="cr") 
     + s(dayinyear, by=as.numeric(Site=="3"), bs="cr"), 
     random=list(Replicate=~ 1), data=Data)


Z<- predict(m2,type="terms",se.fit=T) #Z is the predicted value and SE
Z2<-melt(Z$fit)                     #Z was in wide form, Z2 is long form

Z2$dayinyear<-Data$dayinyear        #add the original values for he predictors 
Z2$SE<-melt(Z$se.fit)$value
Z2$UP<-Z2$value+2*Z2$SE
Z2$Low<-Z2$value-2*Z2$SE

Z2<-Z2[Z2$value!=0,] #gets rid of excess zeroes

Z2<-Z2[order(Z2$dayinyear),]



xyplot(value~dayinyear|X2,data=Z2,type="l",col="black",as.table=T,
     prepanel=function (x,y,...)list(ylim=c(min(Z2$Low),max(Z2$UP))),
     panel=function(x,y,groups,subscripts,...){
       panel.xyplot(x,y,...)
       panel.lines(Z2$UP[subscripts]~Z2$dayinyear[subscripts],lty=2, col="red")
       panel.lines(Z2$Low[subscripts]~Z2$dayinyear[subscripts],lty=2, col="red")
     }
) 

请注意，我将起始 data.frame 的名称从 data 更改为 Data

编辑 - 我添加了两条虚线，显示每个图的 +/- 2 SE