[R-sig-ME] gamm plots in lattice
Ulf Köther
ukoether at uke.de
Mon Jun 19 16:48:31 CEST 2017
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
More information about the R-sig-mixed-models
mailing list