2
$\begingroup$

I am using the mgcv package in R to fit logistic GAMs to survey data. In one of my models I use an interaction between two covariates. I am currently trying to fit this model using RStan and so am constructing the basis functions myself using B-splines. However I am struggling to reproduce the basis functions from the model.matrix output of a GAM including a ti() smooth term of one variable. I illustrate the problem below:

#Set-up
library(mgcv)
set.seed(1)

#Generate random data to fit gam to
y <- rnorm(1000,0,1)
x <- sample(1:30,1000,replace=T) # integers between 1 and 30

#Fit gam with 1D s() smooth term with a B-spline basis (bs = "ps")
gam_s <- gam(y ~ s(x, bs = "ps", k = 6, m=c(2,1)))

#Extract model matrix
X_s <- model.matrix(gam_s)

#Fit gam with 1D ti() smooth term with a B-spline basis (bs = "ps")
gam_ti <- gam(y ~ ti(x, bs = "ps", k = 6, m=list(c(2,1))))

#Extract model matrix
X_ti <- model.matrix(gam_ti)

#Compute vector of indexes of unique elements of x for easy plotting
index <- numeric()
for (i in 1:30)
index[i] <- which(x==i)[1]

#Plot results
par(mfrow=c(1,2))
plot(0,0,xlim=c(1,30),ylim=range(X_s),xlab="x",ylab="B(x)",type="n",main="1D s() basis functions for a B-spline basis")
for (i in 2:6)
  lines(1:30,X_s[index,i],col=i-1)
plot(0,0,xlim=c(1,30),ylim=range(X_ti),xlab="x",ylab="B(x)",type="n",main="1D ti() basis functions for a B-spline basis")
for (i in 2:6)
  lines(1:30,X_ti[index,i],col=i-1)

#Do the same for a cubic regression spline basis (bs = "cr")
gam_s <- gam(y ~ s(x, bs = "cr", k = 6))
X_s <- model.matrix(gam_s)
gam_ti <- gam(y ~ ti(x, bs = "cr", k = 6))
X_ti <- model.matrix(gam_ti)
plot(0,0,xlim=c(1,30),ylim=range(X_s),xlab="x",ylab="B(x)",type="n",main="1D s() basis functions for a cubic regression spline basis")
for (i in 2:6)
  lines(1:30,X_s[index,i],col=i-1)
plot(0,0,xlim=c(1,30),ylim=range(X_ti),xlab="x",ylab="B(x)",type="n",main="1D ti() basis functions for a cubic regression spline basis")
for (i in 2:6)
  lines(1:30,X_ti[index,i],col=i-1)

Running this R code, you should see that the first two plots using the B-spline basis and s()/ti() are different whereas the next two plots using the cubic regression spline basis and s()/ti() are identical. My question is why should the basis make a difference? For the B-spline basis I have worked out how to derive the s() plot - the original basis is transformed using a QR decomposition in order to satisfy the sum-to-zero constraint. But the ti() plot doesn't look like anything you could get by transforming a B-spline basis. The 1D basis functions are important because I use the tensor product of them to model the interaction. Any help would be greatly appreciated!

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.