projection/prediction of NSERC Discovery award dates

nserc_dates.R

library(tidyverse); theme_set(theme_bw(base_size=18))
library(png)
library(mgcv)

## https://stackoverflow.com/questions/9917049/inserting-an-image-to-ggplot2
## dates are approximately when university research offices notify recipients
##  (before ResearchNet posting, dates may vary by university)
this_year <- 2026
(df <- read.table(text="
2025 April 1
2024 March 27
2023 March 29
2022 April 11
2021 April 14
2020 April 23
2019 April 11
2018 April 13
2017 April 12
2016 April 7
2015 April 1
2014 April 4
2013 March 21
")
    %>% transmute(
            year=V1,
            date=lubridate::ymd(paste(this_year,V2,V3)))
)

## FIXME: add manual legend (blue = thin plate spline, purple = cubic spline, red = Gaussian process)
## FIXME: add horizontal grid lines
## FIXME: colour ribbons?

if (!file.exists("nserc_logo.png")) {
    download.file("https://nserc-crsng.canada.ca/sites/default/files/2025-12/NSERC_DIGITAL.zip",
                  dest = "nserc_logo.zip")
    unzip("nserc_logo.zip")
    file.copy("NSERC_DIGITAL/png/NSERC_RGB.png", "nserc_logo.png")
}
gam_tp <- gam(as.numeric(date) ~ s(year, bs = "tp", k=8), data = df)
gam_cs <- gam(as.numeric(date) ~ s(year, bs = "cs", k=8), data = df)
## cheating a bit to get a Gaussian process that's fairly quickly mean-reverting
## m = c(-2, 4, 1.5) -> power-exponential with scale 4 and power 1.5, i.e.
##  correlation \propto (d/4)^(1.5)
gam_gp <- gam(as.numeric(date) ~ s(year, bs = "gp", k=8, m=c(-2, 4, 1.5)), data = df)

pframe0 <- data.frame(year = seq(2013, 2026, by = 0.1))
pfun <- function(x) {
  pp0 <- predict(x, newdata = pframe0, se.fit = TRUE)
  with(pp0, data.frame(pframe0, date = fit, lwr = fit-2*se.fit, upr = fit+2*se.fit)) |>
    mutate(across(c(date, lwr, upr), as.Date))
}

pframe <- purrr::map_dfr(list("thin plate" = gam_tp, "cubic spline" = gam_cs, "GP" = gam_gp),
                         pfun,
                         .id = "model")
gg0 <- (ggplot(df,aes(year,date))
  + expand_limits(x=this_year + 1)
  + geom_vline(xintercept=this_year, lty=2)
  + scale_x_continuous(breaks=seq(2013,this_year,by=2))
  + geom_point()
)
cvec <- palette.colors(4)[-1] ## Okabe-Ito without black
gg1 <- (gg0
  + geom_line(data = pframe, aes(colour = model), linewidth = 2)
  + geom_ribbon(data = pframe, aes(ymin = lwr, ymax = upr, fill = model, colour = model),
                linewidth = 0.5, alpha = 0.3)
  + scale_colour_manual(values = cvec)
  + scale_fill_manual(values = cvec)
)
print(gg1)
gp_pred <- filter(pframe, model == "GP", year == 2026)
fmt_str <- "%B %d"
pred_date <- with(gp_pred, format(c(date, lwr, upr), fmt_str))
pred_str <- sprintf("GP prediction: %s (%s - %s)", pred_date[1], pred_date[2], pred_date[3])

ymin <- lubridate::ymd(paste(this_year, "mar", "20"))
ymax <- lubridate::ymd(paste(this_year, "mar", "23"))
logo <- readPNG("nserc_logo.png")
gg2 <- gg1 + annotation_raster(logo, xmin= 2019.5, xmax = 2021.5,
                               ymin = ymin, ymax = ymax) +
  annotate(geom = "label", x = 2019, y = lubridate::ymd(paste(this_year, "mar", "26")),
           label = pred_str)
                                
print(gg2)

ggsave("nserc_dates.png")

## setup for crazy range of polynomial fits
yvec <- seq(2013, this_year,length=51)
predvals <- purrr::map_dfr(1:7,
                           ~data.frame(year=yvec,
                                       date=predict(lm(date~poly(as.numeric(year),.),data=df),
                                                    newdata=data.frame(year=yvec))),
                           .id="order") %>%
    mutate(across(date, ~as.Date(.,origin=as.Date("1970-01-01"))))

gg1 <- gg0 + geom_line(data=predvals,aes(colour=factor(order)))

print(gg1)

ggsave("nserc_dates_poly.png")
predvals %>% filter(order %in% 1:2,year==2021)