7
$\begingroup$

There are numerous procedures for functional data clustering based on orthonormal basis functions. I have a series of models built with the GAMM models, using the gamm() from the mgcv package in R. For fitting a long-term trend, I use a thin plate regression spline. Next to that, I introduce a CAR1 model in the random component to correct for autocorrelation. For more info, see eg the paper of Simon Wood on thin plate regression splines or his book on GAM models.

Now I'm a bit puzzled in how I get the correct coefficients out of the models. And I'm even less confident that the coefficients I can extract, are the ones I should use to cluster different models.

A simple example, using:

#runnable code
require(mgcv)
require(nlme)
library(RLRsim)
library(RColorBrewer)

x1 <- 1:1000
x2 <- runif(1000,10,500)

fx1 <- -4*sin(x1/50)
fx2 <- -10*(x2)^(1/4)
y <- 60+ fx1 + fx2 + rnorm(1000,0,5)

test <- gamm(y~s(x1)+s(x2))
# end runnable code

Then I can construct the original basis using smoothCon :

#runnable code
um <- smoothCon(s(x1),data=data.frame(x1=x1),
         knots=NULL,absorb.cons=FALSE)
#end runnable code

Now,when I look at the basis functions I can extract using

# runnable code
X <- extract.lmeDesign(test$lme)$X
Z <- extract.lmeDesign(test$lme)$Z

op <- par(mfrow=c(2,5),mar=c(4,4,1,1))
plot(x1,X[,1],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,X[,2],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,8],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,7],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,6],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,5],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,4],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,3],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,2],ylab="Basis function",xlab="X",type="l",lwd=2)
plot(x1,Z[,1],ylab="Basis function",xlab="X",type="l",lwd=2)
par(op)
# end runnable code

they look already quite different. I can get the final coefficients used to build the smoother by

#runnable code
Fcoef <- test$lme$coef$fixed
Rcoef <- unlist(test$lme$coef$random)
#end runnable code

but I'm far from sure these are the coefficients I look for. I fear I can't just use those coefficients as data in a clustering procedure. I would really like to know which coefficients are used to transform the basis functions from the ones I get with smoothCon() to the ones I extract from the lme-part of the gamm-object. And if possible, where I can find them. I've read the related articles, but somehow I fail to figure it out myself. All help is appreciated.

$\endgroup$
4
$\begingroup$

If I understand correctly, I think you want the coefficients from the $gam component:

> coef(test$gam)
(Intercept)     s(x1).1     s(x1).2     s(x1).3     s(x1).4     s(x1).5 
 21.8323526   9.2169405  15.7504889  -3.4709907  16.9314057 -19.4909343 
    s(x1).6     s(x1).7     s(x1).8     s(x1).9     s(x2).1     s(x2).2 
  1.1124505  -3.3807996  21.7637766 -23.5791595   3.2303904  -3.0366406 
    s(x2).3     s(x2).4     s(x2).5     s(x2).6     s(x2).7     s(x2).8 
 -2.0725621  -0.6642467   0.7347857   1.7232242  -0.5078875  -7.7776700 
    s(x2).9 
-12.0056347

Update 1: To get at the basis functions we can use predict(...., type = "lpmatrix") to get $Xp$ the smoothing matrix:

Xp <- predict(test$gam, type = "lpmatrix")

The fitted spline (e.g. for s(x1)) can be recovered then using:

plot(Xp[,2:10] %*% coef(test$gam)[2:10], type = "l")

You can plot this ($Xp$) and see that it is similar to um[[1]]$X

layout(matrix(1:2, ncol = 2))
matplot(um[[1]]$X, type = "l")
matplot(Xp[,1:10], type = "l")
layout(1)

I pondered why these are not exactly the same. Is it because the original basis functions have been subject to the penalised regression during fitting???

Update 2: You can make them the same by including the identifiability constraints into your basis functions in um:

um2 <- smoothCon(s(x1), data=data.frame(x1=x1), knots=NULL, 
                 absorb.cons=TRUE)
layout(matrix(1:2, ncol = 2))
matplot(um2[[1]]$X, type = "l", main = "smoothCon()")
matplot(Xp[,2:10], type = "l", main = "Xp matrix") ##!##
layout(1)

Note I have not got the intercept in the line marked ##!##.

You ought to be able to get $Xp$ directly from function PredictMat(), which is documented on same page as smoothCon().

$\endgroup$
  • $\begingroup$ it's a start, but somehow I can't get the combined spline using these coefficients or a combination. Care to guide me on how I can go from the basis functions constructed with smoothCon to the final smoother, using only coefficient multiplications (basically doing it by hand)? $\endgroup$ – Joris Meys Dec 2 '10 at 16:29
  • $\begingroup$ Thx for the update, but indeed, they're not the same. Actually, some basis functions are turned upside down even. Thin plate regression splines are indeed a penalized analysis, and the mixed in gamm isn't making things easier either. I'm afraid I'll have to get back to the math as well. $\endgroup$ – Joris Meys Dec 2 '10 at 16:50
  • 1
    $\begingroup$ In relation to Joris' 2nd comment above, my Update 2, shows that the difference is due to um in Joris' code not including the identifiability constraints in the basis functions. If you include these constraints into the basis function, smoothCon() returns the same thing as Xp $\endgroup$ – Reinstate Monica - G. Simpson Dec 3 '10 at 11:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.