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

Here's a version using a bit more EpiNow2 functionality and the nonmechanistic model (this is based on the PR version at epiforecasts/EpiNow2#504 which can be installed by using remotes::install_github("epiforecasts/EpiNow2@dist-interface") and making a cup of coffee).

# 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)
library(dplyr)

# 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 distribution (directly from Aloon)
incubation_period <- LogNormal(meanlog = 1.63, sdlog = 0.5, max = 20)
## Set onset to death period (form Verity: mean 17.8, CV 0.45 -> sd 8.01)
onset_to_death_period <- Gamma(mean = 17.8, sd = 8, max = 30)
## combine
infection_to_death <- incubation_period + onset_to_death_period

# Plot delay distributions
plot(infection_to_death)

# Plot convolved delay distribution
plot(collapse(discretise(infection_to_death)))

# Make data.frame of infections for EpiNow2
infections_data <- data.frame(
  date = as.Date("2020-03-01") + 1:n_inf,
  cases = data_infections
)

# Make data.frame of reported cases by applying delay + noise
case_data <- adjust_infection_to_report(
  infections_data,
  delay_defs = infection_to_death,
  reporting_model = function(n) rpois(n, lambda = n)
)
setnames(case_data, "cases", "confirm")

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 <- LogNormal(
  mean = covid_serialint$summary_stats$mean,
  sd = covid_serialint$summary_stats$sd,
  max = covid_serialint_discrete_max
)

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

Thanks, some nice streamlining functions there (adjust_infection_to_report seems handy in particular). Could be good opportunity to surface a few of them for more routine usage elsewhere.

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