4
$\begingroup$

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)

Animal weight through time

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)])

DHARMa plot 1 DHARMa plot 2

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)])

DHARMa plot 3 DHARMa plot 4

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

$\endgroup$
3
  • $\begingroup$ Pat, is your animal weight log-transformed prior to being included in the model? If yes, that would complicate interpretation. Anyway, from your plot of the data, it looks like the 6 occasions when the animals were measured were pre-specified - is that correct? If yes, why not code Time as a factor in your model, with levels corresponding to the 6 occasions? Then your model can include tagged + Time in its fixed effects person and Time in its random effects portion. You might even be able to include tagged*Time in your model. $\endgroup$ Commented Sep 23, 2020 at 2:07
  • $\begingroup$ Thanks Isabella. The 6 occasions when the animals were measured were approximately 2-3 months apart, it was intended that they were approximately evenly spaced time points. I had previously coded time as a factor, but was told by a reviewer that they wanted time as a numeric and something to account for the non-linear relationship between time and weight..... $\endgroup$ Commented Sep 23, 2020 at 2:54
  • $\begingroup$ Yes my data have been log transformed to improve homoscedasticity. Variance was originally very uneven across the 6 time points, very little variance at time point 1 and lots of variance at time point 6, variance increased through time $\endgroup$ Commented Sep 23, 2020 at 2:56

1 Answer 1

3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.