[R] plotting gamm results in lattice
Duncan Mackay
dulcalma at bigpond.com
Mon Jun 12 17:26:27 CEST 2017
Hi Maria
If you have problems just start with a small model with predictions and then plot with xyplot
the same applies to xyplot
Try
library(gamm4)
spring <- dget(file = "G:/1/example.txt")
str(spring)
'data.frame': 11744 obs. of 11 variables:
$ WATERBODY_ID : Factor w/ 1994 levels "GB102021072830",..: 1 1 2 2 2 3 3 3 4 4 ...
$ SITE_ID : int 157166 157166 1636 1636 1636 1635 1635 1635 134261 1631 ...
$ Year : int 2011 2014 2013 2006 2003 2002 2013 2005 2013 2006 ...
$ Q95 : num 100 100 80 80 80 98 98 98 105 105 ...
$ LIFE.OE_spring : num 1.02 1.03 1.02 1.06 1.06 1.07 1.12 1.05 1.14 1.05 ...
$ super.end.group : Factor w/ 6 levels "B","C","D","E",..: 1 1 3 3 3 2 2 2 4 4 ...
$ X.urban.suburban : num 0 0 0.07 0.07 0.07 0.53 0.53 0.53 8.07 8.07 ...
$ X.broadleaved_woodland: num 2.83 2.83 10.39 10.39 10.39 ...
$ X.CapWks : num 0 0 0 0 0 0 0 0 0 0 ...
$ Hms_Poaching : int 0 0 10 10 10 0 0 0 0 0 ...
$ Hms_Rsctned : int 0 0 0 0 0 0 0 0 0 0 ...
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))
Warning message:
Some predictor variables are on very different scales: consider rescaling
# plot(model$gam, page=1, font.lab=2, xlab="Residual Q95")
M <-predict(model$gam,type="response",se.fit=T)
# joining the dataset and the predictions keep it "in house" and on path
springP = cbind(spring, M)
springP$upper = with(springP,fit+1.96*se.fit)
springP$lower = with(springP,fit-1.96*se.fit)
str(springP)
'data.frame': 11744 obs. of 15 variables:
$ WATERBODY_ID : Factor w/ 1994 levels "GB102021072830",..: 1 1 2 2 2 3 3 3 4 4 ...
$ SITE_ID : int 157166 157166 1636 1636 1636 1635 1635 1635 134261 1631 ...
$ Year : int 2011 2014 2013 2006 2003 2002 2013 2005 2013 2006 ...
$ Q95 : num 100 100 80 80 80 98 98 98 105 105 ...
$ LIFE.OE_spring : num 1.02 1.03 1.02 1.06 1.06 1.07 1.12 1.05 1.14 1.05 ...
$ super.end.group : Factor w/ 6 levels "B","C","D","E",..: 1 1 3 3 3 2 2 2 4 4 ...
$ X.urban.suburban : num 0 0 0.07 0.07 0.07 0.53 0.53 0.53 8.07 8.07 ...
$ X.broadleaved_woodland: num 2.83 2.83 10.39 10.39 10.39 ...
$ X.CapWks : num 0 0 0 0 0 0 0 0 0 0 ...
$ Hms_Poaching : int 0 0 10 10 10 0 0 0 0 0 ...
$ Hms_Rsctned : int 0 0 0 0 0 0 0 0 0 0 ...
$ fit : num [1:11744(1d)] 1.03 1.04 1.04 1.02 1.01 ...
..- attr(*, "dimnames")=List of 1
.. ..$ : chr "1" "2" "3" "4" ...
$ se.fit : num [1:11744(1d)] 0.00263 0.00266 0.00408 0.00408 0.00411 ...
..- attr(*, "dimnames")=List of 1
.. ..$ : chr "1" "2" "3" "4" ...
$ upper : num [1:11744(1d)] 1.03 1.04 1.05 1.03 1.02 ...
..- attr(*, "dimnames")=List of 1
.. ..$ : chr "1" "2" "3" "4" ...
$ lower : num [1:11744(1d)] 1.02 1.03 1.03 1.01 1 ...
..- attr(*, "dimnames")=List of 1
.. ..$ : chr "1" "2" "3" "4" ...
# this produces a mess of lines for the upper and lower
xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = springP,
lower = springP$lower,
upper = springP$upper,
subscripts = TRUE,
panel = function(x,y, upper, lower, subscripts, ...){
panel.xyplot(x,y, type="smooth")
panel.xyplot(x[order(x)], upper[subscripts][order(x)], lty=2, col="red")
panel.xyplot(x[order(x)], lower[subscripts][order(x)], lty=2, col="red")
#panel.loess(x,y,...) # have not tried to fix these lines- depends on what you want to do
#panel.rug(x = x[is.na(y)], y = y[is.na(x)])
}
)
# smoothing them produces reasonable lines
xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = springP,
lower = springP$lower,
upper = springP$upper,
subscripts = TRUE,
panel = function(x,y, upper, lower, subscripts, ...){
panel.xyplot(x,y, type="smooth")
panel.xyplot(x[order(x)], upper[subscripts][order(x)], type = "smooth", lty=2, col="red")
panel.xyplot(x[order(x)], lower[subscripts][order(x)], type = "smooth", lty=2, col="red")
#panel.loess(x,y,...)
#panel.rug(x = x[is.na(y)], y = y[is.na(x)])
}
)
# by newdata = ...
# take a cross-section of data - reduces the amount of data to be plotted
mx = aggregate(Q95 ~ super.end.group, spring, max, na.rm=T)
newlst <-
do.call(rbind,
lapply(1:6, function(j) expand.grid(Q95 = seq(0,mx[j,2], length = 50), super.end.group = LETTERS[j])
Year =
Hms_Rsctned =
...)) # fill in yourself
newd <- data.frame(newlst)
and predict using newdata = newd
# starting with a smaller model
model2 <-
gamm4(LIFE.OE_spring~s(Q95, by=super.end.group), data=spring,
random=~(1|WATERBODY_ID/SITE_ID))
newlst <-
do.call(rbind,
lapply(1:6, function(j) expand.grid(Q95 = seq(0,mx[j,2], length = 50), super.end.group = LETTERS[j])) )
newd <- data.frame(Q95 = newlst)
# The following is untested
# I have not used gamm4 for a while; is model$gam correct? Had problems here - have no time to go into the help pages
M3 <- predict(model2$gamm, newdata = newd, type="response",se.fit=T)
# keep together
newdat = cbind(newd, M3)
# If you want CIs
newdat <- within(newdat, {
lower = fit-1.96*se.fit
upper = fit+1.96*se.fit
})
# plot- something like
xyplot(fitted(model$gam) ~ Q95 |super.end.group, data = springP,
Q2 = newddat$Q95,
lower = newdat$lower,
upper = newdat$upper,
subscripts = TRUE,
panel = function(x,y, Q2, upper, lower, subscripts, ...){
panel.xyplot(x,y, type="smooth")
panel.xyplot(Q2, upper, lty=2, col="red")
panel.xyplot(Q2, lower, lty=2, col="red")
#panel.loess(x,y,...)
#panel.rug(x = x[is.na(y)], y = y[is.na(x)])
}
)
If you are putting a key onto the plot see
? xyplot and par.settings and key
Regards
Duncan
Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au
-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Maria Lathouri via R-help
Sent: Monday, 12 June 2017 20:57
To: R-help Mailing List
Subject: [R] plotting gamm results in lattice
Dear all,
I hope that you can help me on this. I have been struggling to figure this out but I haven't found any solution.
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))
plot(model$gam, page=1, font.lab=2, xlab="Residual Q95")
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
More information about the R-help
mailing list