Collin Edwards

#rstats question for mgcv

If I want to test whether there is support for a smooth x factor effect, I can just include the smooth and the smooth x factor as separate terms to create nested models, right? As in

library(mgcv)
set.seed(0)
dat<-gamSim(5,n=200,scale=2)
##make explicit factor
dat$fac = factor(letters[dat$x0])

m = gam(y ~ s(x1) + fac + s(x1, by = fac), data = dat)
m1 = gam(y ~ s(x1) + fac, data = dat)
anova(m, m1, test = "F")

@noamross @gavinsimpson

Jul 26, 2023 at 7:44pmWeb
Show threadCollin EdwardsJul 26, 2023

@noamross @gavinsimpson

It feels like the `by = factor` term is redundant there, and should lead to issues, but I suspect that just means I'm not thinking about shrinkage and wiggliness appropriately. Certainly the code runs that way.

Show threadNoam RossJul 26, 2023

@collinedwards @gavinsimpson always have to look up our own paper for this :)

"group-specific intercepts are not incorporated into factor `by` variable smoothers (as would be the case with a factor smoother or a tensor product random effect)."

https://peerj.com/articles/6876/#p-78

Hierarchical generalized additive models in ecology: an introduction with mgcv

In this paper, we discuss an extension to two popular approaches to modeling complex structures in ecological data: the generalized additive model (GAM) and the hierarchical model (HGLM). The hierarchical GAM (HGAM), allows modeling of nonlinear functional relationships between covariates and outcomes where the shape of the function itself varies between different grouping levels. We describe the theoretical connection between HGAMs, HGLMs, and GAMs, explain how to model different assumptions about the degree of intergroup variability in functional response, and show how HGAMs can be readily fitted using existing GAM software, the mgcv package in R. We also discuss computational and statistical issues with fitting these models, and demonstrate how to fit HGAMs on example data. All code and data used to generate this paper are available at: github.com/eric-pedersen/mixed-effect-gams.

PeerJ
Show threadNoam RossJul 26, 2023
@collinedwards @gavinsimpson
So fac + s(x1, by = fac) are not redundant at all. Whether there are collinearity issues with s(x1) and s(x1, by = fac) depends, which is why in the paper we have an example with an m=1 penalty. That's a bit of matter of what you are testing for and how you are interpreting difference, and if you have collinearity issues exploding your CIs.
Show threadCollin EdwardsJul 27, 2023
@noamross @gavinsimpson Awesome, thanks for all this! I'll play around to see if we're running into collinearity issues, and this paper seems super helpful for a whole host of projects!
Show threadGavin SimpsonJul 27, 2023

@collinedwards @noamross I think in this instance you'd be better off numerically using

y ~ s(x1) + s(x1, fac, bs = "sz")

& compare that with

y ~ s(x1) + fac

as the `sz` basis will set things up so that those difference smooths are orthogonal to the main f(x1) term. This sz basis also includes the group means hence no fac

Fit both models with `method = "ML"` if you are going to attempt a GLRT but do read `?anova.gam` about the multi-model feature of that function

Show threadNoam RossJul 27, 2023
@gavinsimpson @collinedwards 👆I keep forgetting about the newer "sz" basis
Show threadDavid Lawrence MillerJul 28, 2023
@noamross @gavinsimpson @collinedwards same!
Show threadCollin EdwardsJul 28, 2023

@gavinsimpson @noamross This is awesome!! Thanks for all the help!

We're already doing comparisons with AIC, and the differences are quite large. I was just thinking providing a reasonable P value might help with the review process. So I'm not super worried that the P values are biased towards being low.

Show threadChristos Argyropoulos MD, PhDJul 27, 2023
@collinedwards @noamross @gavinsimpson hm, not quite. The smooth terms map to generalized linear mixed models , so if you fit the gam with REML , you will be comparing GLMM with different fixed and random effects with GLR tests which in general is a no-no.
Show threadChristos Argyropoulos MD, PhDJul 27, 2023
@collinedwards @noamross @gavinsimpson one thing you could do through mgcv is to get the design matrices for your model and use the low level functions of lmer/nlme and fit it using regular max likelihood if you want a p value. If you don't want a p-value, you can fit the more general model and have mgcv penalize the terms away using null space penalization from within mgcv.
Or you can use AIC and appeal to your audience's common sense
Show threadCollin EdwardsJul 28, 2023
@ChristosArgyrop @noamross @gavinsimpson Hahaha, definitely we're already using AIC. My thinking is to do double duty by also providing P values (with caveats), just so that we don't have less statistically aware reviewers turn down the analysis.
Show threadChristos Argyropoulos MD, PhDJul 28, 2023
@collinedwards @noamross @gavinsimpson I had to do this once and used the selection smoothers (a mix of L0+L2 penalties in the current formulation). Reviewers did not complain as the interaction was expected in the system under consideration . EVEN REVIEWER 2 WAS OK WITH IT
Show threadCollin EdwardsJul 28, 2023
@ChristosArgyrop @noamross @gavinsimpson EVEN REVIEWER 2??
Show threadChristos Argyropoulos MD, PhDJul 29, 2023
@collinedwards @noamross @gavinsimpson It was a very rewarding experience, yet as healthcare professional I could not shake the feeling they were experiencing a stroke that left them unable to reject.