Plotting gam model results in lattice

I am fitting a model using gam from the mgcv package. I am storing the result in model and so far I have been looking at the smooth components using plot(model). I have recently started using lattice and like its output. So I am wondering, is it possible to plot these graphs using lattice?

This my data set: 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)

How can I do plot this model in lattice package with three panels representing my three sites smoother,please?

You also might noticed that I have used the dayinyear instead of proper month format(the first column in the data). This is because Generalized additive models do not handle categorical variables. However, I would like to represent the time in my graph with the names of months (like in first column), Does any one know the way forward for that in a lattice plot?

Solution

Here is a general way to do it using some fake data. You will need to tweak this to make sure the names are as you like,

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 is the predicted value for each predictor and X2is where the grouping variable is (indicates which data belongs to each predictor). If you are working we these a lot you should rename things to be clearer. The order part just avoids spaghetti plots

You can control the way the x-axis is labeled using the at and labels arguments for the x-axis in the scales argument. For details see ?xyplot

Update - Here is a version that works with this data

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")
     }
)

Note that I changed the name of the starting data.frame from data to Data

EDIT - I have added the two dashed lines that show + /- 2 SE for each plot