Skip to content

Instantly share code, notes, and snippets.

@seabbs
Last active February 7, 2022 16:02
Show Gist options
  • Save seabbs/fb1bc9c79c3dd7117f9314cb97e71615 to your computer and use it in GitHub Desktop.
Save seabbs/fb1bc9c79c3dd7117f9314cb97e71615 to your computer and use it in GitHub Desktop.
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 …
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)
@seabbs
Copy link
Author

seabbs commented Jan 27, 2022

deterministic-deconvolution-sim

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment