Hello,
I am trying to generate the design matrix associated with the tprs smooth in
R mgcv from scratch. I want to construct the matrix X = [UkDkZk | T] which
is produced from the function smoothCon following the steps provided in
Simon Wood''s 2003 Thin Plate Regression Splines paper (Appendix A). I
have
had a crack at it but so far I am unable to replicate the design matrix
exactly. See below for an example.
My questions are twofold. 1) What transformations are being applied to the
matrix T (with absorb.cons=FALSE) behind the scenes in smoothCon? And 2)
what transformations are being applied to UkDkZk in smoothCon (if any)?
Thanks!
Casey
######################### Libraries
library(mgcv)
library(SpatioTemporal) # uses the crossDist() function to calculate
distances (ie. r in eta_{md}(r))
######################### Generate Data
x1 <- sin(rnorm(100))
x2 <- rnorm(100)*x1
x2 <- x2[order(x1)]
x1 <- sort(x1)
######################## Use mgcv functions to get the design matrix [UkDkZk
| T]
um <- smoothCon(s(x1, x2, fx=TRUE, m=2, k=20), data=data.frame(x1=x1,
x2=x2), knots=NULL, absorb.cons=FALSE)
X1 <- um[[1]]$X
######################### Use the steps in Appendix 1 of Wood (2003) to
construct the same matrix
n <- length(x1); m=2; d=2; M <- choose(m+d-1, d); k <- 20
eta.func <- function(m, d, r){
if(d%%2 == 0){
(-1)^(m+1+d/2)/(2^(2*m-1)*pi^(d/2)*factorial(m-1)*factorial(m-d/2))*r^(2*m-d
)*log(r)
} else {
gamma(d/2-m)/(2^(2*m)*pi^(d/2)*factorial(m-1))*r^(2*m-d)
}
}
T <- cbind(rep(1, length(x1)), x1, x2)
E <- eta.func(m=m, d=d, r=crossDist(cbind(x1, x2)))
diag(E) <- 0
Ek <- eigen(E, symmetric=TRUE)
Uk <- Ek$vectors[, 1:k]
Dk <- diag(Ek$values[1:k])
Ek <- Uk%*%Dk%*%t(Uk)
tUkT <- t(Uk)%*%T
QR <- qr(tUkT)
Q <- qr.Q(QR, complete=TRUE)
Zk <- Q[, (M+1):k]
X2 <- cbind(Uk%*%Dk%*%Zk, T)
####################### Compare Basis Functions
par(mfrow=c(1,2))
plot(X1[, 1] ~ x1, type="l", ylim=range(X1))
for(i in 2:ncol(X1)){ lines(X1[, i] ~ x1, col=i) }
plot(X2[, 1] ~ x1, type="l", ylim=range(X1))
for(i in 2:ncol(X2)){ lines(X2[, i] ~ x1, col=i) }
####################### Compare Predictions (unpenalized)
y <- x1+2*x2 + rnorm(100)
mod1 <- lm(y ~ -1 + X1)
pred1 <- X1%*%coef(mod1)
mod2 <- lm(y ~ -1 + X2)
pred2 <- X2%*%coef(mod2)
mod.gam <- gam(y ~ s(x1, x2, fx=TRUE, m=2, k=20))
pred.gam <- predict(mod.gam, data.frame(x1, x2))
par(mfrow=c(1,3))
plot(pred1 ~ y, ylim=c(-3, 3))
abline(0,1)
plot(pred2 ~ y, ylim=c(-3, 3))
abline(0,1)
plot(pred.gam ~ y, ylim=c(-3, 3)) # pred.gam and pred1 are the same
abline(0,1)
[[alternative HTML version deleted]]