gam-tmb-vs-mgcv.R

library(glmmTMB)
library(mgcv)
library(magrittr)

## Continuous covariate
x <- seq(1,10, length=100)

## Set up penalized thin-plate regression spline for x
sm <- mgcv::smoothCon(s(x), data=as.data.frame(x))[[1]]
## null space columns
null.sp <- tail(1:ncol(sm$X), sm$null.space.dim)
## Fixed effects (i.e., spline null space model)
Xf <- sm$X[,null.sp]
## Random effects
S <- sm$S[[1]][-null.sp,-null.sp]
L <- mroot(solve(S))
## TPRS random effects design matrix
Xr <- sm$X[,-null.sp] %*% L
## Make Xr orthogonal to Xf
Xr <- (diag(nrow(Xf))-Xf%*%solve(crossprod(Xf))%*%t(Xf))%*%Xr

## Simulated data set
set.seed(123)
b <- rnorm(8,0,4)
beta <- c(1,0.7)
lambda <- exp(Xf%*%beta + Xr%*%b)
y <- rpois(length(x),lambda)

## Fit with glmmTMB
k <- ncol(Xr)
# make a fake grouping variable
g <- rep(1, length(y)) 
## Set up TMB 'map' argument to have only 1 variance parameter for Xr 
##  This is the TPRS smoothing parameter
tmb_map <- list(theta = factor(c(rep(1, k), rep(NA, k*(k-1)/2))))

ftmb <- glmmTMB(y ~ 0+Xf +(0+Xr|g), family=poisson, map=tmb_map)
ptmb <- predict(ftmb,se=T) %>% as.data.frame()

## Fit with mgcv::gam, method='REML'
fmgcv <- mgcv::gam(y~s(x), family=poisson, method="REML")
pmgcv <- predict(fmgcv, se=TRUE, uncond=TRUE) %>% as.data.frame()

## Compare fits
plot(x, y)
lines(x, fitted(ftmb), col='blue', lwd=7)
lines(x, exp(ptmb[,1]+1.96*ptmb[,2]), col='blue', lty=3, lwd=7)
lines(x, exp(ptmb[,1]-1.96*ptmb[,2]), col='blue', lty=3, lwd=7)
text(1, 10, labels=c("glmmTMB"), col='blue', pos=4, cex=2)

lines(x, fitted(fmgcv), col='red', lwd=3)
lines(x, exp(pmgcv[,1]+1.96*pmgcv[,2]), col='red', lty=3, lwd=3)
lines(x, exp(pmgcv[,1]-1.96*pmgcv[,2]), col='red', lty=3, lwd=3)
text(1, 9, labels=c("mgcv::gam"), col='red', pos=4, cex=2)