Confidence intervals for derivatives from GAM predictions

Question

I have built the following GAM model in mgcv

wt9 <- gam(weight_t ~ 
             tagged + 
             sex_t0 +
             s(age.x, by = tagged, k = 5) +
             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 made population averaged predictions from this model.

# Create new data frame to predict from
pred.dat <- data.frame(tagged = c(rep(0, 752), rep(1, 752)), 
                       sex_t0 = c(rep("f", 376), rep("m", 376), rep("f", 376), rep("m", 376)),
                       age.x = c(rep(seq(9, 384, 1), 4)),
                       scale_id = rep(1, 1504))

# Define factors in new data frame
pred.dat$tagged <- factor(pred.dat$tagged)
pred.dat$sex_t0 <- factor(pred.dat$sex_t0)
pred.dat$scale_id <- factor(pred.dat$scale_id)

# Population averaged predictions from fitted gam wt9
preds <- predict(wt9, 
                 newdata = pred.dat, 
                 exclude = c("s(scale_id)", 
                             "s(age.x,scale_id)"), 
                 se = T)

# Combine predictions to new data frame for plotting
pred.dat <- cbind(pred.dat, fit = preds$fit)
pred.dat <- cbind(pred.dat, se.fit = preds$se.fit)

# Calculate 95% CI for predictions from predicted standard errors
pred.dat$lci <- pred.dat$fit - (1.96*pred.dat$se.fit)
pred.dat$uci <- pred.dat$fit + (1.96*pred.dat$se.fit)

# Plot predicted weight for tagged and untagged, and male and female, animals through time (+/- 95% CI)
mycolours1 <- brewer.pal(4, "Blues")[3:4]
mycolours2 <- brewer.pal(4, "Greens")[3:4]

f2a.1 <- ggplot(pred.dat, aes(x = age.x, y = fit, colour = tagged:sex_t0, fill = tagged:sex_t0)) + 
  geom_line(size = 1.5) +
  geom_ribbon(aes(ymin = lci, ymax = uci), alpha = 0.2, colour = NA) + 
  scale_colour_manual(labels = c("Untagged female", "Untagged male", "Tagged female", "Tagged male"), values = c(mycolours1, mycolours2)) +
  scale_fill_manual(labels = c("Untagged female", "Untagged male", "Tagged female", "Tagged male"), values = c(mycolours1, mycolours2)) +
  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)") + 
  ylab("Body mass (g)"); f2a.1

I would now like to create an equivalent figure for the first derivatives of these curves. I can create and plot the derivatives manually but am struggling to get confidence intervals. I have followed this post to obtain the first derivatives manually.

I understand that there is no function to automatically calculate first derivatives and confidence intervals from model predictions. The derivatives() function from the gratia package will calculate first derivatives and confidence intervals from a fitted GAM but not from model predictions, see answer here.

eps <- 1e-7
X0 <- predict(wt9, 
              newdata = pred.dat, 
              exclude = c("s(scale_id)", 
                          "s(age.x,scale_id)"), 
              se = T, 
              type = 'lpmatrix')

pred.datFeps_p <- pred.dat
pred.datFeps_p$age.x <- pred.datFeps_p$age.x + eps

X1 <- predict(wt9, 
              newdata = pred.datFeps_p, 
              exclude = c("s(scale_id)", 
                          "s(age.x,scale_id)"), 
              se = T, 
              type = 'lpmatrix')

# finite difference approximation of first derivative
# the design matrix
Xp <- (X1 - X0) / eps

# first derivative
fd_d1 <- Xp %*% coef(wt9)

test <- cbind(pred.dat, fd_d1)

ggplot(test, aes(x = age.x, y = fd_d1, colour = tagged:sex_t0, fill = tagged:sex_t0)) + 
geom_line(size = 1.2) + 
scale_colour_manual(labels = c("Untagged female", "Untagged male", "Tagged female", "Tagged male"), values = c(mycolours1, mycolours2))

Q: How can I obtain confidence intervals for these first derivative curves? Noting that these are model predictions and not plotted directly from the fitted GAM.

Answer 1

You can do posterior simulation to draw a large set of samples from the posterior distribution of the model, and then for each sample (which is one set of curves if you are predicting for all your groups over a grid of values in age.x) compute the derivatives and store the values.

This gives you a posterior distribution for the derivative of each curve for each value of age.x. You can then summarise those posterior distributions using some suitable quantiles (0.025 and 0.975 probability quantiles say for a 95% interval) to directly estimate the confidence interval on the derivative.

This is all assuming a Gaussian approximation to the posterior, which works OK a lot of the time but can fail in some others. Simon Wood has implemented INLA and a simple Metropolis Hastings sampler to also generate samples from the posterior distribution of the model, which shouldn't have the same problems that Gaussian approximation has in some situations.

If you want to simplify the process (currently only for Gaussian approximation) you can use fitted_samples() from gratia to get the posterior draws for the input data you specify. Then you can compute the derivatives from those posterior draws, etc.