[R-sig-ME] gamm plots in lattice

[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:
> 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)])       }       )
--

_____________________________________________________________________

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