ABSTRACT
During the coronavirus disease (COVID-19) pandemic, researchers attempted to estimate the number of averted and avertible outcomes due to non-pharmaceutical interventions or vaccination campaigns to quantify public health impact. However, the estimands used in these analyses have not been previously formalized. It is also unclear how these analyses relate to the broader framework of direct, indirect, total, and overall causal effects of an intervention under interference. In this study, using potential outcome notation, we adjust the direct and overall effects to accommodate analyses of averted and avertible outcomes. We use this framework to interrogate the commonly-held assumption in empirical studies that vaccine-averted outcomes via direct impact among vaccinated individuals (or vaccine-avertible outcomes via direct impact among unvaccinated individuals) is a lower bound of vaccine-averted (or avertible) outcomes overall. To do so, we describe a susceptible-infected-recovered-death model stratified by vaccination status. When vaccine efficacies wane, the lower bound fails for vaccine-avertible outcomes. When transmission or fatality parameters increase over time, the lower bound fails for both vaccine-averted and avertible outcomes. Only in the simplest scenario where vaccine efficacies, transmission, and fatality parameters are constant over time, outcomes averted via direct impact among vaccinated individuals (or outcomes avertible via direct impact among unvaccinated individuals) is shown to be a lower bound of overall impact on vaccine-averted (or avertible) outcomes. In conclusion, the lower bound can fail under common violations to assumptions on constant vaccine efficacy, pathogen properties, or behavioral parameters over time. In real data analyses, estimating what seems like a lower bound on overall impact through estimating direct impact may be inadvisable. By classifying estimands for averted and avertible outcomes and examining their relations, this study improves conduct and interpretation of research evaluating impact of infectious disease interventions.
1. INTRODUCTION
Recently, researchers have estimated the number of COVID-19 deaths (or infections) averted by vaccination campaigns in the United States,1,2 Israel,3,4 Chile,5 Brazil,6 and Japan.7 Similarly, other studies, including our own8, have attempted to estimate vaccine-avertible deaths—defined as the number of deaths that could have been averted by vaccination, but were not because of a failure to vaccinate the unvaccinated.8,9 Most empirical studies quantify direct protection conferred by vaccination among vaccinated individuals, but they typically assume the overall impact of vaccination is larger due to indirect protection.4–6 Indeed, these studies often claim that their estimates represent a lower bound on the overall impact of vaccination, although this claim has not been carefully verified.
To support this claim, many of these studies draw, either implicitly or explicitly, on the causal effect framework developed by Halloran and Struchiner,10 in which the overall effect of vaccines can be decomposed into direct and indirect components.11 However, the Halloran and Struchiner framework has not yet been formally extended to cover the specific estimands targeted in vaccine-averted and avertable analyses, which estimate impact. Existing effect estimands are defined by contrasts of individual risk12,13 and have been used to estimate vaccine efficacy in clinical trials.14,15 However, observational studies, post-licensure studies, and policy-makers are often equally or more interested in quantifying public health impact of vaccination in terms of counts16 such as the number of infections, hospitalizations, and deaths in a group of individuals, instead of risk in each individual.
Motivated by recent empirical studies on vaccine-averted and avertible COVID-19 deaths, this article seeks to fill these gaps by: (1) Clearly defining impact estimands as corollaries of the direct and overall effect estimands for averted and avertible outcomes (i.e., counts), and (2) determining the conditions under which direct impact is a lower bound on overall impact. To ground our discussion, we introduce a susceptible-infected-recovered-death (SIRD) model stratified by vaccination status to investigate the direct impact and overall impact of vaccination under different scenarios.
The rest of the article is structured as follows: In Section 2, we summarize the definitions of the direct, indirect, total, and overall effects introduced by Hudgens and Halloran,12 and provide an alternative partitioning of overall effect into components that will align better with the estimands targeted by vaccine-averted and avertible analyses. In Section 3, we propose count outcome corollaries for the direct and overall effects, show how they map onto estimands for vaccine-averted and avertible analyses, and formalize the claim that the direct impact constitutes a lower bound on overall impact. Section 4 outlines a transmission model to simulate vaccine-averted and avertible outcomes. In Section 5, we examine the conditions under which outcomes averted via direct impact among vaccinated individuals (or outcomes avertible via direct impact among unvaccinated individuals) is and is not a lower bound of overall impact on vaccine-averted (or avertible) outcomes.
2. DIRECT, INDIRECT, OVERALL, AND TOTAL EFFECTS
Hudgens and Halloran12 previously defined causal estimands for direct, indirect, overall, and total effects in the two-stage randomized trial, as summarized below.
Setup and notation
Consider a two-stage randomized trial with m groups indexed by i = 1, …, m, such that each group consists of N individuals indexed by j = 1, …, N with a large group size N. All groups are assumed to be of same size N. Partial interference is assumed: Individuals make contacts (leading to interference) within the same group, but individuals in different groups make no contacts. For ease of exposition, assume interest lies in quantifying the effect of vaccination which is a one-time event before the start of an outbreak. Let Aij = 1 if individual j in each group i is vaccinated and Aij = 0 otherwise. Let Ai = (Ai1, Ai2, …, AiN) and Ai,−j = (Ai1, Ai2, …, Aij−1, Aij+1, …, AiN), hereafter referred to as allocation programs.11 Let ai and ai,−j denote possible realizations of Ai and Ai,−j, respectively. Let 𝒜(N) denote the set of all possible 2N vaccine allocations for a group of size N, for which ai ∈ 𝒜(N).
Let Yij(ai) denote the potential binary outcome for individual j in group i with allocation program ai and let Yij(ai,−j, a) denote the potential binary outcome when individual j has vaccination status a and the rest of group i has vaccination status ai,−j.
Individual, group, and population average potential outcomes
Hudgens and Halloran12 define the marginal individual average potential outcome as
and the individual average potential outcome12 as:
where 
is the probability distribution of vaccine allocation program Ai with parameter α ∈ [0,1] representing the proportion vaccinated within group i. Specifically, 
is the probability distribution of Ai conditional on 
. Note that here we use type A parameterization, which gives same effect definitions as type B parameterization suggested by VanderWeele and Tchetgen Tchetgen11 when N is large. Definitions of 
, type A, and type B parameterizations are given in eAppendix 1. Hudgens and Halloran12 further define (marginal) group average potential outcomes 
and 
by averaging over individuals within groups. They also define (marginal) population average potential outcomes12 
and 
.
Population Average Direct, Indirect, Overall, and Total Causal Effects
Hudgens and Halloran define the population average direct casual effect12 as 
, which compares 
when an individual is unvaccinated versus when vaccinated, holding fixed the proportion vaccinated (α). They define the population average indirect casual effect12 as 
, which compares 
for an unvaccinated individual in a group with α proportion vaccinated versus with α′ proportion vaccinated and is hereafter referred to as indirect effect for the unvaccinated. As noted in the literature,17 indirect effect can be analogously defined when an individual is vaccinated: 
, which is hereafter referred to as indirect effect for the vaccinated. They also define the population average total causal effect12 as 
which compares 
when an individual is unvaccinated in a group with α proportion vaccinated versus when the individual is vaccinated in a group with α′ proportion vaccinated. Following from the definition, TE(α, α′) can be decomposed as DE(α′) + IEUnvax(α, α′) or DE(α) + IEvax(α, α′). Finally, they define the population average overall casual effect12 as 
, which compares 
for a typical individual in a group with α proportion vaccinated versus with α′ proportion vaccinated.
Overall Effect Partitioning
Previously, Hudgens and Halloran12 showed that when α = 0, OE(α, α′) is the weighted sum α′ ⋅ TE(α, α′) + (1 − α′) ⋅ IEUnvax(α, α′). In this study, for the purpose of establishing the direct impact as a lower bound of the overall impact that is shown in Section 3, we generalize the partitioning of OE(α, α′) to any α′ > α in Theorem 1.
Theorem 1 (overall effect partitioning)
Theorem 1 is proved in eAppendix 2 and graphically illustrated in Figure 1. Theorem 1 expresses OE(α, α′) as a weighted average of three effects: 1) IEvax(α, α′), 2) TE(α, α′), and 3) IEUnvax(α, α′). Intuitively, if individuals are classified by their vaccination status under a pair of counterfactuals wherein the group has α or α′ proportion vaccinated (α′ > α), then IEvax(α, α′) is in operation for proportion α of individuals who are vaccinated under both counterfactuals; TE(α, α′) is in operation for proportion α′ − α of individuals who are vaccinated under α′ but unvaccinated under α; and IEUnvax(α, α′) is in operation for proportion 1 − α′ of individuals who are unvaccinated under both counterfactuals.
The two rectangles represent a pair of counterfactuals wherein the group has α or α′ proportion vaccinated (α′ > α). Individuals fall into three categories based on their vaccination status under the counterfactuals: 1) The dotted region represents those who are unvaccinated under both counterfactuals (hereafter referred to as “never-vaccinated”) and for whom IEUnvax(α, α′) is in operation; 2) the gridded region represents those who are unvaccinated under α but vaccinated under α′ and for whom TE(α, α′) is in operation (hereafter referred to as “additionally-vaccinated”); and 3) the stripped region represents those who are vaccinated under both counterfactuals (hereafter referred to as “always-vaccinated”) and for whom IEvax(α, α′) is in operation. Theorem 1 shows that OE(α, α′) is a weighted average of three effects: 1) IEUnvax(α, α′), 2) TE(α, α′), and 3) IEvax(α, α′), each weighted by the proportion of individuals for whom the effect is in operation respectively: 1) 1 − α′ for the never-vaccinated, 2) α′ − α for the additionally-vaccinated, and 3) α for the always-vaccinated.
3. ESTIMANDS FOR ANALYSES OF AVERTED AND AVERTIBLE OUTCOMES
Up until now, all effect estimands (except the indirect effect for vaccinated)17 have been previously defined by Hudgens and Halloran in a two-stage randomized trial. To define estimands for analyses of averted and avertible outcomes, we use and expand on the terminology of two-stage randomization. Our goal is to be explicit about the definition of estimands in recent observational studies, to assess the applicability of the lower bound argument, and to place these observational studies in the context of target trial emulation, which is playing a growing role in causal inference.
The original overall effect (OE) and direct effect (DE) are defined for individual risk,10 and their magnitudes are not comparable because OE and DE use different denominators. However, empirical observational studies often assume vaccination prevents more cases overall than directly for a count outcome, by assuming that multiplying OE with total population size is greater than multiplying DE with the number of vaccinated individuals in the presence of indirect protection for unvaccinated individuals.4–9 To verify this assumption, Section 3 defines two estimands for averted and avertible outcomes—direct impact and overall impact—that align with the targeted estimands in the literature on vaccine-averted and avertible outcomes,1–9 by multiplying the original effect estimands with the number of individuals for whom the effect is in operation. Sections 4 and 5 then verify conditions under which direct impact is or is not a lower bound of overall impact.
This paper adds to the existing literature by: (1) Introducing an additional time index t to examine outcomes at multiple timepoints post-vaccination, (2) defining the direct impact and overall impact (i.e., corollaries of the direct and overall effect estimands for count outcomes) to be used in averted and avertible outcome analyses, (3) defining the relationships between direct impact and overall impact, and (4) most importantly, using analytical and simulation approaches to identify conditions under which direct impact may or may not be a lower bound of overall impact.
Notation
To examine the change of direct impact and overall impact over time, we incorporate time index t so that 
denotes the population average cumulative incidence of having developed the outcome by time t, and similarly for 
and all the effect estimands defined in Section 2. t = 0 denotes start of follow-up. In principle, t can be of any timescale, while t is in the unit of days in the simulation in Section 5.
Now let α1 denote a particular vaccination proportion chosen to be implemented in the group, let α2 denote some hypothetical higher vaccination proportion (i.e., α2 > α1), and let α0 denote some hypothetical lower vaccination proportion (i.e., α0 < α1).
Motivating Examples and Causal Questions
During the COVID-19 pandemic, determining the total number of infections (or deaths) averted by vaccination has been of great public health interest.1–7,18–20 Vaccine-averted infections (or deaths) is an impact estimand based on the causal question: How many infections (or deaths) have been averted under the particular proportion vaccinated (α1) compared to the counterfactual in the absence of vaccination (α0 = 0)?
Alternatively, researchers have also shown interest in estimating the vaccine-avertible deaths, that is those that could have been averted by vaccination but were not because of a failure to vaccinate.8,9 The causal question in this case is: How many infections (or deaths) could have been averted under full vaccination (α2 = 1), but were not averted given the particular proportion vaccinated (α1)?
In general, estimands can be defined by comparing the number of infections (or deaths) in a typical group under the particular proportion vaccinated (α1) versus when the group has a lower (α0) or higher (α2) proportion vaccinated.9 We term these as impact estimands because previous literature has referred to vaccine-averted infections (or deaths) as a (population) impact of vaccination.5,7,18,19
Overall Impact
For α0 < α1, overall impact is defined as:
For α2 > α1, overall impact is defined as:
Overall impact (i.e., overall effect multiplied with population size) directly answers the two aforementioned causal questions on vaccine-averted and avertible outcomes. Therefore, overall impact is an estimand for capturing the vaccine averted and avertible outcomes. For reasons that will become clearer later, we force overall impact to have the same sign when α1 is compared to a lower (α0) or higher (α2) hypothetical value. Consequently, arguments of overall impact (δO) are ordered according to the order of the values (e.g., α0 comes before α1 where α0 < α1).
Note that mathematical modelling studies implicitly refer to overall impact when estimating vaccine-averted deaths by simulating the epidemic trajectory under the counterfactual with a hypothetical proportion vaccinated (e.g., α0 = 0) and comparing it with the estimated number under the particular vaccination campaign (α1).
Furthermore, by Theorem 1 and then substituting TE(t, α0, α1) = DE(t, α1) + IEUnvax(t, α0, α1), we have the decomposition:
for α1 > α0. The first term on the right-hand side of the last line of equation (2.1) scales IEvax(t, α0, α1) by number vaccinated under α0. The second term scales DE(t, α1) by the additional number vaccinated under α1 compared to α0. The third term scales IEUnvax(t, α0, α1) by number unvaccinated under α0.
Similarly, by Theorem 1 and then substituting TE(t, α1, α2) = DE(t, α1) + IEvax(t, α1, α2), we also have the decomposition:
for α2 > α1. The first term on the right-hand side of the last line of equation (2.2) scales IEvax(t, α0, α1) by number vaccinated under α2. The second term scales DE(t, α1) by the additional number vaccinated under α2 compared to α1. The third term scales IEUnvax(t, α0, α1) by number unvaccinated under α2.
Direct Impact
Because the overall impact estimand (δO) is defined as a contrast under two different vaccination proportions, at minimum researchers need to observe two noninteracting groups where each is “assigned” to one of the vaccination strategies of interest to estimate overall impact.12 Ideally, vaccination would be assigned via a two-stage randomized trial. However, these trials are rarely conducted because they are expensive and often hard to justify.21 Instead, most empirical studies only observe a single group under one vaccination proportion (α1), in which case it is only feasible to estimate an impact corollary of the direct effect. Examples include observational studies comparing discrete hazards or incidence rates between vaccinated and unvaccinated individuals based on national vaccine data systems.4–9 These studies multiply the direct effect with number vaccinated (or unvaccinated) to estimate outcomes averted (or avertible) via direct impact among vaccinated (or unvaccinated) individuals,4–9 and then generally assume that it is a lower bound of overall impact on averted (or avertible) outcomes in the entire population.4–6,8 In this section, we formalize this definition of direct impact and show how it relates to the overall impact.
Some empirical observational studies have used commonly available data to estimate deaths averted via direct impact among vaccinated individuals4–7 based on formulas similar to N ⋅ (α1 − α0) ⋅ DE(t, α1) by setting α0 = 0, while other studies have estimated deaths avertible via direct impact among unvaccinated individuals8,9 based on formulas similar to N ⋅ (α2 − α1) ⋅ DE(t, α1) by setting α2 to be greater than α1 (Note these studies have also considered increases in proportion vaccinated over time, such as occurred during a vaccination campaign, which here we ignore for simplicity). Following the literature,4–9,18 let the direct impact (δD) for any α, α′ ∈ [0,1] be δD(t, α, α′) = N ⋅ |α′ − α| ⋅ DE(t, α′).
In particular, for α1 > α0, we have
and for α2 > α1,
Note that δD (t, α0, α1) or δD (t, α2, α1) is not a meaningful causal estimand by itself because δD (t, α1) is conditional on α1 only and does not account for the change in direct effect under the counterfactual when the proportion vaccinated is α0 or α2 instead of α1. Importantly, now direct impact can be a lower bound of overall impact. This is because Theorem 1 provides a useful decomposition of OE into TE, IEUnvax, and IEvax, leading to the relationship between OE and |α′ − α| ⋅ DE in equations 2.1 and 2.2, such that we can map direct impact onto N ⋅ |α′ − α| ⋅ DE(t, α′) and show that direct impact can be a lower bound of overall impact. That is, equations (2.1) and (2.2) now represent:
for α1 > α0, and
for α2 > α1. Based on equations 1.1 and 1.2, overall impact has the same sign when α1 is compared to a lower (α0) or higher (α2) hypothetical value. Therefore, we can formalize the assumption often made in the empirical vaccine-averted and avertible death studies—that is, direct impact is a lower bound of overall impact by considering Claim 1:
and
If Claim 1 is true, direct impact, which can be estimated using commonly available data,4–9,18 is a useful lower bound of overall impact that is relevant for policy-making and retrospective policy evaluation requiring samples from a population of groups as in two-stage randomized trials.12
Following the empirical studies estimating vaccine-averted or avertible outcomes, we consider two special cases of Claim 1. In Claim 1a, a particular vaccination proportion α1 is compared to a hypothetical of no vaccination (α0 = 0) to quantify vaccine-averted outcomes.
Claim 1a (vaccine-averted outcomes):
In words, Claim 1a asserts that for α1 > α0 = 0, outcomes averted via direct impact of current vaccination among the vaccinated individuals δD(t, 0, α1) is a lower bound of total vaccine-averted outcomes among both vaccinated and unvaccinated individuals δO(t, 0, α1).
In Claim 1b, α1 is compared to a hypothetical of near-universal vaccination (α2 = 0.9 > α1) to quantify vaccine-avertible outcomes, following the literature.9
Claim 1b (vaccine-avertible outcomes):
In words, Claim 1b asserts that for α2 = 0.9 > α1, outcomes avertible via direct impact of current vaccination among some unvaccinated individuals δD(t, 0.9, α1) is a lower bound of vaccine-avertible outcomes among both vaccinated and unvaccinated individuals δO(t, α1, 0.9) comparing current proportion vaccinated (α1) to a proportion of α2 = 0.9. Note that we do not compare to full vaccination (α2 = 1) because direct effect would be undefined if the group is fully vaccinated, and there would not be an epidemic if the group is fully vaccinated with a highly effective vaccine. In the following sections, we use a susceptible-infected-recovered-death (SIRD) model to verify the conditions under which the Claims may or may not hold.
Based on the effect partitioning results in equation 4.1 or 4.2, Claims 1a and 1b would be true if IEvax and IEUnvax are non-negative. However, it is not immediately intuitive when we might expect that to be the case. Therefore, in the next two sections we describe a transmission model that we use to check Claims 1a and 1b under various scenarios (Section 4), and then describe conditions under which the direct impact is or is not a lower bound for the overall impact, using both analytical and simulation-based approaches (Section 5).
4. TRANSMISSION MODEL
The SIRD Model with Vaccination at Baseline
To model the impact of vaccination on infection and death, an SIRD model is used to represent a well-mixed group in a two-stage randomized trial assuming partial interference.12 To simulate direct impact and overall impact, we simulate a typical group with a large size under a pair of counterfactual vaccination proportions. In this simulation, the (marginal) group average potential outcome is equivalent to (marginal) population average potential outcome because there is only one group in the population. The model consists of four states for a vaccinated or unvaccinated individual—susceptible, infectious, recovered, and death due to infection. We assume that the group has been randomly assigned with a vaccination policy wherein the proportion vaccinated is α, and individuals have been randomly assigned with vaccination status a at baseline (a = 1 for vaccinated and a = 0 for unvaccinated; for equation [5], subscript v denotes vaccinated and subscript u for unvaccinated). The vaccine is “leaky” in protecting against infection and infection-related death—that is, vaccinated individuals have the susceptibility reduced by a factor θ against infection (i.e., vaccine efficacy against infection [VEinfection] is (1 − θ) • 100%) and have the susceptibility reduced by another factor κ against death (i.e., vaccine efficacy against death given infection [VEdeath|infection] is (1 − κ) • 100%). Individuals mix homogeneously such that each vaccinated or unvaccinated susceptible individual is equally likely to contact any infectious individual. Vaccinated and unvaccinated infectious individuals are equally contagious. The transmission dynamics are:
where γ = recovery rate, 
is the hazard rate of infection, with β the number of effective contacts made by a typical infectious individual per unit time, and μ = probability of death due to infection. In equation (5), SU,α(t) and Sv,α(t) denote, respectively, the number of susceptible individuals who are unvaccinated and vaccinated, IU,α(t) and Iv,α(t) for the infectious individuals, RU,α(t) and Rv,α(t) for the recovered individuals who are no longer at risk, and DU,α(t) and Dv,α(t) for those who died due to infection. N(t) denotes the sum of all compartments at time t. eFigure 1 shows the model flowchart, and eTable 1 shows the parameter values used in simulation.
Software
All model simulations and visualization are conducted using R 4.2.2 (R Foundation for Statistical Computing, Vienna, Austria).22 All models are implemented using R package odin.23 Code is available at https://github.com/katjia/population_level_effects.
5. WHEN IS DIRECT IMPACT A LOWER BOUND OF OVERALL IMPACT?
Rationale
The SIRD model has many simplifying assumptions compared to real-world settings. If one can show a counterexample to Claim 1 based on the simplest SIRD model, then Claim 1 is not guaranteed to be true in more general and realistic models. To identify counterexamples of Claim 1, we consider common violations to assumptions on time-invariant parameters so that indirect effects can be negative: (1) Number of effective contacts may increase over time due to meteorological factors,24 lifting of non-pharmaceutical interventions,25 and seasonal variation in social contacts, (2) infection-fatality risk may increase over time,26 (3) waning immunity clearly occurred in the Delta and Omicron waves of the COVID-19 pandemic.27,28 If Claim 1 fails under these common violations, it would be overly optimistic for the empirical studies to assume that Claim 1 holds in reality.
Scenarios
We check the Claims under 5 scenarios (Table). Scenario 1 refers to the SIRD model in equation (5) with time-invariant parameters. eTable 1 lists the model parameters for which one or two parameters vary under each Scenario separately: Scenario 2 increases the number of effective contacts made by a typical infectious individual per day (β) from 0.15 to 0.6 from Day 300 onwards; Scenario 3 increases probability of death due to infection (μ) from 0.01 to 0.1 from Day 300 onwards; Scenario 4 allows both VEinfection and VEdeath|infection to wane linearly after Day 100 reaching 0% at Day 300; and Scenario 5 combines Scenarios 2 (increasing β) and 4 (waning VEs).
Proof of Claim 1 at the end of outbreak under Scenario 1
For Scenario 1 (i.e., time-invariant parameters), in eAppendix 4, we prove that Claim 1 (i.e., direct impact is a lower bound of overall impact for any two proportions vaccinated) holds in the SIRD model at the end of outbreak (i.e., at t → ∞).
Simulations
For t < ∞, we first verify Claims 1a and 1b (special cases of Claim 1) through simulation based on parameters specified in eTable 1 and initial conditions in eTable 2. We specify α0 = 0 versus α1 = 0.7 in the pair of trajectories to verify Claim 1a and α1 = 0.7 versus α2 = 0.9 to verify Claim 1b. If Claims 1a and 1b both hold for the specified parameters, Latin hypercube sampling is conducted to generate alternative sets of proportions vaccinated and model parameters to verify the full Claim 1. If only one of Claim 1a or 1b holds, Latin hypercube sampling is conducted to verify that the Claim holds under alternative proportion vaccinated and model parameters (e.g., trying alternative values for α1 while fixing α0 = 0 for Claim 1a).
Briefly, Claim 1 only holds under Scenario 1 (i.e., time-invariant parameters), but it does not hold under any other Scenarios. The Table summarizes the results. Figures 2 and 3 show trajectories of direct impact and overall impact throughout the epidemic to verify Claims 1a and 1b, respectively. Figures 2 and 3 show that Claims 1a and 1b hold under Scenario 1. Moreover, given alternative sets of proportions vaccinated and model parameters, Latin hypercube sampling verifies that Claim 1 holds under Scenario 1 (eAppendix 5). Under Scenario 2 where β increases, Claims 1a and 1b do not hold: Direct impacts are not lower bounds of overall impacts (Figures 2 and 3) due to negative indirect effects (eFigures 3 to 5). Under Scenario 3 where μ increases, Claims 1a and 1b hold for infections but not deaths due to negative indirect effects for death (eFigures 3 to 5). Under Scenario 4 where VEs wane, only Claim 1a (vaccine-averted outcomes) holds (Figure 2). Given different values for α1 (while holding constant α0 = 0) and alternative sets of model parameters, Latin hypercube sampling verifies that Claim 1a holds under Scenario 4 (eAppendix 5). However, Claim 1b (vaccine-avertible outcomes) does not hold (Figure 3) due to negative indirect effects (eFigures 4 and 5). Finally, under Scenario 5 where β increases and VEs wane, Claims 1a and 1b do not hold.
Scenario 1, all parameters are time-invariant; Scenario 2, the number of effective contacts made by a typical infectious individual per day (β) increases from 0.15 to 0.6 at Day 300; Scenario 3, probability of infection-related death (μ) increases from 0.01 to 0.1 at Day 300; Scenario 4, vaccine efficacies against infection and death start to wane linearly after Day 100 reaching 0% at Day 300; and Scenario 5, the combination of Scenarios 2 and 4.
Scenario 1, all parameters are time-invariant; Scenario 2, the number of effective contacts made by a typical infectious individual per day (β) increases from 0.15 to 0.6 at Day 300; Scenario 3, probability of infection-related death (μ) increases from 0.01 to 0.1 at Day 300; Scenario 4, vaccine efficacies against infection and death start to wane linearly after Day 100 reaching 0% at Day 300; and Scenario 5, the combination of Scenarios 2 and 4.
6. DISCUSSION
Motivated by recent research on estimating deaths averted by COVID-19 vaccination, this study adjusts estimands defined by Hudgens and Halloran to accommodate analyses on averted and avertible outcomes due to infectious disease interventions, thereby enabling researchers to distinguish the estimands when conducting and interpreting related studies. An epidemic model is simulated to verify the commonly held claim that outcomes averted (or avertible) via direct impact among the vaccinated (or unvaccinated) individuals is a lower bound of the overall impact on averted (or avertible) outcomes. Based on the SIRD model, we show that the lower bound fails when transmission or fatality parameter increases, or vaccine efficacies wane, implying that the lower bound is not guaranteed to hold for more general and realistic models. Consequently, it would be overly optimistic for empirical studies4–9 to assume that they have estimated a lower bound of true number of averted (or avertible) outcomes through estimating the direct impact.
When indirect effects for the vaccinated and unvaccinated are both positive (that is, higher vaccination coverage yields positive number of averted or avertible outcomes), the direct impact is a lower bound of overall impact (equations 4.1 and 4.2). In general, the indirect effects may be positive because vaccination reduces the infectious individuals at a given time, such that infection-naïve individuals are less likely to be infected.29 However, as we have shown above, there are scenarios under which direct impact is not guaranteed to be a lower bound of overall impact due to negative indirect effects (Scenarios 2–5).
First, when the number of effective contacts made by a typical infectious individual per day (β) increases over time (Scenario 2), overall impact on infection decreases (and can be negative) because at the early stage of outbreak, the less-vaccinated group has many infected and recovered with sterilizing immunity; while the more-vaccinated group has more susceptible (i.e., infection-naïve) individuals who have escaped the earlier infections and will experience a higher force of infection at the later time (eFigure 6). β is affected by probability of transmission per contact and number of contacts per day. Transmission probability depends on meteorological factors such as absolute humidity for the influenza A virus,24 behavioral factors such as usage of personal protective equipment,30 and biological factors such as changes in host immunity, and evolution of strains (a proper examination of the impact of co-circulating strains on the indirect effect requires modeling the cross immunity).31 The number of contacts can be affected by implementation or discontinuation of non-pharmaceutical interventions and seasonal variation in social contacts.
Second, when the probability of infection-related death (μ) increases over time (Scenario 3), overall impact on death decreases because the extensively vaccinated group(s) has more who escape the earlier infections and then experience the high fatality at the later stage of outbreak (eFigure 6). It is plausible for lethality of pathogens to increase over time: Disease severity was found to increase in the autumn wave of the 1918 flu pandemic compared to the spring-summer wave of the same pathogen in that year. Some evidence from the A/H1N1 pandemic in 2009 also suggests that lethality of a subtype may increase in the second wave compared to the first wave.26 Increasing lethality implies that vaccination at the beginning of outbreak can lead to a negative overall impact by postponing cases. On the other hand, if the fatality rate increases with the infection peak due to the sudden shortage of healthcare resources, the overall impact on death will be more positive because vaccines delay infection and flatten the epidemic curve. Likewise, if infection-fatality rates decline progressively due to improvements in care,32,33 then vaccination that delays the epidemic can have amplified positive overall impact.
Third, the overall impact of increased vaccination coverage decreases and may become negative when vaccine efficacies wane (Scenario 4). In particular, the vaccination proportion of α2 = 0.9 may result in more infections and deaths than the proportion of α1 = 0.7. eFigure 7 shows the epidemic curves given α1 = 0.7 and α2 = 0.9 proportions vaccinated. When α1 = 0.7, the epidemic peaks earlier and slows down due to the build-up of recovered individuals with sterilizing immunity. Consequently, return of full susceptibility (i.e., loss of protection) among vaccinated individuals is too late to rescue the epidemic. However, when α2 = 0.9, epidemic is delayed and a pool of individuals who are at risk of infection is built up, such that return of susceptibility among vaccinated individuals can rescue the epidemic.
The current study has some limitations. First, vaccination is assumed to be a one-time event at baseline before the start of outbreak, but in reality, vaccine rollout is continuous over time and may occur during outbreak. A formalization of the causal effects under a time-varying regime (e.g., by specifying a sequence of proportions vaccinated 
over time) will be left to future work. Second, throughout, we have considered the case where vaccination occurs at random. However, in most empirical settings there may be strong confounding due to staged rollout of vaccines as well as differences in vaccine acceptance by behavioral and health characteristics. Such confounding if uncontrolled threatens the validity of inferences about the magnitude of effects, whether or not a bound is valid in ideal (unconfounded) circumstances. Third, to illustrate the counterexamples to the Claims, the current study uses a stylized model with multiple simplifications. For example, this model does not consider multiple risk groups, which may have had heterogenous susceptibility to adverse outcomes and heterogeneous mixing patterns, leading to a negative indirect effect for the unvaccinated such that the Claims do not hold. Vaccination for a subgroup could cause negative indirect effect in other subgroups by increasing risk for more severe complications. For example, empirical evidence showed that low rubella vaccination coverage in children increased rubella incidence in the 15 years or over and incidence of congenital rubella in newborns.34 Mathematical modelling suggested that infant vaccination against varicella could increase the rate of reactivation (i.e., zoster) in the entire population.35,36 The examples suggest that the potential for negative indirect effects may be greater when subgroups differ in important properties, such as (in the rubella and varicella examples) risk for adverse outcomes.
Another limitation of our model is that in the scenarios considered, changes in lethality and transmission were assumed to occur at a fixed time, whereas in reality they might well occur either in response to pathogen evolution37 or to behavioral changes that are affected by the epidemic trajectory. However, our goal was not to describe the details of a particular epidemic, but to describe qualitatively the conditions under which infectious disease interventions may not prevent more cases overall than directly. Finally, the SIRD model does not consider deaths due to other factors, meaning that simulation experiments based on this model are applicable to studies whose outcome of interest is a consequence of infection. Over a short time frame (e.g., 1 year), it is acceptable to consider deaths due to infection only, when other causes of deaths are negligible. We also did not consider possible adverse events after vaccination, although adverse events have important policy implications. Our focus here is on defining estimands for the protective impact of interventions on averting disease outcomes. Future studies can extend the estimands to investigate adverse events due to interventions.
In conclusion, this study examines the commonly held assumption of empirical vaccine-averted and avertible analyses that direct impact is a lower bound on overall impact due to indirect protection. We show that the use of direct impact as a lower bound of overall impact in averted and avertible outcome studies is reliable only under very strong and often unrealistic assumptions. Therefore, empirical studies should not assume this relation because parameters may vary over time in reality. Alternatively, if researchers want to estimate the averted and avertible outcomes, a transmission model should be used to capture the overall impact, thereby improving clarity about what is estimated despite the expense of making additional modelling assumptions.
Data Availability
This is a simulation study with code available at https://github.com/katjia/population_level_effects.
eAppendix 1. Type A versus type B parameterizations
In the main text, the potential outcomes 
and 
are defined under type A parameterization. eAppendix 1 will define type A and B parameterizations. We assume a large group size N such that results are equivalent under both types of parameterizations.1
Let ρ(Ai = ai; α) denote the probability that group i receive allocation program ai given parameter α.
According to VanderWeele and Tchetgen Tchetgen,1 a parameterization is said to be of type A with parameter α = (N, K) for group i if the allocation program Ai is randomly allocated conditional on 
with probability mass function:
(B) A parameterization is said to be of type B with parameter α if the allocation program Ai is randomly assigned to individuals in group i according to the known Bernoulli probability mass function:
According to VanderWeele and Tchetgen Tchetgen,1 Hudgens and Halloran 2 defined causal effects under Type A parameterization.
The marginal individual average potential outcome of Hudgens and Halloran 2 is:
where 
and ρ(Ai = s; α) follows Definition 1 (A) under type A parameterization. The individual average potential outcome 2 is:
where 
and ρ (A i, −j = s, Aij = a; α) follows Definition 1 (A) under type A parameterization.
Alternatively, VanderWeele and Tchetgen Tchetgen1 proposed the definitions of potential outcomes under Type B parameterization.
Under type B parameterization, the marginal individual average potential utcome1 is:
where 
and ρ(Ai = s; α) follows Definition 1 (B) under type B parameterization. The individual average potential outcome1 is:
where 
and ρ(Ai,−j = s; α) follows Definition 1(B) under type B parameterization.
eAppendix 2. Proof of Theorem 1 (overall effect partitioning)
The first line follows by definition; the second follows from line 1 of proof of Theorem 3 from VanderWeele and Tchetgen Tchetgen1; the third distributes terms, the fourth adds and subtracts 
and 
the fifth and sixth rearrange terms; and the last applies definitions of IEvax, TE, and IEUnvax. The proof can also be shown visually as illustrated in Figure 1.
eAppendix 3. The susceptible-infected-recovered-death (SIRD) model
1. Model structure
eFigure 1 shows the flowchart for the transmission model described in equation (5) of the main text.
2. Model parameters
To verify Claim 1a, a pair of trajectories are simulated with same parameters except for the proportions vaccinated (α0 = 0 versus α1 = 0.7). Similarly, to verify Claim 1b, another pair of trajectories are simulated with α1 = 0.7 versus α2 = 0.9. Epidemic trajectories are simulated based on the SIRD model in equation (5)—first under Scenario 1 (i.e., time-invariant parameters) and then Scenarios 2 to 5 (i.e., time-dependent parameters). eTable 1 specifies parameter values for simulations. All models are implemented using R package odin,3 whereas the actual solution of the differential equations is conducted with the deSolve package using numerical solvers “lsoda” and “ode45.”4
eAppendix 4. Proving Claim 1 under Scenario 1 in the SIRD model at t → ∞
1. Definitions
Let 0 ≤ α0 < α1 < α2 ≤ 1. Recall Claim 1:
and
By expanding equations (4.1) and (4.2) and suppressing the time notation (t), we have:
for α1 > α0, and
for α2 > α1.
We first give an intuitive explanation of δD(α0, α1) being a lower bound of δO(α0, α1). Recall that individuals fall into three categories based on vaccination status under the counterfactuals with α0 or α1 proportion vaccinated: 1) Those who are unvaccinated under both counterfactuals (referred to as “never-vaccinated” and represented by the dotted region in Figure 1); 2) those who are unvaccinated under α0 but vaccinated under α1 (referred to as “additionally-vaccinated” and represented by the gridded region in Figure 1); and 3) those who are vaccinated under both counterfactuals (referred to as “always-vaccinated” and represented by the stripped region in Figure 1). If we update risk of the “always-vaccinated” under α0 from 
to 
and update risk of the “never-vaccinated” under α0 from 
to 
, then we have an updated right-hand side of equation (6.1):
which is reduced to δD(α0, α1) and is lower than the actual overall effect because the original risk 
in the first term of the last line of equation (6.1) is higher than the updated 
(to be proved in Section 3), and the original risk 
in the last term of the last line of equation (6.1) is also higher than the updated 
(to be proved in Section 3). Therefore, δD(α0, α1) in the updated equation is lower than δO(α0, α1) in the original equation (6.1) and Claim 1 holds for α1 > α0.
Similarly, for α2 > α1, if we update risk of “always-vaccinated” under α2 from 
to 
and update risk of “never-vaccinated” under α2 from 
to 
, we have on the right-hand side of equation (6.2):
which is reduced to δD(α2, α1) and is lower than the actual overall effect because the original risk 
in the first term of the last line of equation (6.2) is lower than the updated 
(to be proved in Section 3), and the original risk 
in the last term is also lower than the updated 
(to be proved in Section 3). Therefore, δD(α2, α1) in the updated equation is lower than δO(α1, α2) in the original equation (6.2) and Claim 1 holds for α2 > α1. Next, we will formally prove Claim 1 under Scenario 1 at t → ∞.
2. Sufficient Conditions for Claim 1 to hold
Based on equation (4.1), Claim 1 holds for all α0 < α1 when IEvax(t, α0, α1) and IEUnvax(t, α0, α1) are both non-negative. Similarly, based on equation (4.2), Claim 1 holds for all α2 > α1 when IEvax(t, α1, α2) and IEUnvax(t, α1, α2) are both non-negative.
3. Proof of Claim 1 under Scenario 1 in the Susceptible-infected-recovered (SIR) Model at t → ∞
When probability of death due to infection is zero, the SIRD model in equation (5) reduces to a SIR model stratified by vaccination status. IEUnvax(∞, α0, α1) is positive for α1 > α0 in the SIR model by proving that 
, the final epidemic fraction in the unvaccinated individuals (i.e., the fraction of the unvaccinated individuals that became infected at the end of the outbreak), decreases with α. Similarly, IEvax(∞, α0, α1) is positive by proving that 
, the final epidemic fraction in the vaccinated individuals, also decreases with α.
3.1 Proof that 
anddecrease with α
To prove that 
and 
decrease with the proportion vaccinated α, we prove that 
, the average final epidemic fraction, decreases with α. 
is the weighted average of 
and 
.5 Assuming that 1) the initial proportion susceptible is close to one, 2) the total group size N is near infinite, and 3) transmission events follow a Poisson process, the average final epidemic size is:
where θ = 1 − VEinfection/100% given a leaky vaccine, and R0 is the basic reproduction number, defined as the expected number of new infections caused by a single infected individual during the infectious period in a completely susceptible population.
Define two separate functions with independent arguments, Z and α (Miller5 used a similar type of argument to arrive at the final epidemic size):
Here the equality LHS(⋅) = RHS(⋅) is not assumed; in particular, RHS(⋅) takes two independent arguments, Z and α. Plot LHS(Z) and RHS(Z, α) in one graph, where R0 = 4, θ = 0.2, and α = 0, 0.3 and 0.7.
LHS(Z) = RHS(Z, α) when Z = 0, which is a trivial solution. We will now prove that if a positive (non-trivial) root exists such that LHS(Z) = RHS(Z, α), the value of this root decreases when α increases.
We first show that the first derivative of RHS(Z, α) w.r.t Z at zero is greater than 1, which is the derivative of LHS(Z). That is:
where Re is the effective reproduction number when proportion α of the group is vaccinated with a leaky vaccine.6 The above inequality holds because Re must be greater than 1 in order for any epidemic to occur, meaning that for some ϵ > 0, RHS(ϵ, α) > LHS(ϵ). In other words, RHS(Z, α) exceeds LHS(Z) for small values of Z.
Next, we prove that RHS(Z, α) is a concave function by computing the second derivative w.r.t Z:
for α ∈ [0,1]. Because the second derivative is negative for α ∈ [0,1], RHS(Z, α) is a concave function of Z for any possible values of α. Taken together, the first and second derivatives indicate that RHS(Z, α) is above LHS(Z) when Z is greater than but close to zero, and RHS(Z, α) bends towards LHS(Z) (i.e., the diagonal) and by assumption (that there is a nontrivial root of LHS(Z) = RHS(Z, α)) coincides with the diagonal at some value of Z (See plot above).
Finally, consider two values α1 > α0 giving rise to two curves RHS(Z, α = α0) and RHS(Z, α = α1). Let α1 = 0.7 and α0 = 0.3 as the plot illustrates. Because 0 < θ < 1, [exp(−θR0Z) − exp(−R0Z)] in equation (9) is always positive at any given Z > 0. Therefore, RHS(Z, α = α0) > RHS(Z, α = α1) at any given Z > 0 when α1 > α0. Suppose the upper curve RHS(Z, α0) meets LHS(Z) at Z = z0. Then RHS(z0, α1) < z0 = LHS(z0). Given that all expressions are continuous and that RHS(ϵ, α0) > LHS(ϵ), there must have been some value z1 to the left of z0 where RHS(z1, α1) crossed LHS(z1). In other words, RHS(z1, α1) = LHS(z1) for some value z1 < z0, proving that the root z1 of RHS(z1, α1) = LHS(z1) is less than the root z0 of RHS(z0, α0) = LHS(z0) for α1 > α0, QED. These non-trivial roots are unique because RHS(Z, α) is concave.
Any value Z that satisfies LHS(Z) = RHS(Z, α) will satisfy equation (7). Suppose z0 is a positive root for LHS(Z) = RHS(Z, α = α0) such that 
, and similarly suppose z1 is a positive root for LHS(Z) = RHS(Z, α = α1) such that 
. Therefore, we have 
when α1 > α0.
Given that 
,5 
decreases with α because 
decreases with α. Similarly, given that 
,5 
decreases with α.
3.2 Proof of IEUnvax(∞, α0, α1) > 0 and IEVax(∞, α0, α1) > 0 for α1 > α0
In Section 3.1, we proved that 
decreases with α. Therefore, 
for α1 > α0. We also proved that π1,α (∞) decreases with α, and therefore 
for α1 > α0.
Since IEUnvax(∞, α0, α1) and IEvax(∞, α0, α1) are positive, Claim 1 holds for any α1 > α0 at t → ∞. The proof in Section 3 also applies to any α2 > α1 in showing that IIEUnvax(∞, α1, α2) and IIEvax(∞, α1, α2) are positive, such that Claim 1 holds under Scenario 1 in the Susceptible-infected-recovered (SIR) Model at t → ∞.
4. Proof of Claim 1 under Scenario 1 in the Susceptible-infected-recovered-death (SIRD) Model at t → ∞
Now consider the SIRD model with non-zero probability of death due to infection. Based on the transmission dynamics in equation (5), those who died of infection are fixed proportion of those who became infected at the end of epidemic. We have:
for the unvaccinated individuals and
for the vaccinated individuals.
For the unvaccinated individuals, by substituting the right-hand side of equation (10.1) into the equation of IEUnvax,death, we have:
We show that IEUnvax,death(∞, α0, α1) has the same direction as IEUnvax,infection(∞, α0, α1). If IEUnvax,infection(∞, α0, α1) is positive (Section 3), so is IEUnvax,death(∞, α0, α1).
Similarly, for the vaccinated individuals, by substituting the right-hand side of equation (10.2) into the equation of IEvax,death, we have:
If IEvax,infection(∞, α0, α1) is positive (Section 3), so is IEvax,death(∞, α0, α1).
In sum, IEvax(∞, α0, α1) and IEUnvax(∞, α0, α1) are positive for both death and infection given α1 > α0 (similarly, IEvax(∞, α1, α2) and IEUnvax(∞, α1, α2) are positive for α2 > α1). Therefore, Claim 1 holds under Scenario 1 in the Susceptible-infected-recovered-death (SIRD) model at t → ∞.
eAppendix 5. Latin hypercube sampling to verify Claim 1 under Scenario 1 and Claim 1a under Scenario 4
1. Verifying Claim 1 under Scenario 1
As Figures 2 and 3 show, Claims 1a and 1b hold under Scenario 1 given the particular proportions vaccinated (α0 = 0 versus α1 = 0.7 for Claim 1a and α1 = 0.7 versus α2 = 0.9 for Claim 1b) and parameter values specified in eTable 1. To test the robustness of Claim 1 under alternative parameter combinations, Latin hypercube sampling (LHS) is used to randomly draw 1,000 sets of proportions vaccinated and parameters. Initial conditions are same as in eTable 2. LHS samples are first generated from uniform distributions using the lhs package7 and then scaled to match the range of each parameter as shown in eTable 3.
To exclude results due to roundoff errors during numerical integration of the differential equations, Claim 1 is disproven based on the criteria δO − δD < −10−7at any time t, rather than the conventional criteria of < 0. From the LHS procedure, Claim 1 holds under Scenario 1 for all the 1,000 iterations. eFigure 2 shows direct impact and overall impact for 50 random LHS samples.
2. Verifying Claim 1a under Scenario 4
Figure 2 also shows that Claim 1a holds under Scenario 4 given the parameter values specified in eTable 1. Similarly, LHS is conducted to verify if Claim 1a holds under Scenario 4 with different values for α1 (while holding constant α0 = 0) and alternative sets of model parameters; parameter ranges are specified in eTable 3. From the LHS procedure, Claim 1a holds under Scenario 4 for all the 1,000 iterations (eFigure 2).
eAppendix 6. Trajectories of indirect effects for the vaccinated and unvaccinated.
Based on equation (4.1), Claim 1a does not hold when IEUnvax(t, 0, 0.7) is negative, as shown in eFigure 3. Similarly, based on equation (4.2), Claim 1b may not hold when IEUnvax(t, 0.7,0.9) or IEvax(t, 0.7,0.9) is negative, as shown in eFigures 4 and 5, respectively.
Scenario 1, all parameters are time-invariant; Scenario 2, the number of effective contacts made by a typical infectious individual per day (β) increases from 0.15 to 0.6 at Day 300; Scenario 3, probability of infection-related death (μ) increases from 0.01 to 0.1 at Day 300; Scenario 4, vaccine efficacies against infection and death start to wane linearly after Day 100 reaching 0% at Day 300; and Scenario 5, the combination of Scenarios 2 and 4.
eAppendix 7. Trajectories of proportion susceptible by vaccination status over time
Scenario 1, all parameters are time-invariant; Scenario 2, the number of effective contacts made by a typical infectious individual per day (β) increases from 0.15 to 0.6 at Day 300; Scenario 3, probability of infection-related death (μ) increases from 0.01 to 0.1 at Day 300; Scenario 4, vaccine efficacies against infection and death start to wane linearly after Day 100 reaching 0% at Day 300; and Scenario 5, the combination of Scenarios 2 and 4.
eAppendix 8. Epidemic curves under Scenario 4
The dashed line represents the number of effectively susceptible persons = Su(t) + θ(t) ⋅ Sv(t) where Su(t) is the number of unvaccinated susceptibles and θ(t) ⋅ Sv(t) is the number of vaccinated susceptibles multiplied with θ(t) = 1 − VEinfection(t)/100%. The dashed lines show that more individuals at risk of infection (weighted by that risk) are built up under α2 = 0.9 than α1 = 0.7 before the epidemic peaks.
Footnotes
Methods and results remain the same. Abstract, introduction, and conclusion have been revised to highlight the contributions and update the key takeaways.
REFERENCES
