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