In bbolker/McMasterPandemic: Pandemic Model

knitr::opts_chunk$set(echo=FALSE)

library(tidyverse)
library(glmmTMB)
library(broom.mixed)
library(McMasterPandemic)
library(ggplot2); theme_set(theme_bw())
library(colorspace)
library(TMB)
## make nice defaults
scale_colour_discrete <- function(...) {
  colorspace::scale_colour_discrete_qualitative(...)
}
scale_fill_discrete <- function(...) {
  colorspace::scale_fill_discrete_qualitative(...)
}

goals

explore phenomenological fits to epidemic time-series data
cumulative incidence (wrong?) vs incidence vs
different functional forms: Gaussian (IHME/probit) vs logit vs Richards (vs SIR)?
choice of response distribution (negative binomial 1 vs 2, log-Normal, least-squares, ... ?)
EpiEstim, EpiEstim plus priors?

source("ontario_clean.R")

(ggplot(ont_recent,aes(x=vday,y=value,colour=var))
    + geom_point()
    + scale_y_log10()
    + geom_smooth(method="lm",
                  formula=y~poly(x,2),
                  se=FALSE)
    + scale_x_continuous(limits=c(-2,NA))
)

compile("epigrowthfit_TMB.cpp")
## fit linear, quad mixed? models
## logit, Richards?
## epigrowthfit or mle2 or ... ?
## plot predictions and CIs
dyn.load("epigrowthfit_TMB.so")
## data; include t=-1 as initial time (for cumulated incidence for first step)
tmb_data <- with(filter(ont_recent,
                        vday>=0,var=="Hospitalization"),
                 list(t=c(-1,vday),x=value))
m <- MakeADFun(data=tmb_data
             , parameters=list(log_thalf=5,
                               log_K=5,
                               log_r=-1,
                               log_nb_disp=0,
                               log_p=0)
             , DLL="epigrowthfit_TMB")
m$fn(m$par)
fit1 <- with(m,nlminb(start=par,
                      objective=fn,
                      gradient=gr))
exp(setNames(fit1$par,gsub("log_","",names(fit1$par))))
m$report()$inccurve