glmmTMB with smooth terms and random slopes

I'm a newb using glmmTMB, especially with smooths. The first model (m1) gives me what I expect: parallel smooths for each site since there are random intercepts for site.

My question is: is the second model (m2), where I also include a random slope for DOY, legit? In other words, do those smooths represent the site-specific smooths, which is what I want? The reason I ask is that the fixed effect is a smooth, but the random slope term is not. Or do I need to use a different syntax to get what I want? Or is what I want not possible? Thanks!

library(glmmTMB)
library(marginaleffects)
library(ggplot2)
data("Salamanders")

Salamanders$site <- as.character(Salamanders$site) # To avoid weird error in predict
m1 <- glmmTMB(count ~ s(DOY) + (1|site), family = poisson, data = Salamanders, REML = T)
nd <- datagrid(DOY = seq(-2, 1.5, 0.1), site = unique(m1$frame$site), model = m1)
nd$prediction1 <- predict(m1, newdata = nd)
ggplot(nd, aes(DOY, prediction1, group = site)) + geom_line()

m2 <- glmmTMB(count ~ s(DOY) + (DOY|site), family = poisson, data = Salamanders, REML = T)
nd$prediction2 <- predict(m2, newdata = nd)
ggplot(nd, aes(DOY, prediction2, group = site)) + geom_line()

^{Created on 2025-08-24 with reprex v2.1.0}

tl;dr no, the model you have here doesn't have site-specific smooths; it has site-specific linear trends. You should probably switch to mgcv if possible.

The random slope term ((DOY|site)) is interpreted just as it would be in a model without any smooths, i.e. as a term

+ b_{0,i} + b_{1,i} * DOY

with b ~ MVN(0, Sigma).

What you want (I think) is called a factor smooth; this is easily available in mgcv, and you should read Pedersen et al's excellent 2019 paper that explains all about the different ways you can model group-specific smooths.

library(mgcv)
library(ggplot2)
m1 <- gam(count ~ s(DOY, site, bs = "fs"), family = poisson, data = Salamanders)
nd <- expand.grid(DOY = seq(-2, 1.5, 0.1), site = unique(Salamanders$site))
nd$prediction1 <- predict(m1, newdata = nd)
ggplot(nd, aes(DOY, prediction1, group = site)) + geom_line()

Pedersen, Eric J., David L. Miller, Gavin L. Simpson, and Noam Ross. 2019. “Hierarchical Generalized Additive Models in Ecology: An Introduction with Mgcv.” PeerJ 7 (May): e6876. https://doi.org/10.7717/peerj.6876.