Thanks Simon,
I wasn't aware that "predict.gam" could be used in this situation. I
followed you advice and managed to extract the data I needed.
For people interested, I add the code with comments I used here:
#############################################
# Full code for extracting partial predictions
# Start the package mgcv
library(mgcv)
# Simulate some data (this is from Simon Wood's example in the
# package "mgcv" manual)
n<-200
sig <- 2
dat <- gamSim(1,n=n,scale=sig)
# Check the data
dat
## It looks alright :-)
# Run the GAM
b<-gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat)
# Plot the partial predictions (You may need to rescale your window to
# see the non-linearity
par(mfrow=c(1,3))
plot(b, scale=0)
### Clearly, the variables x0 and x2 are non-linear!
# I would like to extract the "solid line" in the graphs and the
# associated standard errors from the plots. Thus, I use "predict"
# and change to a data.frame:
c<-data.frame(predict(b,type="terms",se.fit=TRUE)$fit)
names(c)<-c("pp_s.x0.","pp_s.I.x1.2..","pp_s.x2.")
d<-data.frame(predict(b,type="terms",se.fit=TRUE)$se.fit)
names(d)<-c("se_s.x0.","se_s.I.x1.2..","se_s.x2.")
# Merge the three datasets:
f=cbind(dat,c,d)
#Restrict the number of variables to the ones I am interested in:
f<-f[,c("x0","pp_s.x0.", "se_s.x0.")]
names(f)
### This data, i.e. "f", can now be exported or used for further
calculations within R
#plot the data
par(mfrow=c(1,1))
plot(pp_s.x0.~x0,type="p",pch=2,cex=0.50)
sequence=order(x0)
lines(x0[sequence],pp_s.x0.[sequence])
lines
rug(jitter(x0))
##########################################
Thanks,
Staffan
--
Staffan Roos, PhD
Research Ecologist
BTO Scotland
Staffan,
Take a look at ?predict.gam. You need to use type="terms" with se=TRUE to
get
the quantities that you need.
Simon
ps. with binary data, method="REML" or method="ML" for the gam fit are
often
somewhat more reliable than the default.
On Monday 14 September 2009 19:30, Staffan Roos wrote:
> Dear package mgcv users,
>
>
>
> I am using package mgcv to describe presence of a migratory bird species
> as
> a function of several variables, including year, day number (i.e.
> day-of-the-year), duration of survey, latitude and longitude. Thus, the
> "global model" is:
>
>
>
> global_model<-gam(present ~ as.factor(year) + s(dayno, k=5) + s(duration,
> k=5) + s(x, k=5) + s(y, k=5), family = binomial(link = logit), data =
> presence, gamma =1.4)
>
>
>
> The model works fine, suggesting that all the variables have strong,
> non-linear, influences on the probability of detecting the species. My
> main
> interest is in the effect "dayno" has on presence, given the inclusion of
> the other explanatory variables. Thus, I would like to extract the values
> of the partial prediction of "dayno" and its associated 2 standard errors
> above and below the main effect, as shown by the code
> "plot(global_model)".
>
>
>
> I have tried to extract the values by using
> "fitted.values(global_model,dayno)", but when plotted, the figure doesn't
> look like the partial prediction for "dayno". Instead, it gives me a very
> scattered figure ("shotgun effect").
>
>
>
> Has anyone tried to extract the partial predictions? If so, please could
> you advise me how to do this?
>
>
>
> Thanks,
>
>
>
> Staffan
>
>
>
> P.S.. For comparison, please have a look at Simon Woods paper in R News,
> 1(2):20-25, June 2001, especially the figures showing the partial
> predictions of Mackerel egg densities. Using those graphs as an example, I
> would like to extract the partial predictions for e.g. "s(b.depth)", given
> the inclusion of the other variables.
>
> --
> Staffan Roos, PhD
> Research Ecologist
> BTO Scotland
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented,
> minimal,
> self-contained, reproducible code.
--
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603 www.maths.bath.ac.uk/~sw283
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
View this message in context:
http://www.nabble.com/How-to-extract-partial-predictions%2C-package-mgcv-tp25441149p25476107.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
