First derivative of fitted GAM changes according to specified model distribution

Question

I have built two GAMs, one with family = gaussian specification and one with family = Gamma(link = "log") specification.

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

gam.gamma <- gam(weight_t ~ 
             tagged + 
             sex_t0 +
             s(age.x, k = 6) + 
             s(scale_id, bs = "re") + 
             s(age.x, scale_id, bs = "re"), 
           data = long, 
           method = "REML", 
           family = Gamma(link = "log"))

Both models produce very similar predictions (below), which show that body mass increases non-linearly with age. The depicted/predicted relationship is however not hugely non-linear.

I would now like to get the first derivative of this curve using finite differences to look at how growth changes with age. For gam.gaussian this is easily done:

fd <- derivatives(gam.gaussian, newdata = pred.dat2,  term = "s(age.x)")

f7 <- ggplot(fd, aes(x = data, y = derivative)) + 
  geom_line(size = 1.2) + 
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, colour = NA) +
  theme_classic() + 
  theme(axis.title.x = element_text(face = "bold", size = 14), 
        axis.title.y = element_text(face = "bold", size = 14), 
        axis.text.x = element_text(size = 12), 
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size = 12)) + 
  xlab("Age (days)") + 
  ylab("BM % daily growth (g)"); f7

However, using the same code for gam.gamma produces a completely different plot.

I could not find any info to say exactly what is being plotted, however I suspect that it is specific growth rate (SGR):

SGR = [log(size_t+1) - log(size_t)] / [time_t+1 - time_t]

I assume that SGR is plotted because by specifying a non-gaussian distribution I am essentially saying that body mass does not increase linearly with age? ....maybe?

Both graphs essentially show the same thing - growth rapidly increases early in life and plateaus later in life. However, as the relationship between body mass and age is not hugely non-linear, and the first growth figure is easier and more intuitive to interpret, I am keen to calculate the first derivative from gam.gamma by the same method used for gam.gaussian.

My two questions are:

1. Why does changing the model distribution change the method by which the derivatives are calculated?
2. How would I calculate derivatives from gam.gamma using the same method as what is used for gam.gaussian so that I can produce the same figure (the first growth figure)? Is this an appropriate/reasonable thing to do (given that the relationship between body mass and age is not hugely non-linear and the first growth figure is easier and more intuative to interpret)?

Answer 1

The derivatives() function was intended for computing derivatives for single smooths on the scale of the link function.

For the Gaussian GAM, the link and response scales are the same thing as the link function is the identity function. For the gamma GAM, the link and response scales are different; the derivatives are computed on the log scale here as you used a log link. We do not backtransform those derivatives to the response scale.

The best plot to visualise what you are getting derivatives for is that shown by draw(), which shows the smooths on the link scale (just like plot.gam() does), not model predictions as derivatives() doesn't yet work for model predictions, only single smooths.

If you compare the plots for the smooth of interest produced by draw() for your two models, you should see why the derivatives look quite different.

If you want derivatives of model predictions, you'll need to either do the finite differences yourself using model predictions for the covariate values you want, or wait until I find some time to add this as an option to gratia. If you go with the faster solution (do it yourself) you can use the functionality I wrote to support derivatives(), but much of it is not exported so preface functions with gratia::: or just copy the code into your own functions (respecting the licence if you ever distribute your code).