Preventing and controlling infectious disease outbreaks requires making strategic use of limited vaccine stockpiles. This challenge is especially acute for diseases such as cholera and typhoid fever, where global supplies remain constrained and timely deployment is crucial. In principle, directing vaccines toward high-risk populations should yield the greatest benefit, but real-world uncertainty—in outbreak timing, risk signals, and operational response—can make simple pro-rata allocation surprisingly competitive.
Existing modeling studies typically assume that outbreaks will occur and therefore evaluate strategies only conditional on an epidemic happening. What remains unclear is how to choose between pre-emptive, reactive, or mixed vaccination strategies when outbreaks are uncertain and when economic constraints matter.
To address this gap, we developed a general analytical framework that identifies the vaccination strategy minimizing total expected societal cost, including both vaccination expenditures and potential outbreak losses. Our approach integrates limited vaccine supply, heterogeneous outbreak risks, targeting accuracy, and economic trade-offs.
Using a combination of closed-form analytical results and Monte Carlo simulations, we characterize optimal allocation rules across a wide range of operational settings. This framework provides practical, quantitative guidance for policymakers designing vaccination strategies under uncertainty.
This component covers the cost of inpatient and outpatient care for all cholera cases:
\[ C_{\text{illness}}^{\text{direct}} = N_{\text{inpatient}}\times C_{\text{inpatient}} + N_{\text{outpatient}}\times C_{\text{outpatient}} . \]
We estimate:
\[ C_{\text{illness}}^{\text{indirect}}= N_{\text{cases}}\times d \times w , \]
where
We compare two approaches:
\[ C_{\text{death}}^{\text{HC}} = D \times \sum_{t=1}^{YLL} \frac{G}{(1+r)^t} \]
\[ C_{\text{death}}^{VSL} = D \times V \]
Where:
The total outbreak cost or the cost from infection is
\[ C_{\text{outbreak}} = C_{\text{I}} = C_{\text{illness}}^{\text{direct}} + C_{\text{illness}}^{\text{indirect}} + C_{\text{death}}. \]
We use:
\[ C_{\text{vaccination}} = C_{\text{V}} =N_{\text{people}}\times (c_{\text{dose}} + c_{\text{delivery}}), \]
where \(c_{\text{dose}}\) = vaccine cost per dose, \(c_{\text{delivery}}\) = delivery cost per dose.
Let:
We consider a single population with outbreak probability \(p\), outbreak cost \(C_{\mathrm{I}}\), vaccination cost \(C_{\mathrm{V}}\), and reactive vaccination effectiveness \(r \in [0,1]\). The outbreak cost \(C_{\mathrm{I}}\) reflects direct medical costs, productivity losses, and mortality-related losses. The vaccination cost \(C_{\mathrm{V}}\) includes vaccine and delivery costs.
Let \(X \sim \mathrm{Bernoulli}(p)\) denote an indicator of an outbreak, where \(X=1\) if an outbreak occurs and \(X=0\) otherwise.
Under pre-emptive vaccination, the population is vaccinated regardless of whether an outbreak occurs. Any potential outbreak is fully averted, so the total cost is deterministic:
\[ C_{\mathrm{pre}}(X) = C_{\mathrm{V}}. \]
Thus,
\[ \mathbb{E}[C_{\mathrm{pre}}] = C_{\mathrm{V}}. \]
Under reactive vaccination, the population is vaccinated only if an outbreak occurs. If \(X=1\), the vaccination cost \(C_{\mathrm{V}}\) is incurred and residual outbreak costs equal \((1-r) C_{\mathrm{I}}\). If \(X=0\), there is no cost. The total cost random variable is:
\[ C_{\mathrm{react}}(X) = X\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr). \]
Taking expectation,
\[ \mathbb{E}[C_{\mathrm{react}}] = p\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr). \]
This calculation assumes no penalty or loss from unused vaccine stock. In practice, vaccine expiry or wastage may incur additional costs, but these are omitted here for analytical clarity.
Define the outbreak-to-vaccination cost ratio:
\[ R = \frac{C_{\mathrm{I}}}{C_{\mathrm{V}}}. \]
Normalizing by \(C_{\mathrm{V}}\):
\[ c_{\mathrm{pre}}^{(1)} = \frac{\mathbb{E}[C_{\mathrm{pre}}]}{C_{\mathrm{V}}} = 1. \]
\[ c_{\mathrm{react}}^{(1)} = \frac{\mathbb{E}[C_{\mathrm{react}}]}{C_{\mathrm{V}}} = p\bigl[1 + (1-r)R\bigr]. \]
These represent the normalized per-population expected costs of the two strategies with the superscript \((1)\) indicating one-population case.
The threshold outbreak probability \(p_{\mathrm{thr}}^{(1)}\) at which the two strategies yield equal expected cost is given by:
\[ c_{\mathrm{pre}}^{(1)} = c_{\mathrm{react}}^{(1)} \quad\Longrightarrow\quad p_{\mathrm{thr}}^{(1)} = \frac{1}{1 + (1-r)R}. \]
Interpretation:
Pre-emptive preferred when \[ p > p_{\mathrm{thr}}^{(1)} \quad\Longleftrightarrow\quad c_{\mathrm{pre}}^{(1)} < c_{\mathrm{react}}^{(1)}. \]
Reactive preferred when \[ p < p_{\mathrm{thr}}^{(1)} \quad\Longleftrightarrow\quad c_{\mathrm{pre}}^{(1)} > c_{\mathrm{react}}^{(1)}. \]
Indifference point \[ p = p_{\mathrm{thr}}^{(1)} \quad\Longleftrightarrow\quad c_{\mathrm{pre}}^{(1)} = c_{\mathrm{react}}^{(1)}. \]
The threshold \(p_{\mathrm{thr}}^{(1)}\) decreases as the reactive effectiveness \(r\) decreases or the outbreak cost ratio \(R\) increases. Thus, pre-emptive vaccination becomes more favorable when reactive vaccination is less effective or when outbreaks are relatively more costly.
# Pre-emptive cost (normalized by C_V)
cost_pre_one <- function(p, R, r) {
# p and r unused but kept for a consistent interface
1
}
# Reactive cost (normalized by C_V)
cost_react_one <- function(p, R, r) {
p * (1 + (1 - r) * R)
}
p_star_one <- function(R, r) {
1 / (1 + (1 - r) * R)
}In the figure below, \(c_s^{(1)}\) denotes the normalized per-population expected cost for strategy \(s \in \{\mathrm{pre}, \mathrm{react}\}\).
Key insights are: 1. The cost remains constant for the pre-emptive strategy whereas it increases with the outbreak probability \(p\) for the reactive strategy. 2. There is a threshold outbreak probability \(p_{\mathrm{thr}}^{(1)}\) at which the costs of the pre-emptive and reactive strategies are equal, below which reactive vaccination is economically preferable. 3. The threshold \(p_{\mathrm{thr}}^{(1)}\) is higher for greater effectiveness of the reactive strategy, \(r\). 4. As long as vaccine-induced protection at the population level is maintained, the pre-emptive strategy becomes preferable over an extended period of duration as the overall probability of an outbreak would increase over time.
# --- Parameters for the illustration ---
p_vals <- seq(0, 1, 0.01) # mean outbreak probability
R_vals <- c(5)
r_vals <- c(0.3, 0.7)
df_cost <- expand.grid(
p = p_vals,
r = r_vals,
R = R_vals
) |>
dplyr::mutate(
c_pre = cost_pre_one(p, R, r),
c_react = cost_react_one(p, R, r),
r_label = paste0("italic(r)==", r)) |>
tidyr::pivot_longer(
cols = c(c_pre, c_react),
names_to = "strategy",
values_to = "cost"
) |>
dplyr::mutate(
strategy = factor(
strategy,
levels = c("c_pre", "c_react"),
labels = c("Pre-emptive", "Reactive")
)
)
# threshold probability
# df_pstar <- expand.grid(r = r_vals, R = R_vals) |>
# dplyr::mutate(
# p_star = p_star_one(R, r),
# r_label = paste0("r = ", r)
# )
df_pstar <- expand.grid(r = r_vals, R = R_vals) |>
dplyr::mutate(
p_star = p_star_one(R, r),
r_label = paste0("r = ", r),
# [NEW] Create the dynamic label string here
# ~ adds a space, italic(r)==r adds the value
label_expr = paste0("italic(p)[thr]^(1) ~ (italic(r) == ", r, ")")
)
# data for braces
cost_react_max1 <- max(filter(df_cost, strategy == "Reactive", r_label == "italic(r)==0.7")$cost)
cost_react_max2 <- max(filter(df_cost, strategy == "Reactive", r_label == "italic(r)==0.3")$cost)
df_brace1 <- data.frame(x=c(0.96,1), y=c(1,cost_react_max1))
df_brace2 <- data.frame(x=c(0.97,1), y=c(1,cost_react_max2))
df_brace3 <- data.frame(x=c(0,0.03), y=c(0,1))
# convenient variables for the overhead-cost text position
text_overhead_x <- max(df_brace3$x) + 0.01 # 0.04
text_overhead_y <- mean(df_brace3$y) + 1 # 1.5
annot_text_size <- 4
ggplot(df_cost, aes(x = p, y = cost, color = strategy,
linetype = r_label)) +
geom_line(linewidth = 1) +
geom_vline(
data = df_pstar,
aes(xintercept = p_star),
linetype = "dashed",
color = "firebrick"
)+
geom_text(
data = df_pstar,
aes(x = p_star, y = Inf, label = label_expr),
parse = TRUE,
inherit.aes = FALSE,
hjust = -0.03,
vjust = 1.1,
size = 3
) +
ggbrace::stat_brace(
data = df_brace1,
mapping = aes(x, y),
outside = FALSE, rotate = 270,
linewidth = 1, inherit.aes = FALSE
)+
annotate("text", x = min(df_brace1$x)-0.01,
y = mean(df_brace1$y),
label = "Cost of delayed\nresponse (r=0.7)",
hjust = 1, size = annot_text_size, lineheight = 0.9) +
ggbrace::stat_brace(
data = df_brace2,
mapping = aes(x, y),
outside = FALSE, rotate = 270,
linewidth = 1, inherit.aes = FALSE
)+
annotate("text", x = min(df_brace2$x)-0.01, y = mean(df_brace2$y),
label = "Cost of delayed\nresponse (r=0.3)",
hjust = 1, size = annot_text_size, lineheight = 0.9) +
ggbrace::stat_brace(
data = df_brace3,
mapping = aes(x, y),
outside = FALSE, rotate = 90,
linewidth = 1, inherit.aes = FALSE
)+
annotate("text", x = max(df_brace3$x)+0.01,
y= mean(df_brace3$y)+1,
label = "Vaccination cost",
hjust = 0, size = annot_text_size, lineheight = 0.9) +
annotate("segment",
x = text_overhead_x + 0.03,
y = text_overhead_y - 0.1,
xend = 0.02, yend = 0.75,
arrow = arrow(length = grid::unit(0.15, "cm")),
colour = "black"
) +
scale_color_manual("", values=c("firebrick","steelblue"))+
scale_linetype_discrete(labels = scales::label_parse())+
labs(
x = expression("Probability of an outbreak "~italic(p)),
y = expression("Normalized per-population expected cost " ~ italic(c)[s]^'(1)'),
color = "",
linetype = "")+
theme_light() +
theme(legend.position = "top")+
guides(
linetype = guide_legend(
override.aes = list(color = "steelblue")
)
)
r_vals <- seq(0.01, 0.99, by = 0.01)
R_vals <- c(1)
data <- expand.grid(r = r_vals, R = R_vals) %>%
dplyr::mutate(pstar = p_star_one(R=R, r=r))
ggplot(data, aes(x = r, y = pstar)) +
geom_line(linewidth = 1) +
geom_abline(slope=1, intercept = 0, linetype="dotted",
linewidth=1, color="firebrick")+
labs(x = expression("Reactive effectiveness "~italic(r)),
y = expression("Threshold outbreak probability "~italic(p)[thr]^(1)))+
theme_light()+
theme(legend.position = "top")+
annotate("text", size=4, x=0,y=1,hjust=0,vjust=1, label="Pre-emptive")+
annotate("text", size=4, x=1,y=0,hjust=1,vjust=0, label="Reactive")
r_vals <- 0.5
R_vals <- seq(0.1, 50, by=0.1)
data <- expand.grid(r = r_vals, R = R_vals) %>%
dplyr::mutate(pstar = p_star_one(R=R, r=r))
ggplot(data, aes(x = R, y = pstar)) +
geom_line(linewidth = 1) +
labs(x = expression("Cost ratio "~italic(R)),
y = expression("Threshold outbreak probability "~italic(p)[thr]^(1)))+
theme_light()+
theme(legend.position = "top")+
annotate("text", size=4, x=10,y=1,hjust=0,vjust=1, label="Pre-emptive")+
annotate("text", size=4, x=10,y=0,hjust=0,vjust=0, label="Reactive")
r_vals <- seq(0.01, 0.99, by = 0.01)
R_vals <- c(0.1, 1, 10)
data <- expand.grid(r = r_vals, R = R_vals) %>%
dplyr::mutate(pstar = p_star_one(r=r, R=R))
ggplot(data, aes(x = r, y = pstar, color = factor(R))) +
geom_line(linewidth = 1) +
geom_abline(slope=1, intercept = 0, linetype="dotted",
linewidth=1, color="firebrick")+
labs(x = expression("Reactive effectiveness "~italic(r)),
y = expression("Threshold outbreak probability "~italic(p)[thr]^(1)),
color = expression(italic(R))) +
theme_light()+
theme(legend.position = "top")+
annotate("text", size=4, x=0,y=1,hjust=0,vjust=1, label="Pre-emptive")+
annotate("text", size=4, x=1,y=0,hjust=1,vjust=0, label="Reactive")
# Grid over r and R
r_vals <- seq(0, 1, length.out = 100)
R_vals <- 10 ^ seq(-2, 2, length.out = 100)
# --- 2. Create Tidy Data for ggplot ---
# Instead of outer(), we use expand.grid() to get x, y, z columns
df_heatmap <- expand.grid(R = R_vals, r = r_vals) %>%
mutate(p_val = p_star_one(R, r))
# --- 3. Create Label Data ---
# Locations to label regions (Pre-emptive vs Reactive)
df_labels <- data.frame(
R = c(10, 0.1), # R_pre, R_reac
r = c(0.2, 0.9), # r_pre, r_reac
lab = c("Pre-emptive\nregion", "Reactive\nregion")
)
# Compute the curve p_thr^(1)=0.5 ---
# For each r, find R such that p_star_one(R,r)=0.5
df_boundary <- lapply(r_vals, function(rr) {
# Define function in R to solve
f_root <- function(R) p_star_one(R, rr) - 0.5
# Solve in log10(R) space for stability
sol <- tryCatch(
uniroot(f_root, interval = c(0.01, 100)),
error = function(e) NULL
)
if (is.null(sol)) return(NULL)
data.frame(
r = rr,
R = sol$root
)
}) %>% bind_rows()
# Pick a midpoint location for annotation
i_mid <- round(nrow(df_boundary) / 2)
annot_R <- df_boundary$R[i_mid]
annot_r <- df_boundary$r[i_mid]
# --- 4. Plotting ---
ggplot(df_heatmap, aes(x = R, y = r)) +
# Use geom_raster for heatmaps (faster/smoother than geom_tile for dense grids)
geom_raster(aes(fill = p_val)) +
# Add red boundary line
geom_line(data = df_boundary, aes(x = R, y = r),
color = "firebrick", linewidth=1) +
# Add annotation for p=0.5
annotate("text",
x = annot_R,
y = annot_r,
label = expression(italic(p) == 0.5),
color = "firebrick",
size = 4,
hjust = -0.2) +
# Add text labels for the regions
geom_text(data = df_labels, aes(label = lab),
color = "white", size = 4) +
# Scales
scale_x_log10(
breaks = c(0.01, 0.1, 1, 10, 100),
labels = c("0.01", "0.1", "1", "10", "100"),
expand = c(0,0) # Removes whitespace at edges
) +
scale_y_continuous(expand = c(0,0)) +
# Colors (Viridis matches your original Plotly choice)
scale_fill_viridis_c(
option = "viridis",
name = expression(italic(p)[thr]^(1)) # Mathematical expression for legend
) +
# Labels and Theme
labs(
x = expression("Cost ratio "~italic(R) == C[I] / C[V]),
y = expression("Reactive effectiveness "~italic(r))
) +
theme_minimal() +
theme(
panel.grid = element_blank(), # Heatmaps usually look better without grids
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12),
legend.key.height = unit(1.5, "cm") # Make legend bar taller
)
In a single population, the normalized per-population expected cost of the pre-emptive strategy, \(c_{\mathrm{pre}}^{(1)}\), remains at 1, while the reactive strategy cost increases linearly with the outbreak probability \(p\). Their intersection defines the threshold \(p_{\mathrm{thr}}^{(1)}\). Pre-emptive vaccination is favored when \(r\) is low and \(R\) is high. This analysis assumes no costs associated with unused vaccines, a simplifying assumption that can be relaxed in more detailed models.
We now extend the model to \(n\) independent populations, each with the same outbreak probability \(p\). The total vaccine stockpile is limited and can cover only a fraction \(0 < f \le 1\) of the populations—equivalently, we can conduct \(fn\) pre-emptive or reactive vaccination campaigns in total.
Importantly, we treat vaccination at the population level: each population is either fully vaccinated or not vaccinated at all. We do not consider scenarios where all populations receive reduced coverage; instead, the capacity constraint applies solely to the number of populations that can be vaccinated.
The superscript \((n)\) in \(c_{\mathrm{pre}}^{(n)}\), \(c_{\mathrm{react}}^{(n)}\), and \(c_{\mathrm{mixed}}^{(n)}\) indicates that these are normalized per-population costs in the \(n\)-population equal-risk setting.
We first consider two benchmark strategies that both respect the capacity constraint. We solve the model for large \(n\) such that the proportion of the population that experience outbreaks is approximately \(p\).
All \(fn\) campaigns are used pre-emptively:
The per-population expected cost is
\[ C_{\mathrm{pre}}^{(n)} = f C_{\mathrm{V}} + (1-f)\,p C_{\mathrm{I}}. \]
Normalizing by \(C_{\mathrm{V}}\) gives
\[ c_{\mathrm{pre}}^{(n)} = \frac{C_{\mathrm{pre}}^{(n)}}{C_{\mathrm{V}}} = f + (1-f)pR. \]
Under the pure reactive strategy, no population is vaccinated pre-emptively, and all \(fn\) campaigns are reserved for outbreak response.
Two regimes arise:
Capacity is sufficient to respond to essentially all outbreaks (\(fn \ge pn\)). For each population, with probability \(p\), an outbreak occurs and is reactively vaccinated. The associated cost is \(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\). With probability \((1-p)\), no outbreak occurs and no cost is incurred.
Thus
\[ C_{\mathrm{react}}^{(n)} = p\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr), \qquad c_{\mathrm{react}}^{(n)} = p\bigl[1 + (1-r)R\bigr], \quad\text{for } f \ge p. \]
On average there are more outbreaks than available campaigns (\(pn > fn\)). Assuming reactive campaigns are allocated uniformly at random among outbreaks, the fraction of outbreaks that receive a reactive campaign is \(f/p\).
For a given population, with probability \(p\), an outbreak occurs and conditional on outbreak, with probability \(f/p\), reactive vaccination occurs. The associated cost is \(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\). Also, conditional on outbreak, with probability \(1 - f/p\), no campaign occurs and the associated cost is simply \(C_{\mathrm{I}}\).
The per-population expected cost is
\[ \begin{aligned} C_{\mathrm{react}}^{(n)} &= p \left[ \frac{f}{p}\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr) + \left(1 - \frac{f}{p}\right) C_{\mathrm{I}} \right] \\ &= f\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr) + (p-f) C_{\mathrm{I}} \\ &= f C_{\mathrm{V}} + \bigl[p - fr\bigr] C_{\mathrm{I}}. \end{aligned} \] Normalizing,
\[ c_{\mathrm{react}}^{(n)} = \frac{C_{\mathrm{react}}^{(n)}}{C_{\mathrm{V}}} = f + \bigl[p - fr\bigr] R, \qquad \text{for } f < p. \]
We can summarize the pure reactive cost as
\[ c_{\mathrm{react}}^{(n)} = \begin{cases} f + (p - fr)R, & f < p,\\[6pt] p\bigl[1 + (1-r)R\bigr], & f \ge p. \end{cases} \]
cost_pre_multi <- function(p, R, r, f) {
f + (1 - f) * p * R
}
cost_react_multi <- function(p, R, r, f) {
ifelse(
f < p,
# reactive-limited
f + (p - f * r) * R,
# reactive-rich
p * (1 + (1 - r) * R)
)
}
cost_diff_multi <- function(p, R, r, f) {
c_pre <- cost_pre_multi(p, R, r, f)
c_react <- cost_react_multi(p, R, r, f)
c_pre - c_react
}
p_star_multi <- function(R, r, f, tol = 1e-10) {
# Candidate 1: reactive-limited regime (p > f): p* = r
p_star_RL <- r
valid_RL <- (p_star_RL > f) && (p_star_RL >= 0) && (p_star_RL <= 1)
if (valid_RL) {
diff_val <- cost_diff_multi(p_star_RL, R = R, r = r, f = f)
if (abs(diff_val) < 1e-6) {
return(p_star_RL)
}
}
# Candidate 2: reactive-rich regime (p <= f)
denom <- R * (r - f) - 1
if (abs(denom) < tol) {
p_star_RR <- NA_real_
} else {
p_star_RR <- -f / denom
}
valid_RR <- !is.na(p_star_RR) &&
(p_star_RR >= 0) && (p_star_RR <= f) && (p_star_RR <= 1)
if (valid_RR) {
diff_val <- cost_diff_multi(p_star_RR, R = R, r = r, f = f)
if (abs(diff_val) < 1e-6) {
return(p_star_RR)
}
}
# Fallback: numeric search
p_grid <- seq(1e-6, 1 - 1e-6, length.out = 2001)
vals <- cost_diff_multi(p_grid, R = R, r = r, f = f)
idx <- which(vals[-1] * vals[-length(vals)] <= 0)
if (length(idx) == 0) {
return(NA_real_)
}
i <- idx[1]
lower <- p_grid[i]
upper <- p_grid[i + 1]
uniroot(cost_diff_multi,
lower = lower, upper = upper,
R = R, r = r, f = f, tol = tol)$root
}The normalzied per-population expected cost linearly increases with \(f\) for the pre-emptive strategy whereas, in case of reactive strategy, it increases only up to a point where \(f=p\) and then remains constant.
\(p = 0.4, R = 1\).
p_val <- 0.4
R_val <- 1
r_vals <- c(0.2, 0.4, 0.6)
f_grid <- seq(0, 1, by = 0.01)
df_cost <- expand.grid(
f = f_grid,
r = r_vals
) |>
dplyr::mutate(
c_pre = cost_pre_multi(p = p_val, R = R_val, r = r, f = f),
c_react = cost_react_multi(p = p_val, R = R_val, r = r, f = f),
r_label = paste0("r = ", r)
) |>
tidyr::pivot_longer(
cols = c(c_pre, c_react),
names_to = "strategy",
values_to = "cost"
) |>
dplyr::mutate(
strategy = factor(
strategy,
levels = c("c_pre", "c_react"),
labels = c("Pre-emptive", "Reactive")
)
)
ggplot(df_cost, aes(x = f, y = cost, color = strategy)) +
geom_line(linewidth = 1) +
geom_vline(
xintercept = p_val,
linetype = "dashed",
linewidth = 0.8,
color = "black"
) +
facet_wrap(~ r_label, nrow = 1) +
labs(
x = expression("Vaccinated fraction"~italic(f)),
y = expression("Normalized per-population expected cost " ~ italic(c)[s]^'(n)'),
color = ""
) +
annotate(
"text",
x = p_val,
y = max(df_cost$cost, na.rm = TRUE),
label = "f == p",
parse = TRUE,
hjust = -0.1,
vjust = 1.2,
size = 4
) +
theme_light() +
theme(
legend.position = "top"
)
\(p = 0.4, R = 5\).
p_val <- 0.4
R_val <- 5
r_vals <- c(0.2, 0.4, 0.6)
f_grid <- seq(0, 1, by = 0.01)
df_cost <- expand.grid(
f = f_grid,
r = r_vals
) |>
dplyr::mutate(
c_pre = cost_pre_multi(p = p_val, R = R_val, r = r, f = f),
c_react = cost_react_multi(p = p_val, R = R_val, r = r, f = f),
r_label = paste0("r = ", r)
) |>
tidyr::pivot_longer(
cols = c(c_pre, c_react),
names_to = "strategy",
values_to = "cost"
) |>
dplyr::mutate(
strategy = factor(
strategy,
levels = c("c_pre", "c_react"),
labels = c("Pre-emptive", "Reactive")
)
)
ggplot(df_cost, aes(x = f, y = cost, color = strategy)) +
geom_line(linewidth = 1) +
geom_vline(
xintercept = p_val,
linetype = "dashed",
linewidth = 0.8,
color = "black"
) +
facet_wrap(~ r_label, nrow = 1) +
labs(
x = expression("Vaccinated fraction"~italic(f)),
y = expression("Normalized per-population expected cost " ~ italic(c)[s]^'(n)'),
color = ""
) +
annotate(
"text",
x = p_val,
y = max(df_cost$cost, na.rm = TRUE),
label = "f == p",
parse = TRUE,
hjust = -0.1,
vjust = 1.2,
size = 4
) +
theme_light() +
theme(
legend.position = "top"
)
r_vals <- seq(0, 1, by = 0.005)
R_vals <- c(0.1, 1, 10)
f_val <- 0.5
df_multi <- expand.grid(r = r_vals, R = R_vals) |>
dplyr::rowwise() |>
dplyr::mutate(p_star = p_star_multi(R = R, r = r, f = f_val)) |>
dplyr::ungroup()
ggplot(df_multi, aes(x = r, y = p_star, color = factor(R))) +
geom_line(linewidth = 1) +
geom_abline(slope = 1, intercept = 0,
linetype = "dotted",
linewidth = 1, color = "firebrick") +
geom_abline(slope = 0, intercept = f_val,
linetype = "dotted",
linewidth = 1, color = "steelblue") +
scale_y_continuous(limits = c(0, 1)) +
labs(
x = expression(italic(r)),
y = expression(italic(p)[thr]^(italic(n))),
color = expression(italic(R))
) +
theme_light() +
theme(legend.position = "top") +
annotate("text", size = 5,
x = 0, y = 1, hjust = 0, vjust = 1,
label = "Pre-emptive") +
annotate("text", size = 4,
x = 1, y = 0, hjust = 1, vjust = 0,
label = "Reactive")+
annotate("text", size = 4,
x = 1, y = f_val, hjust = 1, vjust = -0.3,
label = expression(italic(f)))
f_val <- 0.3
r_vals <- seq(0, 1, length.out = 100)
R_vals <- 10 ^ seq(-2, 2, length.out = 100)
p_mat_n <- outer(
R_vals, r_vals,
Vectorize(function(R, r) p_star_multi(R = R, r = r, f = f_val))
)
p_mat_n_t <- t(p_mat_n)
R_pre <- 10
r_pre <- 0.2
p_pre_star <- p_star_multi(R = R_pre, r = r_pre, f = f_val)
p_pre <- min(1, p_pre_star + 0.6)
R_reac <- 10
r_reac <- 0.8
p_reac_star <- p_star_multi(R = R_reac, r = r_reac, f = f_val)
p_reac <- 0.4
df_labels_n <- data.frame(
R = c(R_pre, R_reac),
r = c(r_pre, r_reac),
p = c(p_pre, p_reac),
lab = c("pre-emptive", "reactive")
)
p_equal_f_mat <-
matrix(f_val, nrow = length(r_vals), ncol = length(R_vals))
camera <- list(
eye = list(x = 0.1, y = 0.08, z = 0.1)
)
titlefontsize <- 34
tickfontsize <- 17
fig_pstar_multi <-
plot_ly() %>%
add_surface(
x = R_vals,
y = r_vals,
z = p_mat_n_t,
colorscale = "Viridis",
showscale = TRUE,
colorbar = list(
title = list(
text = "<i>p</i><sub>thr</sub><sup>(<i>n</i>)</sup>",
font = list(size = titlefontsize)
)
)
) %>%
add_trace(
data = df_labels_n,
x = ~R,
y = ~r,
z = ~p,
type = "scatter3d",
mode = "text",
text = ~lab,
textposition = "middle center",
textfont = list(size = titlefontsize, color = "black"),
showlegend = FALSE
) %>%
add_surface(
x = R_vals,
y = r_vals,
z = p_equal_f_mat,
opacity = 0.5,
showscale = FALSE,
colorscale = list(c(0, "grey50"), c(1, "grey50")),
name = "p=f"
) %>%
add_trace(
x = 100,
y = 1,
z = f_val,
type = "scatter3d",
mode = "text",
text = "<i>p</i> = <i>f</i>",
textfont = list(size = titlefontsize, color = "black"),
showlegend = FALSE
) %>%
layout(
scene = list(
camera = camera,
yaxis = list(
title = "<i>r</i>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)),
xaxis = list(
title = "<i>R</i> =C<sub>I</sub>/C<sub>V</sub>",
type = "log",
tickvals = c(0.01, 0.1, 1, 10, 100),
ticktext = c("0.01", "0.1", "1", "10", "100"),
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)
),
zaxis = list(
title = "<i>p<sub>*</sub></i>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize))
)
)
fig_pstar_multi\(f=0.5\)
# Parameters
f_val <- 0.5 # capacity fraction
p_val <- 0.6 # mean outbreak probability (not directly used here)
r_vals <- seq(0, 1, length.out = 200)
R_vals <- 10 ^ seq(-2, 2, length.out = 200)
# --- 1. Heatmap data: p_thr^(n) across (R, r) ---
df_heatmap_n <- expand.grid(R = R_vals, r = r_vals) %>%
mutate(
p_thr_n = mapply(
function(R_single, r_single) {
p_star_multi(R = R_single, r = r_single, f = f_val)
},
R, r
)
)
# --- 2. Multiple boundary curves for p_thr^(n) ---
p_levels <- seq(0.2, 0.7, by = 0.1)
df_boundary_n <- lapply(p_levels, function(p_target) {
tmp <- lapply(R_vals, function(RR) {
f_root <- function(r) p_star_multi(R = RR, r = r, f = f_val) - p_target
sol <- tryCatch(
uniroot(f_root, interval = c(0, 1)),
error = function(e) NULL
)
if (is.null(sol)) return(NULL)
data.frame(
R = RR,
r = sol$root,
p_target = p_target
)
}) %>% bind_rows()
tmp
}) %>% bind_rows()
# --- 3. Midpoints along each boundary curve for annotation ---
df_boundary_labels <- df_boundary_n %>%
group_by(p_target) %>%
slice(round(n() / 2)) %>%
ungroup() %>%
mutate(
label = sprintf("p[thr]^{(n)}==%.1f", p_target)
)
# --- 4. Region labels (optional) ---
R_pre <- 10
r_pre <- 0.2
R_reac <- 10
r_reac <- 0.8
df_labels_n <- data.frame(
R = c(R_pre, R_reac),
r = c(r_pre, r_reac),
lab = c("pre-emptive", "reactive")
)
# --- 5. Heatmap with multiple boundary lines (all same style) ---
ggplot(df_heatmap_n, aes(x = R, y = r)) +
# Heatmap
geom_raster(aes(fill = p_thr_n)) +
# Boundary curves (same color, no legend)
geom_line(
data = df_boundary_n,
aes(x = R, y = r, group = p_target),
color = "firebrick",
linewidth = 0.8
) +
# Line labels
geom_text(
data = df_boundary_labels,
aes(x = R, y = r, label = label),
color = "firebrick",
size = 3.5,
vjust = -0.3,
parse = TRUE
) +
# Region labels
geom_text(
data = df_labels_n,
aes(label = lab),
color = "white",
size = 4
) +
# Axes
scale_x_log10(
breaks = c(0.01, 0.1, 1, 10, 100),
labels = c("0.01", "0.1", "1", "10", "100"),
expand = c(0, 0)
) +
scale_y_continuous(expand = c(0, 0)) +
# Colorbar
scale_fill_viridis_c(
option = "viridis",
name = expression(italic(p)[thr]^{(n)})
) +
# Axis labels
labs(
x = expression("Cost ratio " * italic(R) == italic(C)[I] / italic(C)[V]),
y = expression("Reactive effectiveness " * italic(r))
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12),
legend.key.height = unit(1.5, "cm")
)
\(f=0.3\)
# Parameters
f_val <- 0.3 # capacity fraction
p_val <- 0.6 # mean outbreak probability (not directly used here)
r_vals <- seq(0, 1, length.out = 200)
R_vals <- 10 ^ seq(-2, 2, length.out = 200)
# --- 1. Heatmap data: p_thr^(n) across (R, r) ---
df_heatmap_n <- expand.grid(R = R_vals, r = r_vals) %>%
mutate(
p_thr_n = mapply(
function(R_single, r_single) {
p_star_multi(R = R_single, r = r_single, f = f_val)
},
R, r
)
)
# --- 2. Multiple boundary curves for p_thr^(n) ---
p_levels <- seq(0.2, 0.7, by = 0.1)
df_boundary_n <- lapply(p_levels, function(p_target) {
tmp <- lapply(R_vals, function(RR) {
f_root <- function(r) p_star_multi(R = RR, r = r, f = f_val) - p_target
sol <- tryCatch(
uniroot(f_root, interval = c(0, 1)),
error = function(e) NULL
)
if (is.null(sol)) return(NULL)
data.frame(
R = RR,
r = sol$root,
p_target = p_target
)
}) %>% bind_rows()
tmp
}) %>% bind_rows()
# --- 3. Midpoints along each boundary curve for annotation ---
df_boundary_labels <- df_boundary_n %>%
group_by(p_target) %>%
slice(round(n() / 2)) %>%
ungroup() %>%
mutate(
label = sprintf("p[thr]^{(n)}==%.1f", p_target)
)
# --- 4. Region labels (optional) ---
R_pre <- 10
r_pre <- 0.2
R_reac <- 10
r_reac <- 0.8
df_labels_n <- data.frame(
R = c(R_pre, R_reac),
r = c(r_pre, r_reac),
lab = c("pre-emptive", "reactive")
)
# --- 5. Heatmap with multiple boundary lines (all same style) ---
ggplot(df_heatmap_n, aes(x = R, y = r)) +
# Heatmap
geom_raster(aes(fill = p_thr_n)) +
# Boundary curves (same color, no legend)
geom_line(
data = df_boundary_n,
aes(x = R, y = r, group = p_target),
color = "firebrick",
linewidth = 0.8
) +
# Line labels
geom_text(
data = df_boundary_labels,
aes(x = R, y = r, label = label),
color = "firebrick",
size = 3.5,
vjust = -0.3,
parse = TRUE
) +
# Region labels
geom_text(
data = df_labels_n,
aes(label = lab),
color = "white",
size = 4
) +
# Axes
scale_x_log10(
breaks = c(0.01, 0.1, 1, 10, 100),
labels = c("0.01", "0.1", "1", "10", "100"),
expand = c(0, 0)
) +
scale_y_continuous(expand = c(0, 0)) +
# Colorbar
scale_fill_viridis_c(
option = "viridis",
name = expression(italic(p)[thr]^{(n)})
) +
# Axis labels
labs(
x = expression("Cost ratio " * italic(R) == italic(C)[I] / italic(C)[V]),
y = expression("Reactive effectiveness " * italic(r))
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12),
legend.key.height = unit(1.5, "cm")
)
Because \(c_{\mathrm{react}}^{(n)}\) is piecewise, \(p_{\mathrm{thr}}^{(n)}\) also has a piecewise form.
Here the equality occurs in the reactive-limited regime, so we set \[ c_{\mathrm{pre}}^{(n)} = f + (1-f)pR, \qquad c_{\mathrm{react}}^{(n)} = f + (p - fr)R. \]
Setting them equal and simplifying: \[ p_{\mathrm{thr}}^{(n)} = r. \] In this regime, the threshold, \(p_{\mathrm{thr}}^{(n)}\) , is independent of \(R\): only the relative size of \(f\) and \(p\) matters.
Here the equality occurs in the reactive-rich regime, so we set \[ c_{\mathrm{pre}}^{(n)} = f + (1-f)pR, \qquad c_{\mathrm{react}}^{(n)} = p\bigl[1 + (1-r)R\bigr]. \]
Equating and solving for \(p\):
\[ p_{\mathrm{thr}}^{(n)} = \frac{f}{1 + R(f - r)}. \]
Consistency with the reactive-rich regime requires \(p_{\mathrm{thr}}^{(n)} \le f\), which holds whenever \(f \ge r\) and \(R > 0\).
Putting the two regimes together: \[ p_{\mathrm{thr}}^{(n)}(f,r,R) = \begin{cases} r, & f \le r,\\[8pt] \dfrac{f}{1 + R(f - r)}, & f > r. \end{cases} \]
Interpretation:
We now allow the fixed capacity \(f\) to be split between pre-emptive and reactive use.
Let
subject to the capacity constraint
\[ f_{\mathrm{pre}} + f_{\mathrm{react}} = f. \]
We parameterize by a mixing parameter \(\alpha \in [0,1]\):
\(\alpha = 1\) means the pure pre-emptive strategy whereas \(\alpha = 0\) means the pure reactive strategy.
The pre-emptive component is straightforward:
\[ C_{\mathrm{pre\text{-}part}}^{(n)} = f_{\mathrm{pre}} C_{\mathrm{V}} = \alpha f C_{\mathrm{V}}. \]
The remaining fraction \((1 - f_{\mathrm{pre}}) = (1 - \alpha f)\) is not pre-emptively vaccinated. Among these:
Two regimes arise, depending on whether reactive capacity is rich or limited among non-pre-emptive populations:
\[ f_{\mathrm{react}} \ge p(1 - f_{\mathrm{pre}}) \quad\Longleftrightarrow\quad (1-\alpha)f \ge p(1-\alpha f). \]
\[ f_{\mathrm{react}} < p(1 - f_{\mathrm{pre}}) \quad\Longleftrightarrow\quad (1-\alpha)f < p(1-\alpha f). \]
If \((1-\alpha)f \ge p(1-\alpha f)\), then every outbreak in the non-pre-emptive group can receive a reactive campaign.
For non-pre-emptive populations:
Thus the reactive component is
\[ C_{\mathrm{react\text{-}part}}^{(n)} = (1-\alpha f)\, p\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr). \]
Total per-population cost for the mixed strategy is:
\[ C_{\mathrm{mixed}}^{(n)} = \alpha f C_{\mathrm{V}} + (1-\alpha f)\, p\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr). \]
Normalizing,
\[ c_{\mathrm{mixed}}^{(n)}(\alpha) = \frac{C_{\mathrm{mixed}}^{(n)}}{C_{\mathrm{V}}} = \alpha f + (1-\alpha f)\, p\bigl[1 + (1-r)R\bigr], \quad \text{if } (1-\alpha)f \ge p(1-\alpha f). \]
If \((1-\alpha)f < p(1-\alpha f)\), there are more outbreaks than available reactive campaigns in the non-pre-emptive group. The fraction of such outbreaks that receive a campaign is
\[ q = \frac{f_{\mathrm{react}}}{p(1-f_{\mathrm{pre}})} = \frac{(1-\alpha)f}{p(1-\alpha f)}. \]
Among non-pre-emptive populations that experience an outbreak:
The reactive component is
\[ \begin{aligned} C_{\mathrm{react\text{-}part}}^{(n)} &= (1-f_{\mathrm{pre}})p \Bigl[ q\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr) + (1-q)C_{\mathrm{I}} \Bigr] \\ &= f_{\mathrm{react}}\bigl(C_{\mathrm{V}} + (1-r)C_{\mathrm{I}}\bigr) + \bigl[(1-f_{\mathrm{pre}})p - f_{\mathrm{react}}r \bigr] C_{\mathrm{I}} \\ &= f_{\mathrm{react}} C_{\mathrm{V}} + \bigl[(1-f_{\mathrm{pre}})p - f_{\mathrm{react}}r \bigr] C_{\mathrm{I}}. \end{aligned} \]
Total per-population cost:
\[ \begin{aligned} C_{\mathrm{mixed}}^{(n)} &= f_{\mathrm{pre}} C_{\mathrm{V}} + f_{\mathrm{react}} C_{\mathrm{V}} + \bigl[(1-f_{\mathrm{pre}})p - f_{\mathrm{react}}r \bigr] C_{\mathrm{I}}\\ &= f C_{\mathrm{V}} + \bigl[(1-f_{\mathrm{pre}})p - f_{\mathrm{react}}r \bigr] C_{\mathrm{I}}. \end{aligned} \]
Substituting \(f_{\mathrm{pre}} = \alpha f\), \(f_{\mathrm{react}} = (1-\alpha)f\) and normalizing,
\[ c_{\mathrm{mixed}}^{(n)}(\alpha) = \frac{C_{\mathrm{mixed}}^{(n)}}{C_{\mathrm{V}}} = f + R\Bigl[(1-\alpha f)p - (1-\alpha)f r \Bigr], \quad \text{if } (1-\alpha)f < p(1-\alpha f). \]
An equivalent form is
\[ c_{\mathrm{mixed}}^{(n)}(\alpha) = f + R\Bigl[p - f r + \alpha f(r-p)\Bigr], \quad \text{if } (1-\alpha)f < p(1-\alpha f), \]
which makes the linear dependence on \(\alpha\) explicit.
Let \(\alpha^*\) denote the optimal fraction of vaccine capacity \(f\) allocated to pre-emptive vaccination. The optimal strategy depends on whether vaccine capacity is scarce (\(f \le p\)) or abundant (\(f > p\)).
When capacity is insufficient to cover the expected outbreak size, the strategy is a simple corner solution determined entirely by the reactive effectiveness \(r\):
\[ \alpha^* = \begin{cases} 0 \quad (\text{pure reactive}), & \text{if } r > p \\ 1 \quad (\text{pure pre-emptive}), & \text{if } r < p \\ \text{any } \alpha \in [0,1], & \text{if } r = p \end{cases} \]
When capacity exceeds the expected outbreak size, the optimal strategy depends on the cost ratio \(R\) relative to the threshold \(R_{\mathrm{thr}}\):
\[ \alpha^* = \begin{cases} 1 \quad (\text{pure pre-emptive}), & \text{if } r < p \\ 0 \quad (\text{pure reactive}), & \text{if } r > p \text{ and } R < R_{\mathrm{thr}} \\ \frac{f - p}{f(1 - p)} \quad (\text{mixed strategy}), & \text{if } r > p \text{ and } R \ge R_{\mathrm{thr}} \end{cases} \]
Note: In the abundant regime with high outbreak costs (\(R \ge R_{\mathrm{thr}}\)), the mixed strategy \(\alpha^* = \frac{f - p}{f(1 - p)}\) ensures that pre-emptive vaccination is maximized just enough so that the remaining capacity can fully cover the remaining risk reactively.
# Mixed strategy cost c_mix^{(n)}(alpha)
cost_mix_multi <- function(alpha, p, R, r, f) {
alpha <- pmax(pmin(alpha, 1), 0)
f_pre <- alpha * f
f_react <- (1 - alpha) * f
cond_rich <- (1 - alpha) * f >= p * (1 - f_pre)
c_rich <- f_pre + (1 - f_pre) * p * (1 + (1 - r) * R)
c_lim <- f + R * ((1 - f_pre) * p - f_react * r)
ifelse(cond_rich, c_rich, c_lim)
}
opt_alpha_multi <- function(p, R, r, f, grid_len = 1001) {
alpha_grid <- seq(0, 1, length.out = grid_len)
costs <- cost_mix_multi(alpha = alpha_grid,
p = p, R = R, r = r, f = f)
idx_min <- which.min(costs)
list(
alpha_star = alpha_grid[idx_min],
cost_star = costs[idx_min],
alpha_grid = alpha_grid,
cost_grid = costs
)
}
## ============================================================
## 1. Analytic solution for alpha* (equal-risk, large-n)
## ============================================================
alpha_star_equal_analytic <- function(p, R, r, f, tol = 1e-6) {
stopifnot(p > 0, p < 1, f > 0, f <= 1, R >= 0, r >= 0, r <= 1)
regime_capacity <- if (f <= p + tol) "scarce" else "abundant"
# Handle degenerate case r = p explicitly
if (abs(r - p) < tol) {
# Many alphas can be optimal; we pick a convenient representative
return(list(
alpha_star = 0, # arbitrary choice in [0,1]
regime_capacity = regime_capacity,
regime_detail = "degenerate_r_eq_p",
R_star = NA_real_,
alpha_bar = if (regime_capacity == "abundant") (f - p) / (f * (1 - p)) else NA_real_
))
}
if (regime_capacity == "scarce") {
# f <= p
if (r > p) {
alpha_star <- 0 # pure reactive
regime_detail <- "scarce_pure_reactive"
} else { # r < p
alpha_star <- 1 # pure pre-emptive
regime_detail <- "scarce_pure_preemptive"
}
return(list(
alpha_star = alpha_star,
regime_capacity = regime_capacity,
regime_detail = regime_detail,
R_star = NA_real_,
alpha_bar = NA_real_
))
}
# Abundant capacity: f > p
alpha_bar <- (f - p) / (f * (1 - p)) # regime boundary
alpha_bar <- max(min(alpha_bar, 1), 0) # safety clamp
if (r < p) {
# Pure pre-emptive
alpha_star <- 1
regime_detail <- "abundant_pure_preemptive"
R_star <- NA_real_
} else {
# r > p
R_star <- (1 - p) / (p * (1 - r))
if (R < R_star) {
alpha_star <- 0
regime_detail <- "abundant_pure_reactive"
} else {
alpha_star <- alpha_bar
regime_detail <- "abundant_mixed_at_boundary"
}
}
list(
alpha_star = alpha_star,
regime_capacity = regime_capacity,
regime_detail = regime_detail,
R_star = if (exists("R_star")) R_star else NA_real_,
alpha_bar = alpha_bar
)
}p <- 0.3
R <- 10
r <- 0.6
f <- 0.5
opt_res <- opt_alpha_multi(p = p, R = R, r = r, f = f,
grid_len = 1001)
df <- tibble::tibble(
alpha = opt_res$alpha_grid,
cost = opt_res$cost_grid
)
alpha_star <- opt_res$alpha_star
alpha_star_anal <- alpha_star_equal_analytic(p = p, R = R, r = r, f = f, tol = 1e-6)
ggplot(df, aes(x = alpha, y = cost)) +
geom_line(linewidth = 1) +
geom_vline(
xintercept = alpha_star,
linetype = "dashed",
color = "firebrick",
linewidth = 0.9
) +
labs(
title = bquote(italic(p) == .(p) ~ "," ~
italic(r) == .(r) ~ "," ~
italic(R) == .(R) ~ "," ~
italic(f) == .(f)),
x = expression("Fraction of capacity pre-emptive " ~ italic(alpha)),
y = expression("Normalized per-population expected cost " ~ italic(c)[mixed]^(italic(n)))
) +
geom_vline(
xintercept = alpha_star_anal$alpha_star,
linetype = "dotted",
color = "firebrick",
linewidth = 0.9
) +
annotate(
"text",
x = alpha_star_anal$alpha_star,
y = max(df$cost, na.rm = TRUE),
label = "alpha^'*'",
parse = TRUE,
hjust = -0.1,
vjust = 1.2,
color = "firebrick",
size = 4
) +
theme_light()
3D surface of \(\alpha^*(R,r)\):
titlefontsize <- 34
tickfontsize <- 17
p_val <- 0.3
f_val <- 0.6
r_vals <- seq(0, 1, length.out = 100)
R_vals <- 10 ^ seq(-2, 2, length.out = 100)
alpha_mat <- outer(
R_vals, r_vals,
Vectorize(function(R, r) {
opt_alpha_multi(
p = p_val,
R = R,
r = r,
f = f_val,
grid_len = 301
)$alpha_star
})
)
alpha_mat_t <- t(alpha_mat)
camera <- list(
eye = list(x = 0.1, y = 0.08, z = 0.1)
)
fig_alpha_surface <-
plot_ly() %>%
add_surface(
x = R_vals,
y = r_vals,
z = alpha_mat_t,
colorscale = "Viridis",
showscale = TRUE,
colorbar = list(
title = list(
text = "<i>α</i><sup>*</sup>",
font = list(size = titlefontsize)
)
)
) %>%
layout(
scene = list(
camera = camera,
yaxis = list(
title = "<i>r</i>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)),
xaxis = list(
title = "<i>R</i> = C<sub>I</sub> / C<sub>V</sub>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize),
type = "log",
tickvals = c(0.01, 0.1, 1, 10, 100),
ticktext = c("0.01", "0.1", "1", "10", "100")
),
zaxis = list(
title = "<i>α</i><sup>*</sup>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)
)
)
)
fig_alpha_surfaceWe determine the optimal fraction of capacity \(\alpha^*\) to allocate pre-emptively. This depends on the total capacity regime (\(f\) vs mean risk \(\bar{p}\)) and the specific parameters (\(r, R, p_H, p_L\)).
Let \(\bar{p} = f_H p_H + f_L p_L\) be the mean outbreak probability.
In this regime, capacity is insufficient to cover all expected outbreaks. The system is Reactive-Limited. Decisions are driven by comparing reactive effectiveness (\(r\)) directly to risk probabilities (\(p_H, p_L\)).
\[ \alpha^* = \begin{cases} 1 \quad (\text{pure pre-emptive}), & \text{if } r < p_L \\ \min\left(1, \frac{f_H}{f}\right) \quad (\text{focused on high-risk}), & \text{if } p_L \le r \le p_H \\ 0 \quad (\text{pure reactive}), & \text{if } r > p_H \end{cases} \]
Interpretation: \(r < p_L\): Reactive vaccination is weak. Pre-empt everyone possible. \(p_L \le r \le p_H\): Reactive is better than pre-empting low-risk, but worse than pre-empting high-risk. Prioritize High-Risk. \(r > p_H\): Reactive is highly effective. Wait and respond to outbreaks.
In this regime, capacity exceeds expected outbreaks. The system is initially reactive-rich. Pre-emptive vaccination is only justified if the cost of outbreaks (\(R\)) is high enough to outweigh the efficiency of waiting.
We define two cost thresholds: \[ R_H^* = \frac{1 - p_H}{p_H(1-r)}, \quad R_L^* = \frac{1 - p_L}{p_L(1-r)} \]
(Note: since \(p_H > p_L\), it implies \(R_H^* < R_L^*\). It is “cheaper” to justify protecting high-risk groups.)
\[ \alpha^* = \begin{cases} 0 \quad (\text{pure reactive}), & \text{if } R < R_H^* \\ \min\left(1, \frac{f_H}{f}\right) \quad (\text{focused on high-risk}), & \text{if } R_H^* \le R < R_L^* \\ 1 \quad (\text{pure pre-emptive}), & \text{if } R \ge R_L^* \end{cases} \]
Interpretation: \(R < R_H^*\): Outbreaks are cheap. Use pure reactive strategy. \(R_H^* \le R < R_L^*\): Outbreaks are costly enough to pre-empt high-risk, but low-risk is still effectively managed reactively. \(R \ge R_L^*\): Outbreaks are catastrophic. Pre-empt everyone possible.
The intermediate strategy \(\alpha = \min(1, f_H/f)\) is the “Focused” strategy where high-risk populations are prioritized. The total expected cost \(c(\alpha)\) varies based on whether the specific capacity constraints for high-risk and low-risk groups are met.
| Capacity Scenario | Condition | Optimal \(\alpha\) | Total Cost \(c(\alpha)\) |
|---|---|---|---|
| Insufficient for High-Risk | \(f < f_H\) | \(1\) | \(f + R \big[ (f_H - f)p_H + f_L p_L \big]\) |
| Sufficient for H, Limited for L | \(f_H \le f < f_H + f_L p_L\) | \(\frac{f_H}{f}\) | \(f + R \big[ f_L p_L - (f - f_H)r \big]\) |
| Sufficient for H, Abundant for L | \(f \ge f_H + f_L p_L\) | \(\frac{f_H}{f}\) | \(f_H + f_L p_L \big[ 1 + (1-r)R \big]\) |
Breakdown: 1. Insufficient: We use all \(f\) on High-Risk. Zero reactive capacity. Cost is vaccination + all unprevented outbreaks in both groups. 2. Limited for L: We saturate high-risk (\(f_H\)). Remaining capacity (\(f-f_H\)) is used reactively on low-risk but is exhausted. Cost is total capacity \(f\) + net unprevented outbreaks in Low-Risk. 3. Abundant for L: We saturate high-risk (\(f_H\)). Remaining capacity is sufficient to cover all Low-Risk outbreaks. Cost is pre-emptive (\(f_H\)) + reactive campaigns used (\(f_L p_L\)) + failed reactive attempts.
We extend the equal-risk \(n\)-population model to allow heterogeneous outbreak probabilities across populations and to incorporate targeting accuracy for pre-emptive vaccination. Costs are normalized by \(C_{\mathrm{V}}\), and the superscript \((n)\) denotes the \(n\)-population setting.
The pre-emptive and reactive fractions satisfy
\[
f_{\mathrm{pre}} + f_{\mathrm{react}} = f.
\]
Assume each population has a latent outbreak hazard
\[ \Lambda_i \sim \mathrm{Exponential}(\theta), \qquad \theta > 0. \]
The outbreak probability for population \(i\) is
\[
p_i = 1 - e^{-\Lambda_i}, \qquad p_i \in (0,1).
\]
Using the transformation \(p = 1 - e^{-\lambda}\), we obtain the density
\[
f_P(p) = \theta (1-p)^{\theta - 1}, \qquad 0 < p < 1,
\] so the outbreak-probability random variable
\[
P \sim \mathrm{Beta}(1,\theta).
\]
Key properties:
Mean outbreak probability
\[
p_{\mathrm{mean}} = \mathbb{E}[P] = \frac{1}{1+\theta}.
\]
If \(\theta < 1\): mass is concentrated near \(p = 1\) (many high-risk populations).
If \(\theta = 1\): \(P \sim \mathrm{Uniform}(0,1)\).
If \(\theta > 1\): risks are skewed toward 0.
We treat \(p_i\) as i.i.d. draws from \(\mathrm{Beta}(1,\theta)\).
Under perfect ranking, pre-emptive vaccination targets the populations with the highest outbreak probabilities.
Let
\[
q = f_{\mathrm{pre}} = \alpha f
\] be the pre-emptively vaccinated population fraction.
Let \(p_{\mathrm{cut}}(q)\) satisfy
\[
\Pr(P \ge p_{\mathrm{cut}}(q)) = q.
\]
The CDF of \(P \sim \mathrm{Beta}(1,\theta)\) is
\[
F_P(p) = 1 - (1-p)^\theta,
\] so
\[
\Pr(P \ge t) = (1-t)^\theta.
\]
Thus,
\[
q = (1 - p_{\mathrm{cut}}(q))^\theta
\quad\Longrightarrow\quad
p_{\mathrm{cut}}(q) = 1 - q^{1/\theta}.
\]
Define
\[
p_{\mathrm{pre}}(q) = \mathbb{E}[P \mid P \ge p_{\mathrm{cut}}(q)].
\]
Direct calculation yields
\[
p_{\mathrm{pre}}(q)
= 1 - \frac{\theta}{\theta + 1}\, q^{1/\theta},
\qquad 0 < q \le 1.
\]
Checks:
If \(q = 1\):
\[
p_{\mathrm{pre}}(1) = \frac{1}{1+\theta} = p_{\mathrm{mean}}.
\]
As \(q \to 0\):
\[
p_{\mathrm{pre}}(q) \to 1.
\]
Let
\[
p_{\mathrm{rem}}(q) = \mathbb{E}[P \mid P < p_{\mathrm{cut}}(q)].
\]
Using the law of total expectation:
\[
p_{\mathrm{mean}}
= q\,p_{\mathrm{pre}}(q)
+ (1-q)\,p_{\mathrm{rem}}(q),
\] so
\[
p_{\mathrm{rem}}(q)
= \frac{p_{\mathrm{mean}} - q\,p_{\mathrm{pre}}(q)}{1-q}.
\]
Substituting \(p_{\mathrm{mean}} = 1/(\theta+1)\) and the expression above for \(p_{\mathrm{pre}}(q)\) gives
\[ p_{\mathrm{rem}}(q) = \frac{ \dfrac{1}{\theta+1} - q + \dfrac{\theta}{\theta+1} q^{1 + 1/\theta} }{1 - q}. \]
Using \(q = f_{\mathrm{pre}} = \alpha f\), we write
\[
p_{\mathrm{pre}}(\alpha f)
= 1 - \frac{\theta}{\theta+1}(\alpha f)^{1/\theta},
\] \[
p_{\mathrm{rem}}(\alpha f)
=
\frac{
\dfrac{1}{\theta+1}
- \alpha f
+ \dfrac{\theta}{\theta+1} (\alpha f)^{1 + 1/\theta}
}
{1 - \alpha f}.
\]
Normalized expected cost:
\[ c_{\mathrm{pre}}^{(n)}(\alpha) = f_{\mathrm{pre}} = \alpha f. \]
Fraction:
\[ 1 - f_{\mathrm{pre}} = 1 - \alpha f. \]
Mean outbreak probability:
\[ p' = p_{\mathrm{rem}}(\alpha f). \]
Reactive capacity:
\[
f_{\mathrm{react}} = (1-\alpha)f.
\]
Effective capacity per remaining population:
\[
f' = \frac{f_{\mathrm{react}}}{1 - f_{\mathrm{pre}}}
= \frac{(1-\alpha)f}{1 - \alpha f}.
\]
\[ c_{\mathrm{rem}}^{(n)}(p',f') = f' + \bigl(p' - fr' \bigr)R. \]
Total cost:
\[
c_{\mathrm{mixed}}^{(n)}(\alpha)
= f + R\Bigl[(1-\alpha f)p' - (1-\alpha)f r \Bigr].
\]
\[ c_{\mathrm{rem}}^{(n)}(p',f') = p'\,\bigl[1 + (1-r)R\bigr]. \]
\[ c_{\mathrm{mixed}}^{(n)}(\alpha) = \alpha f + (1-\alpha f)\,p'\,\bigl[1 + (1-r)R\bigr]. \]
Unlike the equal-risk case, \(p'\) is nonlinear, so \(c_{\mathrm{mixed}}^{(n)}(\alpha)\) is not piecewise-linear, and no closed-form optimal \(\alpha^*\) exists.
We therefore compute \(\alpha^*\) by 1D minimization: \[ \alpha^* = \arg\min_{\alpha\in[0,1]} c_{\mathrm{mixed}}^{(n)}(\alpha). \]
For incomplete ranking, we interpolate between perfect ranking (\(\rho = 1\)) and random targeting (\(\rho = 0\)):
\[ p_{\mathrm{pre}}(\rho, q) = p_{\mathrm{mean}} + \rho\bigl(p_{\mathrm{pre}}^{*}(q) - p_{\mathrm{mean}}\bigr), \]
$$
p_{\mathrm{rem}}(\rho, q)= p_{} + (p_{}^{*}(q) - p_{}), $$
where $p_{\mathrm{pre}}^{*}(q)$ and $p_{\mathrm{rem}}^{*}(q)$ correspond to $\rho = 1$.
- $\rho = 1$: perfect targeting
- $\rho = 0$: random targeting
- $0 < \rho < 1$: partial enrichment
Mixed-cost formulas remain unchanged, replacing $p'$ with $p_{\mathrm{rem}}(\rho,\alpha f)$.As \(\rho\) decreases, enrichment weakens and the optimal pre-emptive share \(\alpha^*(\rho)\) shifts toward 0.
The functions below implement:
mu_pre_star_beta1theta <- function(q, theta) {
1 - theta / (theta + 1) * q^(1/theta)
}
mu_rem_star_beta1theta <- function(q, theta) {
mu <- 1 / (1 + theta)
if (q <= 0) return(mu)
if (q >= 1) return(NA_real_)
mu_pre <- mu_pre_star_beta1theta(q, theta)
(mu - q * mu_pre) / (1 - q)
}
mu_pre_beta1theta_rho <- function(q, theta, rho) {
mu <- 1 / (1 + theta)
mu_star <- mu_pre_star_beta1theta(q, theta)
mu + rho * (mu_star - mu)
}
mu_rem_beta1theta_rho <- function(q, theta, rho) {
mu <- 1 / (1 + theta)
mu_star <- mu_rem_star_beta1theta(q, theta)
mu + rho * (mu_star - mu)
}
c_mix_beta1theta <- function(alpha, f, r, R, p_mean, rho) {
theta <- 1 / p_mean - 1
q <- alpha * f
if (q < 0 || q >= 1) return(Inf)
frac_rem <- 1 - q
f_react <- (1 - alpha) * f
mu_rem <- mu_rem_beta1theta_rho(q, theta, rho)
f_prime <- f_react / frac_rem
p_prime <- mu_rem
c_pre <- q
if (f_prime < p_prime) {
c_rem_per_pop <- f_prime + (p_prime - f_prime * r) * R
} else {
c_rem_per_pop <- p_prime * (1 + (1 - r) * R)
}
c_rem_total <- frac_rem * c_rem_per_pop
c_pre + c_rem_total
}
opt_alpha_beta1theta <- function(f, r, R, p_mean, rho,
grid_len = 201) {
alpha_grid <- seq(0, 1, length.out = grid_len)
costs <- sapply(alpha_grid, function(a)
c_mix_beta1theta(a, f = f, r = r, R = R,
p_mean = p_mean, rho = rho)
)
idx_min <- which.min(costs)
list(
alpha_star = alpha_grid[idx_min],
alpha_grid = alpha_grid,
cost_grid = costs
)
}
cost_profile_df <- function(f, r, R, p_mean, rho_vals,
alpha_grid = seq(0, 1, by = 0.01)) {
purrr::map_dfr(rho_vals, function(rho_i) {
costs <- sapply(alpha_grid, function(a)
c_mix_beta1theta(
alpha = a,
f = f,
r = r,
R = R,
p_mean = p_mean,
rho = rho_i
)
)
alpha_star_i <- opt_alpha_beta1theta(
f = f,
r = r,
R = R,
p_mean = p_mean,
rho = rho_i
)$alpha_star
tibble::tibble(
alpha = alpha_grid,
cost = costs,
rho = rho_i,
alpha_star = alpha_star_i
)
})
}Example scan across parameters
R_vals <- c(5, 20, 50)
r_vals <- c(0.3, 0.6, 0.9)
f_vals <- c(0.1, 0.3, 0.6)
p_vals <- c(0.1, 0.3, 0.6)
rho_vals <- c(0, 0.5, 1.0)
param_grid <- expand.grid(
R = R_vals,
r = r_vals,
f = f_vals,
p_mean = p_vals,
rho = rho_vals,
KEEP.OUT.ATTRS = FALSE
)
results <- param_grid %>%
dplyr::mutate(
opt = purrr::pmap(
list(f, r, R, p_mean, rho),
~ opt_alpha_beta1theta(
f = ..1, r = ..2, R = ..3,
p_mean = ..4, rho = ..5,
grid_len = 201
)
),
alpha_star = purrr::map_dbl(opt, "alpha_star")
)Cost profiles and the \(\alpha^*\) by \(\rho\):
f_val <- 0.3
r_val <- 0.6
R_val <- 20
p_mean_val <- 0.4
# Title with italic symbols
title_expr <- bquote(
"Mixed strategy cost vs " * alpha ~
"(" *
italic(f) == .(f_val) * "," ~
italic(r) == .(r_val) * "," ~
italic(R) == .(R_val) * "," ~
italic(p)[mean] == .(p_mean_val) *
")"
)
# title_expr <- bquote(
# "Mixed strategy cost vs " * alpha ~
# "(" *
# italic(f) == .(f_val) * "," ~
# italic(r) == .(r_val) * "," ~
# italic(R) == .(R_val) * "," ~
# italic(p)[mean] == .(p_mean_val) *
# ")"
# )
rho_vals <- c(1, 0.5, 0)
df_cost <-
cost_profile_df(
f = f_val,
r = r_val,
R = R_val,
p_mean = p_mean_val,
rho_vals = rho_vals
) |>
dplyr::mutate(
rho_lab = factor(
rho,
levels = rho_vals,
labels = paste0("italic(rho)==", rho_vals)
)
)
# Vertical-line positions (one per facet)
df_vline <- df_cost |>
dplyr::distinct(rho_lab, alpha_star)
# Label positions for α*
df_label <- df_cost |>
dplyr::group_by(rho_lab) |>
dplyr::summarise(
alpha_star = dplyr::first(alpha_star),
y = max(cost, na.rm = TRUE),
.groups = "drop"
)
axis_title_size <- 18 # axis label size
axis_text_size <- 14 # tick label size
facet_title_size <- 16 # facet label size
plt <- ggplot(df_cost, aes(x = alpha, y = cost)) +
geom_line(linewidth = 1) +
geom_vline(
data = df_vline,
aes(xintercept = alpha_star),
linetype = "dashed",
color = "firebrick",
linewidth = 0.9
) +
geom_text(
data = df_label,
aes(x = alpha_star, y = max(y)),
label = "alpha^'*'",
parse = TRUE,
hjust = -0.1, # nudges label right of the line
vjust = 1.2, # slightly above the top
color = "firebrick",
size = 4
) +
facet_wrap(~ rho_lab, labeller = label_parsed) +
labs(
title = title_expr,
x = expression("Fraction of capacity pre-emptive " * alpha),
y = expression(italic(c)[mix]^"(n)")
) +
scale_x_continuous(
breaks = seq(0, 1, by = 0.25),
labels = function(x) sub(".00$", "", sprintf("%.2f", x))
)+
theme_light(base_size=14)+
theme(
axis.title.x = element_text(size = axis_title_size),
axis.title.y = element_text(size = axis_title_size),
axis.text.x = element_text(size = axis_text_size),
axis.text.y = element_text(size = axis_text_size),
strip.text = element_text(size = facet_title_size)
)
plt
3D surface of \(\alpha^*\)
f_val <- 0.3 # vaccination capacity fraction
theta <- 4 # Beta(1, theta) heterogeneity (choose as desired)
rho_val <- 0.6 # targeting accuracy rho
p_mean <- 1 / (1 + theta)
r_vals <- seq(0, 1, length.out = 100)
R_vals <- 10 ^ seq(-2, 2, length.out = 100)
alpha_star_multi <- function(R, r, f, p_mean, rho) {
opt_alpha_beta1theta(
f = f,
r = r,
R = R,
p_mean = p_mean,
rho = rho,
grid_len = 301
)$alpha_star
}
alpha_mat <- outer(
R_vals, r_vals,
Vectorize(function(R, r) {
alpha_star_multi(
R = R,
r = r,
f = f_val,
p_mean = p_mean,
rho = rho_val
)
})
)
alpha_mat_t <- t(alpha_mat)
camera <- list(eye = list(x = 0.1, y = 0.08, z = 0.1))
titlefontsize <- 34
tickfontsize <- 17
fig_alpha_star_multi <-
plot_ly() %>%
add_surface(
x = R_vals,
y = r_vals,
z = alpha_mat_t,
colorscale = "Viridis",
showscale = TRUE,
colorbar = list(
title = list(
text = "<i>α</i><sup>*</sup>",
font = list(size = titlefontsize)
)
)
) %>%
layout(
scene = list(
camera = camera,
yaxis = list(
title = "<i>r</i>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)
),
xaxis = list(
title = "<i>R</i> = C<sub>I</sub>/C<sub>V</sub>",
type = "log",
tickvals = c(0.01, 0.1, 1, 10, 100),
ticktext = c("0.01", "0.1", "1", "10", "100"),
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)
),
zaxis = list(
title = "<i>α</i><sup>*</sup>",
titlefont = list(size = titlefontsize),
tickfont = list(size = tickfontsize)
)
)
)
fig_alpha_star_multiOptimal \(\alpha^*\) satisfies: \[ \alpha^* = \arg\min_{\alpha \in [0,1]} c_{\mathrm{mixed}}^{(n)}(\alpha;\rho), \] considering both regime boundaries and interior solutions.
Let \[ p'(\alpha) := \mu_{\mathrm{rem}}(\rho,\alpha f), \qquad p'_{\alpha}(\alpha) := \frac{d}{d\alpha} p'(\alpha). \]
The derivative is: \[ \frac{\partial}{\partial \alpha} c_{\mathrm{mix}}^{(n)}(\alpha;\rho) = R \Bigl[ f r - f p'(\alpha) + (1-\alpha f)\, p'_{\alpha}(\alpha) \Bigr]. \]
Setting to zero: $$
f r - f p’(_*) + (1-** f), p’*{}(_*) = 0.
$$
Equivalently, define the effective risk threshold \[ p_{\mathrm{eff}}^{\mathrm{RL}}(\alpha;\rho,\theta) = p'(\alpha) - \frac{1-\alpha f}{f}\, p'_{\alpha}(\alpha). \]
Then interior optimality implies: $$
r = p_{}^{}(_*;,).
$$
Derivative: \[ \frac{\partial}{\partial \alpha} c_{\mathrm{mix}}^{(n)}(\alpha;\rho) = f + \Bigl[-f p'(\alpha) + (1-\alpha f) p'_{\alpha}(\alpha)\Bigr] \bigl[1 + (1-r)R\bigr]. \]
Interior optimum satisfies: \[ f + \Bigl[-f p'(\alpha_*) + (1-\alpha_* f) p'_{\alpha}(\alpha_*)\Bigr] \bigl[1 + (1-r)R\bigr] = 0. \]
This may be written as: \[ r = p_{\mathrm{eff}}^{\mathrm{RR}}(\alpha_*;\rho,\theta,R). \]
As in the equal-risk case, optimal \(\alpha_*\) often occurs at boundaries:
If \[ f'(\alpha) < p'(\alpha) \quad \forall \alpha, \] the entire domain is reactive-limited, the derivative has the form
\[ \frac{\partial c_{\mathrm{mix}}^{(n)}}{\partial \alpha} = f R \bigl[r - p_{\mathrm{eff}}^{\mathrm{RL}}(\alpha)\bigr]. \]
Thus: - If \(r > p_{\mathrm{eff}}^{\mathrm{RL}}(\alpha)\) for all \(\alpha\): pure reactive, \(\alpha_* = 0\). - If \(r < p_{\mathrm{eff}}^{\mathrm{RL}}(\alpha)\) for all \(\alpha\): pure pre-emptive, \(\alpha_* = 1\). - Otherwise: unique interior \(\alpha_*\) solving \(r = p_{\mathrm{eff}}^{\mathrm{RL}}(\alpha_*;\rho,\theta)\).
If \[ f'(\alpha) \ge p'(\alpha) \quad \forall \alpha, \] then the derivative equality \[ \frac{\partial c_{\mathrm{mix}}^{(n)}}{\partial \alpha} = 0 \] gives the interior condition \[ r = p_{\mathrm{eff}}^{\mathrm{RR}}(\alpha_*;\rho,\theta,R). \]
Corner solutions occur if the derivative’s sign is constant.
If \(f'(\alpha)\) crosses \(p'(\alpha)\), candidates for \(\alpha^*\) include:
The optimal solution is whichever yields the smallest \(c_{\mathrm{mixed}}^{(n)}\).
Then \(p'(\alpha) = \mu\) and \(p'_\alpha(\alpha) = 0\). The heterogeneous model collapses to the equal-risk model with \(p=\mu\).
In the reactive-limited regime: \[ \frac{\partial c_{\mathrm{mix}}^{(n)}}{\partial \alpha} = fR(r - \mu), \] a constant.
Thus: \[ \alpha_* = \begin{cases} 0, & r > \mu,\\ 1, & r < \mu,\\ \text{any } \alpha, & r = \mu. \end{cases} \]
In the reactive-rich regime, the cost is linear in \(\alpha\) and the solution depends on \(R\) exactly as in the equal-risk abundant-capacity case.
For \(\theta = 2\), \[ \mu_{\mathrm{rem}}^*(q) = \frac{\frac{1}{3} - q + \frac{2}{3} q^{3/2}}{1 - q}. \]
In the reactive-limited regime, the mixed cost becomes: \[ c_{\mathrm{mixed}}^{(n)}(q) = f + \frac{R}{3}(2q^{3/2} - 3q - 3r(f-q) + 1), \qquad q = \alpha f. \]
Differentiating: \[ \frac{d c}{d q} = 0 \quad\Longrightarrow\quad \sqrt{q_*} + r - 1 = 0. \]
Thus: $$ q_*^{} = (1-r)^2, _*^{} = ,
$$
subject to regime consistency (\(f'(\alpha_*) < p'(\alpha_*)\)) and \(0 \le \alpha_* \le 1\).
This illustrates how strong heterogeneity and perfect targeting yield explicit closed forms. For general \(\theta\) and \(\rho\), \(\alpha_*\) satisfies the implicit analytic FOCs above.
Heterogeneity and targeting modify the effective mean risk in the remaining group: \[ p'(\alpha;\rho,\theta) = \mu_{\mathrm{rem}}(\rho,\alpha f). \]
The mixed cost remains piecewise analytic; \(\alpha^*\) must be one of:
No targeting (\(\rho=0\)): collapses exactly to the equal-risk case.
Perfect targeting and fatter tails: increase the value of pre-emptive use; interior optimum may be explicit (e.g., \((1-r)^2/f\) for \(\mathrm{Beta}(1,2)\)).
To allow fatter right tails while remaining analytically tractable, we model the outbreak probability \(p_i\) as a mixture of two Beta components:
Formally, let \[ p_i \sim (1-\omega)\,\mathrm{Beta}(1,\theta_{\mathrm{L}}) + \omega\,\mathrm{Beta}(1,\theta_{\mathrm{H}}), \qquad 0 \le \omega \le 1,\quad 0 < \theta_{\mathrm{H}} < \theta_{\mathrm{L}}. \]
Here: - \(\omega\) is the mixture weight for the high–risk component, - \(\theta_{\mathrm{L}}\) controls the baseline risk distribution (mostly lower risks), - \(\theta_{\mathrm{H}}\) controls the high–risk tail (smaller \(\theta_{\mathrm{H}}\) \(\Rightarrow\) more mass near 1).
The resulting density on \((0,1)\) is \[ f(p) = (1-\omega)\,\theta_{\mathrm{L}}(1-p)^{\theta_{\mathrm{L}}-1} + \omega\,\theta_{\mathrm{H}}(1-p)^{\theta_{\mathrm{H}}-1}, \qquad 0 < p < 1. \]
The mean outbreak probability remains explicit: \[ \mu := \mathbb{E}[p_i] = (1-\omega)\,\frac{1}{1+\theta_{\mathrm{L}}} + \omega\,\frac{1}{1+\theta_{\mathrm{H}}}. \]
If we wish to parameterize the model by a target mean \(p_{\mathrm{mean}}\) and a “tail–heaviness” configuration \((\omega,\theta_{\mathrm{H}})\), we can set \[ p_{\mathrm{mean}} = \mu \quad\Longrightarrow\quad \frac{1}{1+\theta_{\mathrm{L}}} = \frac{p_{\mathrm{mean}} - \omega/(1+\theta_{\mathrm{H}})}{1-\omega}, \] which yields \[ \theta_{\mathrm{L}} = \frac{1-\omega}{p_{\mathrm{mean}} - \dfrac{\omega}{1+\theta_{\mathrm{H}}}} - 1, \] provided \(p_{\mathrm{mean}} > \omega/(1+\theta_{\mathrm{H}})\) so that the baseline component is well–defined.
Under perfect targeting (\(\rho = 1\)), we again vaccinate the highest–risk fraction \[ q = f_{\mathrm{pre}} = \alpha f \] of populations. There exists a threshold \(t(q)\) such that \[ \Pr(p_i \ge t(q)) = q. \] In the mixture model, \[ \Pr(p_i \ge t) = (1-\omega)(1-t)^{\theta_{\mathrm{L}}} + \omega(1-t)^{\theta_{\mathrm{H}}}, \] so \(t(q)\) is defined implicitly by \[ q = (1-\omega)\bigl(1-t(q)\bigr)^{\theta_{\mathrm{L}}} + \omega\bigl(1-t(q)\bigr)^{\theta_{\mathrm{H}}}. \]
There is generally no closed–form for \(t(q)\), but it is one–dimensional and smooth, so \(t(q)\) can be obtained numerically for any given \((\omega,\theta_{\mathrm{L}},\theta_{\mathrm{H}})\).
Conditional means in the pre–emptive and remaining groups are still analytically expressible because each component is Beta:
\[ \mu_{\mathrm{pre}}(q) = \mathbb{E}[p_i \mid p_i \ge t(q)] = \frac{\displaystyle \int_{t(q)}^{1} p\,f(p)\,\mathrm{d}p}{q}, \] - Remaining group:
\[ \mu_{\mathrm{rem}}(q) = \mathbb{E}[p_i \mid p_i < t(q)] = \frac{\displaystyle \int_{0}^{t(q)} p\,f(p)\,\mathrm{d}p}{1-q}. \]
Each integral decomposes into a weighted sum of Beta–type integrals such as \[ \int p(1-p)^{\theta_k-1}\,\mathrm{d}p, \qquad k \in \{\mathrm{L},\mathrm{H}\}, \] which have simple closed forms in terms of powers of \((1-p)\).
Setting \(q = \alpha f\) as before, we obtain \[ \mu_{\mathrm{pre}}(\alpha f), \qquad \mu_{\mathrm{rem}}(\alpha f), \] and plug these into the same mixed–strategy cost expressions as in the \(\mathrm{Beta}(1,\theta)\) case.
Let \(f_{\mathrm{pre}} = \alpha f\) and \(f_{\mathrm{react}} = (1-\alpha)f\) as before. In the heterogeneous mixture setting:
\[ c_{\mathrm{pre\text{-}part}}^{(n)}(\alpha) = f_{\mathrm{pre}} = \alpha f, \]
\[ f' = \frac{f_{\mathrm{react}}}{1-f_{\mathrm{pre}}} = \frac{(1-\alpha)f}{1-\alpha f}. \]
We then reuse the equal–risk formulas on the remaining group:
\[ c_{\mathrm{mixed}}^{(n)}(\alpha) = f + R\Bigl[(1-\alpha f)\,\mu_{\mathrm{rem}}(\alpha f) - (1-\alpha)f r\Bigr]. \] - If \(f' \ge p'\) (reactive–rich among remaining populations),
\[ c_{\mathrm{mixed}}^{(n)}(\alpha) = \alpha f + (1-\alpha f)\,\mu_{\mathrm{rem}}(\alpha f)\bigl[1 + (1-r)R\bigr]. \]
As in the single–Beta case, the optimal pre–emptive fraction
\[ \alpha^* = \arg\min_{0 \le \alpha \le 1} c_{\mathrm{mix}}^{(n)}(\alpha) \]
characterizes the best split of capacity between pre–emptive and reactive use. The mixture mainly changes the mapping
\[ \alpha \;\mapsto\; \mu_{\mathrm{rem}}(\alpha f) \]
by enriching the high–risk tail through \((\omega,\theta_{\mathrm{H}})\), typically pushing \(\alpha^*\) upward (i.e., making pre–emptive vaccination more attractive) when a small subset of populations carries substantially higher outbreak risk.