Add support for random slopes in `nbGAM()`

[jr-leary7]

• use when comparing GAM results to the GLMM backend in scLANE
• include random intercept per-subject & random slope on basis spline of pseudotime
• e.g., gamlss(Y ~ bs(PT) + re(random = ~ 1 + bs(PT) | subject), data = sc_data, family = "NBI")