0
$\begingroup$

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

enter image description here

$\endgroup$
4
  • 1
    $\begingroup$ 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? $\endgroup$ – Gavin Simpson Dec 3 '19 at 17:39
  • $\begingroup$ 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. $\endgroup$ – Jason Dec 3 '19 at 18:30
  • $\begingroup$ 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. $\endgroup$ – Gavin Simpson Dec 3 '19 at 19:56
  • $\begingroup$ 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). $\endgroup$ – Jason Dec 3 '19 at 21:50
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.