[R] plotting gam curve against predictors

[Gavin Simpson]

Suresh Krishna wrote:
> 
> 
> Sonya Ku wrote:
> 
>> Hi I am just beginning to learn R and have fitted several GAM to my
>> species presense/absence data.
>>
>> I have used plot(x,y) using fitted.values as a y variable against
>> predictors. However, it is hard to see general relationships where
>> there is wide spread in predicted values for any x.
>>
>> So I'd like to 1) plot general curves, not individual fitted values 
>> against predictors
> 
> 
> What exactly do you mean by "general curves" or "general relationships" ?
> 
> -s.

I think Sonya meant plotting the response curve of the species over the 
full environmental gradient in question rather then plotting a line that 
goes through the fitted values. If you have a small number of samples or 
they are clumped along that gradient, just plotting a line through the 
fitted points can result in a strange response curve.

The simplest way is to create a new set of predictors variable(s) that 
are evenly spread across the range of your environmental data. Then use 
the predict method to predict the abundance of each species for each 
value of the predictor you created.

e.g. assuming you used mgcv for GAMs (not stated), using dummy code:

# fit the model
gam.mod <- gam(spp ~ s(pH), family = binomial(link = "logit"))
# create some new predictor data aross range of your measured values
# Here we create 100 evenly spaced values
new.dat <- data.frame(pH = seq(min(pH), max(pH), length = 100))
# predict abundance for these new values
gam.pred <- predict(gam.mod, newdata = new.dat, type = "response")
# plot the response curve
plot(new.dat, gam.pred, type = "l")

if you want the standard errors as well then:

# note se.fit = TRUE, gam.pred contains to items in a list $fit and $se.fit
gam.pred <- predict(gam.mod, newdata = new.dat, type = "response",
	se.fit = TRUE)
# so we have to use pam.pred$fit for the predicted values
plot(new.dat, gam.pred$fit, type = "l")
# and +/- the se.fit for the standard errors
lines(new.dat, gam.pred$fit + gam.pred$se.fit, lty = "dotted")
lines(new.dat, gam.pred$fit, type = "l", lty = "dotted")

See ?predict.gam for more information.

I assume the procedure will be similar if you use package gam instead, 
the details of the predict methods may differ though.

HTH

Gav
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Gavin Simpson                     [T] +44 (0)20 7679 5522
ENSIS Research Fellow             [F] +44 (0)20 7679 7565
ENSIS Ltd. & ECRC                 [E] gavin.simpsonATNOSPAMucl.ac.uk
UCL Department of Geography       [W] http://www.ucl.ac.uk/~ucfagls/cv/
26 Bedford Way                    [W] http://www.ucl.ac.uk/~ucfagls/
London.  WC1H 0AP.
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%