Interpretation of GAM diagnostics

Question

I am hoping to get some advice on my below gam diagnostics and whether my model needs further refining or is adequate as is?

I have conducted a randomised controlled trial to test if tagging impacts the growth of pythons. Pythons were allocated to two groups (tagged vs untagged) and their weight measured at six approximately evenly spaced time points starting from hatching and ending almost 400 days later.

I have built a series of GAMs and used AIC for variable selection. My final GAM is below. It includes:

• a discrete fixed effect for tagged vs untagged pythons
• a discrete fixed effect for sexes
• separate smooths for sexes through time/across python age
• a random intercept for individual python
• a random slope for individual python through time/across age

wt9 <- gam(weight_t ~ 
             tagged + 
             sex_t0 +
             s(age.x, by = sex_t0, k = 5) + 
             s(scale_id, bs = "re") + 
             s(age.x, scale_id, bs = "re"), 
           data = long, 
           method = "REML")

I then used gam.check() to assess diagnostics for my model.

Note that my model is being used to describe trends and not for prediction. I have also included (or tried and then removed) all of the possible explanatory variables I have. In additional I have tried log transforming my response, due to small weights at start of study and large weights at end of study, but this did not improve my residual plots.

gam.check(wt9, rep = 500) 

Method: REML   Optimizer: outer newton
full convergence after 6 iterations.
Gradient range [-0.001751545,0.00288209]
(score 8980.895 & scale 508965.4).
Hessian positive definite, eigenvalue range [0.001751487,559.9235].
Model rank =  411 / 411 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                        k'      edf k-index p-value
s(age.x):sex_t0f  4.00e+00 3.89e+00       1    0.52
s(age.x):sex_t0m  4.00e+00 3.86e+00       1    0.52
s(scale_id)       2.00e+02 9.29e-03      NA      NA
s(age.x,scale_id) 2.00e+02 1.74e+02      NA      NA

I am happy with most diagnostics, but in particular I would appreciate any second opinions on the QQ plot and Resids vs. linear pred. plot? I am unsure if these are substantial enough deviations from the ideal/optimal plots to warrant further model refining (although I am not sure where I would go to from here).

Answer 1

Some of those diagnostics plots look pretty ropey.

First though, it looks like your response can't be negative and yet your model allows it to be negative. Is this the case (that the response is non-negative)?

If so, you likely want to switch to a distribution that has support on the positive (or non-negative) reals, such as the Gamma, (Tweedie), rather than use the Gaussian distribution.

Even though the check for a large enough basis suggests everything is OK, the EDFs of the by smooths are very close to 4. I would increase k to 8 or 10 even, and refit the model to be sure that you really aren't being too restrictive in terms of the basis.

Then I would start decomposing the residuals to see what might be leading to some of the structure we see; try plotting the residuals against the terms in the model and the other explanatory variables you said you;d include but then removed.

But you want to do this after you've switched to a more appropriate distribution if the response can;t be negative.

Instead of smooths of age by sex plus linear effects of age per subject, you might consider having random smooths of age by sex:

wt9 <- gam(weight_t ~ 
             tagged + 
             sex_t0 +
             s(age.x, by = sex_t0, k = 8) +
             s(age.x, scale_id, bs = 'fs'), # <---
           data = long, 
           method = "REML")

where in the indicated line I have added a random factor smooth term that include random intercepts plus random smooths of age by subject.

You could even think of doing that random factor smooth by sex

wt9 <- gam(weight_t ~ 
             tagged + 
             sex_t0 +
             s(age.x, scale_id, bs = 'fs', by = sex_t0), # <---
           data = long, 
           method = "REML")

where this would mean that the sexes could have different smoothness parameters but the individuals of each sex all shared the same smoothness parameter for their sex.