Demo of testing smooth interaction terms in an additive model

demo_smooth_interactions.R

library(tidyverse)
library(mgcv)
library(gratia) #For nicer gam plots
library(viridis) #For nicer colors

#if you don't have above libraries run
# install.packages("tidyverse")
# install.packages("mgcv")
# install.packages("gratia")
# install.packages("viridis")

#In this example with the built-in dataset mtcars, mpg is the outcome (analagous to shock), 
#and hp and disp are continuous covariates (analagous to age and running experience)

#Exploratory data analysis
mtcars %>%
  ggplot(aes(x=disp, y=mpg)) + 
  geom_point(size=3)

mtcars %>%
  ggplot(aes(x=hp, y=mpg)) + 
  geom_point(size=3)

#Use color to visualize outcome
mtcars %>% 
  ggplot(aes(x=disp, y=hp, color=mpg)) + 
  geom_point(size=4) + 
  scale_color_viridis()
#Coverage of the disp-hp space is decent so we can try an interaction term if we want


#Fit additive model (smooth effects, no interaction)
additive_mod <- gam(mpg ~ s(disp, bs="cr", k=5) + s(hp, bs="cr", k=5),
                    data = mtcars,
                    method = "REML")
draw(additive_mod, scales = "fixed") #can use mgcv::plot(additive_mod) but gratia::draw() makes nicer plots
summary(additive_mod)

#Few things to notice:
# you should always use method="REML" unless you really know what you're doing
# bs = "cr" is cubic regression splines, good default choice
# knots: k = 5 based on Jarek's recommendation: use min(35, sqrt(n)) where n is # unique observations (27 & 22 here)
# Notice how both effects are significantly nonlinear - would not be modeled well by a linear model! 
# the p-values here are *approximate*, so don't over-intepret them. 


#Is there evidence for an interaction between disp and hp? Or is an additive model (above) sufficient?
interaction_mod <- gam(mpg ~ s(disp, bs="cr", k=5) + s(hp, bs="cr", k=5) +
                         ti(disp, hp, bs="cr", k=c(5,5)),
                       data = mtcars,
                       method = "REML")
summary(interaction_mod)

#Additional notes
# ti(...) sets up a tensor product basis specifically for testing interaction terms
# p = 0.9878 for the interaction term, so no evidence for interaction

# you can use gratia::draw(interaction_mod, scales="fixed")
# but it is a little misleading until you notice how small the partial effects  
# are for the interaction term (-0.6 to 0.6 vs -10 to 10 for the additive effects)

# the p values for the additive terms shouldn't be interpreted independently; 
# this model's purpose is only to test the significance of the interaction term
# Same cavaeats as above apply re: approximate p-values

# In this example we'd conclude that the additive model (additive_mod) is sufficient
# to explain our data, and we have no evidence for interactions.