# First derivative of fitted GAM changes according to specified model distribution

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",


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)?
• This is apparently a question about an R function called derivatives but you don't tell us where you got this function from or what it is documented to do. There is no function by that name in the core R packages or in the gam package or in the mgcv package. Nov 10, 2020 at 9:20
• @GordonSmyth It's from my gratia package Nov 10, 2020 at 16:23

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

• Is there a way to calculate the total derivative of a variable used in interaction(s)? Example: quantity ~ age + s(price) + s(price, by=age) The reason I'm asking is that this would be useful for demand elasticity modelling, where elasticity is defined as: d(quantity) / d(price) * price / quantity Note: quantity would be the predicted quantity, and price would be the observed price. Jun 18, 2021 at 17:10
• @Gavin Simpson how would you calculate derivatives of smoothers and the corresponding CIs on the response scale? Aug 26, 2022 at 16:17
• @statmerkur The easiest way would be to create two new data frames; one with the values of the covariate you want derivatives for and a second with those values shifted a little. Then a forward difference (the finite difference one usually thinks of). Call these fit_original & fit_forward respectively and eps is the amount you shifted the original values. The first forward finite difference is then: (fit_forward - fit_original) / eps. Getting an uncertainty interval is more involved, requiring some math that I'm not familiar with or using posterior simulation. Sep 5, 2022 at 8:36
• ... so repeat the above steps for a large number of draws from the posterior of the model. This is a little trickier to do as it requires knowledge of some of {mgcv}'s lower-level functions. This is pretty close to the top of my TODO list for {gratia} so it might be implemented in the next couple of weeks. If you can't wait, perhaps ask your question as its own question on this site and I'll write some rough code to show how it would be done. Sep 5, 2022 at 8:38