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?
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: example.txt
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20170612/3a9fedae/attachment-0001.txt>
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
Hi Duncan
Thank you very much for your help.?
The function model$gam is correct; it applies to gamm4 as well.?
I believe that it is the smooth term the reason that the plots from lattice are
not the same as the plots from the plot(model$gam) but at least I can get the
smooth line and the confidence intervals.?
Kind regards,Maria
???? 4:26 ?.?. ???????, 12 ??????? 2017, ?/? Duncan Mackay <dulcalma at
bigpond.com> ??????:
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
[[alternative HTML version deleted]]