projection of young adult full vax rates in ON

youngadult_ON_vaxproj.R

##https://www.dataquest.io/blog/r-api-tutorial/
## https://tonyelhabr.rbind.io/post/nested-json-to-tidy-data-frame-r/
## https://api.covid19tracker.ca/docs/1.0/vaccine-age-groups
library(httr)
library(jsonlite)
library(tibble)
library(purrr)
library(ggplot2)
library(dplyr)
library(nlme)


url <- "https://api.covid19tracker.ca/vaccines/age-groups/province/on"
h <- httr::GET(url)
ugly <- (h$content
    |> rawToChar()
    |> fromJSON(simplifyVector=FALSE)
)

## wrangle data
getf <- function(x) {
    tibble(date = x$date,
           vax = fromJSON(x$data)[["18-29"]][["full"]])
}
res <- (purrr::map_dfr(ugly$data, getf)
    %>% mutate(across(date, as.Date))
    %>% mutate(vaxprop = vax/tail(vax,1)*38.24)
)

## chop off initial time window (poor fit to 4-param logistic)
res2 <- (res
    %>% filter(date>as.Date("2021-04-01"))
    %>% mutate(t=as.numeric(date-min(date)))
)


## fit four-parameter logistic
n1 <- nls(vaxprop ~ SSfpl(t, v_start, v_end, d_mid, scal),
    data = res2)
pp <- data.frame(t=0:200)
pp$vaxprop <- predict(n1, newdata=pp)
pp$date <- res2$date[1] + pp$t


## compute confidence intervals by pseudo-bootstrap sampling
## posterior prediction intervals: assume sampling dist of params is MVN
psamp <- MASS::mvrnorm(1000, mu = coef(n1), Sigma = vcov(n1))
pvals <- apply(psamp,1,
               function(x) with(as.list(x),
          v_start + (v_end-v_start)*(1/(1+exp(-(pp$t-d_mid)/scal)))))
penv <- t(apply(pvals, 1, quantile, c(0.025, 0.975)))
pp$lwr <- penv[,1]
pp$upr <- penv[,2]

## plot
theme_set(theme_bw())
ggplot(res2, aes(date, vaxprop)) + geom_point() +
    geom_line(data=pp, colour = "red") +
    geom_ribbon(data=pp, aes(ymin =lwr, ymax = upr),
                fill = "red", alpha=0.2, colour = NA) +
    labs(y="predicted/observed proportion fully vax\n(18-29 year olds)")