First derivative and confidence intervals from fitted gam

Question

I have created mock data which are available here

My data relate to the growth of animals through time and have the below variables:

read.csv(test.csv, header = T)
test$sex_t0 <- factor(test$sex_t0)
test$tagged <- factor(test$tagged)
test$scale_id <- factor(test$scale_id)

colnames(test)
[1] "weight_t" "age.x"    "sex_t0"   "tagged"   "scale_id"

test$scale_id is a unique identifier

test$tagged indicates if the animal had an ear tag (1) or not (0)

test$sex_t0 indicates if the animal was a male (m) or female (f)

test$age.x is the age of the animal

test$weight_t represents weight measurements at six different time point for each animal, NA shows that the animals had been removed form the study

Using this data I fit the following gam in mgcv

gam1 <- 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 = test, 
+            method = "REML")

This gam includes discrete fixed effect for tagged and sex_t0, separate smooth functions for sexes through time s(age.x, by = sex_t0, k = 5), a random intercept for individual s(scale_id, bs = "re") and a random slope for individuals through time s(age.x, scale_id, bs = "re").

I would now like to calculate and plot the rate of change in growth for tagged and untagged animals and male and female animals through time/against age. This will show if animals grew more rapidly at certain ages or not. To do this we can calculate the first derivatives of the spline functions from our fitted gam. Gavin Simpson has two great blog post on this here and here. There is also an example of calculating the first derivative from a very simple fitted gam here. However, I struggle to follow these examples and can't seem to find an example where someone has calculated the first derivative from a more complex gam that also includes random effects - any help would be very much appreciated?

Edit: I managed to load Gavin Simpson derivSimulCI() function from his GitHub page here, which 'generates posterior simulations for the first derivatives of the spline terms in an additive model'. However, this gives me an error I am yet to solve.

library(devtools)
tmpf <- tempfile()
download.file("https://gist.githubusercontent.com/gavinsimpson/ca18c9c789ef5237dbc6/raw/295fc5cf7366c831ab166efaee42093a80622fa8/derivSimulCI.R",
              tmpf, method = "auto")
source(tmpf)

fd <- derivSimulCI(gam1, samples = 10000)
Error in Summary.factor(c(2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,  : 
  ‘min’ not meaningful for factors

traceback() gives the below

7.
stop(gettextf("%s not meaningful for factors", sQuote(.Generic))) 
6.
Summary.factor(structure(c(2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,  ... at C:\Users\taggarp\AppData\Local\Temp\Rtmp0u2MEM\file8b14519e163a#8
5.
seq(min(x), max(x) - (2 * eps), length = n) at C:\Users\taggarp\AppData\Local\Temp\Rtmp0u2MEM\file8b14519e163a#8
4.
FUN(X[[i]], ...) 
3.
lapply(X = X, FUN = FUN, ...) 
2.
sapply(model.frame(mod)[, m.terms, drop = FALSE], function(x) seq(min(x), 
    max(x) - (2 * eps), length = n)) at C:\Users\taggarp\AppData\Local\Temp\Rtmp0u2MEM\file8b14519e163a#8
1.
derivSimulCI(wt9, samples = 10000)

score 0 · Answer 1 · answered Nov 06 '20 at 01:41

I managed to find an answer using the gratia package...

library(gratia)

# Calculate and plot growth rate for male and female pythons across ages
# Extract first derivatives for smooths of age.x by sex_t0
fd <- fderiv(gam1, newdata = pred.dat,  term = "age.x")

# Calculate 95% confidence intervals for deriviatives
ci <- confint(fd, type = "confidence")

# Attach deriviatives and confidence intervals to pred.dat for plotting
fd.plot1 <- cbind(pred.dat, ci)

ggplot(fd.plot1, aes(x = age.x, y = est, colour = sex_t0, fill = sex_t0)) + 
  geom_line(size = 1.2) + 
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, colour = NA) + 
  scale_colour_manual(labels = c("Female", "Male"), values = c("#F8766D", "#00BFC4")) +
  scale_fill_manual(labels = c("Female", "Male"), values = c("#F8766D", "#00BFC4")) +
  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), legend.title = element_blank()) + 
  xlab("Age (days since hatch)") + 
  ylab("Growth rate (g/day)")

I don't think it is possible to produce a similar figure for tagged and untagged animals as there is no such smooth in gam1

You should use the `derivatives()` function instead of `fderiv()`. These derivatives are all on the scale of the link function at the smooth level. Several people have requested to have derivatives for predictions from the model, be they on the link or response scale. I’ll look into adding this to *gratia* — Gavin Simpson, Nov 07 '20 at 20:27
Thanks Gavin. If I understand correctly this will still give me an approximation of the rate at which animals grew through time over my study? I just need to back-transform the results of ```derivatives()``` to get plot on response and not link scale? — Pat Taggart, Nov 08 '20 at 22:02
Yes, but you could just do the finite difference predictions on the response scale and compute the derivatives from there. — Gavin Simpson, Nov 09 '20 at 23:16

First derivative and confidence intervals from fitted gam

1 Answers1