[R] plotting gamm results in lattice

Duncan Mackay dulcalma at bigpond.com
Mon Jun 12 17:26:27 CEST 2017


Hi Maria

If you have problems just start with a small model with predictions and then plot with xyplot
the same applies to xyplot
Try 

library(gamm4)

spring <- dget(file = "G:/1/example.txt")
 str(spring)
'data.frame':   11744 obs. of  11 variables:
 $ WATERBODY_ID          : Factor w/ 1994 levels "GB102021072830",..: 1 1 2 2 2 3 3 3 4 4 ...
 $ SITE_ID               : int  157166 157166 1636 1636 1636 1635 1635 1635 134261 1631 ...
 $ Year                  : int  2011 2014 2013 2006 2003 2002 2013 2005 2013 2006 ...
 $ Q95                   : num  100 100 80 80 80 98 98 98 105 105 ...
 $ LIFE.OE_spring        : num  1.02 1.03 1.02 1.06 1.06 1.07 1.12 1.05 1.14 1.05 ...
 $ super.end.group       : Factor w/ 6 levels "B","C","D","E",..: 1 1 3 3 3 2 2 2 4 4 ...
 $ X.urban.suburban      : num  0 0 0.07 0.07 0.07 0.53 0.53 0.53 8.07 8.07 ...
 $ X.broadleaved_woodland: num  2.83 2.83 10.39 10.39 10.39 ...
 $ X.CapWks              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Hms_Poaching          : int  0 0 10 10 10 0 0 0 0 0 ...
 $ Hms_Rsctned           : int  0 0 0 0 0 0 0 0 0 0 ...

model<-
gamm4(LIFE.OE_spring~s(Q95, by=super.end.group)+Year+Hms_Rsctned+Hms_Poaching+X.broadleaved_woodland +X.urban.suburban+X.CapWks, data=spring,
      random=~(1|WATERBODY_ID/SITE_ID))
Warning message:
Some predictor variables are on very different scales: consider rescaling

# plot(model$gam, page=1, font.lab=2, xlab="Residual Q95")

M <-predict(model$gam,type="response",se.fit=T)

# joining the dataset and the predictions keep it "in house" and on path
springP = cbind(spring, M)
springP$upper = with(springP,fit+1.96*se.fit)
springP$lower = with(springP,fit-1.96*se.fit)
str(springP)
'data.frame':   11744 obs. of  15 variables:
 $ WATERBODY_ID          : Factor w/ 1994 levels "GB102021072830",..: 1 1 2 2 2 3 3 3 4 4 ...
 $ SITE_ID               : int  157166 157166 1636 1636 1636 1635 1635 1635 134261 1631 ...
 $ Year                  : int  2011 2014 2013 2006 2003 2002 2013 2005 2013 2006 ...
 $ Q95                   : num  100 100 80 80 80 98 98 98 105 105 ...
 $ LIFE.OE_spring        : num  1.02 1.03 1.02 1.06 1.06 1.07 1.12 1.05 1.14 1.05 ...
 $ super.end.group       : Factor w/ 6 levels "B","C","D","E",..: 1 1 3 3 3 2 2 2 4 4 ...
 $ X.urban.suburban      : num  0 0 0.07 0.07 0.07 0.53 0.53 0.53 8.07 8.07 ...
 $ X.broadleaved_woodland: num  2.83 2.83 10.39 10.39 10.39 ...
 $ X.CapWks              : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Hms_Poaching          : int  0 0 10 10 10 0 0 0 0 0 ...
 $ Hms_Rsctned           : int  0 0 0 0 0 0 0 0 0 0 ...
 $ fit                   : num [1:11744(1d)] 1.03 1.04 1.04 1.02 1.01 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr  "1" "2" "3" "4" ...
 $ se.fit                : num [1:11744(1d)] 0.00263 0.00266 0.00408 0.00408 0.00411 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr  "1" "2" "3" "4" ...
 $ upper                 : num [1:11744(1d)] 1.03 1.04 1.05 1.03 1.02 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr  "1" "2" "3" "4" ...
 $ lower                 : num [1:11744(1d)] 1.02 1.03 1.03 1.01 1 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr  "1" "2" "3" "4" ...

# this produces a mess of lines for the upper and lower
xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = springP,
       lower = springP$lower,
       upper = springP$upper,
       subscripts = TRUE,
       panel = function(x,y, upper, lower, subscripts, ...){
                  panel.xyplot(x,y, type="smooth")
                  panel.xyplot(x[order(x)], upper[subscripts][order(x)], lty=2, col="red")
                  panel.xyplot(x[order(x)], lower[subscripts][order(x)], lty=2, col="red")
                  #panel.loess(x,y,...) #  have not tried to fix these lines- depends on what you want to do
                  #panel.rug(x = x[is.na(y)], y = y[is.na(x)])
                }
)

# smoothing them produces reasonable lines
xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = springP,
       lower = springP$lower,
       upper = springP$upper,
       subscripts = TRUE,
       panel = function(x,y, upper, lower, subscripts, ...){
                  panel.xyplot(x,y, type="smooth")
                  panel.xyplot(x[order(x)], upper[subscripts][order(x)], type = "smooth", lty=2, col="red")
                  panel.xyplot(x[order(x)], lower[subscripts][order(x)], type = "smooth", lty=2, col="red")
                  #panel.loess(x,y,...)
                  #panel.rug(x = x[is.na(y)], y = y[is.na(x)])
                }
)

# by newdata = ...
# take a cross-section of data - reduces the amount of data to be plotted
mx = aggregate(Q95 ~ super.end.group, spring, max, na.rm=T)
newlst <-
do.call(rbind,
        lapply(1:6, function(j) expand.grid(Q95 = seq(0,mx[j,2], length = 50), super.end.group = LETTERS[j])
                                            Year =
                                            Hms_Rsctned =
                                            ...))   # fill in yourself
newd <- data.frame(newlst)

and predict using newdata = newd

# starting with a smaller model
model2 <-
gamm4(LIFE.OE_spring~s(Q95, by=super.end.group), data=spring,
      random=~(1|WATERBODY_ID/SITE_ID))

newlst <-
do.call(rbind,
        lapply(1:6, function(j) expand.grid(Q95 = seq(0,mx[j,2], length = 50), super.end.group = LETTERS[j])) )
newd <- data.frame(Q95 = newlst)
# The following is untested
# I have not used gamm4 for a while; is  model$gam correct? Had problems here - have no time to go into the help pages
M3 <- predict(model2$gamm, newdata = newd, type="response",se.fit=T)

# keep together 
newdat = cbind(newd, M3)

# If you want CIs
newdat <- within(newdat, {
    lower = fit-1.96*se.fit
    upper = fit+1.96*se.fit
})

# plot- something like
xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = springP,
       Q2    = newddat$Q95,
       lower = newdat$lower,
       upper = newdat$upper,
       subscripts = TRUE,
       panel = function(x,y, Q2, upper, lower, subscripts, ...){
                  panel.xyplot(x,y, type="smooth")
                  panel.xyplot(Q2, upper, lty=2, col="red")
                  panel.xyplot(Q2, lower, lty=2, col="red")
                  #panel.loess(x,y,...)
                  #panel.rug(x = x[is.na(y)], y = y[is.na(x)])
                }
)

If you are putting a key onto the plot see
? xyplot and par.settings and key 

Regards

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au

-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Maria Lathouri via R-help
Sent: Monday, 12 June 2017 20:57
To: R-help Mailing List
Subject: [R] plotting gamm results in lattice

Dear all, 
I hope that you can help me on this. I have been struggling to figure this out but I haven't found any solution.
I am running a generalised mixed effect model, gamm4, for an ecology project. Below is the code for the model:
model<-gamm4(LIFE.OE_spring~s(Q95, by=super.end.group)+Year+Hms_Rsctned+Hms_Poaching+X.broadleaved_woodland             +X.urban.suburban+X.CapWks, data=spring, random=~(1|WATERBODY_ID/SITE_ID))
plot(model$gam, page=1, font.lab=2, xlab="Residual Q95")

I am trying to plot the results in lattice for publication purposes so I need to figure this out. I have been struggling but I think I have reached a dead end. 

Here is what I have been able to code:
M<-predict(model$gam,type="response",se.fit=T)
upr<- M$fit + (1.96 * M$se.fit)lwr<- M$fit - (1.96 * M$se.fit)
library(lattice)xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = spring, gm=model,       prepanel=function (x,y,...)list(ylim=c(min(upr),max(lwr))),       panel = function(x,y, gm, ...){              panel.xyplot(x,y, type="smooth")         panel.lines(upr,lty=2, col="red")         panel.lines(lwr,lty=2, col="red")         panel.loess(x,y,...)         panel.rug(x = x[is.na(y)],                   y = y[is.na(x)])       }       )
But, unfortunately, this is not what I get when I have the simple plot(model$gam). 

I have also attached a reproducible example in case you want to see for yourself. I hope that someone here has come up with a similar problem and can help me on this.
Thank you very much for your time.
Kind regards,Maria 



More information about the R-help mailing list