A script containing the functions required to simulation infections from a Poisson process, calculate prevalence assuming a deterministic convolution with a vector of detection probabilities, deconvolve prevalence to infections by numerically estimating the inverse of the convolution matrix, and then plotting a simulation of this process. *Note …

deterministic-deconvolution.R

library(data.table)
library(snakecase)
library(ggplot2)

# make a convolution matrix
convolution_matrix <- function(vec1, vec2) {
  lvec1 <- length(vec1)
  lvec2 <- length(vec2)
  conv <- matrix(0, nrow = lvec1, ncol = lvec1)
  for (s in 1:lvec1) {
    l <- max(1, s - lvec2 + 1)
    conv[s, l:s] <- tail(rev(vec2), min(s, lvec2))
  }
  return(conv)
}

# simulate infections from a Poisson process
simulate_infections <- function(
  init = 100,
  growth = c(
    rep(0.05, 50), rep(0.01, 50), rep(-0.05, 100), rep(0.2, 25),
    rep(0.15, 25), rep(-0.1, 25), rep(-0.15, 25)
  )
) {
    cases <- rep(0, length(growth))
    cases[1] <- rpois(1, init * exp(growth[1]))
    for (i in seq_along(growth[-1])) {
      cases[i + 1] <- rpois(1, cases[i] * exp(growth[i + 1]))
    }
  return(cases)
}

# simulate prevalence
simulate_prevalence <- function(infs, dist = c(1, 0.6, 0.4, 0.2, 0.1)) {
  conv_mat <- convolution_matrix(infs, dist)
  prev <- as.vector(conv_mat %*% infs)
  return(prev)
}

# deconvolve prevalence by solving for the inverse
# this is unstable when initial dist entries are small and so
# here conditioning on postivity
deconvolve_prevalence <- function(prev, dist = c(1, 0.6, 0.4, 0.2, 0.1)) {
  conv_mat <- convolution_matrix(prev, dist)
  est_infs <- solve(conv_mat, prev)
  return(est_infs)
}

# Simulate infections -> convolve to prevalence
# solve for the inverse of the convolution matrix + prevalence
infections_to_deconvolved_infections <- function(
  init = 100,
  growth = c(
    rep(0.05, 50), rep(0.01, 50), rep(-0.05, 100), rep(0.2, 25),
    rep(0.15, 25), rep(-0.1, 25), rep(-0.15, 25)
  ),
  dist = c(1, 0.6, 0.4, 0.2, 0.1, 0.05, 0.01)
) {
  infs <- simulate_infections(init = init, growth = growth)
  conv_mat <- convolution_matrix(infs, dist)
  prev <- as.vector(conv_mat %*% infs)
  est_infs <- solve(conv_mat, prev)
  sims <- data.table::data.table(
      time = 1:length(infs),
      infections = infs,
      prevalence = prev,
      estimated_infections = est_infs
  )
  return(sims[])
}

# melt simulations and give nice names
melt_simulations <- function(sims) {
  sims <- melt(
    sims,
    id.vars = "time",
    value.name = "value",
    variable.name = "metric"
  )
  sims[,
   metric := snakecase::to_any_case(as.character(metric), case = "sentence")]
  return(sims[])
}

# Plot simulations
plot_infection_simulations <- function(sims) {
  ggplot(sims) +
    aes(x = time, y = value, col = metric) +
    geom_point(alpha = 0.5) +
    geom_line() +
    scale_y_continuous(trans = scales::log_trans()) +
    scale_color_brewer(palette = "Dark2") +
    theme_bw() +
    theme(legend.position = "bottom") +
    labs(color = "Metric", y = "Value (log)", x = "Time")
}

infections_to_deconvolved_infections() |>
  melt_simulations() |>
  plot_infection_simulations()

ggsave("deterministic-deconvolution-sim.png", width = 8, height = 8)