AIC to determine optimal degrees of freedom for natural spline in GLMM?

Question

Is it appropriate to use AIC to determine the optimal degrees of freedom for a natural spline?

I have measured 200 animals at six points in time. My data look like below.

plot(long$t_days, long$lweight_t)

To capture the non-linear relationship between animal weight and time I am using a natural spline implemented through the ns() function in the splines package.

I first built the below model, in which my spline has two degrees of freedom. This model has an AICc value of 519.

w8 <- glmmTMB(lweight_t ~ 
+                   tagged + ns(t_days, df = 2) + (t_days | scale_id), 
+                 data = long, family = gaussian) 

AICcmodavg::AICc(w8)
[1] 519.0959

When I check the residuals of this model using the DHARMa package it appears that my model does not adequately capture variability in animal weight through time.

sim_resid_w8 <- simulateResiduals(fittedModel = w8, n = 250)
plot(sim_resid_w8)
plotResiduals(sim_resid_w8, 
+               form = long$t_days[!is.na(long$lweight_t)])

If I increase the number of degrees of freedom in my spline my models AICc values and residuals improve substantially. A spline with 5 degrees of freedom gives me the best AICc score - 276. Beyond 5 degrees of freedom my models AICc score and residuals improve little.

w8 <- glmmTMB(lweight_t ~ 
+                   tagged + ns(t_days, df = 5) + (t_days | scale_id), 
+                 data = long, family = gaussian) 
AICcmodavg::AICc(w8)
[1] 276.1959

sim_resid_w8 <- simulateResiduals(fittedModel = w8, n = 250)
plot(sim_resid_w8)
plotResiduals(sim_resid_w8, 
+               form = long$t_days[!is.na(long$lweight_t)])

Is it appropriate to use AIC in this way to determine the optimal degrees of freedom for a natural spline?

This is the first time I have used splines in GLMMs, but I am aware that overfitting is a problem when using splines and acknowledge that I only have six data points for each of the 200 animals in my dataset. I was unsure if it would be appropriate to have a spline with five degrees of freedom with so few points per individual.

Thanks

Edit: I see this post suggests that using AIC to select optimal knots/degrees of freedom for splines might be ok

Answer 1

Pat, I still see some evidence of heteroscedasticity even after log-transforming the data.

If your interest is in describing the overall shape of the temporal trend in log weight for a typical animal in your study, why not use the bam() function in the mgcv package of R? Then you can consider three different models:

# Model 1: random intercepts model
m1 <- bam(lweight_t ~ tagged + 
                      s(t_days) + 
                      s(scale_id, bs = "re"), 
          data = long)

# Model 2: random intercepts and slopes model 
m2 <- bam(lweight_t ~ tagged + 
                      s(t_days) + 
                      s(scale_id, bs = "re") + 
                      s(scale_id, t_days, bs="re"), 
          data = long)

# Model 3: random smooths model 
m3 <- bam(lweight ~ tagged + 
                    s(t_days) + 
                    s(t_days, scale_id, bs="fs", m=1), 
          data = long)

The paper GENERALISED ADDITIVE MIXED MODELS FOR DYNAMIC ANALYSIS IN LINGUISTICS: A PRACTICAL INTRODUCTION by Márton Sóskuthy (https://arxiv.org/pdf/1703.05339.pdf) does a nice job at explaining the difference between these 3 models.

Notice that you do not have to specify the degree of smoothness of the smooth s(t_days) in your model - it will automatically be estimated.

You can compare the 3 models in terms of their AIC values (as well as adjusted R-squared values and deviance explained) to see which one seems most appropriate for your data.

The link http://jacolienvanrij.com/Tutorials/GAMM.html explains how you can visualize the results produced by each model using the itsadug package in R.

The suggested models are referred to as hierarchical generalized additive models and have also been explored in the paper Hierarchical generalized additive models in ecology: an introduction with mgcv by Eric Pedersen et al., which is available here: https://peerj.com/articles/6876.pdf.