Skip to content

gamm plots in lattice

3 messages · Maria Lathouri, Ulf Köther

Original
Maria Lathouri
# 
Dear all,?
I am new in list. I hope that you can help me on this.?
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))
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?
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: example.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20170608/ac10b235/attachment-0001.txt>
11 days later
Ulf Köther
# 
Dear Maria,

since it appears that no one has answered your question until now, I
will give you some hints how to proceed:

Caveat: I have not used lattice for a long time and therefore I will
give you an ggplot2-answer because I have no time to figure out the
details for lattice. But this is only the plotting side - at the end,
both graphic-systems should provide similar plots if you follow some
basic rules.

If you do not want to use ggplot2 but lattice anyway, maybe this post
will get you going:

https://www.r-bloggers.com/confidence-bands-with-lattice-and-r/

Good luck, Ulf

---------
R-Code:
---------

# Read your data:
dat <-  dget("D:/example.txt")
dat$SITE_ID <- factor(dat$SITE_ID)

library(gamm4)
library(ggplot2)

# You should include "super.end.group" also as a factor because
# your model has 6 smoothers, and each smoother is automatically centred
# around 0. The extra main term "super.end.group" allows for a vertical
# shift for the other 5 smoothers (Period 2).

m1 <- gamm4(LIFE.OE_spring ~ super.end.group + s(Q95, by = 	
            super.end.group) + Year + Hms_Rsctned + Hms_Poaching +
            X.broadleaved_woodland + X.urban.suburban + X.CapWks,
            data = dat, random = ~(1|WATERBODY_ID/SITE_ID))

# You want to reproduce this one, right?
plot(m1$gam, pages = 1)

# 1. You need new data to be predicted, not the old ones. Here
# Every variable in the model must be present. Which values you choose
# depends on what you want to present. Here I chose the first year and
# zero for everything else, but more often the mean of the variables is
# the smarter choice. The values of Q95 are chosen from min to max with
# 100 values in between for plotting:

newDat <- expand.grid(super.end.group = levels(dat$super.end.group),
                      Q95 = seq(from = min(dat$Q95, na.rm = TRUE),
                                to = max(dat$Q95, na.rm = TRUE),
                                length = 100),
                      Year = 2002,
                      Hms_Rsctned = 0,
                      Hms_Poaching = 0,
                      X.broadleaved_woodland = 0,
                      X.urban.suburban = 0,
                      X.CapWks = 0,
                      WATERBODY_ID = "GB102021072830",
                      SITE_ID = "157166")

# Then you predict with the new data:
datM <- predict(m1$gam, type = "response",
                 se.fit = TRUE, newdata = newDat)

# If you use a different family like "poisson" or any other than
# the gaussian, you need to use type = "link", and after
# calculating the lower and upper limits, you have to
# manually apply the inverse link function yourself on the fit and
# on the upper and lower limit. With a gaussian distribution, this is
# not necessary:
#
# datM2 <- predict(m1$gam, type = "link",
#                 se.fit = TRUE, newdata = newDat)
# all.equal(datM$fit, datM2$fit)

# Put the fit and the limits in the new data frame from which you
# predicted the response to get them in order with the variable
# "super.end.group":

newDat$fit <- datM$fit
newDat$upr <- datM$fit + (1.96 * datM$se.fit)
newDat$lwr <- datM$fit - (1.96 * datM$se.fit)

# Now some simple plotting, with the limits on the y-axis chosen to your
# data. Here you see that the smoothers are not centred around zero but
# on the point predicted by the model (smoother plus an individual
# intercept for each level of "super.end.group"):

ggplot(newDat, aes(x = Q95, y = fit, group = super.end.group)) +
	theme_bw() +
    geom_rug(data = dat, aes(x = Q95, y = 0.85), sides = "b") +
    ylim(0.85, NA) +
    geom_ribbon(aes(ymin = lwr, ymax = upr), col = NA, fill = "grey",
	alpha = 0.3) +
    geom_line(size = 1.2) +
    facet_wrap(~ super.end.group)





Am 08.06.2017 um 12:15 schrieb Maria Lathouri via R-sig-mixed-models:
--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Rainer Schoppik
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING
Ulf Köther
# 
Dear Maria,

here is an amendment to your gamm plots in lattice-question:

Since I had to procrastinate a little, here is your plot done in
lattice, first defining the positions of the grid lines, then plotting.
You can either choose the polygon-variant which is now commented out in
the code (same as ggplot2-version) or use the calls to panel.xyplot to
produce lines for the CI. The call to panel.abline with x.at and y.at is
a workaround to produce the grid lines aligned with the tick marks:

library(lattice)
library(latticeExtra)

y.at <- pretty(range(c(0.85, 1.1)), 10)
x.at <- pretty(newDat$Q95, 10)

xyplot(fit ~ Q95 | super.end.group, type = "l",
       xlab = "Q95", ylab = "LIFE OE Spring",
       data = newDat, ylim = c(0.85, 1.1),
       scales = list(x = list(at = x.at),
                     y = list(at = y.at)),
       par.settings = list(strip.background = list(col = "lightgrey")),
       panel = function(x, y, subscripts, ...){
           panel.abline(v = x.at,
                        h = y.at, col = "lightgrey")
           panel.xyplot(newDat$Q95[subscripts], newDat$upr[subscripts],
                        type = "l", col = "black", lwd = 2, lty = 2)
           panel.xyplot(newDat$Q95[subscripts], newDat$lwr[subscripts],
                        type = "l", col = "black", lwd = 2, lty = 2)
           # panel.polygon(c(newDat$Q95[subscripts],
	   #		   rev(newDat$Q95[subscripts])),
           #               c(newDat$upr[subscripts],
	   #		   rev(newDat$lwr[subscripts])),
           #               col = "grey", border = NA, ...)
           panel.xyplot(x, y, col = "black", lwd = 2, ...)
	   panel.rug(x = dat$Q95[subscripts], col = 1, end = ...)
})

Have fun..!



Am 19.06.2017 um 15:10 schrieb Ulf K?ther:
--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Rainer Schoppik
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING