[R] package mgcv inconsistency in help files? cyclic P-spline "cs" not cyclic?

Joris Meys jorismeys at gmail.com
Fri Jun 4 15:56:24 CEST 2010


Dear all,

I'm a bit stunned by the behaviour of a gam model using cyclic
P-spline smoothers.  I cannot provide the data, as I have about 61.000
observations from a time series.

I use the following model :
testgam <- gam(NO~s(x)+s(y,bs="cs")+s(DD,bs="cs")+s(TT),data=Final)

The problem lies with the cyclic smoother I use for seasonal trends.
The variable Final$y is a numerical variable, going from 1 to 366,
representing the day of the year. I have hourly data from 2003 until
2009, so each day is represented 168 times in the dataset (apart from
366, that one only 48). DD is the wind direction, going from 1 to
3600, and is also modeled with the same cyclic smoother. Yet, if I
check the predictions, the smoother for y is far from cyclic.

I checked the help files ?smooth.terms, and found about 10 lines apart :

bs="cs" specifies a shrinkage version of "cr".

bs="cs" gives a cyclic version of a P-spline.

When I use the (bs="cc") option, I get the results as I want them, so
I keep with the cyclic cubic splines for now. Yet, I find the
behaviour of bs="cs" puzzling, and I'm wondering whether I missed
something, or if this really is an inconsistency in the package.

I currently run mgcv 1.6-1 on R 2.10.1

A small example showing what I experience. Mind you that here x is in
fact NOT cyclic, whereas in my data I'm sure it has to be :
y <- rep(1:20,200)
x <- 1:4000
DD <- sample(1:360,4000,replace=T)
TT <- sample(-10:10,4000,replace=T)
NO <- TT^2 + (10-y+2)^2 + 10*sin(DD*2*pi/360) - 0.002*sqrt(x) +rnorm(4000,0,100)

model <- gam(NO~s(x)+s(y,bs="cs")+s(DD,bs="cs")+s(TT))
plot(model)

model <- gam(NO~s(x)+s(y,bs="cc")+s(DD,bs="cc")+s(TT))
plot(model)


Cheers
Joris

-- 
Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
Joris.Meys at Ugent.be
-------------------------------
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php



More information about the R-help mailing list