Solved – Use coefficients of thin plate regression splines in a clustering method

clusteringmgcvr

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.

Best Answer

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

Related Question