Fair and efficient vaccine allocation: A generalized Gini index approach

Abstract

The paper proposes an optimization model for the allocation of vaccines to a heterogeneous population composed of several subpopulations with different sizes and epidemiological disease transmission parameters. As the objective, an aggregated function combining a standard utilitarian efficiency criterion with a Gini index–related penalty term is considered. Contrary to previous work, we adopt an outcome equity view: The inequity measure is not based on vaccination fractions or other input factors, but on the fractions of individuals escaping infection, as predicted by an susceptible‐infectious‐removed (SIR) model. An adjusted pro rata (APR) policy of vaccine allocation minimizing inequity in this outcome view is introduced, and a numerical procedure for its determination is presented. The concepts are developed both for the case of segregated subpopulations and for that of interactions between the subpopulations. Interestingly, in a large number of instances, the optimal solution under the aggregated objective function turns out to be identical to APR. Whether APR is locally or even globally optimal in a concrete case depends on the relation of an inequity aversion parameter to certain threshold values. While the local optimality threshold can be determined by linear programming, the determination of the global optimality threshold, as the vaccine allocation problem itself, is a problem of nonconvex optimization. We suggest an exact optimization approach for smaller instances, and propose algorithms building on particle swarm optimization for threshold determination and allocation optimization at larger instances. Extensions to alternative outcome measures such as the number of fatalities are presented as well. In addition to the investigation of randomly generated instances, two test cases from the literature are revisited in the context of the present work. Moreover, a new case study based on data from the COVID‐19 outbreak in Austria in 2020 is introduced and analyzed.

INTRODUCTION

The recent worldwide COVID‐19 crisis showed not only the importance of vaccination as a main instrument to mitigate the effect of a pandemic, but also the difficulties in deploying this instrument in an effective way. A main issue in this context is vaccine availability. Due to the limitations in vaccine production and also due to prices that may overwhelm the financial resources of poorer countries, the case where a high demand of vaccine meets a scarce supply is, unfortunately, not rare.

Some articles (for an example, see Keeling & Shattock, 2012) addressed the problem of minimizing the overall size of an epidemic by vaccine allocation under insufficient supply, but admitted that the resulting solution can be highly inequitable and may therefore not be accepted by the society. It has been stated that in the vaccine allocation literature, a reconciliation between equity and efficiency is missing (Duijzer et al., 2018). This paper aims at overcoming this gap by applying, as the objective function for a proposed optimization model, a social welfare function that combines efficiency with equity, a so‐called “inequity‐averse aggregation function” (IAAF) (see Karsu & Morton, 2015).

The envisaged scenario is prophylactic vaccination of a heterogeneous overall population. A model with multiple subpopulations (patches) is used. We start with the situation where the subpopulations are noninteracting (e.g., spatially segregated) in order to obtain a basic insight into the influence of vaccination policies, and relax this assumption then to the case of interaction.

Contrary to some previous works that considered fairness in their vaccine allocation models, usually in the form of equity constraints, this paper aims at health outcome equity rather than at input equity. As to the distinction between these two notions, see Pressman et al. (2021). In our context, the focus on outcome equity means that the ideal state with respect to fairness is not that every person has the same probability of getting vaccinated, no matter whether in their environment, the infection risk is vast or minimal, but rather that everybody has the same chance of not being infected (or, alternatively, of not needing hospitalization, or of not dying from the disease, etc.).

To make the motivation for using outcome equity quite clear, let us take the case of breast cancer as an example. It is well known that this disease primarily affects women, with a lifetime incidence rate of roughly 1/8. In rare cases, also men can be affected; however, their lifetime incidence rate is only about 1/800. Suppose now that a vaccine that permanently protects against breast cancer will be developed in the future, and 10 million doses of this vaccine will be available in a country with 20 million inhabitants, half of whom are women. The solution maximizing input equity would allocate 5 million doses to women (subpopulation 1) and 5 million doses to men (subpopulation 2). As a result, ≈625, 000 women, but only ≈6250 men would still get the disease.

It may be doubted that this is a perfectly fair allocation: Why should it be fair to spend half of the available scarce resource for further improvements in a subpopulation that is much better off anyway with respect to the considered disease? Most healthcare decision makers (DMs) would probably tend in this situation to focus the resource allocation on that subpopulation that more urgently needs it. Short calculation shows that by spending a share of 100/101 on the female group and only a share of 1/101 on the male group, perfect outcome equity is achieved: In both subpopulations, the number of people still getting the disease would be about the same, namely ≈12, 400. This paper adopts the described outcome equity view without claiming that it is the ultimate solution to the difficult fairness quantification problem.

Our basic model follows that in Duijzer et al. (2018), but extends it by the outcome equity feature and by a more explicit consideration of the case of interaction between the patches. Moreover, we assume that the given stockpile of the vaccine is available already before or immediately after the onset of the epidemic outbreak, and that no later replenishment by new supply will take place. This assumption, made in several publications before, opens ample opportunities for generalizations in future research.

As the social welfare function, we employ a linear combination of average benefit and Gini's mean absolute difference of benefits. This measure, termed Gini aggregation function (GAF) in the sequel, is well founded in theory: It is characterized by a set of axioms and has been shown to be a natural extension of the Gini index, the possibly most prominent equity measure in the economic literature, to the simultaneous consideration of equity and efficiency (see Argyris et al., 2022; Porath & Gilboa, 1994). The weight given to the Gini term is an inequity aversion parameter, which can be interpreted and elicited.

In addition to the mentioned Gini‐based measure, we shall also investigate the simpler Rawlsian IAAF, which takes the outcome of the worst‐off individual as its evaluation criterion. It will be shown that the optimal solution under the Rawlsian IAAF typically coincides with a suitable modification of the simple pro rata (PR) vaccine allocation, the allocation applying the same vaccination fraction in each subpopulation. A main topic of our investigation will be when this modification, called adjusted pro rata (APR) allocation in the sequel, is locally or globally optimal under the GAF measure. We answer the question of local optimality by linear programming approach. The question of global optimality is more difficult to address, because it amounts to the solution of a nonconvex optimization problem, which we tackle by using the particle swarm optimization (PSO) metaheuristic.

Surprisingly, it turns out that for regional partition into subpopulations, without or with interaction, APR is typically even globally optimal for a large range of the inequity aversion parameter. Thus, for applications to geographical clusters, a recommendation of using APR can be given independently of the precise degree of inequity aversion, which may vary among DMs, as long as the DMs at least agree on a certain problem‐specific (usually small) minimum weight they are willing to give to the equity objective. For partition into age groups, this conclusion cannot be drawn anymore. We shall also address the case where instead of the number of infected individuals, the number of fatalities is the considered outcome quantity. Again, differences between regional and age‐based stratification can be observed.

Summarizing, the main contributions of this work are the following:

A fairness model for vaccine allocation in a pandemic is developed, based on the notion of outcome equity instead of input equity.

An objective function aggregating fairness and efficiency by means of the generalized Gini index is proposed, and the resulting optimization problem is analyzed within the context of a multipatch epidemiological model, both without and with interactions between the patches.

An APR allocation of vaccine is defined and shown to be a suitable choice especially in typical situations of a geographical segregation of the patches (subpopulations).

A new case study referring to the first phase of the COVID‐19 pandemic in Austria is introduced.

The paper is organized as follows: Section 2 shortly reviews related literature. Section 3 introduces the investigated vaccine allocation problem by deriving it from the epidemiological SIR model on the one hand and quantitative approaches to handle the equity issue in optimization on the other hand. In this context, also the APR vaccine allocation is defined. In Section 4, some analytical insights into the structure of the problem are presented. Section 5 deals with numerical methods for solving the problem under consideration. Section 6 extends the multipatch model to the case of interactions between the patches, while Section 7 generalizes it to limited vaccine efficacy and to patch‐specific weights. In Section 8, three case studies are provided. Section 9 discusses the results and limitations of the paper from the viewpoint of application, and Section 10 concludes the paper.

RELATED LITERATURE

There is a broad literature on quantitative approaches to optimal vaccine allocation based on suitable epidemiological models. For general introductions to epidemic modeling, see, for example, Brauer et al. (2019) and Daley and Gani (2001). Let us focus on those papers that are most relevant in the context of this paper, and outline their relations to the latter.

Hill and Longini Jr (2003) deal with the question, which vaccination fractions in a heterogeneous population composed of a certain number of groups are suitable in order to eliminate the possibility of an epidemic in the whole population (in other words, in order to keep the reproduction number below 1). After describing the threshold surface of critical vaccine allocations, they address the optimization problem of minimizing the total required number of vaccinations while remaining at threshold level.

A starting point for this investigation is Keeling and Shattock (2012). The authors ask for the optimal distribution of vaccine over spatially segregated subpopulations without or with interaction in the case of limited supply. They build their model on the classical final epidemic size formula by Kermack and McKendrick (1927), extended to the situation where the total population consists of several subgroups. It is shown that the distribution minimizing the total number of cases depends heavily on the size of the given stockpile. Some of these distribution policies are rather extreme insofar as they drastically advantage certain subgroups at the cost of other subgroups, which, as the authors admit, might be considered unethical.

Duijzer et al. (2016) define two vaccine allocation problems for multiple populations: (i) Find the vaccine allocations to subpopulations that maximizes the herd effect, that is, the number of individuals who escape infection without being vaccinated themselves. (ii) Find that critical vaccination coverage, that is, that vaccine allocation entailing a reproduction number of 1 in the whole population, that needs the least amount of vaccines. The authors show that problems (i) and (ii) are actually equivalent to each other, and provide an efficient algorithm (and for special cases even analytic solutions) for them.

For a sudden epidemic outbreak in a population partitioned in regional subgroups, Duijzer et al. (2018) address the problem of determining the optimal allocation of a vaccine stockpile that is insufficient to vaccinate the entire population. The focus of the work is on a situation without interaction between the subgroups, but also the case of interaction is touched. The authors identify a convex–concave property and other characteristics of the function quantifying the herd effect, elaborate on the role of the dose‐optimal vaccination fraction, and derive a simple heuristic for the vaccine allocation problem. Similarly as Keeling and Shattock (2012), the authors observe that mathematically optimal solutions are frequently inequitable, such that one has to choose between equity and efficiency. Nevertheless, the authors fully acknowledge the importance of the equity aspect in vaccine allocation and conclude that “new ideas are needed to reconcile equity and efficiency” in this area. This paper takes up the last‐mentioned research issue and exploits some useful analytical techniques from Duijzer et al. (2018) for investigating it. However, not only the extension to inequity aversion, but also the full treatment of the case of interaction between subgroups, require new analytical results that are provided in this paper.

Westerink‐Duijzer et al. (2020) consider a variant of the multipopulation situation where a central DM is lacking, but health agencies negotiate on sharing their available doses. The concept of the “core” from cooperative game theory is used to model this case. The authors use a case study on the redistribution of influenza vaccines to illustrate their approach.

As Duijzer et al. (2016) and some other articles, Enayati and Özaltın (2020) capture the vaccine allocation problem in the form of minimizing the number of vaccine doses needed to extinguish an emerging outbreak in its early stages, that is, to ensure that the reproduction number is less than or equal to unity. Their work achieves an essential progress on the one hand by proposing an exact solution approach, based on discretization with multiparametric disaggregation, on the other hand by incorporating equity constraints expressed in terms of the Gini coefficient. There are close relations to this paper, but three main differences have to be mentioned: (i) We consider outcome equity instead of input equity, which significantly increases the computational complexity. (ii) Instead of minimizing the number of doses required to extinguish the outbreak, we minimize the size of the epidemic for a given stockpile of doses. (iii) We do not use an equity constraint, but rather combine equity and efficiency to a compact (albeit parameterized) objective function satisfying well‐defined axioms.

In the context of the COVID‐19 pandemic, the trade‐off between efficiency and equity has also been addressed in some recent preprints as Bertsimas et al. (2020) or Chen et al. (2020). While these works use more involved models of the disease dynamics (adapted to COVID‐19) than this paper and also include the practically important aspects of timing and age groups, they put less emphasis on analytical insights into optimal vaccine allocation decisions. Insofar, this paper is rather in the line of articles as Duijzer et al. (2016) or Enayati and Özaltın (2020). Again in the context of COVID‐19, Munguía‐López and Ponce‐Ortega (2021) compare different fairness‐related objectives, as the Rawlsian criterion or the Nash solution from cooperative game theory, to the utilitarian objective where the sum of all individual utilities is maximized. The model does not specify an epidemiological model for representing the infection dynamics.

Pressman et al. (2021) aim at the development of an empirically measurable vaccine equity index for COVID‐19 vaccinations. In particular, the index should be able to quantify inequities between different racial or ethnic groups. For the purposes of this paper, Pressman et al. (2021) is especially relevant by the clear distinction the authors make between input equity and outcome equity (the authors use the terms “equality” and “equity” for these two notions, respectively), their proposed measure being of outcome equity type. Rastegar et al. (2021) build on Enayati and Özaltın (2020), but they modify the objective function definition by maximizing the minimum delivery‐to‐demand ratio per group and time period, a definition of Rawlsian type that captures efficiency and equity in one single function. This paper adopts a Gini index–based IAAF with the same aim; however, it also investigates a Rawlsian version of the problem, both within the framework of outcome equity rather than input equity. A further recent investigation on fair vaccine distribution can be found in Balcik et al. (2022).

For a review on models, methods, and applications of optimization under equity considerations, see Karsu and Morton (2015). Therein, a central notion of this paper, IAAFs, is discussed in detail, based on works as Kostreva et al. (2004). A special case of IAAFs are ordered weighted averages (OWAs), which occur under different names already early in the literature. For example, they are called generalized Gini indices in Weymark (1981) or symmetric comonotonically linear functionals in Ben‐Porath et al. (1997). In this paper, a particular class of OWAs will be used as social welfare functions allowing a simultaneous consideration of efficiency and equity: the class of linear combinations of the utilitarian measure (total or average benefit) with Gini's mean absolute difference (for the latter, see, e.g., Yitzhaki & Schechtman, 2013). This class is a parameterized set with a single parameter expressing the weight given to the Gini component. We call these functions GAFs. They are not new; for example, they have been investigated in detail in Porath and Gilboa (1994). Therein, besides a direct definition, also an axiomatic characterization in the form of seven axioms has been given. Recent applications of GAFs can be found in several papers as Gutjahr and Fischer (2018), Mostajabdaveh et al. (2019), or Gutjahr (2021). Moreover, also Eisenhandler and Tzur (2019a, 2019b) choose in their investigations a special GAF as their evaluation measure for inequity‐averse optimization. Another measure of OWA type, but different from GAFs, is used in Filippi et al. (2021).

Compared to the above‐cited literature, the new contribution of this paper to optimal vaccine allocation under fairness considerations can be seen mainly in the following issues: (a) The presented model is based on outcome equity rather than on input equity. (b) An IAAF of Gini type is applied as the objective function of the optimization approach in order to have a compact representation of both efficiency and equity. (c) An APR allocation with close connections to the Rawlsian solution is defined and shown to be optimal for a large range of preference profiles also under the Gini‐oriented evaluation perspective. (d) Numerical solution techniques for the proposed model are developed. (e) The model is extended to the practically relevant case of interactions between subgroups. (f) A new case study is presented, and two case studies from the literature are reanalyzed from a viewpoint including inequity aversion.

PROBLEM FORMULATION

Epidemiological model

The backbone of the epidemiological component of the proposed model is the well‐known SIR dynamics. As Duijzer et al. (2018) (cf. also the literature cited in Brauer et al., 2019, chapter 5), we apply it in a multipopulation context: It is assumed that the overall population decomposes into m subpopulations (patches), each of which is further subdivided into the three compartments “susceptible,” “infected,” and “removed.” The size of patch j is

$N_j$

, and

$N=N_1+\cdots +N_m$

denotes the size of the overall population. By

$s_j(t)$

,

$i_j(t)$

, and

$r_j(t)$

, the fractions of susceptible, infected, and removed individuals in patch j, respectively, are denoted. At every point t in time,

$s_j(t)+i_j(t)+r_j(t)=1$

. The SIR model is defined by the following system of differential equations:

1

Therein,

${\beta }_{jk}$

is the transmission rate from patch k to patch j, and γ is the recovery rate. (As some other authors, e.g., Ma & Earn, 2006, we assume here that the recovery rate does not depend on the patch, but the model can easily be generalized to patch‐dependent values of γ.) Furthermore, it is assumed that boundary conditions

$s_j(0)=s_{0j}$

,

$i_j(0)=i_{0j}$

, and

$r_j(0)=r_{0j}$

$(j=1,\text{…},m)$

are given.

As Duijzer et al. (2018), we suppose a situation where at a time τ in the early stage of an epidemic, a limited stockpile V of doses of a vaccine is available. For simplicity, we assume

$\tau =0$

. The DM has to find a distribution of the V available doses over the m subpopulations in a both fair and efficient way. It is assumed that vaccination leads to full immunity of the vaccinated person against the disease under consideration, that this immunity persists over the planning horizon, and that no repeated or booster vaccinations are necessary. A vaccinated individual is considered as sterile, that is, unable to transmit the disease to other individuals. Formally, the decision is described by the vector

$f=(f_1,\text{…},f_m)$

, where

$f_j$

denotes the fraction of people from subpopulation j that are to be vaccinated.

The basic output variables on which our health outcome evaluation rests are the variables

$y_j$

, where

$y_j=y_j(f)$

denotes the fraction of individuals from patch j who never get infected. This fraction is the sum of two components: (i) the fraction 

$f_j$

of vaccinated individuals in patch j and (ii) the fraction

$G_j(f)$

of individuals who, although not being vaccinated (i.e., remaining susceptible), will escape infection by the herd effect. The argument f in

$G_j(f)$

takes into account that the herd effect depends on the chosen vaccination distribution

$f=(f_1,\text{…},f_m)$

. Observe that

$G_j(f)$

is equal to the limiting fraction of susceptibles

${\lim }_{t\to \infty }{s}_j(t)$

resulting from vaccination distribution f. We get

$y_j(f)=f_j+G_j(f)$

.

From a utilitarian point of view, that is, when disregarding equity aspects, the resulting optimization problem is the following:

2

This is just the problem investigated in Duijzer et al. (2018), except that we divided the total number of individuals escaping infection by the constant N to get an average, which does not affect the optimal solution. We shall extend (2) by incorporating the issue of equity.

Before turning to the conceptual extension of (2), let us address the question how

$G_j(f)={\lim }_{t\to \infty }{s}_j(t)$

can be determined. In the literature, the determination of the final number of susceptibles has found much consideration (see, e.g., Ma & Earn, 2006), so we can be short about this point. The classical final size formulas describe the limiting fractions of susceptible and removed individuals (the number of infected individuals tends to zero), starting from an initial distribution

$(s_0,i_0,\mathbf{1}-s_0-i_0)$

for the relative sizes of the three compartments, where s ₀ and i ₀ are the vectors of the values

$s_{0j}$

and

$i_{0j}$

, respectively, and

$\mathbf{1}=(1,\text{…},1)$

. If at time

$\tau =0$

, vaccinations according to the vector f take place (we neglect the time needed for vaccination), the initial distribution above has to be replaced by

$(s(0),i(0),r(0))=(s_0-f,i_0,\mathbf{1}-s_0-i_0+f)$

, where

$(s_0,i_0,\mathbf{1}-s_0-i_0)$

is the distribution before vaccination. Following the lines of the derivation of the classical final size formulas from (1), one obtains the system of nonlinear equations

3

for

$j=1,\text{…},m$

, from which the values

$G_j(f)$

can be computed. Therein,

${\sigma }_j={\beta }_{jj}/\gamma$

. We will focus on the case

$s_{0j}+i_{0j}=1$

for all j of a “fresh” epidemic with no removed individuals in the beginning.

A frequently studied special case in the literature is that of no interaction between subpopulations. In this special case,

${\beta }_{jk}=0$

for

$j\ne k$

, and the function

$G_j(f)$

depends only on component

$f_j$

of the vector f, so formula (3) reduces to the mutually independent implicit equations

4

for the determination of

$G_j(f_j)$

(cf. Duijzer et al., 2018, online supplement S1.1).

In the case of no interaction between the subpopulations, the quantity

${\sigma }_j={\beta }_{jj}/\gamma$

equals the basic reproduction number within subpopulation j, the number of new infections caused by a single infectious individual in a completely susceptible subpopulation j. We always assume

${\sigma }_j>1$

for all j.

An alternative representation of

$G_j(f_j)$

is possible through the Lambert W function, which is defined as the solution

$W=W(x)$

of the equation 

$W{e}^W=x$

$(x\ge -1/e)$

:

5

Including equity

It has been recognized in the vaccine allocation literature that applying the utilitarian objective function “sum of the benefits of all individuals,” as expressed formally in the optimization model (2), can lead to highly inequitable and perhaps socially unacceptable vaccine distributions (cf. Duijzer et al., 2018; Enayati & Özaltın, 2020; Keeling & Shattock, 2012). This paper follows Enayati and Özaltın (2020) and other works in the choice to express inequity by a Gini‐related measure, but it deploys this concept in a quite different way. First of all, we rather strive for equity in health outcomes

$y_j$

(with the goal that similar fractions of individuals escape infection) instead of striving for equity in prevention allotments

$f_j$

(with the goal that similar fractions of vaccine are administered). The rationale behind the concept of outcome equity is the idea that a fair distribution of medical resources does not consist in providing everybody with the same amount of treatment, but in giving every patient the amount of treatment needed to cope with their individual degree of disease. For the specific case of vaccine allocation, this principle implies that subpopulations that are at a higher risk for infection should get larger fractions 

$f_j$

of the vaccine than subpopulations at a lower risk. For this reason, the presented model builds on the Gini index of the outcome quantities

$y_j$

rather than on that of the input quantities

$f_j$

. It has to be noted that because of the nonlinear dependence of

$y_j$

on

$f_j$

(see Equation 4), this will considerably increase the computational burden.

Second, instead of splitting the efficiency and the equity issue over objective function and constraints, respectively, we use a single compact IAAF as the overall aggregated assessment of a solution 

$f=(f_1,\text{…},f_m)$

. This is possible since, as mentioned above, the proposed model expresses now the equity component of the evaluation in the same (health outcome–related) measurement units as the utilitarian component. By choosing an IAAF that satisfies well‐established and plausible axioms, consistency and theoretical soundness of the model can be ensured, while a separated treatment of the two evaluation components may produce solutions that lack credibility, as observed by Matl et al. (2017) in the context of vehicle routing.

Gini aggregation functions

To take equity into account, we replace in (2) the objective function by a weighted sum of the utilitarian objective and a 1‐norm deviation of the outcomes across subpopulations:

6

Thus, the objective is now a weighted aggregation of the average benefit and the average of pairwise differences of benefits, with weights of 1 and (

$-\lambda$

), respectively, where the benefit

$b_{\nu }$

, representing the probability that individual ν will escape infection, is given as

$b_{\nu }=y_j(f)$

for individual ν in subpopulation j

$(j=1\text{…},m)$

. For reasons outlined below, we call the objective function in (6) a GAF. Objectives of this type are well founded in the theory of fair decisions: Our GAF can be rewritten as

7

Therein,

$\mu =(1/N){\sum }_{\nu =1}^N{b}_{\nu }$

is the average benefit,

$\lambda \ge 0$

is a parameter expressing the degree of inequity aversion,

$\mathrm{\Delta }=2\mu \mathcal{G}$

is Gini's mean absolute difference of the benefits

$b_{\nu }$

(see, e.g., Yitzhaki & Schechtman, 2013), and

$\mathcal{G}$

is the ordinary Gini index of the benefits

$b_{\nu }$

, which attains the extreme values 0 and 1 in the cases of full equity and full inequity, respectively. By using Gini's mean absolute difference Δ instead of the Gini index 

$\mathcal{G}$

in the weighted difference

$\mu -\lambda \mathrm{\Delta }$

, the quantities μ and Δ are measured in the same units, such that the weight λ can be dimensionless. From Porath and Gilboa (1994), Theorem D, it follows that the only IAAFs satisfying seven axioms on inequity‐averse preferences (see Supporting Information EC.4) are those given by the parameterized family (7), with

$0<\lambda <(1/2)·N/(N-1)$

. Choosing

$\lambda =1/2$

in (7) gives a special IAAF introduced in Eisenhandler and Tzur (2019a, 2019b). Since for large N, this is in a certain sense the most inequity‐averse measure satisfying the mentioned seven axioms, it is of particular interest. For

$\lambda =0$

, the other extreme, we get the utilitarian measure.1The measures of the family (7) are special cases of OWA functions (Argyris et al., 2022; Ben‐Porath et al., 1997).

While the choice of the inequity aversion parameter λ is up to the subjective preferences of the DM, it can be made understandable to the DM what a specific value of λ means by giving simple numerical examples. A way to elicit a rough estimate of λ from the DM is the following: Consider the case of two persons to whom the DM wants to provide benefits of some kind. Let y ₁ and y ₂ be the benefits allotted to persons 1 and 2, respectively, under distribution plan A, and let

$\stackrel{\sim }{y_1}$

and

$\stackrel{\sim }{y_2}$

be the corresponding benefits under distribution plan B. We write

$(y_1,y_2)\cong (\stackrel{\sim }{y_1},\stackrel{\sim }{y_2})$

if the DM is indifferent between plans A and B based on their GAF evaluations. The average value and Gini's mean absolute difference of the benefits under plan A are

$\mu =(1/2)(y_1+y_2)$

and

$\mathrm{\Delta }=(1/2)|y_1-y_2|$

, respectively, and analogously for plan B. By (7), we have

$(y_1,y_2)\cong (\stackrel{\sim }{y_1},\stackrel{\sim }{y_2})$

if

$\lambda =(\mu -\stackrel{\sim }{\mu })/(\mathrm{\Delta }-\stackrel{\sim }{\mathrm{\Delta }})$

.

For example, we may ask the DM how large the benefit x would have to be so that she/he is indifferent between the two distributions (60, 40) and

$(x,x)$

. The estimate of λ results then as

$(50-x)/(\frac{1}{2}·20)=(50-x)/10$

. For (60, 40)≅(47, 47), for example, one obtains

$\lambda =0.30$

. It can be expected that a social DM will hardly prefer (60, 40) over the completely fair and only marginally less efficient distribution (49,49). Therefore, in practice, typical values for λ will presumably be larger than 0.10.

In addition to the use of a GAF as the objective function, which gives Problem (6), it will also be interesting to investigate the use of the simplest IAAF, the measure aggregating the individual benefits 

$b_{\nu }$

by the function

$\min (b_1,\text{…},b_N)$

. This corresponds to the Rawlsian approach to distributional equity, which consists in maximizing the benefit of the worst‐off individual. Following Rawls, one obtains the optimization problem

8

An APR allocation

The well‐known pro rata (PR) allocation in vaccine allocation consists in allotting the available vaccine in proportion to the population sizes 

$N_j$

, yielding equal fractions in each subpopulation. In formal terms, this is the solution

$f=(f_1,\text{…},f_m)$

of problem (2) for which

$f_1=\cdots =f_m$

and

${\sum }_{j=1}^m{N}_j{f}_j=V$

, that is, the solution

$f_j=V/N$

for all

$j=1,\text{…},m$

. Of course, this solution can also be used as a feasible solution of (6) or (8). However, it is obvious that the motivation behind the PR allocation is rather inspired by the idea of input equity than by that of outcome equity. To establish an analogous notion for the context of outcome equity, one has to look for a solution that leads to equal outcome variables (fractions of never infected people) 

$y_j$

in all subpopulations

$j=1,\text{…},m$

: Definition 1

An APR allocation is a vector

$(f_1,\text{…},f_n)$

of vaccination fractions with

${\sum }_{j=1}^m{N}_j{f}_j=V$

that satisfies

$y_1(f)=\cdots =y_m(f)$

.

An APR allocation lets the second term in the objective function of (6) vanish and provides in this way an allocation that is perfect from the equity point of view. In the next section, it will be shown that at least in the more specific scenario without interaction, under mild conditions, an APR allocation exists and is unique. The relations between the APR allocation, the Rawlsian solution, and the GAF solution need a deeper investigation, which will be one of the topics of the following sections. An obvious other topic is of course the question how the three types of solutions can be determined.

ANALYTICAL RESULTS

This section derives analytical results on vaccine allocation under different solution concepts, still on the assumption of no interaction between the subpopulations. For easier orientation, Figure 1 gives an overview.

FIGURE 1

Schematic overview of the analytical results (case without interaction).

Existence and uniqueness of the APR allocation

For abbreviation, set

$$G_0^{max}=\underset{k=1,\text{…},m}{\max }{G}_k(0)\text{and}{s}_0^{min}=\underset{k=1,\text{…},m}{\min }{s}_{0k}.$$

9

$G_0^{max}$

is the maximum herd effect without vaccination in any subpopulation, and

$s_0^{min}$

is the minimum initial fraction of susceptibles in any subpopulation. Since we suppose that vaccination takes place in the very first phase of the epidemic,

$s_0^{min}$

is close to one, so we can safely assume

$G_0^{max}<s_0^{min}$

. Proposition 1

The following conditions are necessary conditions for the existence of an APR solution:

$$(\mathrm{a})V\ge \sum _{j=1}^m{N}_j{y}_j^{-1}\left(G_0^{max}\right),(\mathrm{b})V\le \sum _{j=1}^m{N}_j{y}_j^{-1}\left(s_0^{min}\right).$$

10

The proofs of all mathematical results can be found in Supporting Information EC.2.

Note that in the special case where all

$G_j$

are identical,

$y_j^{-1}(G_0^{max})=y_j^{-1}(G_j(0))=y_j^{-1}(y_j(0))=0$

. Thus, if the epidemiological dynamics in the single subpopulations do not differ too much, the right‐hand hand (r.h.s.) of (10)(a) will be small, that is, the bound is not restrictive. Also Equation (10)(b) is not very restrictive for applications: If the values

$s_{0j}$

are of a similar size, as it is the case if all

$s_{0j}$

are near to one, then the r.h.s. of (10)(b) is approximately equal to

${\sum }_{j=1}^m{N}_j{y}_j^{-1}(s_{0j})={\sum }_{j=1}^m{N}_j{s}_{0j}$

, and the last expression gives the total number of susceptibles at time 0, which cannot be exceeded by the number of vaccinations anyway.

The conditions above, together with

$G_0^{max}<s_0^{min}$

, are also sufficient: Proposition 2

If (10) and

$G_0^{max}<s_0^{min}$

hold, then the APR allocation exists and is unique.

In the following, we shall impose condition (10)(b) in its inessentially strengthened form

$$V<\sum _{j=1}^m{N}_j{y}_j^{-1}\left(s_0^{min}\right).$$

11

Whenever it exists, the APR allocation can be computed within very short runtime by Newton's method, see Supporting Information EC.3. In Section 6, the numerical computation of the APR solution will be explained for the more general case where there is interaction between the patches.

Rawlsian solution and APR allocation

In the context of vaccine allocation optimization disregarding equity, Duijzer et al. (2018) showed that it is optimal to always use the complete vaccine stockpile, in other words, that the first constraint in (2) is satisfied with equality in an optimal solution. It turns out that on the not very restrictive condition (11), also in the Rawlsian context, an optimal solution exhausts the available stockpile: Lemma 1

Suppose that the stockpile size V satisfies (11). Then, in an optimal solution of the Ralwsian problem (8), the constraint

${\sum }_{j=1}^m{N}_j{f}_j\le V$

is met with equality.

With the help of this lemma, we can show that on mild conditions, Rawlsian solution and APR allocation are equivalent: Theorem 1

Suppose that the stockpile size V satisfies (10)(a) and (11). Then, the Rawlsian optimization problem (8) possesses a unique optimal solution, namely, the APR allocation.

The GAF solution

We turn now to the use of a GAF as the objective function, as described by problem formulation (6). Let us start by the observation that also the GAF solution, as the Rawlsian solution, completely exhausts the available stockpile: Lemma 2

Suppose that the stockpile size V satisfies (11). Then, in an optimal solution of the GAF problem (6), the constraint

${\sum }_{j=1}^m{N}_j{f}_j\le V$

is met with equality.

Next, we discuss the relation of the GAF solution to the APR allocation. By definition, the APR allocation reduces the second term of the objective of (6) to the minimal possible value of zero. Thus, the APR allocation corresponds in a certain sense to the limiting case

$\lambda \to \infty$

of (6). (Note, however, that we cannot simply drop the first term of the objective function to describe this limiting case, as long as we do not explicitly require that the first constraint is satisfied with equality, because otherwise we would also get extremely inefficient solutions such as

$f_j=0$

$\forall j$

as optimal solutions.)

Intuitively, one would expect that as λ is reduced from ∞ to a large finite value and further on through a sequence of decreasing numbers, such that gradually more and more weight is shifted from equity to efficiency, the optimal solution vector f should, gradually again, move from the APR allocation to solution vectors taking the efficiency aspect more and more into account. This intuition turns out to be wrong, as we will show. Actually, mixing an increasing concern for efficiency into the objective function of (6) by reducing λ will at first not impair the (exact) optimality of the APR allocation at all. Only for a typically already rather low value of λ, the solution will jump to a vector that pays attention to the efficiency criterion (and change then gradually in the case of a further reduction of λ). Example 1

Let

$m=3$

,

$(N_1,N_2,N_3)=(1000,2000,3000)$

,

$\sigma =(3,2,4)$

,

$V=1800$

, and

$s_{0j}=0.99$

$(j=1,\text{…},3)$

. Figure 2 shows the plot. By Lemma 2, an optimal solution has

$N_1{f}_1+\cdots +N_3{f}_3=V$

, which enables us to reduce the dimension of the search space to the horizontal plane with axes f ₁ and f ₂. We determine f ₃ accordingly, and plot

${\mathrm{\Psi }}_{\lambda }(f)$

. Objective function values below 0.35 have been cut. The two pictures correspond to the values

$\lambda =0.05$

and

$\lambda =0.20$

. The APR allocation results as

$f^{\ast }=(0.3165,0.1759,0.3772)$

, it is the same for both cases, as well as its function value η. It can be seen that while the APR allocation point in the left part of the plot is still dominated for

$\lambda =0.05$

, it has become clearly superior for

$\lambda =0.20$

.

FIGURE 2

Plot of the Gini aggregation function (GAF) objective function 

${\mathrm{\Psi }}_{\lambda }(f)$

from (6) for Example 1 for

$\lambda =0.05$

(left) and

$\lambda =0.20$

(right).

Local and global optimality of the APR allocation

Motivated by the last examples, it is tempting to ask under which conditions the APR allocation already solves the GAF problem to optimality. The next lemma confirms that a threshold value for λ decides on this issue. This is valid both for the question of local and for that of global optimality (with different threshold values). In the sequel,

$f^{\ast }$

will always denote the (unique) APR allocation. Lemma 3

If for some

$\lambda \ge 0$

, the APR allocation

$f^{\ast }$

is a global (resp. local) maximizer of 

${\mathrm{\Psi }}_{\lambda }(f)$

, then

$f^{\ast }$

is also a global (resp. local) maximizer of 

${\mathrm{\Psi }}_{{\lambda }^{\prime }}(f)$

for any

${\lambda }^{\prime }>\lambda$

.

While the results in this section generally refer to the case without interaction, a look at the proof of Lemma 3 shows that the lemma holds in the case of interaction as well.

The next result shows that by choosing λ sufficiently large, the APR allocation can always be made globally (and hence also locally) optimal. Although this property may be expected, as for large λ, the second term of the objective function of (6) dominates the first term, proving the result rigorously for the general case is surprisingly intricate.2 Theorem 2

To each instance of the vaccine allocation problem (6), there is a

$\lambda \ge 0$

such that the APR allocation

$f^{\ast }$

is a global maximizer of 

${\mathrm{\Psi }}_{\lambda }(f)$

.

Theorem 2, together with Lemma 3, suggests the following definition: Definition 2

Let

$f^{\ast }$

be the APR allocation.

• (a)

  The local APR optimality threshold is the smallest λ for which APR is locally optimal:

  ${\lambda }^{loc}=inf\{\lambda \ge 0:f^{\ast }\text{local maximizer of}{\mathrm{\Psi }}_{\lambda }\}$
  .

• (b)

  The global APR optimality threshold is the smallest λ for which APR is globally optimal:

  ${\lambda }^{\ast }=inf\{\lambda \ge 0:f^{\ast }\text{global maximizer of}{\mathrm{\Psi }}_{\lambda }\}$
  .

Thus, by Lemma 3, the APR allocation is locally optimal w.r.t. the GAF measure for

$\lambda >{\lambda }^{loc}$

, and it is not locally optimal w.r.t. the GAF measure for

$\lambda <{\lambda }^{loc}$

. Analogously for global optimality. Obviously,

${\lambda }^{loc}\le {\lambda }^{\ast }$

, since each global optimum is also a local optimum.

Definition 2 extends in a formally identical way to the case where interaction is present, with the only difference that in this case, it can happen that the set for which

$f^{\ast }$

is a (local) maximizer can become empty, such that the infimum attains the value ∞. In the case without interaction, this is excluded by Theorem 2. If in the presence of interaction,

${\lambda }^{loc}$

is finite, the property that the APR allocation is locally optimal w.r.t. the GAF measure for

$\lambda >{\lambda }^{loc}$

and not locally optimal w.r.t. the GAF measure for

$\lambda <{\lambda }^{loc}$

remains valid; analogously for 

${\lambda }^{\ast }$

.

For numerical purposes, a slight extension of Definition 2(b) is necessary: With

$y_1(f_1^{\ast })=\cdots =y_m(f_m^{\ast })=\eta$

, we have also

${\mathrm{\Psi }}_{\lambda }(f^{\ast })=\eta$

for each λ. Then according to Definition 2(b), with S denoting the feasible set,

${\lambda }^{\ast }=inf\{\lambda :{\max }_{f\in S}{\mathrm{\Psi }}_{\lambda }(f)=\eta \}$

. However, since in a fixed point or floating point arithmetic, each numerical function evaluation of

${\mathrm{\Psi }}_{\lambda }(f)$

has only limited computational accuracy, the validity of

${\max }_{f\in S}{\mathrm{\Psi }}_{\lambda }(f)=\eta$

cannot be ascertained on the computational level. As usual in numerical analysis, the statement has to be replaced by

$|{\max }_{f\in S}{\mathrm{\Psi }}_{\lambda }(f)-\eta |\le \overline{\epsilon }$

or, equivalently,

${\max }_{f\in S}{\mathrm{\Psi }}_{\lambda }(f)\le \eta +\overline{\epsilon }$

, with some small tolerance parameter

$\overline{\epsilon }>0$

. In this way, we get the following extension of Definition 2(b): Let

$f^{\ast }$

be the APR allocation,

$\eta ={\mathrm{\Psi }}_{\lambda }(f^{\ast })$

, and

$\overline{\epsilon }\ge 0$

. Then the

$\overline{\epsilon }$

‐tolerant global APR optimality threshold is the number

$${\lambda }_{\overline{\epsilon }}^{\ast }=inf\{\lambda \ge 0:\underset{f\in S}{\max }{\mathrm{\Psi }}_{\lambda }(f)\le \eta +\overline{\epsilon }\}.$$

12

It should be noted that

${\lambda }_{\overline{\epsilon }}^{\ast }\le {\lambda }^{\ast }$

. In particular, it can happen that

${\lambda }_{\overline{\epsilon }}^{\ast }<{\lambda }^{loc}$

.

For the computation of

${\lambda }^{loc}$

, the following result will be crucial: Lemma 4

Let V satisfy (11), and let 

$f^{\ast }$

satisfy the constraints

$0\le f_j^{\ast }\le s_{0j}$

$\forall j$

with strict inequalities. Then,

$\lambda \ge {\lambda }^{loc}$

if and only if the system of linear constraints

13

for the variables

${\nu }_{k\ell }$

$(k,\ell =1,\text{…},m)$

is consistent. Therein,

$$\begin{array}{c}\hfill \xi =\left(N·\prod _{j=1}^m{\rho }_j\right)/\left(\sum _{j=1}^m{N}_j\prod _{\ell \ne j}{\rho }_{\ell }\right)\text{with}{\rho }_j=y_j^{\prime }\left(f_j^{\ast }\right).\end{array}$$

NUMERICAL SOLUTION

In this section, we indicate numerical solution methods for (a) the problem (6) of finding the optimal vaccine allocation w.r.t. a GAF, (b) the problem of determining the local and the global APR optimality thresholds,

${\lambda }^{loc}$

and

${\lambda }^{\ast }$

. We confine us here to a rough outline; the reader can find all details in the Supporting Information.

It has to be mentioned that larger instances of problem (a) cannot be expected to be solvable to optimality because of its nonlinear and even nonconvex nature: Note that already the simpler problem of optimizing vaccine allocation without including a nondifferentiable equity term in the objective function has been recognized as computationally too challenging for exact optimization approaches (see Duijzer et al., 2018; Westerink‐Duijzer et al., 2020).3 Therefore, a two‐track approach has been followed here: First, an enumerative algorithm was developed. This algorithm is based on a suitable grid‐type discretization and enables the solution of smaller problem instances (up to ≈5 patches) to a desired degree of accuracy, providing rigorous lower and upper bounds for the solution value of the maximization problem. In order to identify the bounds, we need to know to what extent the function values in the continuous domain can deviate from those in the grid points. For this purpose, we determine a Lipschitz constant of the objective function and use it to calibrate the mesh size of the grid. The details are described in Supporting Information EC.5.

Second, a PSO metaheuristic (Poli et al., 2007) for the numerical solution of larger instances was implemented and assessed by comparison to the results of the enumerative algorithm on small instances. The PSO algorithm builds on an appropriate encoding of the solution space, together with a repair mechanism for infeasible solutions. The description of the algorithm as well as results on its quality evaluation are provided in Supporting Information EC.6.

From an applied point of view, one would wish to be able to specify upper bounds

$t_j$

on the fraction 

$f_j$

of people to be vaccinated that may be smaller than the initial fraction 

$s_{0j}$

of susceptibles, considering the case where a fraction

$1-t_j$

of the population refuses vaccination or cannot be vaccinated because of medical reasons. Both the discretization approach and the PSO can easily be modified to the case where the constraints

$0\le f_j\le s_{0j}$

are tightened to constraints

$0\le f_j\le t_j$

with

$t_j\le s_{0j}$

. The APR allocation 

$f^{\ast }$

is not affected by bounds

$t_j<s_{0j}$

, as long as 

$f^{\ast }$

satisfies them itself.

Contrary to problem (a), the first of the problems listed under (b), namely, the determination of the local optimality threshold 

${\lambda }^{loc}$

, is still amenable to an exact solution by linear programming. The key tool is here Lemma 4: Obviously, consistency of (13) can efficiently be judged by applying an linear program (LP) solver to (13) with an arbitrary objective function. The local optimality threshold 

${\lambda }^{loc}$

can then easily be determined by bisection search, using Lemma 3 and repeated calls of the solver.

Finally, we addressed the harder problem of the estimation of the global optimality threshold, where we relied again on the discretization approach as well as on the PSO metaheuristic. This led to Algorithm 1 and Algorithm 2, respectively, which are presented in Supporting Information EC.7. Algorithm 1 determines a lower and an upper bound on 

${\lambda }_{\overline{\epsilon }}^{\ast }$

based on the discretized version of the problem and improves then these bounds by narrowing the interval between them. Algorithm 2 repeats attempts to find a tight lower bound on

${\lambda }_{\overline{\epsilon }}^{\ast }$

by applying the PSO to the problem with a λ value that is chosen as near as possible to the supposed optimality threshold. As for each heuristic, the produced solution value is not an unbiased estimate of the exact solution value, but a bound—in our case, a lower bound, so we denote it by LB

$({\lambda }_{\overline{\epsilon }}^{\ast })$

. The comparisons to the intervals determined by Algorithm 1, however, show that the bound LB

$({\lambda }_{\overline{\epsilon }}^{\ast })$

is typically very good (see Supporting Information EC.7) and can therefore be used as a guideline in applications, where it does not make a big difference whether a solution better than the APR allocation cannot be found or does not exist. Example 2

Teytelman and Larson (2013) investigate a vaccine allocation model for the 2009 H1N1 influenza epidemic in the United States. In Table EC.11 in the Supporting Information, the numbers

$N_j$

of inhabitants of

$m=10$

different regions as well as the estimated basic reproduction numbers 

${\sigma }_j={\sigma }_j^{(0)}$

in the 10 regions are given. In order to investigate also possible other pathogens that are more infectious, we consider four additional variants by setting

$\sigma =p·{\sigma }^{(0)}$

, where σ⁽⁰⁾ are the original parameters in Table EC.11, and

$p=1,\text{…},5$

represent the cases of increasingly infectious epidemics. Moreover, let

$V=N/4$

and

$s_{0j}=0.99$

for all j. An application of the above‐described algorithm yields the values 

${\lambda }^{loc}$

for the local optimality thresholds, as well as heuristic estimates (precisely: lower bounds) LB

$({\lambda }_{\overline{\epsilon }}^{\ast })$

for the

$\overline{\epsilon }$

‐tolerant global optimality thresholds (

$\overline{\epsilon }=0.005$

) as indicated below:

14

The resulting optimality threshold values are all distinctly below 0.10, in particular: far below the value of

$\lambda =0.5$

used in Eisenhandler and Tzur (2019a), which allows the conclusion that in each of these five application instances, only a DM whose inequity aversion is comparably weak would have an incentive to deviate from the APR allocation.

Also for 12 randomly generated test instances investigated in Supporting Information EC.6, a majority of 75% produced PSO‐based lower bounds (and even upper bounds) below 0.10, as seen there.

An investigation of the influence of varying stockpile size V on the optimal allocation, analyzed at an example taken from Keeling and Shattock (2012), can be found in Supporting Information EC.8.

INTERACTION

In this section, the concepts and procedures of the previous sections will be extended to the case where there is interaction (positive infection rates) between the subpopulations. Two types of application may be of specific interest: a first one where subpopulations are defined by geographical regions or political districts and a second one where subpopulations are formed by age groups. (Also combinations may be envisaged, see Enayati & Özaltın, 2020.) It will be seen that the results for these two types show marked differences, even after the application of the same concepts. The section will focus on a generalization of the APR allocation and of the solution optimizing the GAF. The relations between APR allocation and the Rawlsian solution in this more general context is a topic of future research.

Computation of the APR allocation in the presence of interactions

An APR allocation has been defined by Definition 1 already in a context including the possibility of interactions. We are now interested in a way how to compute the APR allocation also in the presence of interactions. For this purpose, we build on Equation (3). Slight reformulation and use of the assumption

$s_{0j}+i_{0j}=1$

$(j=1,\text{…},m)$

transforms Equation (3) to the system of equations

15

Now suppose that the characteristic property

$y_1(f)=\cdots =y_n(f)$

of an APR allocation is satisfied. Denote the common value of

$y_j(f)$

by η. Then, the herd effect satisfies

$G_j(f)=\eta -f_j$

for all j. With

$$b_j=(1/\gamma )·\sum _{k=1}^m{\beta }_{jk}(j=1,\text{…},m),$$

16

this turns Equation (15) into

$\log ((\eta -f_j)/(s_{0j}-f_j))=(\eta -1)b_j$

$(j=1,\text{…},m)$

. Resolving for

$f_j$

, we get

$$f_j={\stackrel{\sim }{f}}_j(\eta )=\frac{\eta -s_{0j}\exp (-b_j(1-\eta ))}{1-\exp (-b_j(1-\eta ))}(j=1,\text{…},m).$$

17

By the formula above, an intended fraction 

$\eta \in [0,1[$

of never infected people determines a feasible vaccine allocation

$f=(f_1,\text{…},f_m)$

that achieves this η if and only if

${\sum }_{j=1}^m{N}_j{f}_j\le V$

and

$0\le f_j\le s_{0j}$

$\forall j$

. We look for the largest fraction

$\eta ={\eta }^{\ast }$

that can be achieved by a feasible allocation. This value

${\eta }^{\ast }$

is the solution of

$$\begin{array}{c}\hfill \max \{\eta |\sum _{j=1}^m{N}_j{\stackrel{\sim }{f}}_j(\eta )\le V,0\le {\stackrel{\sim }{f}}_j(\eta )\le s_{0j}(j=1,\text{…},m)\}.\end{array}$$

18

Lemma 5

The following properties hold:

• (a)

The function 

$\eta \mapsto {\stackrel{\sim }{f}}_j(\eta )$

is strictly increasing in

$\eta \in [0,1]$

.

(b)

The function

$\eta \mapsto \eta -s_{0j}\exp (-b_j(1-\eta ))$

is concave in

$\eta \in [0,1]$

.

Lemma 5 can be used as follows for efficiently solving (18): First, observe that constraint

${\stackrel{\sim }{f}}_j(\eta )\le s_{0j}$

reduces to

$\eta \le s_{0j}$

. Next, constraint

${\stackrel{\sim }{f}}_j(\eta )\ge 0$

is equivalent to

$\eta -s_{0j}\exp (-b_j(1-\eta ))\ge 0$

. Because of Lemma 5(b), the latter constraint is satisfied for the η values in an interval

$[{\eta }_j^{min},{\eta }_j^{max}]$

, the bounds of which can easily be determined by bisection search. Finally, by Lemma 5(a), the left‐hand side (l.h.s.) of the first constraint of (18) is strictly increasing in 

$\eta \in [0,1[$

, so the constraint restricts the feasible values to an interval 

$[0,{\eta }^V]$

, the upper bound 

${\eta }^V$

of which can again be immediately computed by bisection search. In total, the feasible set for η results as the intersection of all the indicated intervals, and the solution value

${\eta }^{\ast }$

is the maximum value of this intersection interval. The optimal allocation results as

$f^{\ast }=\stackrel{\sim }{f}({\eta }^{\ast })$

.

If the intersection of the

$m+1$

intervals [0, s ₀] and

$[{\eta }_j^{min},{\eta }_j^{max}]$

$(j=1,\text{…},m)$

contains 

${\eta }^V$

, then

${\eta }^{\ast }={\eta }^V$

, and the allocation 

$f^{\ast }$

determined as above satisfies not only the condition

$y_1(f^{\ast })=\cdots =y_n(f^{\ast })$

, but it also exhausts the stockpile V, so it is an APR allocation. Typically this is the case. However, for very small or very large V, it may also happen that full equality of the

$y_j(f)$

cannot be reached by vaccination, for the same reasons as discussed in Section 4.1 in the more specific interaction‐free context. In practice, one may start with determining

${\eta }^V$

by bisection search, and check then whether the corresponding f fulfills the constraints

$0\le f_j\le s_{0j}$

for all j.

The proposition below derives monotonicity properties concerning the influence of transmission rates and of initial fractions of susceptibles on the APR allocation. Regarding its second part, recall that as a consequence of the assumption of no removed individuals at the beginning, a smaller number of susceptibles at time 0 entails a larger number of infectious persons already at time 0. Proposition 3

• (a)

With

$b_j$

as in (16), suppose that j and

$j^{\prime }$

are two subpopulations with

$N_j=N_{j^{\prime }}$

,

$s_{0j}=s_{0j^{\prime }}$

, and

$b_j<b_{j^{\prime }}$

. Then, in an APR allocation,

$f_j<f_{j^{\prime }}$

, that is, subpopulation 

$j^{\prime }$

gets a larger per‐person fraction of vaccine than subpopulation j. In particular, in the special case without interaction, a larger reproduction number

${\sigma }_j$

implies, ceteris paribus, a larger vaccine allocation 

$f_j$

.

(b)

Suppose that j and

$j^{\prime }$

are two subpopulations with

$N_j=N_{j^{\prime }}$

,

$b_j=b_{j^{\prime }}$

, and

$s_{0j}<s_{0j^{\prime }}$

. Then, in an APR allocation,

$f_j>f_{j^{\prime }}$

, that is, subpopulation 

$j^{\prime }$

gets a smaller per‐person fraction of vaccine than subpopulation j.

It can be conjectured that the GAF solutions converge to the APR allocation 

$f^{\ast }$

as

$\lambda \to \infty$

. In the absence of interaction, this follows from Theorem 2. A weaker version of this statement, based on the key property that

$f^{\ast }$

satisfies

$y_k(f^{\ast })-y_{\ell }(f^{\ast })=0$

for all k and ℓ, is valid also in the presence of interaction: Proposition 4

If the APR allocation exists, then as

$\lambda \to \infty$

, the maximizers

$f_{\lambda }^{opt}$

of

${\mathrm{\Psi }}_{\lambda }$

converge to equal outcomes 

$y_j$

in the following sense:

$y_k(f_{\lambda }^{opt})-y_{\ell }(f_{\lambda }^{opt})\to 0(\lambda \to \infty )$

for all

$k=1,\text{…},m,\ell =1,\text{…},m$

.

Computation of the GAF solution in the presence of interactions

If there is interaction, the solution of the optimization problem (6) becomes especially challenging because already the determination of the objective function value

${\mathrm{\Psi }}_{\lambda }(f)$

for some given vaccine allocation f may become time‐consuming: Note that now, the values

$G_j(f)$

(from which the

${\mathrm{\Psi }}_{\lambda }(f)$

are finally obtained via the

$y_j(f)$

) result as the solutions of the system of nonlinear equations (3). In the literature, different ways how to get the herd effect sizes

$G_j(f)$

have been suggested, for example, to approximate them by direct numerical simulation of the trajectory of the differential equations (1), or to reformulate (3) as an iterative vector equation producing a series of solutions that converge to the fixed point of the vector equation. We tested both approaches and found them of moderate efficiency. To speed up the computation, the following approach was used: In order to get a good initial solution, numerical simulation of the ordinary differential equation (ODE) trajectory with a medium number of iterations was applied. Then, this solution was improved by the multidimensional Newton–Raphson method. Details are presented in Supporting Information EC.9. By this approach, very accurate solutions of the nonlinear system (3) were obtained within comparably short time. Furnished with the described procedure for determining the objective function values 

${\mathrm{\Psi }}_{\lambda }(f)$

as a subroutine, the PSO algorithm outlined in Section 5 worked well also in the presence of interactions.

OTHER GENERALIZATIONS

Vaccine efficacy

Up to now, it was assumed that vaccination is 100% effective and leaves complete immunity. In particular, this assumption entails also that a vaccinated person cannot transmit the disease to any other individual. For many infectious diseases, this is an oversimplification. Therefore, we relax the assumption along the lines of Hill and Longini Jr (2003) (see also Duijzer et al., 2016, online appendix D) by assuming that a vaccination achieves its goal only in a certain fraction

${\phi }_j$

of cases, where the efficacy factor

${\phi }_j\in [0,1]$

can possibly depend on the subpopulation j. The vaccination at time

$t=0$

changes then the initial distribution

$(s_j(0),i_j(0),r_j(0))=(s_{0j},i_{0j},1-s_{0j}-i_{0j})$

in patch j to

$(s_{0j}-{\phi }_j{f}_j,i_0,1-s_{0j}-i_{0j}+{\phi }_j{f}_j)$

. Again, it is supposed that

$r_j(0)=0$

, such that

$s_{0j}+i_{0j}=1$

.

As Duijzer et al. (2016) remark, the introduction of efficacy factors simply rescales the parameters 

$f_j$

, replacing them by values

${\phi }_j{f}_j$

in all formulas except the resource constraint

${\sum }_j{N}_j{f}_j\le V$

, which has to be kept unchanged. In particular, in the optimization problem (6), the expressions

$y_j(f)$

in the objective function have to be replaced by

$y_j({\phi }_1{f}_1,\text{…},{\phi }_m{f}_m)$

, which still allows the application of the PSO‐based solution algorithm. Concerning the APR allocation, the only effect of a replacement of

$f_j$

by

${\phi }_j{f}_j$

is that the r.h.s. of (17) has still to be divided by

${\phi }_j$

. Based on the new values

${\stackrel{\sim }{f}}_j(\eta )$

obtained in this way, program (18) remains formally identical, and it is easily seen that also Lemma 5 extends to the more general situation.

Patch‐specific weights

One might be interested in being able to assign different weights to the numbers of infected cases in the single subpopulations. Let

$w_j$

denote the weight factor assigned to the case that a person from subpopulation j is infected at some time. The factor 

$w_j$

can mean different things, for example, it can stand for the probability of an infected person to be hospitalized, the probability that an infected person needs intensive care, or the infection fatality rate (IFR), that is, the probability of an infected person to die from the disease. To simplify the diction, let us use in the following the interpretation of 

$w_j$

as the probability of hospitalization of an infected person.

Instead of the outcome variables 

$y_j$

, the weighted model builds on the outcome variables

$z_j=w_j(1-y_j)$

representing, for example, the fractions of hospitalized persons in patch j

$(j=1,\text{…},m)$

. The absolute number of hospitalized persons in patch j is then

$N_j{z}_j(f)$

. The variables

$z_j$

are used consistently in both terms of the objective function. One obtains the following problem:

19

Program (19) is a generalization of program (6), where the latter results as the special case 

$w_j=1$

for all j, rewritten as a maximization problem. Observe that (19) minimizes the number of hospitalizations instead of maximizing the number of infection‐free individuals; consequently, the equity penalty term has to be added instead of being subtracted. The case of different hospitalization rates or IFRs in the subpopulations makes the advantage of the outcome equity concept especially clear: Input equity would not provide a criterion for giving more vulnerable groups of individuals a larger chance of being vaccinated.

An APR allocation can be defined for the described weighted case analogously as for the original unweighted case: Generalizing Definition 1, we call a vector

$f=(f_1,\text{…},f_m)$

of vaccination fractions with

${\sum }_{j=1}^m{N}_j{f}_j=V$

an APR allocation for the weighted model, if it satisfies

$z_1(f)=\cdots =z_m(f)$

. The special choice

$w_j=1$

for all j leads again back to the original definition.

To compute the APR allocation for the weighted model, we set now

$z_1=\cdots =z_m=\zeta$

, such that

$z_j=w_j(1-y_j)$

entails

$G_j=y_j-f_j=1-\zeta /w_j-f_j$

$(j=1,\text{…},m)$

. Inserted into (15), this yields

20

By resolving (20) for

$f_j$

, one obtains

$$\begin{array}{c}\hfill f_j={\widehat{f}}_j(\zeta )=\frac{1-\zeta /w_j-s_{0j}\exp (-{\overline{b}}_j\zeta )}{1-\exp (-{\overline{b}}_j\zeta )}(j=1,\text{…},m).\end{array}$$

21

As before, the last formula can be generalized to vaccine efficacy values

${\phi }_j\le 1$

by replacing the denominator on the r.h.s. of (21) with

$[1-\exp (-{\overline{b}}_j\zeta )]·{\phi }_j$

, analogously as in the generalization of (17).

Finally, the optimal ζ is the determined as the solution of the following counterpart of (18):

$$\begin{array}{c}\hfill \min \left\{\zeta \mid \sum _{j=1}^m{N}_j{\widehat{f}}_j(\zeta )\le V,0\le {\widehat{f}}_j(\zeta )\le s_{0j}(j=1,\text{…},m)\right\}.\end{array}$$

22

Susceptible‐exposed‐infectious‐removed (SEIR) and SI ⁿ R models

Our approach can also be extended from SIR models to SEIR models, or, more generally, to SI ⁿ R models. SEIR models address latency by including an “exposed” state E, while SI ⁿ R models assume n different infectious stages with different transmission parameters. As remarked in Ma and Earn (2006), SEIR is a special case of SI ⁿ R, so it suffices to consider SI ⁿ R. Final size formulas for SI ⁿ R in the form of systems of equations are available for particular situations; see, for example, Theorem 9.2 in Ma and Earn (2006), which gives conditions ensuring that the final size formula for a multipatch SI ⁿ R model is of the same form as that for a multipatch SIR model. In such situations, our approach can be applied unchanged. In the general case, an (approximate) APR allocation can be determined as follows: As an initial solution, compute (along the lines of Section 6.1) the APR allocation of the corresponding SIR model based on the final size formula, and fine‐tune this solution then by minimizing the inequity term of (6) on the constraint of stockpile exhaustion, using PSO, where

${\mathrm{\Psi }}_{\lambda }(f)$

is obtained by direct numerical simulation of the SI ⁿ R differential equations. The GAF solution can be computed analogously.

CASE STUDIES

Case Austria

In this subsection, we use publicly available data4 from the first wave of the COVID‐19 pandemic in Spring 2020 in the nine provinces (“federal states”) of Austria. We take this wave as an instance of a severe and highly contagious epidemic and estimate all parameters needed for our model from the real‐world data, though we do not include all features of COVID‐19 in our case description. In particular, we assume an early availability of a vaccine, which was actually not yet available for COVID‐19 in 2020, and we assume that vaccination needs not to be repeated, but leaves long‐term immunity at least for a certain fraction of vaccinated individuals, an assumption that is now well known to be wrong for COVID‐19, but holds for some other infectious diseases. Thus, the goal of this case study is not to provide a fully adequate model for the very complex COVID‐19 pandemic, but to analyze a hypothetical situation based on parameter values as they could be valid for some future epidemic.

First of all, from the time series of daily incidences, daily estimates of the reproduction numbers for March 15 to April 3, 2020, in each of the nine provinces have been estimated (see Table 1). Their average values over the 20 days lie between between

$\mathcal{R}=1.59$

and

$\mathcal{R}=1.85$

. Let us denote these numbers by

${\sigma }_j^{emp}$

$(j=1,\text{…},9)$

. The overall estimated average reproduction number for Austria is

$\mathcal{R}=1.64$

.

TABLE 1

Test cases Austria 1 and 2: The nine provinces (“federal states”) of Austria with their population sizes, estimated reproduction numbers, and case fatality rates (CFRs; in %) in the first wave of COVID‐19 (2020).

Index j	Province	$N_j$	${\sigma }_j^{emp}$	CFR (in %)
1	Burgenland	294,436	1.64	0.425
2	Carinthia	561,293	1.85	0.651
3	Lower Austria	1,684,287	1.63	0.415
4	Upper Austria	1,490,279	1.63	0.346
5	Salzburg	558,410	1.76	0.361
6	Styria	1,246,395	1.66	0.572
7	Tyrol	757,634	1.70	0.245
8	Vorarlberg	397,139	1.72	0.249
9	Vienna	1,911,191	1.59	0.387

To obtain estimates of the transmission rates 

${\beta }_{jk}$

, a generalizable procedure was used, which is described in Supporting Information EC.10. In addition to the empirical reproduction numbers, it exploits the geographical distances between the provinces (cf. Figure 3) by assuming that the contact rates between individuals from different provinces are inversely proportional to these distances. Together with a calibration of the resulting model to the empirically observed reproduction numbers, this allows an estimation of the matrix 

$\beta =({\beta }_{jk})$

.

FIGURE 3

Map of the provinces (“federal states”) of Austria.

The reproduction number derived from the observed incidence data during March 15 to April 3, 2020, cannot be interpreted as the basic reproduction number 

${\mathcal{R}}_0$

of COVID‐19 in Austria, since during most of the time of the considered 20 days, nationwide contact restrictions had already come into effect, such that 

${\mathcal{R}}_0$

is expected to be larger. We did a sensitivity analysis by generating two versions of the case: Instance Austria 1 sticks to the assumption of an epidemic with

$\mathcal{R}=1.64$

, while in instance Austria 2, all parameters 

${\beta }_{jk}$

were increased by 50%, such that

$\mathcal{R}$

equals now 2.46.

A vaccine efficacy of

${\phi }_j=0.55$

for all j is assumed. Moreover, to make the test cases more realistic with respect to the willingness of people to be vaccinated, we imposed an upper bound of 

$t_j=0.90$

on all vaccination fractions 

$f_j$

. Concerning vaccine availability, we assume that the stockpile V is sufficient for 50% of the population, that is,

$V/N=0.5$

.

A. Results for the unweighted model. Table 2 shows the quantities below. Note that the PR allocation is

$f^{PR}=(0.5,\text{…},0.5)$

in the assumed case of

$V/N=0.5$

.

• $f^{util}=arg{\max }_f{\mathrm{\Psi }}_0(f)$
  , the utilitarian allocation, that is, the optimal GAF solution for
  $\lambda =0$
  ,
  TABLE 2

  Test cases Austria 1 and 2, results. Upper tables: utilitarian allocation (left) and adjusted pro rata (APR) allocation (right). Lower table: objective function values.

  Allocation	Province	Austria 1	Austria 2
  $f^{util}$	1	0.5272	0.8588
  	2	0.6241	0.0000
  	3	0.4921	0.8415
  	4	0.4771	0.8206
  	5	0.5703	0.0000
  	6	0.4950	0.0000
  	7	0.5142	0.0000
  	8	0.5256	0.0000
  	9	0.4560	0.8148

  Allocation	Province	Austria 1	Austria 2
  $f^{\ast }$	1	0.5207	0.5130
  	2	0.6078	0.5671
  	3	0.4923	0.4952
  	4	0.4803	0.4877
  	5	0.5606	0.5379
  	6	0.4957	0.4974
  	7	0.5126	0.5080
  	8	0.5227	0.5143
  	9	0.4627	0.4766

  Quantity	Austria 1	Austria 2
  ${\mathrm{\Psi }}_0(f^{util})$	0.8175443	0.5209979
  ${\mathrm{\Psi }}_{0.1}(f^{util})$	0.8172	0.4900
  ${\mathrm{\Psi }}_{0.2}(f^{util})$	0.8168	0.4590
  ${\mathrm{\Psi }}_{\lambda }(f^{\ast })(={\eta }^{\ast })$	0.8175139	0.5104511
  ${\mathrm{\Psi }}_0(f^{PR})$	0.8162	0.5106
  ${\mathrm{\Psi }}_{0.1}(f^{PR})$	0.8136	0.5088
  ${\mathrm{\Psi }}_{0.2}(f^{PR})$	0.8111	0.5070

• ${\mathrm{\Psi }}_{\lambda }(f^{util})$
  , the GAF objective value of the utilitarian allocation, for
  $\lambda \in \{0,0.1,0.2\}$
  ,

• $f^{\ast }$
  , the APR allocation,

• ${\mathrm{\Psi }}_{\lambda }(f^{\ast })={\eta }^{\ast }$
  , the GAF objective value of the APR allocation (identical for all λ),

• ${\mathrm{\Psi }}_{\lambda }(f^{PR})$
  , the GAF objective value of the ordinary PR allocation 
  $f^{PR}$
  , for
  $\lambda \in \{0,0.1,0.2\}$
  .

We can derive the following observations from Table 2:

Test case Austria 1:

• (i)

APR almost reaches the efficiency (measured by the utilitarian objective, the fraction of infection‐free individuals) of the utilitarian solution 

$f^{util}$

. The difference is only

$3.04·{10}^{-5}$

, that is, three more individuals from 100,000 get infected under APR. The PR allocation

$f^{PR}$

is distinctly less efficient at Austria 1, with a difference to the utilitarian solution of 130 more individuals from 100,000 getting infected under PR.

(ii)

For moderate inequity aversion of

$\lambda =0.1$

, the GAF objective function value of

$f^{util}$

is slightly worse than that of APR, for

$\lambda =0.2$

already distinctly worse. Similar statements hold for PR compared to APR at

$\lambda =0.1$

and

$\lambda =0.2$

, but the differences are here more pronounced.

Test case Austria 2:

• (i)

At this test instance, where contact rates have been increased by 50%, the difference of the fraction of infection‐free individuals between

$f^{util}$

and APR is larger than at Austria 1, it amounts to 1055 individuals from 100,000. For PR, the difference to

$f^{util}$

is almost the same as for APR.

(ii)

The quality of 

$f^{util}$

, measured by the GAF 

${\mathrm{\Psi }}_{\lambda }$

, decreases here steeply when moving from

$\lambda =0.1$

to

$\lambda =0.2$

. The PR solution, on the other hand, is only moderately worse than APR for these increased λ values.

(iii)

It is striking that

$f^{util}$

looks here entirely different from

$f^{\ast }$

. In the utilitarian solution, only the inhabitants of four provinces (Burgenland, Lower Austria, Upper Austria, and Vienna) would be vaccinated. From Figure 3, it is seen that this is a connected cluster of provinces in the north‐eastern part of Austria. Vaccine would be completely refused to the western and the southern provinces.

The left graph of Figure 4 shows the dependence of

${\max }_f{\mathrm{\Psi }}_{\lambda }(f)$

on the inequity aversion parameter λ for test case Austria 1. It can be observed that this function converges rather fast to the GAF value

${\eta }^{\ast }$

of the APR allocation. An optimality threshold

${\lambda }^{\ast }\approx 0.025$

can be read off from the figure. The corresponding graph for Austria 2 (which is not shown here) behaves similarly, with an optimality threshold of

${\lambda }^{\ast }\approx 0.035$

.

FIGURE 4

Test case Austria 1. Left graph: unweighted case, plot of the optimal Gini aggregation function (GAF) objective function value 

${\mathrm{\Psi }}_{\lambda }(f)$

in dependence of λ (solid blue curve), and objective function value of the adjusted pro rata (APR) allocation (height of the dashed red line). Right graph: weighted case, plot of the optimal GAF loss function value 

${\widehat{\mathrm{\Psi }}}_{\lambda }(f)$

in dependence of λ (solid blue curve), and loss function value of the APR allocation (height of the dashed red line). Values on the second axes in %.

B. Results for the weighted model. Let us now analyze the behavior of the model for Austria 1 under patch‐specific weights, as described in Section 7. The weights 

$w_j$

are interpreted in the following as IFRs and are approximated by case fatality rates (CFRs). Table 1 contains CFRs as observed for COVID‐19 for the nine Austrian provinces, computed as ratios of cumulated fatalities to cumulated cases (Wikipedia, 2022, data from August 1, 2022). In this interpretation, the variables 

$z_j=z_j(f)$

represent province‐specific long‐term death rates by the disease under a chosen vaccine allocation f. It should be noted that in view of symptom‐free or very mild disease processes, the CFR systematically overestimates the IFR, but the former is nonetheless frequently used as a proxy for the latter since in the first stage of an epidemic, it is almost impossible to reliably estimate the IFR (see Wong et al., 2013). In our model, a long as the proportionality factor between IFR and CFR is constant across the patches, the use of the CFRs as weights instead of the IFRs does not change the resulting allocations

$f^{\ast }$

,

$f^{util}$

etc., nor the resulting 

${\lambda }^{\ast }$

. Only the interpretation of

$N_j{z}_j$

as the number of fatalities in patch j is then of course a pessimistic estimate.

For test case Austria 1, one obtains the following solutions, where

$f^{util}=arg\min {\widehat{\mathrm{\Psi }}}_0(f)$

denotes again the utilitarian allocation:

Notably, the utilitarian solution allocates no vaccine at all to two provinces (Tyrol and Vorarlberg). Furthermore, one gets

${\widehat{\mathrm{\Psi }}}_0(f^{util})=0.000690$

, and

${\zeta }^{\ast }={\widehat{\mathrm{\Psi }}}_0(f^{\ast })=0.000710$

. Recall that

${\widehat{\mathrm{\Psi }}}_0$

represents loss, which explains why

${\zeta }^{\ast }>{\widehat{\mathrm{\Psi }}}_0(f^{util})$

. The ordinary PR allocation

$f^{PR}$

produces a loss function value of

${\widehat{\mathrm{\Psi }}}_0(f^{PR})=0.000757$

.

The right graph of Figure 4 plots the loss function

${\min }_f{\widehat{\mathrm{\Psi }}}_{\lambda }(f)$

in dependence of λ. Again, we see fast convergence to

${\zeta }^{\ast }$

, the objective function value of the APR allocation. An optimality threshold of

${\lambda }^{\ast }\approx 0.16$

can be observed.

The numerical results above allow it to estimate the price of fairness (Bertsimas et al., 2011), represented here by the difference between

${\widehat{\mathrm{\Psi }}}_0(f^{util})$

and

${\widehat{\mathrm{\Psi }}}_0(f^{\ast })$

, which is

$\approx 2·{10}^{-5}$

for case Austria 1. Thus, the APR allocation entails ≈2 additional fatalities per 100,000 inhabitants compared to the utilitarian allocation 

$f^{util}$

in a conservative estimate equalizing CFR and IFR (see the remark above). The ordinary PR allocation turns out as distinctly worse than APR from a utilitarian point of view: From the indicated numbers, it computes as ≈6.7 additional fatalities per 100,000 inhabitants, compared to

$f^{util}$

.

The downside of the utilitarian allocation, on the other hand, is that fatalities would be distributed very unjustly over the provinces: Tyrol, for example, would have to bear

$N_7·{\zeta }^{\ast }\approx 534$

fatalities under APR, but

$N_7·z_7(f^{util})\approx 1117$

fatalities under allocation 

$f^{util}$

, which grants this province no vaccine at all, that is, 583 fatalities more. It is hard to imagine that such an allocation decision would be accepted by the population, so it is not a realistic option (as to this aspect, cf. also Bertsimas et al., 2020). Contrary to that, the APR allocation can be explained to the public as a means to satisfy basic equity demands.

Considering the actual practice in many countries during the COVID‐19 pandemic, it rather seems that political DMs shied away from a utilitarian policy minimizing the total number of infections or of deaths, a policy that would necessitate the mentioned harsh inequalities in vaccine allocation, and opted instead for a broad distribution of vaccine, usually guided by the ordinary PR rule. For example, the European Union committed itself to “ensure that each country receives doses based on a PR population distribution key” (Bolcato et al., 2021). As our results for Austria 1 show, however, the PR policy may entail a much larger death toll than the utilitarian policy. Compared to PR, the GAF solution (which is in the case

$\lambda \ge {\lambda }^{\ast }$

identical to the APR allocation) saves lives. In Section 8.4, it will be shown that this observation can be generalized to some extent from the case Austria 1 to other instances.

Case Teytelman and Larson 2013

A further application case, obtained by extending the case from Teytelman and Larson (2013) (which we already addressed in Example 2) to the situation of interactions between patches, is described in Supporting Information EC.11. The results are similar to those for Austria 1 and 2.

Case Wallinga 2010

In order to investigate also the situation where the population is not stratified into spatial clusters but rather into age groups, let us reanalyze an application case presented in Wallinga et al. (2010) and taken up again in Duijzer et al. (2016). Six different age groups are distinguished, and age group sizes 

$N_j$

and transmission rates 

${\beta }_{jk}$

have been estimated. We reproduce the data in Supporting Information EC.12. The recovery rate is assumed as 

$\gamma =0.286$

. The reproduction number without vaccination results as 

$\mathcal{R}=2.0$

.

As with the test case referring to Austria, we investigate two different versions: a version Wallinga 1 where the original transmission rates 

${\beta }_{jk}$

are used, and a “stress test variant” Wallinga 2, where all

${\beta }_{jk}$

have been increased by 50% in order to model a more infectious strain of the pathogen. As before, a vaccine efficacy of

${\phi }_j=0.55$

for all j is assumed, the upper bound to all vaccination fractions 

$f_j$

is set to 0.9, and a vaccine availability V covering 50% of the population size N is considered.

A. Results for the unweighted model. The following solutions are obtained for test case Wallinga 1:

On can see that the utilitarian solution does not vaccinate the youngest and the oldest age group (age 0–5 and 60+, respectively), and focuses on the age groups 13–19 and 20–39, where the proposed allocation goes to the upper limit 0.9 of the vaccination rate. The APR solution is more balanced, but also emphasizes older children and young adults.5Moreover, one obtains

Evidently, APR is more efficient under the utilitarian evaluation criterion than the ordinary PR allocation, but APR is now outperformed by 

$f^{util}$

under objectives Ψ_0.1 or Ψ_0.2. The optimality threshold 

${\lambda }^{\ast }$

is comparably large ( > 1.1) for Wallinga 1. For Wallinga 2, the threshold

${\lambda }^{\ast }$

is smaller, but still distinctly larger than 0.2. This behavior completely differs from what one observes in Austria 1 or Austria 2 (and in TL 1 and TL 2 in Supporting Information EC.11). We have to conclude that for age‐based stratification, APR is not as favorable as it is for regional stratification.

In view of the large

${\lambda }^{\ast }$

value for Wallinga 1, one may be interested in the trade‐off between efficiency and equity for intermediate values of λ, that is, values with

$0<\lambda <{\lambda }^{\ast }$

. By varying the weight λ within this range and computing the optimizers of the corresponding functions

${\mathrm{\Psi }}_{\lambda }(f)$

, one obtains the supported efficient solutions of the biobjective optimization problem

$\max (\mu (f),-\mathrm{\Delta }(f))$

with μ and Δ as in (7), where a “supported” solution is defined in multiobjective optimization as a Pareto‐efficient solution that can be obtained by the optimization of a weighted sum of the objectives. It is well known that in general, not all efficient solutions are supported. We get a heuristic approximation to the set of supported efficient solutions by applying the PSO‐based optimization algorithm to diverse values of λ within the interesting range. The result is shown in Figure 5. In addition to the found efficient points 

$(\mu (f),\mathrm{\Delta }(f))$

, where inequity is measured by Gini's mean absolute difference, we also plot the corresponding points 

$(\mu (f),\mathcal{G}(f))$

where the Gini index 

$\mathcal{G}(f)=\mathrm{\Delta }(f)/(2\mu (f))$

is used to quantify inequity.

FIGURE 5

Test case Wallinga 1, unweighted case. Blue upward triangles: supported efficient solutions for the biobjective problem with objective functions μ and Δ. Red downward triangles: same solutions with the second coordinate representing the Gini index

$\mathcal{G}$

instead of Gini's mean absolute difference Δ.

Starting from the rightmost point on the frontier, which corresponds to the utilitarian solution, one can see that already a small drop in efficiency by less than 0.5% suffices to achieve a reduction of the Gini index by almost 20%. Near the leftmost point of the frontier, which corresponds to the APR allocation, the marginal gain in equity bought by a reduction of efficiency is weaker. This effect is consistent with related observations in other publications (cf. Chen et al., 2020).

B. Results for the weighted model. As the Wallinga et al. test case supposes an influenza epidemic, weights representing IFRs should reflect the mortality rates of influenza. These rates typically grow fast with the age (an exception was the Spanish flu 1918–1920), with values between 1 and 1000 fatalities per 10⁵ persons (see Wong et al., 2013). The IFRs of different influenza outbreaks vary to a large extent and the partitions into age groups in the available statistics usually differ from that in Wallinga et al. (2010). In order to get a representative scenario of special values for age‐specific OFRs, we assume an IFR of

$3·{10}^{-5}$

for the youngest age group in the Wallinga case, a geometric growth of the IFR from each age group to the next, and a final IFR of

$300·{10}^{-5}$

for the oldest age group. This choice, which is realistic compared to the empirical results cited in Wong et al. (2013), defines the weight vector in our model as

$w=(3.0,7.5,18.9,47.5,119.4,300.0)$

. Let us confine ourselves to the discussion of case Wallinga 1.

It turns out that an APR allocation does not exist anymore for Wallinga 1, weighted as described above. This is not quite surprising: The huge mortality differences between age group 1 and age group 6 cannot be equalized anymore by a suitable vaccination strategy relying on a limited stockpile of vaccine. One may ask whether the goal of providing each age group with the same excess mortality rate from the epidemic is meaningful at all. We cannot expect that the health‐related disadvantage of the elderly can be completely overcome by medical interventions. This consideration questions the application of established fairness measures to age groups, while such measures may be completely appropriate and justifiable when applied to regional subpopulations.

Whereas

$f^{\ast }$

does not exist in the described case, it is possible (using the proposed PSO procedure) to compute

$f^{util}=arg{\min }_f{\widehat{\mathrm{\Psi }}}_0(f)$

or

$arg{\min }_f{\widehat{\mathrm{\Psi }}}_{\lambda }(f)$

for some

$\lambda >0$

. It turns out that already

$f^{util}$

, and the more the allocations

$arg{\min }_f{\widehat{\mathrm{\Psi }}}_{\lambda }(f)$

, shifts the available vaccine as far as possible to the older age groups:

$f^{util}=(0,0,0,0.4316,0.9,0.9)$

(recall the assumed upper bound 0.9 on the vaccination rate).

Utilitarian comparison of PR to APR

Although the purpose of APR is not to yield a good utilitarian evaluation but to maximize outcome fairness, some results in the case studies in Sections 8.1–8.3 suggest that the APR allocation provided by the weighted model might possibly be superior at least to the ordinary PR allocation also w.r.t. the utilitarian assessment by

${\widehat{\mathrm{\Psi }}}_0$

. The current subsection investigates at numerical instances whether this observation can be generalized. In other words, under which circumstances can we expect that replacing PR by APR reduces the number of fatalities? Formally, this amounts to the question when the loss difference

${\widehat{\mathrm{\Psi }}}_0(f^{PR})-{\widehat{\mathrm{\Psi }}}_0(f^{\ast })$

is larger than zero (then PR is worse than APR), and when it is smaller than zero (then PR is better). We focus on the weighted model, which allows it to predict fatalities, but add also some results on the unweighted model that predicts infections.

Comparison results for special cases with only

$m=2$

patches are presented in Supporting Information EC.13. It turns out that APR is usually, but not always better than PR. In order to verify this result also for larger test cases, a total number of 28, 000 test instances with

$m=10$

patches were generated randomly. Different parameter combinations were studied, with an attempt to choose realistic parameter values by keeping them in the magnitude of the values observed in the three case studies. We generated the instances as follows: (i) For

$V/N$

, the two alternatives

$V/N=0.25$

and

$V/N=0.50$

were investigated. (ii) Parameters

${\sigma }_j$

were taken uniformly at random from intervals

$[{\sigma }_{min},{\sigma }_{max}]$

, with the alternatives [1.5, 2.0], [1.7, 1.8], [2.0, 2.5], and [2.2, 2.3] for the intervals. Interaction was modeled by defining a transmission matrix β with

${\sigma }_j-c$

in the main diagonal (

$j=1,\text{…},10$

) and values

$c/9$

in the other entries, where c was set to

$c=0.03$

similarly as in the Teytelman and Larson test case (cf. EC.11 in the Supporting Information). (iii) Parameters

$w_j$

were selected uniformly at random from intervals

$[0.5-a,0.5+a]$

with

$a\in \{0.000,0.025,0.050,0.075,0.100,0.125,0.150\}$

. This produces

$2·4·7=56$

combinations. For each combination, 10 random choices of

$(N_1,\text{…},N_{10})$

were made by drawing each

$N_j$

uniformly at random and normalizing their sum to the fixed value 1000. For each of the resulting distributions

$(N_1,\text{…},N_{10})$

, 50 test instances were generated by drawing the values for

${\sigma }_j$

and

$w_j$

, so that each of the 56 combinations obtained 500 test instances.

Table 3 summarizes the results. We see that only in 3 out of the 56 combinations, PR was superior to APR. These are combinations with small V and a small variation of the weights

$w_j$

. As soon as the interval length for the weights reaches 10% or more of the mean, APR is, in each combination, superior in more than 90% of cases.

TABLE 3

Results for 500 random test instances per entry, each instance with

$m=10$

. The entries show the fractions of cases in which adjusted pro rata (APR) is better than pro rata (PR) w.r.t. the utilitarian measure (weighted model). The parameters

${\sigma }_j$

are taken from the intervals indicated in the second column; the parameters

$w_j$

are taken from

$[0.5-a,0.5+a]$

.

		a
	${\sigma }_j$ in	0.000	0.025	0.050	0.075	0.100	0.125	0.150
$V/N=0.25$	[1.5, 2.0]	0.856	0.862	0.902	0.932	0.972	0.980	0.988
	[1.7, 1.8]	0.578	0.916	0.992	1.000	1.000	1.000	1.000
	[2.0, 2.5]	0.006	0.394	0.910	0.986	1.000	1.000	1.000
	[2.2, 2.3]	0.158	1.000	1.000	1.000	1.000	1.000	1.000
$V/N=0.50$	[1.5, 2.0]	1.000	1.000	1.000	1.000	1.000	1.000	1.000
	[1.7, 1.8]	0.996	0.994	0.984	0.956	0.934	0.932	0.936
	[2.0, 2.5]	0.998	1.000	1.000	0.996	0.996	0.996	0.996
	[2.2, 2.3]	0.894	0.954	0.990	0.996	1.000	1.000	1.000

Similarly, the average sizes of the difference

${\widehat{\mathrm{\Psi }}}_0(f^{PR})-{\widehat{\mathrm{\Psi }}}_0(f^{\ast })$

over the 500 test instances have been computed for each of the 56 combinations in Table 3. For the sake of brevity, we omit the detailed results. The smallest entry was

$-0.248·{10}^{-3}$

for the combination

$V/N=0.25$

,

${\sigma }_j\in [2.0,2.5]$

and

$a=0$

, while the largest entry was

$7.784·{10}^{-3}$

for the combination

$V/N=0.25$

,

${\sigma }_j\in [2.0,2.5]$

and

$a=0.150$

. Thus, the results for

$m=10$

confirm the effect (already seen in the results for

$m=2$

) that an occasional dominance of PR over APR is only very weak, while the dominance of APR over PR is usually more accentuated.

The general tendency to a superiority of APR observed in Table 3 is weakest in the case

$a=0$

, that is,

$w_j=0.5$

for all j. Obviously, this case is equivalent to the unweighted model where the outcome variables are defined by the fractions of (not) infected individuals. It may be interesting to check whether in this case a significant superiority of APR can be stated at all. To clarify this issue, a third series of numerical experiments was carried out. For each of three values of the variation of 

${\sigma }_j$

, a sample of 1000 test instances was generated as follows: (a) N was set to 1000, and 50 different random vectors

$(N_1,\text{…},N_{10})$

with sum equal to N were generated as described under (2) above. (b) For each resulting

$(N_1,\text{…},N_{10})$

, 20 test instances were randomly generated according to

$V\in [250,500]$

uniformly at random,

$\overline{\sigma }\in [1.5,2.5]$

uniformly at random, and

${\sigma }_j\sim \text{normal}(\overline{\sigma },{\text{std}}^2)$

$(j=1,\text{…},m)$

, where std = 0.1, 0.2, or 0.3. (c) The transmission matrix β was derived from the values

${\sigma }_j$

as under (2), again with

$c=0.03$

. For the number

$n_{APR}$

of cases (out of 1000 cases) where APR was better than PR, the following results were obtained:

$n_{APR}=629$

, 720, and 766, for std = 0.1, 0.2, and 0.3, respectively. This shows a not very strong, but statistically highly significant superiority of APR over PR even for the unweighted model.

DISCUSSION AND POLICY IMPLICATIONS

The main message for practice derived from our numerical investigations above is that for epidemics of the considered type, already a comparably small weight given to the equity term in an aggregated evaluation of the vaccine allocation to geographical clusters typically tends to make an APR policy optimal. The necessary adjustment of PR is prompted not by efficiency considerations, but by a replacement of input equity with outcome equity: The goal is equal chances of escaping infection rather than equal chances of getting vaccinated. The proposed APR allocation can be computed by a fast numerical approach. It turns out that under regional stratification with comparably low contact rates between the regions, both the investigated artificially generated instances and the real world–based case studies usually reveal the APR solution as already best possible also in the GAF evaluation, except for DMs with an uncommonly low commitment to the equity goal.

While these results suggest that a healthcare DM faced with geographical subgroups can expect to be on the safe side by a suitably APR policy, this does not hold anymore for age‐based subgroup partitioning of the overall population. In the latter case, an APR allocation needs not to be optimal or may even not exist at all. By estimating the indices proposed in this paper, namely, local and global APR optimality thresholds, the DM can recognize such situations and, if necessary, apply the proposed heuristic procedure to get a tailored solution to the specific allocation problem.

As shown in Section 7, the presented model also allows it to take account of different weights of an infection case in different subpopulations, for example reflecting different rates of severe illness or of mortality. Empirically, marked regional differences in mortality rates can often be observed. Possible reasons might be differing speed of access to hospitals in urban and rural areas, or differing ways of dealing with diseases in different socioeconomic groups. If numbers of fatalities instead of numbers of infections are chosen as the outcome variables, the question of which price of fairness the society should accept becomes especially controversial. This question is first and foremost an ethical one and has to be discussed on the level of values. However, it has also a pragmatic aspect: In our case study Austria 1, it could be seen that the utilitarian solution (minimizing the number of fatalities) would completely deprive two provinces of vaccine. No matter how this is evaluated from an ethical point of view, it seems rather clear that it is not politically implementable in view of the resistance it would provoke. On the other hand, the ordinary PR solution is frequently rather inefficient, as indicated by the results in Section 8.4. The APR allocation may be a good compromise solution between utilitarian and PR allocation in such situations.

Some limitations of this work have to be mentioned. First, the extension from the classical SIR model to SEIR or SI ⁿ R would deserve an experimental investigation. Second, the assumption of vaccine supply at one time rather than by a series of deliveries should be relaxed by suitable model extensions. A natural way to deal with the latter case is to apply our approach repeatedly at different times, using at each time the observed state as the new initial condition (after relaxing the assumption that initially, no individuals are in the state “removed”). Of course, this does not guarantee optimality, and the obtained algorithm should be evaluated by comparison to other dynamic allocation algorithms. Similar remarks are in place for the related case where (as with COVID‐19) repeated vaccinations are necessary to achieve a sufficient degree of immunity.

A practical limitation of the presented approach concerns the possible difficulty to explain the solution aiming at outcome equity to the population. Admittedly, outcome equity is not quite as easy to defend by political DMs as input equity, since the prediction of outcome equity measures is model based and needs the involvement of epidemiological experts. Nevertheless, the difficulty regarding political communication would still be distinctly aggravated by a model‐based approach disregarding equity.

CONCLUSIONS

We extended models for optimal vaccine allocation to subpopulations by using social welfare functions that integrate the utilitarian criterion (“sum of all benefits”) with a fairness criterion (“equality of benefits”), based on the probabilities of being spared from infection. Three concepts were investigated: (a) an outcome‐based adjustment of the popular PR policy of vaccine allocation, (b) the Rawlsian solution (“best outcome for the worst‐off”), and (c) the allocation optimizing a linear combination of the utilitarian objective and Gini's mean absolute difference, called GAF. On the theoretical level, we showed that (a) and (b) are basically equivalent. Moreover, also (a) and (c) tend to produce the same solutions, provided that the degree of inequity aversion, the weight parameter used in the linear combination defining the GAF, is not too low. We defined threshold values for this parameter, above which the APR allocation simultaneously optimizes the GAF, and proposed procedures for computing good estimates for these thresholds. Moreover, exact as well as heuristic solution algorithms for computing the optimal vaccine allocation also in the cases where it is different from APR have been provided. Furthermore, a generalization to subpopulation‐specific weights of infections was investigated, which allows it to analyze questions about fractions of hospitalization or of fatalities.

Numerical results differed for the cases of spatial and of age‐based stratification. In the first case, the optimality threshold is typically low, which recommends the APR allocation as a both fair and efficient policy. In the second case, the APR allocation only optimizes the combined objective function if the DM's inequity aversion is comparably high, so that losses in terms of fairness have to be taken into account in order to achieve an optimal total health benefit. However, the question how “fairness between age groups” should really be measured may remain subject to discussion.

Footnotes

1

For convenience, we shall include in the definition of the GAF the boundary case

$\lambda =0$

of the utilitarian measure, though for

$\lambda =0$

, the GAF is not an IAAF anymore.

2

There may possibly be connections of this result to the theory of optimization with sparsity‐inducing norms, as it has been studied for parameter fitting applications (see, e.g., Bach et al., 2011). The aim is there to get sparse models by letting only a limited subset of the model parameters deviate from zero. This is managed by extending the objective function with a penalty term containing the 1‐norm of the parameter vector. Our problem formulation is not the same as that in sparse optimization, as we do not apply the 1‐norm to the coefficients

$f_j$

themselves, but to differences of nonlinear functions of the coefficients 

$f_j$

. Moreover, these differences are obviously not independent from each other, which further complicates the situation. Nevertheless, in some sense, the APR allocation can be viewed as the case of perfect “sparsity” where all the mentioned differences vanish, and Theorem  guarantees that this can be achieved by giving the penalty term a sufficiently large weight λ. Considering the parallels with sparse optimization, we conjecture that a replacement of the 1‐norm in the penalty term of the GAF definition by some other

${\ell }_p$

norm only preserves the validity of Theorem 2 if

$p\le 1$

.

3

Because of the convex–concave nature of the objective function of (2), already the special case

$\lambda =0$

of our problem belongs to the class of knapsack problems with S‐shaped return function, for which Ağralı and Geunes () have shown NP‐hardness.

4

Federal Ministery of Social Affairs, Health, Care and Consumer Protection, Republic of Austria; see also Wikipedia ().

5

The reader should be aware that this refers to infection rates as the outcome variables; for the consideration of severe courses of the disease, see the results for the weighted model below.