needed_packages <- c(
"deSolve",
"tibble",
"purrr",
"dplyr",
"magrittr",
"RColorBrewer"
)TPT model
All the data and source code files are here.
You can report errors or ask for additional analysis here.
1 Introduction
2 Packages
Required packages:
Installing the missing packages:
to_install <- needed_packages[! needed_packages %in% installed.packages()[,"Package"]]
if (length(to_install)) install.packages(to_install)Loading some packages:
library(tibble)
library(purrr)
library(dplyr)
library(magrittr)3 Utilitary functions
Default line width:
lwd_val <- 2Tuning some functions:
ode2 <- function(...) tibble::as_tibble(as.data.frame(deSolve::ode(...)))
seq2 <- function(...) seq(..., le = 512)
plotl <- function(...) plot(..., type = "l", lwd = 2)
lines2 <- function(...) lines(..., lwd = 2)
text2 <- function(...) text(..., adj = 1)
mtext2 <- function(...) mtext(..., side = 2, line = 2.25)
segments2 <- function(...) segments(..., lwd = 2)
abline4 <- function(...) abline(..., lwd = lwd_val)
abline2 <- function(..., col = 2) abline4(..., col = col)
abline3 <- function(...) abline(..., col = "grey")
plot2 <- function(..., col = 4) plot(..., type = "l", lwd = lwd_val, col = col)
plot3 <- function(...) plot2(..., xlab = "actual duration of treatment (days)")
legend2 <- function(...) legend(..., lty = 1, lwd = lwd_val, bty = "n")
legend3 <- function(x, y, legend, ...) {
tmp <- legend2(x, y, legend = rep(" ", length(legend)),
text.width = max(strwidth(legend)), xjust = 1, yjust = 1, ...)
text(tmp$rect$left + tmp$rect$w, tmp$text$y, legend, pos = 2)
}
legend4 <- function(x, y, legend, ...) {
cex_val <- .8
tmp <- legend(x, y, legend = rep(" ", length(legend)), bty = "n", cex = cex_val,
text.width = max(strwidth(legend)), xjust = 1, yjust = 1, lwd = 2, ...)
text(tmp$rect$left + tmp$rect$w, tmp$text$y, legend, pos = 2, cex = cex_val)
}
legend5 <- function(...) legend4(..., lty = 1)
mclapply2 <- function(...) {
parallel::mclapply(..., mc.cores = parallel::detectCores() - 1)
}
polygon2 <- function(x, y, col = 4, alpha = .2, ...) {
polygon(c(x[1], x, tail(x, 1)), c(0, y, 0), col = adjustcolor(col, alpha),
border = NA)
}
polygon3 <- function(x, y1, y2, col = 4, ...) {
polygon(c(x, rev(x)), c(y1, rev(y2)), col = col, border = NA, ...)
}Color palette functions:
rgb2 <- function(...) rgb(..., maxColorValue = 255)
spectral <- function(n) {
colorRampPalette(rev(RColorBrewer::brewer.pal(9, "Spectral")))(n)
}
yellows <- map_chr(list(
list(255, 241, 195),
list(246, 220, 138),
list(235, 200, 79),
list(204, 159, 44),
list(156, 120, 24),
list(107, 85, 19)), ~ do.call(rgb2, .x))
reds <- map_chr(list(
list(250, 208, 206),
list(237, 164, 167),
list(222, 104, 102),
list(201, 62, 63),
list(156, 39, 38),
list(116, 24, 21)), ~ do.call(rgb2, .x))
blues <- map_chr(list(
list(200, 231, 251),
list(157, 202, 235),
list( 87, 153, 209),
list( 1, 113, 178),
list( 5, 78, 145),
list( 3, 44, 92)), ~ do.call(rgb2, .x))
greens <- map_chr(list(
list(217, 232, 197),
list(160, 201, 121),
list( 98, 179, 70),
list( 69, 148, 51),
list( 35, 113, 48),
list( 21, 59, 28)), ~ do.call(rgb2, .x))
tones <- map_chr(list(
list(247, 245, 239),
list(228, 224, 204),
list(198, 194, 166),
list(168, 163, 134),
list(136, 131, 101),
list( 98, 94, 74)), ~ do.call(rgb2, .x))
greys <- map_chr(list(
list(230, 230, 230),
list(200, 206, 218),
list(153, 163, 181),
list(111, 123, 144),
list( 73, 86, 109),
list( 35, 50, 72)), ~ do.call(rgb2, .x))
skin_tones <- map_chr(list(
list(248, 230, 215),
list(221, 189, 163),
list(191, 153, 125),
list(147, 105, 87),
list(117, 79, 64),
list( 68, 45, 31)), ~ do.call(rgb2, .x))
# template:
#<- map_chr(list(
# list(, , ),
# list(, , ),
# list(, , ),
# list(, , ),
# list(, , ),
# list(, , )), ~ do.call(rgb2, .x))A function that converts a percentage decrease per year into a rate:
find_rate <- function(p = .1, b1 = .5 * p, b2 = 1.5 * p) {
optimise(function(x) abs(exp(- x) - (1 - p)), c(b1, b2))$min
}4 Prophylaxis
4.1 The Hill equation
In what follows we will model many phenomena phenomenologically using the 3-parameter monotonically increasing Hill equation, the general equation of which is
\[ y = Y\frac{x^h}{X^h + x^h} \]
where \(Y \in \mathbb{R}^+\) is the maximum value that \(y\) can take, \(X \in \mathbb{R}^+\) is the value of \(x\) at which \(y\) reaches half of its maximum value \(Y\) and \(h \in \mathbb{R}^+\) is the Hill coefficient that controls the shape of the relationship with S-shape for \(0 < h < 1\) and a simple saturating shape when \(h \ge 1\).
hill <- function(x, X, Y, h) {
x2h <- x^h
Y * x2h / (X^h + x2h)
}Let’s illustrate the properties of this equation:
xs <- seq2(0, 40)
X <- 15
Y <- 1.4
plot(NA, xlim = c(0, max(xs)), ylim = c(0, 1.5), xlab = "x values", ylab = "y values")
abline2(v = X, col = "grey")
abline2(h = Y, col = "grey")
hs <- exp(seq(-3, 3))
cols <- spectral(length(hs))
walk2(hs, cols, ~ lines2(xs, hill(xs, X, Y, h = .x), col = .y))
legend3("bottomright", col = cols, title = "Hill coefficient h",
legend = formatC(round(hs, 2), digits = 2, format = "f"))
axis(1, X, "X")
axis(2, Y, "Y", las = 1)
4.2 Treatment uptake
We can model the probability \(\upsilon(d)\) of treatment uptake as a function of the proposed treatment duration \(d\) using a Hill equation as so:
\[ \upsilon(d | X_\upsilon, Y_\upsilon, Z_\upsilon, h_\upsilon) = Z_\upsilon - Y_\upsilon\frac{d^{h_\upsilon}}{X^{h_\upsilon} + d^{h_\upsilon}} \]
where \((Y_\upsilon, Z_\upsilon) \in [0,1]^2\). In R:
uptake <- function(d, Xv, Yv, Zv, hv) {
Zv - hill(d, Xv, Yv, hv)
}Two examples:
ds <- seq2(0, 30)
plot2(ds, uptake(ds, 15, 1, 1, 4), ylim = 0:1,
xlab = "duration of treatment (days)", ylab = "uptake probability")
lines2(ds, uptake(ds, 10, .7, .9, 1.75), col = 3)
4.3 Treatment adherence
4.3.1 Conditional dropping probability
Call \(p(t)\) the probability that a person who, at day \(t\), has not stopped her treatment yet, decides to stop it. This probability could be constant and independent of the duration the person has been in the treatment. On the contrary, we could imagine that this probability depends on the duration the person has been in the treatment. And this dependency could result from 2 effects:
- a fatigue effect whereby \(p(t)\) would increase with \(t\);
- a finish-line effect whereby \(p(t)\) would decreases nearby the end of treatment.
In order to account for all these possibilities, we model the fatigue effect as a probability \(f(t)\) increasing with time \(t\) being in the treatment according to a Hill equation:
\[ f(t|X_f, Y_f, h_f) = Y_f\frac{t^{h_f}}{X_f^{h_f} + t^{h_f}} \]
with \(0 \leq Y_f \leq 1\). Note too here that \(f(0)=0\). The finish-line effect is modelled as a decreasing mitigation weight \(m(t) \in [0, 1]\) of \(f(t)\), and it is also modelled with a Hill equation as so:
\[ m(t|X_m, Y_m, h_m) = 1 - Y_m\frac{t^{h_m}}{X_m^{h_m} + t^{h_m}} \]
where \(0 \leq Y_m \leq 1\). Note too that \(m(0) = 1\). From here the conditional probability \(p(t)\) reads
\[ p(t|\theta_m, \theta_f) = m(t|\theta_m) \times f(t|\theta_f) \]
Let’s look at some examples using the following function that plots these probabilities:
plot_probabilities <- function(
pd = 40, # proposed duration of treatment
Xf = 25, Yf = .8, hf = 7, # fatigue
Xmp = .875, Ym = 1, hm = 50, # mitigation of fatigue
lwd = 2,
ylim = 0:1,
add = FALSE,
lgnd = TRUE) { # whether to add the legend to the plot or not
Xm <- Xmp * pd
ts <- seq2(0, pd)
fatigue <- hill(ts, Xf, Yf, hf)
mitigation <- 1 - hill(ts, Xm, Ym, hm)
if (! add) plot(ts, fatigue, ylim = ylim, col = 3, type = "l", lwd = lwd,
xlab = "time in treatment (days)",
ylab = "probability or mitigation factor")
lines2(ts, mitigation, col = 2)
lines2(ts, mitigation * fatigue, col = 4)
if (lgnd & ! add)
legend2("left", col = c(3, 2, 4),
legend = c("fatigue", "mitigation of fatigue", "dropping probability"))
}Note here the pd parameter that is the proposed duration of the treatment. Note too that the \(X_m\) parameter is expressed as a proportion of pd. Executing the function:
plot_probabilities()
Note that this model accounts for the special case where \(p(t)\) is a constant independent of time:
plot_probabilities(hf = 0, hm = 0, Ym = .05, Yf = .4)
4.3.2 Duration of treatment distribution
From \(p(t)\) we can express the probability that somebody stops her treatment exactly at time \(t\) as
\[ P(t|\theta_m, \theta_f) = p(t|\theta_m, \theta_f)\prod_{x=0}^{{}^{-}t}(1 - p(x|\theta_m, \theta_f))^{dx} \]
The density of probability of the actual duration of treatment (aka effective duration of treatment) can then be expressed by:
\[ \varphi(t|d, \theta_m, \theta_f) = \frac{P(t|\theta_m, \theta_f)} {\int_{0}^d P(x|\theta_m, \theta_f)dx} \]
where \(d\) is the proposed duration of treatment. Let’s see this in R. First a function that computes the dropping probability:
dropping_probability <- function(
pd = 40, # proposed duration of treatment
Xf = 25, Yf = .8, hf = 7, # fatigue
Xmp = .875, Ym = 1, hm = 50, # mitigation of fatigue
by = pd / (le - 1), le = 512) { # sampling day in the treatment
ts <- seq(0, pd, by)
tibble(x = ts,
y = hill(ts, Xf, Yf, hf) * (1 - hill(ts, Xmp * pd, Ym, hm)))
}Note the parameterization that is similar to the parameterization of the plot_probabilities() function. This is what it gives:
dropping_probability()# A tibble: 512 × 2
x y
<dbl> <dbl>
1 0 0
2 0.0783 2.36e-18
3 0.157 3.02e-16
4 0.235 5.16e-15
5 0.313 3.87e-14
6 0.391 1.84e-13
7 0.470 6.61e-13
8 0.548 1.94e-12
9 0.626 4.95e-12
10 0.705 1.13e-11
# ℹ 502 more rowswhere the first column is the time \(t\) in treatment and the second column the conditional density of probability \(p(t)\) of dropping. This function is used by the following one to compute the probability distribution of the effective duration of treatment:
effective_duration <- function(
pd = 40, # proposed duration of treatment
Xf = 25, Yf = .8, hf = 7, # fatigue
Xmp = .875, Ym = 1, hm = 50, # mitigation of fatigue
by = pd / (le - 1), le = 512) { # sampling effective durations of treatment
dp <- dropping_probability(pd, Xf, Yf, hf, Xmp, Ym, hm, by, le)
y <- dp$y
dens <- c(y[1], y[-1] * cumprod(1 - y[-length(y)]))
tibble(x = dp$x,
y = dens / sum(dens))
}Note again the same parameterization as for the plot_probabilities() and dropping_probability() functions. Let’s try it:
with(effective_duration(), {
plot3(x, y, ylab = "probability density")
polygon2(x, y)
})
The following function (with the same parameterization as the 3 previous ones) helps visualizing the relationship between the proposed duration of treatment (vertical grey lines) and the distribution of the actual duration of treatment (in blue):
explore_dropping <- function(
Xf = 200, Yf = .15, hf = 7,
Xmp = .875, Ym = 1, hm = 50,
by = .1, ylim = c(0, .2),
proposed_durations = c(20, 40, seq(60, 100, 10))) {
proposed_durations <- sort(proposed_durations)
eds <- map(proposed_durations,
~ effective_duration(.x, Xf, Yf, hf, Xmp, Ym, hm, by))
opar <- par(mfrow = 2:1, plt = c(.1, .96, .21, .93))
plot_probabilities(tail(proposed_durations, 1),
Xf, Yf, hf, Xmp, Ym, hm, lwd = 4, ylim = ylim, lgnd = FALSE)
walk(head(proposed_durations, -1),
~ plot_probabilities(.x, Xf, Yf, hf, Xmp, Ym, hm, add = TRUE))
plot(NA, xlim = c(0, max(proposed_durations)),
ylim = c(0, max(unlist(map(eds, ~ .x$y)))),
xlab = "effective duration of treatment (days)",
ylab = "probability density")
walk(eds, ~ with(.x, {
lines2(.x$x, .x$y, col = 4)
polygon2(x, y)
}))
abline2(v = proposed_durations, col = "grey")
par(opar)
}Let’s play with it:
explore_dropping(Xmp = .95, hm = 50, ylim = c(0, .0007))
A function to draw durations from the distribution:
draw_durations <- function(
pd = 40, # proposed duration of treatment
Xf = 25, Yf = .8, hf = 7, # fatigue
Xmp = .875, Ym = 1, hm = 50, # mitigation of fatigue
by = pd / (le - 1), le = 512, # sampling effective durations of treatment
n = 1000) { # number to draw
with(effective_duration(pd, Xf, Yf, hf, Xmp, Ym, hm, by, le), sample(x, n, TRUE, y))
} Let’s try it:
draw_durations(n = 16) [1] 15.73386 16.82975 17.29941 12.83757 18.70841 16.51663 15.57730 12.68102
[9] 14.09002 14.87280 14.16830 13.30724 16.98630 17.92564 12.28963 15.342474.4 Treatment efficacy
We can assume that the treatment efficacy \(\varepsilon(\phi)\) as a function of the actual duration \(\phi\) of treatment also follows a Hill equation:
\[ \varepsilon(\phi|X_\varepsilon, Y_\varepsilon, h_\varepsilon) = Y_\varepsilon\frac{\phi^{h_\varepsilon}} {X_\varepsilon^{h_\varepsilon} + \phi^{h_\varepsilon}} \]
Let’s see 2 examples:
ts <- seq2(0, 30)
plot3(ts, hill(ts, 5, 1, .8), ylim = 0:1, ylab = "treatment efficacy")
lines2(ts, hill(ts, 20, 1, 5), col = 3)
legend2("topleft", legend = c("optimistic", "pessimistic"), col = 4:3)
4.5 Full model
The proportion \(\pi\) of recent infections that can be sterilized then reads:
\[ \pi(d|\tau, \theta_\upsilon, \theta_m, \theta_f, \theta_\varepsilon) = \tau \times \upsilon(d| \theta_\upsilon) \times \varepsilon(\phi(d|\theta_m,\theta_f)| \theta_\varepsilon) \]
where \(\tau\) is the proportion of recent infections that can be identified through contact tracing and \(d\) is the proposed duration of treatment. In R:
pi_value <- function(
tau = .99,
pd = 10,
Xv = 15, Yv = 1, hv = 4, Zv = 1,
Xf = 200, Yf = .15, hf = 7,
Xmp = .875, Ym = 1, hm = 50,
Xe = 10, Ye = 1, he = 1,
by = pd / (le - 1), le = 512,
n = 1000) {
tau *
uptake(pd, Xv, Yv, Zv, hv) *
hill(draw_durations(pd, Xf, Yf, hf, Xmp, Ym, hm, by, le, n), Xe, Ye, he)
}Note the similar parameterization as for functions above. Let’s use this function to look at the effect of \(\tau\) on \(\pi\):
tau_vals <- c(.2, .5, .75, .9, 1)
out1 <- tau_vals |>
map(~ pi_value(tau = .x, pd = 10, n = 1e4, le = 1e6)) |>
map(density)
plot2(NA, xlim = 0:1, ylim = c(0, max(unlist(map(out1, extract2, "y")))),
xlab = expression(paste("proportion ", pi, " of new infections sterilized")),
ylab = "density of probability")
walk2(out1, seq_along(tau_vals), ~ with(.x, {
polygon2(x, y, col = .y)
lines2(x, y, col = .y)
}))
add_legend_sterilization <- function(...) {
legend2("topright", legend = paste0(rev(100 * tau_vals), "%"),
col = rev(seq_along(tau_vals)),
title = expression(paste("perentage ", tau,
" of recent infections identified")))
}
add_legend_sterilization()
A function that plots a quantile distribution along the \(y\)-axis:
quantplot <- function(x, y, probs = c(.025, seq(.05, .95, .05), .975),
col = 4, alpha = .07, smooth = TRUE, ...) {
lowess2 <- function(y) lowess(x, y, ...)$y
y <- map_dfr(y, quantile, probs)
if (smooth) y <- map_dfc(y, lowess2)
nbcol <- ncol(y)
test <- nbcol %% 2
nb <- (nbcol - test) / 2 - 1
walk(0:nb,
~ polygon3(x, y[[1 + .x]], y[[nbcol - .x]], col = adjustcolor(col, alpha)))
if (test) lines2(x, y[[ceiling(nbcol / 2)]], col = col)
}Let’s generate the distribution of the proportion of new infections sterilized for various proposed durations of treatment (about 1’40”):
pd_vals <- seq(1, 40, le = 512)
out <- map(tau_vals,
~ mclapply2(pd_vals,
function(x) pi_value(pd = x, tau = .x, n = 1e4, le = 1e6)))This gives:
plot(NA, xlim = c(0, 45), ylim = 0:1,
xlab = "proposed duration d of treatment (days)",
ylab = expression(paste("proportion ", pi, " of new infections sterilized")))
walk2(out,
seq_along(out),
~ quantplot(pd_vals, .x, probs = seq(.025, .975, le = 75), col = .y,
alpha = .03, smooth = FALSE))
add_legend_sterilization()
A dashboard that gathers everything regarding prophylaxis:
dashboard_prt <- function(
max_dur = 40,
tau = .99,
Xv = 15, Yv = 1, hv = 4, Zv = 1,
Xf = 200, Yf = .15, hf = 7,
Xmp = .875, Ym = 1, hm = 50,
Xe = 10, Ye = 1, he = 1,
le = 512, # nb of pts along the x axes
npop = 1000, # nb of pts along the distribution to draw from
nspl = 1000, # nb of pts to draw from the distribution (nspl >> npop is better)
probs = seq(.025, .975, le = 75),
ylim1 = 0:1, ylim3 = 0:1, ylim4 = 0:1,
col1 = 4, col2 = 4, col3 = 4, col4 = 4, col5 = 4,
alpha = .1, smooth = TRUE, f1 = 2 / 3, f2 = .08, ...) {
min_dur <- .1
ds <- seq(min_dur, max_dur, le = le)
# 1. treatment uptake
plot2(ds, uptake(ds, Xv, Yv, Zv, hv), ylim = ylim1, col = col1,
xlab = "proposed duration of treatment (days)", ylab = "uptake probability")
# 2. distribution of actual durations of treatment
proposed_durations <- seq(10, 40, 10) #pretty(2:max_dur)
eds <- map(proposed_durations,
~ effective_duration(.x, Xf, Yf, hf, Xmp, Ym, hm, by = .1))
plot(NA, xlim = c(0, max_dur),
ylim = c(0, max(unlist(map(eds, ~ .x$y)))),
xlab = "actual duration of treatment (days)",
ylab = "probability density")
walk(eds, ~ with(.x, {
lines2(.x$x, .x$y, col = col2)
polygon2(x, y)
}))
abline2(v = proposed_durations, col = "grey")
# 3. durations of treatment vs proposed duration of treatment
plot(NA, xlim = c(0, max_dur),
ylim = c(0, max_dur),
xlab = "proposed duration of treatment (days)",
ylab = "of treatment (days)")
mtext("actual duration", 2, 2.25)
eds <- mclapply2(ds,
function(x)
draw_durations(x, Xf, Yf, hf, Xmp, Ym, hm, le = npop, n = nspl))
quantplot(ds, eds, probs, col = col3, alpha, smooth, f = f1, ...)
abline2(0, 1, col = "grey")
# 4. pharmaco-dynamics treatment efficacy
plot3(ds, hill(ds, Xe, Ye, he), ylim = ylim3, col = col4, ylab = "treatment efficacy")
# 5. intervention efficacy
plot(NA, xlim = c(0, 45), ylim = ylim4,
xlab = "proposed duration of treatment (days)",
ylab = "infections sterilized")
mtext("proportion of new", 2, 2.25)
out <- mclapply2(ds,
function(x) pi_value(tau, x, Xv, Yv, hv, Zv, Xf, Yf, hf, Xmp, Ym,
hm, Xe, Ye, he, le = npop, n = nspl))
quantplot(ds, out, probs, col = col5, alpha, smooth, f = f2, ...)
lines2(ds, 1 - map_dbl(out, ~ mean(1 - .x)), col = 2, lty = 3)
invisible(list(d = ds, pi = out))
}Let’s try it:
opar <- par(mfrow = c(3, 2), cex = 1, plt = c(.25, .95, .25, .9))
dashboard_prt(probs = seq(.001, .999, le = 9), alpha = .2, npop = 1e4, nspl = 1e3,
col2 = 4, col4 = 4, ylim4 = c(0, .5))
par(opar)
5 Epidemiology
5.1 A minimalist SID
5.1.1 Dynamics
\[ \begin{align} \frac{dS}{dt} &= \gamma I + \alpha D - (1 - \pi) \beta DS \\ \frac{dI}{dt} &= (1 - \pi) \beta DS - (\delta + \gamma) I \\ \frac{dD}{dt} &= \delta I - \alpha D \end{align} \]
R code:
sid_dyn <- function(N, I0 = 0, D0, beta, gamma, pi, delta, alpha, times) {
ode2(c(S = N - I0 - D0, I = I0, D = D0),
times,
function(time, state, pars) {
with(as.list(c(state, pars)), {
infection <- (1 - pi) * beta * D * S
dS <- gamma * I + alpha * D - infection
dI <- infection - (delta + gamma) * I
dD <- delta * I - alpha * D
list(c(dS, dI, dD))
})
},
c(beta = beta, gamma = gamma, pi = pi, delta = delta, alpha = alpha))
}5.1.2 Equilibrium
\[ \begin{align} S^* &= N - I^* - D^* \\ I^* &= \frac{\alpha}{\delta} D^* \\ D^* &= \frac{\delta}{\alpha + \delta} N - \frac{\alpha}{\alpha + \delta} \frac{\delta + \gamma}{(1 - \pi)\beta} \end{align} \]
R code:
sid_equ <- function(N, beta, gamma, pi, delta, alpha) {
Dstar <- delta * N / (alpha + delta) - alpha * (delta + gamma) /
((alpha + delta) * (1 - pi) * beta)
Istar <- alpha * Dstar / delta
c(S = N - Istar - Dstar, I = Istar, D = Dstar)
}A function that plots the dynamics and equilibrium values of the infection prevalence and the disease incidence and prevalence:
plot_sid_dyn_equ <- function(
N = 1e5,
D0 = 350,
beta = .000001,
gamma = 0,
pi = 0,
delta = .0001,
alpha = 1 / 70,
times = seq2(0, 25000)) {
sims <- sid_dyn(N = N, I0 = 0, D0 = D0, beta = beta, gamma = gamma, pi = pi,
delta = delta, alpha = alpha, times = times)
equs <- sid_equ(N = N, beta = beta, gamma = gamma, pi = pi,
delta = delta, alpha = alpha)
plotl2 <- function(...) plotl(..., xlab = "time (year)")
opar <- par(mfrow = c(1, 3), cex = 1, plt = c(.25, .95, .25, .9))
with(sims, {
plotl2(time, I, col = 2, ylab = "infection prevalence")
abline(h = equs["I"], col = 2)
plotl2(time, D, col = 3, ylab = "disease prevalence")
abline(h = equs["D"], col = 3)
plotl2(time, delta * I, col = 4, ylab = "yearly disease incidence")
abline(h = delta * equs["I"], col = 4)
})
par(opar)
}Let’s do some verification:
plot_sid_dyn_equ()
5.1.3 Calibration
In absence of TPT (i.e. \(\pi = 0\)), let’s further assume that \(\gamma = 0\). Let’s call \(d\) the yearly disease incidence. Then,
\[ \begin{align} \delta &= \frac{d^*}{I^*} \\ \alpha &= \frac{D^*}{d^*} \\ \beta &= \frac{N}{\alpha} - \frac{\alpha + \delta}{\alpha \delta} D^* \end{align} \]
Corresponding R code:
sid_param_val <- function(inf_prev, dis_prev, dis_incd, N = 1e5) {
delta <- dis_incd / inf_prev
alpha <- dis_incd / dis_prev
c(alpha = alpha,
beta = alpha * delta / (delta * N - (alpha + delta) * dis_prev),
delta = delta)
}In Vietnam:
- prevalence of infection: somewhere between 10 and 30% of the population?
- prevalence of disease: 322 / 100,000
- yearly incidence of disease: 176 / 100,000
(sid_p_val <- sid_param_val(2e4, 322, 176)) alpha beta delta
5.465839e-01 6.859909e-06 8.800000e-03 Let’s see:
plot_sid_dyn_equ(D0 = 322,
beta = sid_p_val[["beta"]],
delta = sid_p_val[["delta"]],
alpha = sid_p_val[["alpha"]],
times = seq2(0, 6000))
5.2 A two-speed SIID
5.2.1 Dynamics
\[ \begin{align} \frac{dS}{dt} &= \gamma_s I_s + \gamma_f I_f + \alpha D - (1 - p\pi)\beta D S\\ \frac{dI_s}{dt} &= (1 - p) \beta D S - (\gamma_s + \delta_s) I_s \\ \frac{dI_f}{dt} &= (1-\pi) p \beta D S - (\gamma_f + \delta_f) I_f \\ \frac{dD}{dt} &= \delta_s I_s + \delta_f I_f - \alpha D \end{align} \]
with \(N = S + I_s + I_f + D\). R code:
siid_dyn <- function(N = 1e5, Is0 = 0, If0 = 0, D0,
pi, p, beta, gammas, gammaf, deltas, deltaf, alpha, times) {
ode2(c(S = N - Is0 - If0 - D0, Is = Is0, If = If0, D = D0),
times,
function(time, state, pars) {
with(as.list(c(state, pars)), {
infection <- beta * D * S
dS <- gammas * Is + gammaf * If + alpha * D - (1 - p * pi) * infection
dIs <- (1 - p) * infection - (gammas + deltas) * Is
dIf <- (1 - pi) * p * infection - (gammaf + deltaf) * If
dD <- deltas * Is + deltaf * If - alpha * D
list(c(dS, dIs, dIf, dD))
})
},
c(pi = pi, p = p, beta = beta, gammas = gammas, gammaf = gammaf,
deltas = deltas, deltaf = deltaf, alpha = alpha))
}5.2.2 Equilibrium
\[ \begin{align} A &= (\alpha\delta_f + [1 + \delta_s]\gamma_f)[1 - p] + (\alpha\delta_s + [\alpha + \delta_f]\gamma_s)[1 - \pi]p + [1 - p\pi]\delta_s\delta_f \\ D^* &= \frac{[1 - p]\delta_s\gamma_f + [1 - \pi]p\gamma_s\delta_f + [1 - p\pi]\delta_s\delta_f} {A} N - \frac{\delta_s\gamma_f + \gamma_s\delta_f + \delta_s\delta_f}{A} \frac{\alpha}{\beta} \\ I_f^* &= \frac{(N - [\alpha / \delta_s + 1]D^*)(1 - \pi)p\beta D^*} {\gamma_f + \delta_f + (1 - \delta_f / \delta_s)(1 - \pi)p\beta D^*} \\ I_s^* &= \frac{\alpha D^* - \delta_f I_f^*}{\delta_s} \\ S^* &= N - I_s^* - I_f^* - D^* \end{align} \]
R code:
siid_equ <- function(N = 1e5, pi, p, beta, gammas, gammaf, deltas, deltaf, alpha) {
A <- (alpha * deltaf + (1 + deltas) * gammaf) * (1 - p) +
(alpha * deltas + (alpha + deltaf) * gammas) * (1 - pi) * p +
(1 - p * pi) * deltas * deltaf
D <- (((1 - p) * deltas * gammaf +
(1 - pi) * p * gammas * deltaf +
(1 - p * pi) * deltas * deltaf) * N -
(deltas * gammaf + gammas * deltaf + deltas * deltaf) * alpha / beta) / A
If <- (N - (alpha / deltas + 1) * D) * (1 - pi) * p * beta * D /
(gammaf + deltaf + (1 - deltaf / deltas) * (1 - pi) * p * beta * D)
Is <- (alpha * D - deltaf * If) / deltas
c(S = N - Is - If - D,
Is = Is,
If = If,
D = D)
}A function that plots the dynamics and equilibrium values of the infection prevalence and the disease incidence and prevalence:
plot_siid_dyn_equ <- function(
N = 1e5,
Is0 = 0,
If0 = 0,
D0 = 350,
pi = 0,
p = .9,
beta = 1e-6,
gammas = 0,
gammaf = 0,
deltas = 1e-6,
deltaf = 1e-4,
alpha = 1 / 70,
times = seq2(0, 700000)) {
sims <- siid_dyn(N, Is0, If0, D0,
pi, p, beta, gammas, gammaf, deltas, deltaf, alpha, times)
equs <- siid_equ(N, pi, p, beta, gammas, gammaf, deltas, deltaf, alpha)
plotl2 <- function(...) plotl(..., xlab = "time (year)")
opar <- par(mfrow = c(2, 2), cex = 1, plt = c(.25, .95, .25, .9))
with(sims, {
incidence_dyn <- deltas * Is + deltaf * If
incidence_equ <- deltas * equs["Is"] + deltaf * equs["If"]
plotl2(time, Is + If, col = 1, ylab = "infection prevalence")
abline(h = equs["Is"] + equs["If"], col = 1)
plotl2(time, D, col = 2, ylab = "disease prevalence")
abline(h = equs["D"], col = 2)
plotl2(time, incidence_dyn, col = 3, ylab = "yearly disease incidence")
abline(h = incidence_equ, col = 3)
plotl2(time, deltaf * If / incidence_dyn, col = 4,
ylab = "proportion of fast infections")
abline(h = deltaf * equs["If"] / incidence_equ, col = 4)
})
par(opar)
}Let’s see:
plot_siid_dyn_equ()
Let’s verify that we get the same result as the SID when \(p = 1\):
N <- 1e5
beta <- .000001
gamma <- gammas <- gammaf <- 0
pi <- 0
delta <- deltas <- deltaf <- .0001
alpha <- 1 / 70
p <- 1
sid_equ(N, beta, gamma, pi, delta, alpha) S I D
14285.7143 85118.4565 595.8292 siid_equ(N, pi, p, beta, gammas, gammaf, deltas, deltaf, alpha) S Is If D
1.428571e+04 -1.243450e-10 8.511846e+04 5.958292e+02 5.2.3 Calibration
When \(\pi = 0\), we can re-arrange the equilibrium equations so that we have:
\[ \begin{align} D^* &= \frac{ ([\gamma_f + \delta_f]\delta_s + [\gamma_s \delta_f - \delta_s \gamma_f]p)\beta N - (\gamma_s \delta_f + [\gamma_f + \delta_f]\delta_s)\alpha}{ (\alpha + \delta_s)\delta_f + (1 - p)(1 + \delta_s)\gamma_f + ([\delta_f + \alpha]\gamma_s - [\delta_f - \delta_s]\alpha)p}\frac{1}{\beta} \\ I^* &= \frac{(N - D^*)(\delta_f - \delta_s)p\beta - (\gamma_f + \delta_f)\alpha} {(\delta_f - \delta_s)p\beta D^* - (\gamma_f + \delta_f)\delta_s}D^* \\ d_f^* &= \frac{(\delta_s N - (\alpha + \delta_s)D^*)\delta_f p\beta D^*} {(\gamma_f + \delta_f)\delta_s - (\delta_f - \delta_s)p\beta D^*} \\ d_s^* &= \left( \alpha + \frac{(\delta_s N - [\alpha + \delta_s] D^*)p\beta\delta_f} {(\delta_f - \delta_s)p\beta D^* - (\gamma_f + \delta_f)\delta_s} \right) D^* \end{align} \]
where \(I^*\) and \(D^*\) are the infection and disease prevalences and \(d_f^*\) and \(d_s^*\) are the disease incidences from recent and ancient infections respectively. In R, this is what it looks like:
siid_equ2 <- function(N = 1e5, p, beta, gammas, gammaf, deltas, deltaf, alpha) {
D <- (((gammaf + deltaf) * deltas +
(gammas * deltaf - deltas * gammaf) * p) * beta * N -
(gammas * deltaf + (gammaf + deltaf) * deltas) * alpha) /
(((alpha + deltas) * deltaf + (1 - p) * (1 + deltas) * gammaf +
((deltaf + alpha) * gammas - (deltaf - deltas) * alpha) * p) * beta)
c(D = D,
I = ((N - D) * (deltaf - deltas) * p * beta - (gammaf + deltaf) * alpha) * D /
((deltaf - deltas) * p * beta * D - (gammaf + deltaf) * deltas),
df = (deltas * N - (alpha + deltas) * D) * deltaf * p * beta * D /
((gammaf + deltaf) * deltas - (deltaf - deltas) * p * beta * D),
ds = D * (alpha + (deltas * N - (alpha + deltas) * D) * p * beta * deltaf /
((deltaf - deltas) * p * beta * D - (gammaf + deltaf) * deltas))
)
}Verification:
siid_equ(N = 1e5, pi = 0, p = .9, beta = .00001, gammas = 1 / 70, gammaf = 0,
deltas = 1 / 50, deltaf = 1 / 5, alpha = 1 / 70) S Is If D
1490.683 3856.629 5950.227 88702.461 siid_equ2(N = 1e5, p = .9, beta = .00001, gammas = 1 / 70, gammaf = 0,
deltas = 1 / 50, deltaf = 1 / 5, alpha = 1 / 70) D I df ds
88702.46085 9806.85592 1190.04544 77.13257 Another verification:
plot_siid_dyn_equ2 <- function(
N = 1e5,
Is0 = 0,
If0 = 0,
D0 = 350,
pi = 0,
p = .9,
beta = 1e-6,
gammas = 0,
gammaf = 0,
deltas = 1e-6,
deltaf = 1e-4,
alpha = 1 / 70,
times = seq2(0, 500000)) {
sims <- siid_dyn(N, Is0, If0, D0,
pi, p, beta, gammas, gammaf, deltas, deltaf, alpha, times)
equs <- siid_equ2(N, p, beta, gammas, gammaf, deltas, deltaf, alpha)
plotl2 <- function(...) plotl(..., xlab = "time (year)")
opar <- par(mfrow = c(2, 2), cex = 1, plt = c(.25, .95, .25, .9))
with(sims, {
plotl2(time, Is + If, col = 1, ylab = "infection prevalence")
abline(h = equs["I"], col = 1)
plotl2(time, D, col = 2, ylab = "disease prevalence")
abline(h = equs["D"], col = 2)
plotl2(time, deltaf * If, col = 3, ylab = "yearly fast disease incidence")
abline(h = equs["df"], col = 3)
plotl2(time, deltas * Is, col = 4, ylab = "yearly slow disease incidence")
abline(h = equs["ds"], col = 4)
})
par(opar)
}plot_siid_dyn_equ2()
If, additionally, \(\gamma_s = \gamma_f = 0\), then the above equilibrium becomes:
\[ \begin{align} D^* &= \frac{ \delta_f\delta_s(\beta N - \alpha)}{ (\alpha + \delta_s)\delta_f - (\delta_f - \delta_s)\alpha p}\frac{1}{\beta} \\ I^* &= \frac{(N - D^*)(\delta_f - \delta_s)p\beta - \delta_f\alpha} {(\delta_f - \delta_s)p\beta D^* - \delta_f\delta_s}D^* \\ d_f^* &= \frac{(\delta_s N - (\alpha + \delta_s)D^*)\delta_f p\beta D^*} {\delta_f\delta_s - (\delta_f - \delta_s)p\beta D^*} \\ d_s^* &= \left( \alpha + \frac{(\delta_s N - [\alpha + \delta_s] D^*)p\beta\delta_f} {(\delta_f - \delta_s)p\beta D^* - \delta_f\delta_s} \right) D^* \end{align} \]
Numerical verification:
siid_equ3 <- function(N = 1e5, p, beta, deltas, deltaf, alpha) {
D <- deltaf * deltas * (beta * N - alpha) /
(((alpha + deltas) * deltaf - (deltaf - deltas) * alpha * p) * beta)
c(D = D,
I = ((N - D) * (deltaf - deltas) * p * beta - deltaf * alpha) * D /
((deltaf - deltas) * p * beta * D - deltaf * deltas),
df = (deltas * N - (alpha + deltas) * D) * deltaf * p * beta * D /
(deltaf * deltas - (deltaf - deltas) * p * beta * D),
ds = (alpha + (deltas * N - (alpha + deltas) * D) * p * beta * deltaf /
((deltaf - deltas) * p * beta * D - deltaf * deltas)) * D
)
}
siid_equ2(N = 1e5, p = .9, beta = .00001, gammas = 0, gammaf = 0, deltas = 1 / 50,
deltaf = 1 / 5, alpha = 1 / 70) D I df ds
86792.4528 11778.9757 1115.9030 123.9892 siid_equ3(N = 1e5, p = .9, beta = .00001, deltas = 1 / 50, deltaf = 1 / 5,
alpha = 1 / 70) D I df ds
86792.4528 11778.9757 1115.9030 123.9892 Now, let’s re-arrange this equilibrium to express 4 of the model’s parameters values as a function of infection and disease prevalences as well as disease incidences:
\[ \begin{align} \alpha &= \frac{d_f^* + d_s^*} {D^*} \\ \beta &= \frac{d_f^* + d_s^*} {(N - I^* - D^*)D^*} \\ \delta_s &= \frac{\delta_f d_s^*}{\delta_f I^* - d_f^*} \\ p &= \frac{d_f^*} {d_f^* + d_s^*} \end{align} \]
In R it reads:
siid_param_val <- function(I, D, ds, df, N, deltaf) {
c(alpha = (df + ds) / D,
beta = (df + ds) / ((N - I - D) * D),
deltas = deltaf * ds / (deltaf * I - df),
p = df / (df + ds))
}Let’s look at some examples:
x <- .9
di <- 176
D <- 322
df <- round(x * di)
siid_param_val(I = 1e4, D = D, ds = di - df, df = df, N = 1e5, deltaf = 1 / 2) alpha beta deltas p
5.465839e-01 6.094960e-06 1.858736e-03 8.977273e-01 siid_param_val(I = 3e4, D = D, ds = di - df, df = df, N = 1e5, deltaf = 1 / 2) alpha beta deltas p
5.465839e-01 7.844425e-06 6.063873e-04 8.977273e-01 Let’s do a numerical verification:
deltaf <- 1 / 2
N <- 1e5
x <- .9
di <- 176
D0 <- 322
df <- x * di
siid_p_val <- siid_param_val(I = 2e4, D = D0, ds = di - df, df = df, N = N,
deltaf = deltaf)
a <- siid_equ(N = N, pi = 0, p = siid_p_val[["p"]], beta = siid_p_val[["beta"]],
gammas = 0, gammaf = 0, deltas = siid_p_val[["deltas"]],
deltaf = deltaf, alpha = siid_p_val[["alpha"]])
b <- siid_equ2(N = N, p = siid_p_val[["p"]], beta = siid_p_val[["beta"]],
gammas = 0, gammaf = 0, deltas = siid_p_val[["deltas"]],
deltaf = deltaf, alpha = siid_p_val[["alpha"]])
a["D"] D
322 sum(a[c("Is", "If")])[1] 20000b["I"] I
20000 (d <- sum(b[c("df", "ds")]))[1] 176b[["df"]] / d[1] 0.95.2.4 Prophylaxis
pi2dis_inc <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
deltaf = 1 / 2,
N = 1e5,
rates = c(.02, .1, .2),
by = .1,
times = 2025 + seq2(0, 30),
eps = 5,
lgnd = TRUE,
show_rates = TRUE) {
df <- x * di
siid_p_val <- siid_param_val(I = I0, D = D0, ds = di - df, df = df, N = N,
deltaf = deltaf)
start_val <- siid_equ(N = N, pi = 0, p = siid_p_val[["p"]],
beta = siid_p_val[["beta"]],
gammas = 0, gammaf = 0, deltas = siid_p_val[["deltas"]],
deltaf = deltaf, alpha = siid_p_val[["alpha"]])
pi_effect <- function(pi) {
siid_dyn(N = 1e5, Is0 = start_val[["Is"]], If0 = start_val[["If"]], D0, pi,
p = siid_p_val[["p"]],
beta = siid_p_val[["beta"]], gammas = 0, gammaf = 0,
deltas = siid_p_val[["deltas"]], deltaf = deltaf,
alpha = siid_p_val[["alpha"]], times) |>
mutate(d = deltaf * If + siid_p_val[["deltas"]] * Is)
}
pi_vals <- seq(0, 1, .1)
maxt <- max(times)
mint <- min(times)
minx <- mint - eps
plot(NA, xlim = c(minx, maxt), ylim = c(0, di),
xlab = "year", ylab = "(/year/100,000)")
mtext("disease incidence", 2, 2.25)
cols <- spectral(length(pi_vals))
pi_vals |>
map(pi_effect) |>
walk2(cols, ~ with(.x, lines2(time, d, col = .y)))
segments2(minx, di, mint, di, col = cols[1])
if (show_rates) {
xs <- seq(mint, maxt, by)
walk2(rates, 2:4,
~ lines2(xs, di * exp(- find_rate(.x) * (xs - mint)), lty = .y))
}
if (lgnd) {
legend2("bottomleft", col = cols, title = expression(pi),
legend = formatC(round(pi_vals, 2), digits = 2, format = "f"))
legend(2048, 175, lty = 2:4, title = "annual decrease", bty = "n", lwd = 2,
legend = paste0(formatC(round(rates, 2), digits = 2, format = "f"), "%"))
}
}Let’s see:
pi2dis_inc()
5.3 Serial SIID
5.3.1 Dynamics
\[ \begin{align} \frac{dS}{dt} &= \gamma_1 I_1 + \gamma_2 I_2 + \alpha D - p\beta DS \\ \frac{dI_1}{dt} &= p\beta DS - (\gamma_1 + \delta_1 + \sigma)I_1 \\ \frac{dI_2}{dt} &= \sigma I_1 - (\delta_2 + \gamma_2) I_2 \\ \frac{dD}{dt} &= \delta_1 I_1 + \delta_2 I_2 - \alpha D \end{align} \]
siid2_dyn <- function(N = 1e5, I10 = 0, I20 = 0, D0,
p, beta, sigma, alpha, gamma1, gamma2, delta1, delta2, times) {
ode2(c(S = N - I10 - I20 - D0, I1 = I10, I2 = I20, D = D0),
times,
function(time, state, pars) {
with(as.list(c(state, pars)), {
infection <- p * beta * D * S
dS <- gamma1 * I1 + gamma2 * I2 + alpha * D - infection
dI1 <- infection - (gamma1 + delta1 + sigma) * I1
dI2 <- sigma * I1 - (delta2 + gamma2) * I2
dD <- delta1 * I1 + delta2 * I2 - alpha * D
list(c(dS, dI1, dI2, dD))
})
},
c(p = p, beta = beta, sigma = sigma, alpha = alpha, gamma1 = gamma1,
gamma2 = gamma2, delta1 = delta1, delta2 = delta2))
}delta1 <- 1 / 2
delta2 <- 1 / 20
siid2_dyn(N = 1e5,
I10 = 10000,
I20 = 0,
D0 = 300,
p = .1,
beta = .0001,
sigma = 1 / 2,
alpha = 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
delta1,
delta2,
times = seq2(0, 100)) |>
with(plotl(time, delta1 * I1 + delta2 * I2))
abline(h = 0)
delta1 <- 1 / 2
delta2 <- .0008941636
siid2_dyn(N = 1e5,
I10 = 316.8,
I20 = 19683.2,
D0 = 322,
p = 0,
beta = .00001985216,
sigma = 1.091077,
alpha = .5465839,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
delta1,
delta2,
times = seq2(0, 30)) |>
with(plotl(time, delta1 * I1 + delta2 * I2, ylim = c(0, 176)))
abline(h = 0)
5.3.2 Equilibrium
\[ \begin{align} S^* &= N - I_1^* - I_2^* -D^* \\ I_1^* &= \frac{[(\gamma_2\delta_1 + [\delta_1 + \sigma]\delta_2)p\beta N - (\gamma_1 + \delta_1 + \sigma) (\gamma_2 + \delta_2)\alpha](\gamma_2 + \delta_2)\alpha} {[(\alpha + \delta_1)(\gamma_2 + \delta_2) + (\alpha + \delta_2)\sigma] [\gamma_2\delta_1 + (\delta_1 + \sigma)\delta_2]p\beta} \\ I_2^* &= \frac{\sigma}{\delta_2 + \gamma_2} I_1^*\\ D^* &= \frac{\delta_1 I_1^* + \delta_2 I_2^*}{\alpha} \end{align} \]
siid2_equ <- function(N, p, beta, sigma, alpha, gamma1, gamma2, delta1, delta2) {
I1 <- ((gamma2 * delta1 + (delta1 + sigma) * delta2) * p * beta * N -
(gamma1 + delta1 + sigma) * (gamma2 + delta2) * alpha) *
(gamma2 + delta2) * alpha /
(((alpha + delta1) * (gamma2 + delta2) + (alpha + delta2) * sigma) *
(gamma2 * delta1 + (delta1 + sigma) * delta2) * p * beta)
I2 <- sigma * I1 / (delta2 + gamma2)
D <- (delta1 * I1 + delta2 * I2) / alpha
c(S = N - I1 - I2 - D, I1 = I1, I2 = I2, D = D)
}5.3.3 Calibration
\[ \begin{align} \alpha &= \frac{d_1^* + d_2^*}{D^*} \\ \beta &= \frac{(\delta_1 + \gamma_1 - \gamma_2)d_1^* + \delta_1d_2^* + \gamma_2\delta_1 I^*} {(N - I^* - D^*)\delta_1 D^*} \\ \sigma &= \frac{\delta_1 d_2^* - \gamma_2 d_1^* + \delta_1\gamma_2 I^*}{d_1^*} \\ \delta_2 &= \frac{\delta_1 d_2^*}{\delta_1 I^* - d_1^*} \end{align} \]
siid2_param_val <- function(I, D, d1, d2, delta1, gamma1, gamma2, N) {
c(
alpha = (d1 + d2) / D,
beta = ((delta1 + gamma1 - gamma2) * d1 + delta1 * d2 + gamma2 * delta1 * I) /
((N - I - D) * delta1 * D),
sigma = (delta1 * d2 - gamma2 * d1 + delta1 * gamma2 * I) / d1,
delta2 = delta1 * d2 / (delta1 * I - d1)
)
}The next 2 functions are functions that generate diagnostic plots following calibration. The first one shows the dynamics of leaving the compartment of recent infections:
leaving_recent_infections <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5) {
pexp <- .99
xtext <- .98
d1 <- x * di
p_val <- siid2_param_val(I0, D0, d1, di - d1, delta1, gamma1, gamma2, N)
lambda <- delta1 + p_val[["sigma"]] + gamma1
xmax <- ceiling(qexp(pexp, lambda))
xs <- seq2(0, xmax)
ys <- pexp(xs, lambda)
plot(NA, xlim = c(0, xmax), ylim = 0:1,
xlab = "leaving recent infections (years)",
ylab = "cumulative probability")
polygon3(xs, ys * gamma1 / lambda, rep(0, length(xs)), col = tones[1])
polygon3(xs, ys * (gamma1 + delta1) / lambda, ys * gamma1 / lambda,
col = skin_tones[1])
polygon3(xs, ys * (gamma1 + delta1) / lambda, ys, col = greys[1])
abline3(v = pretty(0:xmax))
abline3(h = seq(0, 1, .1))
lines2(xs, ys * gamma1 / lambda, col = tones[3])
lines2(xs, ys * (gamma1 + delta1) / lambda, col = skin_tones[3])
lines2(xs, ys, col = greys[3])
text2(xtext * xmax, .5, "long-lasting infections", col = greys[3], font = 2)
text2(xtext * xmax, .15, "disease", col = skin_tones[3], font = 2)
}The second one shows the distribution of the duration of disease together with R\(_0\) and the probability of long-lasting infections to develop disease.
disease_duration_distribution <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5) {
pexp <- .99
x1 <- .9
x2 <- .95
y <- seq(.93, .48, -.15)
eps <- .025
ps <- seq(0, .9, .1)
d1 <- x * di
p_val <- siid2_param_val(I0, D0, d1, di - d1, delta1, gamma1, gamma2, N)
alpha <- p_val[["alpha"]]
delta2 <- p_val[["delta2"]]
xmax <- ceiling(qexp(pexp, alpha))
xs <- seq2(0, xmax)
ys <- dexp(xs, alpha)
ymax <- max(ys)
qtls <- qexp(ps, alpha)
cols <- colorRampPalette(yellows)(length(qtls))
plot(NA, xlim = c(0, xmax), ylim = c(0, ymax),
xlab = "duration of disease (years)",
ylab = "density of probability")
walk2(qtls, cols,
function(q, col) {
sel <- q <= xs
polygon3(xs[sel], ys[sel], rep(0, length(xs[sel])), col = col)
})
abline4(v = 1 / alpha)
abline4(v = log(2) / alpha, lty = 2)
text2(x1 * xmax, y[1] * ymax,
bquote(paste("R"[0], " = ", .(round(N * p_val[["beta"]] / alpha, 2)))))
text2(x1 * xmax, y[2] * ymax,
paste("Pr[d|oi] =", round(delta2 / (delta2 + gamma2), 2)))
text2(x1 * xmax, y[3] * ymax, paste("median =", round(log(2) / alpha, 2)))
text2(x1 * xmax, y[4] * ymax, paste("mean =", round(1 / alpha, 2)))
text(x2 * xmax, (y[3] - eps) * ymax, "yr")
text(x2 * xmax, (y[4] - eps) * ymax, "yr")
}5.3.4 Prophylaxis
Here we use the epidemiological model to investigate the epidemiological consequences of prophylaxis. The first function generates the function used to explore the effect of prophylaxis:
make_p_effect <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5,
times = seq2(2025, 2030)) {
d1 <- x * di
p_val <- siid2_param_val(I0, D0, d1, di - d1, delta1, gamma1, gamma2, N)
alpha <- p_val[["alpha"]]
beta <- p_val[["beta"]]
sigma <- p_val[["sigma"]]
delta2 <- p_val[["delta2"]]
start_val <- siid2_equ(N, 1, beta, sigma, alpha, gamma1, gamma2, delta1, delta2)
function(p) {
siid2_dyn(N, start_val[["I1"]], start_val[["I2"]], start_val[["D"]],
p, beta, sigma, alpha, gamma1, gamma2, delta1, delta2, times) |>
mutate(d = delta1 * I1 + delta2 * I2)
}
}The above function is used by the next 4 functions used to vizualize the effect of prophylaxis under various angles. The first one looks at the effect of prophylaxis as a function of time for various levels of prophylaxis.
p2dis_inc <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5,
p_vals = seq(0, 1, .1),
rates = c(.02, .1, .2),
by = .05,
times = 2025 + seq2(0, 30),
eps = 5, xl = 2048, yl = 175) {
p_effect <- make_p_effect(I0, D0, di, x, delta1, gamma1, gamma2, N, times)
maxt <- max(times)
mint <- min(times)
minx <- mint - eps
plot(NA, xlim = c(minx, maxt), ylim = c(0, di),
xlab = "year", ylab = "(/year/100,000)")
mtext2("disease incidence")
abline3(v = pretty(minx:maxt))
abline3(h = pretty(0:di))
cols <- colorRampPalette(reds)(length(p_vals))
p_vals |>
map(p_effect) |>
walk2(rev(cols), ~ with(.x, lines2(time, d, col = .y)))
segments2(minx, di, mint, di, col = cols[1])
xs <- seq(mint, maxt, by)
walk2(rates, 2:4,
~ lines2(xs, di * exp(- find_rate(.x) * (xs - mint)), lty = .y))
}The second one looks at the yearly rate of disease incidence decrease (in percentage) as a function of time for various levels of prophylaxis:
p2dec_rate <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5,
p_vals = seq(0, 1, .1),
rates = c(.02, .1, .2),
times = seq2(2025, 2030),
maxt = NULL,
xl = 2048, yl = 175) {
p_effect <- make_p_effect(I0, D0, di, x, delta1, gamma1, gamma2, N, times)
decrease_rate <- function(x) {
dt <- diff(x$time)[1]
x |>
mutate(time = time + dt,
d = c(diff(d) / head(d, -1), NA)) |>
filter(row_number() <= n() - 1) |>
mutate(d = 100 * (1 - exp(d / dt)))
}
cols <- colorRampPalette(reds)(length(p_vals))
if (is.null(maxt)) maxt <- max(times)
mint <- min(times)
tmp <- p_vals |>
map(p_effect) |>
map(decrease_rate)
ymax <- max(map_dbl(tmp, ~ max(.x$d)))
plot(NA, xlim = c(mint, maxt), ylim = c(0, ymax),
xlab = "year", ylab = "% decrease / year")
abline3(v = pretty(mint:maxt))
abline3(h = pretty(c(0, ymax)))
walk2(tmp, rev(cols), ~ with(.x, lines2(time, d, col = .y)))
abline4(h = 100 * rates, lty = 2:4)
}The third one looks at the year a targeted level of disease incidence is reached as a function of the efficacy of prophylaxis and this for various levels of targeted disease incidence:
time2val <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5,
target_vals = NULL,
times = seq2(2025, 2055),
d_vals = head(seq2(0, 1), -1)) {
if (is.null(target_vals)) target_vals <- seq(.1, .9, .1) * di
cols <- rev(colorRampPalette(blues)(length(target_vals)))
p_effect <- make_p_effect(I0, D0, di, x, delta1, gamma1, gamma2, N, times)
compute_data <- function(target) {
time_vals <- unlist(
mclapply2(d_vals, function(x) with(p_effect(x), approx(d, time, target)$y)))
sel <- ! is.na(time_vals)
tibble(d = 1 - d_vals[sel], t = time_vals[sel])
}
plot(NA, xlim = 0:1, ylim = range(times),
xlab = "efficacy of prophylaxis",
ylab = "year reaching target")
abline3(v = pretty(0:1))
abline3(h = pretty(times))
target_vals |>
map(compute_data) |>
walk2(cols, ~ with(.x, lines2(d, t, col = .y)))
}The last one looks at the level of disease incidence that is reached as a function of the efficacy of prophylaxis, and this for various time horizons:
inc_at_time <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5,
horizons = NULL,
times = seq2(2025, 2055),
d_vals = head(seq2(0, 1), -1)) {
if (is.null(horizons)) horizons <- pretty(times)
cols <- colorRampPalette(greens)(length(horizons))
p_effect <- make_p_effect(I0, D0, di, x, delta1, gamma1, gamma2, N, times)
compute_data <- function(hor) {
inc_vals <- unlist(
mclapply2(d_vals, function(x) with(p_effect(x), approx(time, d, hor)$y)))
sel <- ! is.na(inc_vals)
tibble(d = 1 - d_vals[sel], i = inc_vals[sel])
}
plot(NA, xlim = 0:1, ylim = c(0, di),
xlab = "efficacy of prophylaxis",
ylab = "(/year/100,000)")
mtext("disease incidence", 2, 2.25)
abline3(v = pretty(0:1))
abline3(h = pretty(0:di))
horizons |>
map(compute_data) |>
walk2(cols, ~ with(.x, lines2(d, i, col = .y)))
}5.3.5 Dashboard
The following function puts together the diagnostic plots of the calibration and the plots of the exploration of the prophylaxis:
dashboard_epi <- function(
I0 = 2e4,
D0 = 322,
di = 176,
x = .9,
delta1 = 1 / 2,
gamma1 = 1 / 60,
gamma2 = 1 / 60,
N = 1e5,
p_vals = seq(0, 1, .1),
target_vals = NULL,
rates = c(.02, .1, .2),
horizons = seq(2025, 2055, 5),
by = .1,
times1 = seq2(2025, 2055),
times2 = seq2(2025, 2030),
d_vals = head(seq2(0, 1), -1),
eps = 2) {
# 1. duration of "recent infections"
leaving_recent_infections(I0, D0, di, x, delta1, gamma1, gamma2, N)
# 2. duration of disease
disease_duration_distribution(I0, D0, di, x, delta1, gamma1, gamma2, N)
# 3. decrease in disease incidence
p2dis_inc(I0, D0, di, x, delta1, gamma1, gamma2, N, p_vals, rates, by, times1, eps)
legend4("topright", lty = 2:4,
legend = paste0(formatC(round(100 * c(.02, .1, .2)),
digits = 0, width = 2, format = "f"), "%/yr"))
# 4. disease incidence as function of prophylaxis efficacy and year
inc_at_time(I0, D0, di, x, delta1, gamma1, gamma2, N, horizons, times1, d_vals)
cols <- colorRampPalette(greens)(7)
cols <- c(cols[1:2], "white", cols[7])
legend5("topright", legend = c("2025", "2030", "...", "2055"), col = cols)
# 5. rate of decrease of disease incidence
p2dec_rate(I0, D0, di, x, delta1, gamma1, gamma2, N, p_vals, rates, times2)
cols <- colorRampPalette(reds)(10)
cols <- c(cols[1:2], "white", cols[10])
legend5("topright", title = "1 - p", legend = c("0.00", "0.10", "...", "1.00"),
col = cols)
# 6. time to target
if (is.null(target_vals)) target_vals <- seq(25, di, 25)
time2val(I0, D0, di, x, delta1, gamma1, gamma2, N, target_vals, times1, d_vals)
ll <- length(target_vals)
cols <- rev(colorRampPalette(blues)(ll))
cols <- c(cols[1:2], "white", cols[ll])
legend5("topright", title = "incidence", col = cols,
legend = c(target_vals[1:2], "...", target_vals[length(target_vals)]))
}Let’s play with it:
opar <- par(mfrow = c(3, 2), cex = 1, plt = c(.25, .95, .25, .9))
system.time(
dashboard_epi(times1 = seq(2025, 2055, le = 512),
times2 = seq(2025, 2030, le = 512),
d_vals = head(seq(0, 1, le = 512), -1))
)
par(opar)opar <- par(mfrow = c(3, 2), cex = 1, plt = c(.25, .95, .25, .9))
system.time(
dashboard_epi(times1 = seq(2025, 2055, le = 100),
times2 = seq(2025, 2030, le = 100),
d_vals = head(seq(0, 1, le = 100), -1))
)
user system elapsed
20.147 14.115 6.181 par(opar)