Algorithm for selecting interactions without very many degrees of freedom (mgcv, gam)

Question

I've got a semi-parametric model that I'm fitting with GAM's (mgcv in R). It is of the form

$$y = \theta + X'\beta + f(Z) + \text{Interactions!} + \epsilon$$

I've got 289 observations.

Many of the interactions that I'd like to specify are tensor (or Kronecker) product interactions. Unpenalized, each interaction takes up $k^n$ degrees of freedom, where k is the number of knots and n is the number of continuous variables in the interaction. Of course, penalization reduces the effective degrees of freedom, but GAM won't fit a model with more to-be-penalized parameters than degrees of freedom.

I've already coerced the models to specify 4 knots per Z variable and specified their location based on subject knowledge -- this isn't optimal in the sense that Simon Wood's thin plate regression splines are optimal, but it seems like a not-bad compromise in order to keep k down.

Anyway, I have a long candidate list of interactions that might be important. The problem is that a model that specifies all of them has >1000 degrees of freedom. I can't use a model selection algorithm that starts saturated and eliminates. I need to use one that starts with only the main effects (all of which I knot that I want to include a priori), and then adds terms.

The problem is that adding one interaction might change the relevance of others, or of a set of others. There are probably hundreds of combinations of interactions that are accessible given my df constraints.

I'm looking for ideas for model selection besides trial and error.

(I'm comparing models by "lowest AIC wins")

Anybody have any ideas or experience with this sort of problem?

Answer 1

One option is to use a boosted GAM. "Boosting" is a generic method to extend many models (linear, GAM, tree-based, etc). In practice, boosting works a lot like LASSO, in that it performs variable selection during the cross-validation tuning process. Boosting builds the model slowly and sequentially, so it can also handle "p>n" models.

Here's code that performs boosted GAMs (with 1st order interaction terms) using the mtcars dataset:

# Boosted GAM Example

# Citations
citation('mboost')
citation('gamboostLSS')

# Libraries
library(data.table)    # For data manipulation
library(mgcv)          # gam()
library(gamboostLSS)   # gamboostLSS() # Loads mboost
# library(gamlss.dist) # Optional. For additional distributions, if needed.
library(ggplot2)       # ggplot()
theme_set(theme_classic())

# Clear all
rm(list=ls(all=TRUE))

# Example data
data(mtcars)
df=copy(mtcars)
rm(mtcars)
setDT(df)
dim(df) # 32 xx - VERY small sample size!

# Explore data
# Let's use "mpg" as our response, and focus on these 4 predictors.
pairs(df[,.(mpg,disp,hp,drat,wt)])

# Modeling functions notes
# ========================.
# Use mboost() or glmboost() for boosted linear regression.
# Use gamboost() for boosted GAMs.
# Use gamboostLSS() if you wish to conditionally model the scale parameter.

# Boost modeling settings
nu=0.01 # Shrinkage paramter. Default= 0.1. Try values between 0.01 and 0.1.
mstopMax=1000L # Number of boosting iterations. Will be trimmed down later using CV. Default = 100. Setting "mstopMax" and "nu" values is a balancing act.
K=5L # Set user global default for smoothing prior. Default = 4.

# Fit
rm(fit)
fit=gamboost(mpg~
                    # Main effects:
                    bbs(disp,df=K)+ 
                    bbs(hp,df=K)+
                    bbs(drat,df=K)+
                    bbs(wt,df=K)+
                    # 1st order interaction terms. Must include main effects separately.
                    bbs(disp,by=hp,df=K)+
                    bbs(disp,by=drat,df=K)+
                    bbs(disp,by=wt,df=K)+
                    bbs(hp,by=drat,df=K)+
                    bbs(hp,by=wt,df=K)+
                    bbs(drat,by=wt,df=K),
                 data=df,
                 control=boost_control(mstop=mstopMax,nu=nu,trace=F))

# Browse methods for mboost and gamboost objects
class(fit) # gamboost mboost
?predict.gamboost

# CV settings
tmp=10 # Calculate CV error every "tmp" boosting iterations
length.out=mstopMax/tmp # Used in make.grid()
# Note: cvrisk() can take a while to run. Set "B" low at first, to see if your "nu" and "mstopMax" are in the ballpark.

# CV tuning - Takes a while
rm(cv)
set.seed(123)
cv=cvrisk(fit,
             folds=cv(model.weights(fit),type='bootstrap',B=50), # Default is 'bootstrap' with B=25. Can also specify 'kfold' with B=10.
             grid=make.grid(mstopMax,length.out=length.out,min=1,log=F))
plot(cv)
mstop(cv) # 223

# All boosted models will converge to a non-boosted version, after enough iterations.
# For this model, CV error was minimized at ~223 iterations.
# Use mstop() to identify this value.

# Determine optimal number of boosting iterations
rm(mstop)
mstop=mstop(cv)
mstop # 223 optimal iterations

# More notes:
# Unlike bagging (eg, random forest), boosting is sequential.
# "fit" object stores ALL boosted iterations.
# Example:
fit           # 1000 iterations, as initially specified
fit[10]       # Set to 10 iterations
fit           # still 10 iterations
fit[30]       # However, the other iterations are not lost :)
fit           # still 30 iterations
fit[mstopMax] # Reset to max iterations

# Plot partial effects
# See how shapes change, if we change the number of iterations
# Overfit model: All 4 main effects and 3 interaction terms were selected.
par(mfrow=c(3,3),mar=c(4.5,4.5,1,1))
plot(fit[500],type='o')

# Underfit model: 2 main effects and 1 interaction terms were selected. Shape of 2 main effects is smoother.
par(mfrow=c(2,3),mar=c(4.5,4.5,1,1))
plot(fit[80])

# CV tuned model: 3 main effects and 3 interaction terms were selected.
par(mfrow=c(3,3),mar=c(4.5,4.5,1,1))
plot(fit[mstop])
par(mfrow=c(1,1))

# Coefs
coef(fit[mstop])
coef(fit[mstop])$`bbs(disp, df = K)`
# coef() is more useful for boosted linear regression...

# CIs - Takes a while to run
rm(ci)
ci=confint(fit[mstop],level=0.95,B=30,B.mstop=10) # Number of iterations set way too low. Actual run will need to set "B" and "B.mstop" much larger.
warnings()

# Plot CIs
par(mfrow=c(2,3),mar=c(4.5,4.5,1,1))
# plot(fit[mstop])
plot(ci,which=1)
plot(ci,which=2)
plot(ci,which=4)
par(mfrow=c(1,1))

# Calculate yhat and pseudo-R2
df[,yhat:=NA_real_]
df[,yhat:=predict(fit[mstop],type='response')]
R2=cor(df$mpg,df$yhat,method='pearson')^2
R2 # 0.9050343

# Plot 1:1
ggplot(df,aes(y=mpg,x=yhat))+
    geom_point()+
    geom_abline(intercept=0,slope=1)

# Generate newdata
xgrid=seq(min(df$disp),max(df$disp),length=100)
rm(newdata)
newdata=expand.grid(disp=xgrid,hp=mean(df$hp),drat=mean(df$drat),wt=mean(df$wt),
                          KEEP.OUT.ATTRS=F,stringsAsFactors=F)
setDT(newdata)
newdata[,yhat:=predict(fit[mstop],newdata=newdata,type='response')]

# Plot fit, assuming other predictors set to their means
ggplot(df,aes(y=mpg,x=disp))+
    geom_point()+
    geom_line(data=newdata,mapping=aes(x=disp,y=yhat),color=4,size=1)

See mboost and gamboostLSS for more information.