[R-sig-ME] gamm plots in lattice

[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:
> 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