[R] How to extract partial predictions, package mgcv

[David Winsemius]

It is nice to see worked examples, but there were some errors that
relate to not including dataframe references. I've paste in code that
runs without error on my machine:

On Sep 16, 2009, at 12:39 PM, Staffan Roos wrote:

>
> 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, data=f)
sequence=order(f$x0)
lines(f$x0[sequence],f$pp_s.x0.[sequence])
lines
rug(jitter(f$x0))

> ##########################################
>
>
> 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?
>>

>> 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.
>>

David Winsemius, MD
Heritage Laboratories
West Hartford, CT