Calculating percentiles (quantiles) from GAM predictions in R

I'm working with a bird migration dataset, and exploring some different metrics to quantify changes in migration phenology (timing) over several decades. There are many different approaches to do this, including fitting a generalized additive model to each year of bird counts, deriving percentiles from the predicted counts over a season, and analyzing trends for the different percentiles by year.

I'm curious what the approach is in R to calculate different percentiles GAM predictions. Using an example dataset from this paper (Fig. 1F) I've fit a negative binomial GAM:

# 2012 birds counts
df_2012 <- structure(list(year = c(2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012), doy = c(205,
206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218,
219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231,
232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244,
245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257,
258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270,
271, 272, 273, 274, 275), birds_per_day = c(NA, NA, NA, NA, NA,
4, 1, 2, 1, 0, 0, 0, 0, 1, 0, 1, 1, 4, 6, 5, 4, 7, 4, 10, 3,
NA, 22, 4, 4, 16, 6, 11, 14, 23, 14, 5, 2, 43, 2, 0, 1, 6, 4,
19, 0, 1, 9, 1, 5, 3, 0, NA, 12, 0, 1, 0, 0, 3, 2, 0, 2, 2, NA,
1, 2, 2, 2, 1, 2, 1, 1)), class = c("spec_tbl_df", "tbl_df",
"tbl", "data.frame"), row.names = c(NA, -71L), spec = structure(list(
cols = list(Species = structure(list(), class = c("collector_character",
"collector")), Season = structure(list(), class = c("collector_character",
"collector")), Year = structure(list(), class = c("collector_double",
"collector")), DOY = structure(list(), class = c("collector_double",
"collector")), X_sp = structure(list(), class = c("collector_double",
"collector")), N = structure(list(), class = c("collector_double",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector")), skip = 1), class = "col_spec"))

# fit gam (number of birds counted per day ~ day of year)
m1 <- gam(
birds_per_day ~ s(doy),
family = nb(),
data = df_2012,
method = "REML"
)
m1

# create new data frame for prediction
newdf <- tibble(
doy = seq(min(df_2012$$doy), max(df_2012$$doy), by = 1)
)

# predictions from GAM
newdf <- cbind(newdf, as.data.frame(predict(m1, newdf, type = "link", se.fit = TRUE)))
newdf <- transform(newdf, fitted = exp(fit), upper_ci = exp(fit + (2 * se.fit)),
lower_ci = exp(fit - (2 * se.fit)))


Here's the prediction plot, with the 10th ('onset' of migration) and 90th ('end' of migration) percentiles added as in Fig. 1F. The percentiles from the fitted GAM should roughly correspond with day 224 (onset of migration) and day 257 (end of migration).

ggplot() +
geom_ribbon(data = newdf, mapping = aes(ymin = lower_ci, ymax = upper_ci, x = doy), alpha = 0.2) +
geom_line(data = newdf, mapping = aes(doy, fitted), size = 1) +
geom_point(data = df_2012, aes(doy, birds_per_day)) +
geom_errorbarh(aes(xmin = 224, xmax = 257, y = 25), color = 'blue') +
annotate("text", x = 224, y = 29, label = "10th \npercentile") +
annotate("text", x = 257, y = 29, label = "90th \npercentile") +
theme_bw() +
labs(
x = 'Day of year',
y = 'Birds per day'
)


• In this example, are you interested in i) predicting the expected count of birds per day from the model for each day in a specific interval, ii) computing the cummulative number of birds from the daily predictions, and then ii) taking the 10th and 90th quantiles of those cummulative birds per day, and returning on what days those quantiles were observed? Dec 3, 2019 at 17:39
• Thanks, clarifying that helped to answer my question...I was confused about what exactly was happening with the quantiles. I was interested in extracting the day of year corresponding to each quantile, which would come from estimating the cumulative number of birds over the course of the season from the daily model predictions. Dec 3, 2019 at 18:30
• Do you need help doing the above with the mgcv::gam() model or are you OK now? (I wasn't clear how/what the quantiles referred to; I may be wrong.) The extra trick will be computing the uncertainty in that interval... which could be done by simulating from the posterior distribution of the model and repeating the above steps on each posterior simulation. Dec 3, 2019 at 19:56
• Hi Gavin - I posted a solution for calculating the day of year (doy) associated with the 10th and 90th percentiles, but I agree it would be more informative to include uncertainty in the doy estimates for the percentiles by repeating the process for posterior simulations (if you don't mind adding that). Dec 3, 2019 at 21:50

Here's the approach I used to calculate the onset (10th percentile) and end (90th percentile) of migration. I created a new column with the cumulative number of birds each day from the predicted daily values, and found the day of year where each percentile was reached:

newdf <- newdf %>%
mutate(
cumulative_birds = cumsum(fitted), # cumulative number of birds
cumulative_perc = cumulative_birds/max(cumulative_birds), # percentage each day
onset = doy[which.max(cumulative_perc >= 0.1)], # 10th percentile, 'onset' of migration
end = doy[which.max(cumulative_perc >= 0.9)] # 90th percentile, 'end' of migration
) %>%
distinct(onset, end)

> newdf
onset end
1   224 257

• Quantiles from GAMs are completely Gaussian-assumption-dependent. Consider semiparametric models that will estimate quantiles, means, etc. in a way that doesn't assume a set distribution. There is a case study in RMS in the chapter on ordinal analysis of continuous Y. Aug 5, 2021 at 12:04