[R] MGCV: overlay fitted (marginal) curves over a plot of the original data

Christoph Scherber Christoph.Scherber at agr.uni-goettingen.de
Tue Jul 16 11:11:28 CEST 2013


Thanks, the sequence of x0 values was clearly too short.

However, is there a way to overlay the (marginal) curve from plot.gam() over a plot of (x,y) values?

Best wishes
Christoph



Am 16/07/2013 11:04, schrieb Simon Wood:
> Probably you didn't want to set x0=0:1? Here is some code, to do what you want.
> (The CI shape is not identical to the plot(b) version as the uncertainty includes
> the uncertainty in the other smooths and the intercept now.)
> 
> library(mgcv)
> set.seed(2)
> dat <- gamSim(1,n=400,dist="normal",scale=2)
> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
> plot(b,select=1)
> 
> plot(y~x0,dat)
> mydata=data.frame(x0=0:200/200,x1=mean(dat$x1),x2=mean(dat$x2),x3=mean(dat$x3))
> pv <- predict(b,mydata,type="response",se=TRUE)
> lines(mydata$x0,pv$fit)
> lines(mydata$x0,pv$fit+2*pv$se.fit,lty=2)
> lines(mydata$x0,pv$fit-2*pv$se.fit,lty=2)
> 
> 
>   
> 
> On 16/07/13 09:52, Christoph Scherber wrote:
>> Dear R users,
>>
>> I´ve stumbled over a problem that can be easily seen from the R code below:
>>
>> - When I use plot.gam() on a fitted model object, I get a nice and well-looking smooth curve for all
>> terms in the model.
>>
>> - However, when I use predict(model) for a given predictor, with values of all other predictors set
>> to their means, the resulting curve doesn´t fit well at all.
>>
>> Is there a way to "overlay" the curve produced by plot.gam() over a plot of the original data?
>>
>> Here´s some reproducible code with mgcv version 1.7-22 on R3.0.1 (Windows 7):
>>
>> ##
>>
>> library(mgcv)
>> set.seed(2)
>> dat <- gamSim(1,n=400,dist="normal",scale=2)
>> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat)
>> plot(b,select=1)
>>
>> plot(y~x0,dat)
>> mydata=data.frame(x0=0:1,x1=mean(dat$x1),x2=mean(dat$x2),x3=mean(dat$x3))
>> lines(0:1,predict(b,mydata,type="response"))
>>
>> ##
>>
>> Best wishes,
>> Christoph
>>
>>
> 
> 
> 



More information about the R-help mailing list