Simulated infection recovery

epinow2_infection.R

# Example code for simulation recovery of infection dynamics from deaths for COVID -----------------------------------------

# Load libraries:

# install.packages("EpiNow2", repos = "https://epiforecasts.r-universe.dev")
library(EpiNow2)

# Simulate data and delays -----------------------------------------

# Simulate infections with a sharp drop
data_infections <- c(seq(220,1000,20),rev(seq(200,1000,200)),rep(200,25))
n_inf <- length(data_infections) # number of days to consider

# Set infection-to-death function pmf (mean from McAloon + Verity)
mean_p <- 5 + 18
scale_p <- 3
shift_p <- 0
p_by_day <- function(x){dgamma(x,shape=mean_p/scale_p,scale=scale_p)}

# Convert distribution parameters:
mean1 <- convert_to_logmean(mean_p,sqrt(scale_p^2*mean_p/scale_p))
sd1 <- convert_to_logsd(mean_p,sqrt(scale_p^2*mean_p/scale_p))

# Plot infection-to-death delay function
max_days <- 70
plot(0:max_days,p_by_day(0:max_days))

# OBSERVATION: would be more elegant to define these via epiparameter::epi_dist() - have noted issue with current quickstart

# Define transition matrix to simulate outcome data
f_matrix <- matrix(0,nrow=n_inf,ncol=n_inf)
n_delay_days <- 30 # maximum delay period to consider

# Convolution function
for(ii in 1:n_inf){
  i_max <- min(ii+n_delay_days-1,n_inf)
  j_max <- min(n_inf-ii+1,n_delay_days)
  f_matrix[ii:i_max,ii] <- p_by_day(0:(j_max-1)) # fill matrix entries
  
}

data_outcomes <- f_matrix %*% data_infections
data_outcomes_noise <- rpois(length(data_outcomes),lambda = data_outcomes) # add some noise

# OBSERVATION: probably a quicker way to do above step using Epinow2::forecast_secondary()

# Make data.frame of cases for EpiNow2
case_data <- data.frame(date = as.Date("2020-03-01")+1:length(data_outcomes_noise), confirm = data_outcomes_noise)

# Version with EpiNow --------------------------------------------------------------

covid_serialint <- 
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = T
  )

covid_serialint_discrete <- 
  epiparameter::discretise(covid_serialint)

covid_serialint_discrete_max <- 
  covid_serialint_discrete$prob_dist$q(p = 0.999)

serial_interval_covid <- 
  dist_spec(
    mean = covid_serialint$summary_stats$mean,
    sd = covid_serialint$summary_stats$sd,
    max = covid_serialint_discrete_max,
    distribution = "lognormal"
  )

# Delay infection to death
incubation_time_covid <- dist_spec(
  mean = mean1,
  sd = sd1,
  max = 50,
  distribution = "lognormal"
)

# OBSERVATION: would be handy to have a wrapper than can format these from 
# epiparameter::epi_dist() stored distributions in a line or two of code

epinow_estimates <- epinow(
  # cases
  reported_cases = case_data,
  # delays
  generation_time = generation_time_opts(serial_interval_covid),
  # delays
  delays = delay_opts(incubation_time_covid),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

base::plot(epinow_estimates)

infection_estimates <- epinow_estimates$estimates$summarised |> filter(variable=="infections")

# Plot comparison
par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,1))
plot(case_data$date,data_infections,ylim=c(0,1e3),xlab="time",ylab="events",type="l",lwd=2)
lines(case_data$date,case_data$confirm,col="red",lwd=2)
lines(case_data$date,data_outcomes,col="red",lwd=1,lty=2)

lines(infection_estimates$date,infection_estimates$median,lwd=2,col="dark blue")
polygon(c(infection_estimates$date,rev(infection_estimates$date)),c(infection_estimates$lower_90,rev(infection_estimates$upper_90)),
        col=rgb(0,0,1,0.1),border=NA)

dev.copy(png,paste0("delay_recover.png"),units="cm",width=15,height=10,res=150)
dev.off()

Yes, good idea. This is a nice case study.

noting that in these lines of this reprex https://gist.github.com/adamkucharski/11227d7ccca4b8e88e3aa9d1f4cd6963#file-epinow2_infection-r-L67-L72 we propagated an misspecification of dist_spect() for "lognormal" distributions from epiverse-trace/howto#34 which the @dist-interface prevents to and is most clearly described in epiforecasts/EpiNow2#581