Modeling the relative risk of incidence and mortality of select vaccine-preventable diseases by wealth group and geographic region in Ethiopia

Immunization is one of the most effective public health interventions, saving millions of lives every year. Ethiopia has seen gradual improvements in immunization coverage and access to child health care services; however, inequalities in child mortality across wealth quintiles and regions remain persistent. We model the relative distributional incidence and mortality of four vaccine-preventable diseases (VPDs) (rotavirus diarrhea, human papillomavirus, measles, and pneumonia) by wealth quintile and geographic region in Ethiopia. Our approach significantly extends an earlier methodology, which utilizes the population attributable fraction and differences in the prevalence of risk and prognostic factors by population subgroup to estimate the relative distribution of VPD incidence and mortality. We use a linear system of equations to estimate the joint distribution of risk and prognostic factors in population subgroups, treating each possible combination of risk or prognostic factors as computationally distinct, thereby allowing us to account for individuals with multiple risk factors. Across all modeling scenarios, our analysis found that the poor and those living in rural and primarily pastoralist or agrarian regions have a greater risk than the rich and those living in urban regions of becoming infected with or dying from a VPD. While in absolute terms all population subgroups benefit from health interventions (e.g., vaccination and treatment), current unequal levels and pro-rich gradients of vaccination and treatment-seeking patterns should be redressed so to significantly improve health equity across wealth quintiles and geographic regions in Ethiopia.

Introduction comprehensive multi-year immunization plans (cMYP), which encompass all components of immunization services (e.g., service delivery, vaccine supply, quality and logistics, disease surveillance, advocacy, social mobilization and communication, and program management) [12,13]. However, effectively translating these strategies and plans into action has been challenging and community engagement has declined in recent years, with waning capacity and acceptability among WDA leaders and low acceptance by community members [12,13].
Since Ethiopia's introduction of the EPI in 1980, immunization coverage has progressively increased; however, there remains a substantial coverage gap and inequality in coverage by socioeconomic stratum and by geographic location [6,7,9,14,15]. Nationally, the proportion of children who were fully immunized was 41% in 2019, but this varies drastically [7]. Only 20% of children in the Afar region (a predominately pastoralist population) are fully vaccinated, compared to 83% in Addis Ababa (an entirely urban population) [7]. Similarly, disparities in immunization coverage are also observed by wealth quintile [7].
Chang and colleagues [10,16] have proposed a first modeling approach to assess the distributional impact of vaccines on VPD incidence and mortality; however, their methodology did not yet account for the joint distribution of risk factors (i.e., the overlap between risk factors, that is accounting for individuals with multiple risk factors) and did not include variation by geographic region, neither did it examine the country of Ethiopia specifically. In this paper, we build on this preliminary work and significantly augment it in assessing the relative distribution of VPD incidence and mortality across socioeconomic strata and by geographic region in Ethiopia, taking into account the joint distribution of risk factors for these diseases.

Methods
We examined four major VPDs in Ethiopia: rotavirus diarrhea, human papillomavirus (HPV), measles, and pneumonia (pneumococcal and Haemophilus influenzae type b, Hib), as well as their five corresponding vaccines: rotavirus vaccine (RV, against rotavirus diarrhea), HPV vaccine, measles-containing vaccine (MCV), pneumococcal conjugate vaccine (PCV, against pneumococcal pneumonia), and penta-3 vaccine (against Hib pneumonia). We evaluated all vaccines in Ethiopia's EPI for inclusion in the analysis; however, only rotavirus diarrhea, measles, pneumonia, and HPV were selected as VPDs for inclusion due to data availability. Data on risk and prognostic factors for these diseases (with the exception of HPV) is available from the Demographic and Health Survey (DHS) disaggregated by wealth quintile and region [9].

Risk and prognostic factors
For each VPD, we selected a number of risk and prognostic factors. Diarrhea, measles, and pneumonia have well-studied and established risk and prognostic factors, but there is more uncertainty surrounding these factors for HPV [17][18][19]. A thorough literature review was conducted to gather evidence on the risk and prognostic factors for HPV, with factors evaluated for use in the analysis based on the strength of evidence and the availability of data (disaggregated by wealth quintile and region) in the Ethiopian DHS [9]. Prognostic factors for HPVrelated disease outcomes, notably cervical cancer, are unavailable in disaggregate and we therefore did not analyze HPV-related deaths. Sections 2 and 3 in S1 Text provide detailed information on our methods and summarize the risk and prognostic factors selected for each VPD (see Table 1 for a brief summary).

Distributional modeling
Our approach substantially extends a previous methodology, which utilizes the population attributable fraction (PAF) and differences in the prevalence of risk and prognostic factors across wealth groups to estimate the distribution of VPD incidence and mortality across those wealth quintiles [16]. In this paper, we augment this methodology by using a linear system of equations to estimate the joint distribution of risk and prognostic factors in population subgroups (i.e. wealth quintile and geographic region). To compute the joint distribution of risk and prognostic factors, we treat each possible combination of risk or prognostic factors as computationally distinct, thereby allowing us to account for individuals with multiple risk factors.
For a given VPD, a simulated population (derived from DHS data) is divided according to all the possible risk factor combinations. The incidence and mortality of the VPD are estimated for each risk factor combination by minimizing a constrained linear system. The parameters of the linear systems are functions of the relative risks, the proportion of the risk factor combinations within the population, and the proportion of cases or deaths attributable to the disease. Optimization of VPD incidence for each risk factor combination was performed starting from a random admissible probability corresponding to either cases or attributable deaths. From this random starting point, we identified a nearby probability combination that provided a better fit to the parameters of the system. Then, from this second point, an improved nearby fit was identified, providing a third point that provided a better fit than the first two points. The process was repeated until a minimizer was identified (further details are provided in Section 2 of S1 Text).
After optimizing the multidimensional risk factor incidences, each individual in the simulated population was assigned a probability of morbidity (case) or mortality (death) according to their risk factor profile, giving rise to a probability distribution that can be further analyzed by wealth quintile, geographic region, or in totality. In order to estimate the impact of vaccination across quintiles and regions, the expected number of cases predicted by the probability distribution was reduced according to vaccination coverage rates and efficacies ( Table 2).
To estimate the distribution of mortality across quintiles and regions we employed three computational strategies. First, we applied estimates of quintile-and region-specific treatment-seeking behavior (e.g., care-seeking for diarrhea and lower respiratory infections, as reported in the DHS) by reducing the number of cases proportionally to estimated rates of Table 1. Selected risk and prognostic factors for diarrhea, measles, pneumonia, and HPV. See Section 3 of S1 Text for more detailed methodology on risk factor selection.

Results
The proportional risk distribution of VPD incidence and mortality across wealth quintiles and geographic regions was estimated for Ethiopia for: a "counterfactual" scenario, a baseline scenario, and a maximized vaccination scenario. Across all VPDs and scenarios, the risk gradients across wealth quintiles show an inverse relationship with wealth, such that the poor have a greater risk of becoming infected with or dying from a VPD (Fig 1). Likewise, by geographic region, regions with more rural and primarily pastoralist and agrarian populations have a greater risk of VPD than regions with richer and predominantly urban populations (Figs 2 and 3).
In the "counterfactual" scenario, VPD incidence and mortality is estimated based on the distribution of risk and prognostic factors and does not include health interventions (i.e., vaccination coverage and treatment-seeking behavior). Compared to the richest quintile, the poorest quintile has 1.78 (95% uncertainty range (UR): 1.65-1.93) times greater risk of rotavirus diarrhea incidence, 1.83 (1.52-2.23) times greater risk of measles incidence, 2.07 (1.97-2.16) times greater risk of pneumococcal and Hib pneumonia incidence, and 1.13 (1.09-1.17) times greater risk of HPV incidence (Table 3). Compared to the gradients observed by wealth quintile, the relative risk gradients by region are much larger. For example, Afar, which has a  The models for all three scenarios include the distribution of risk and prognostic factors. The "counterfactual" scenario presents a scenario without vaccination or treatment-seeking behavior; the baseline scenario presents a scenario with current vaccination coverage and treatment-seeking behavior; and the maximized predominantly rural and pastoralist population, has 3.79 (3.42-4.28) times greater risk of rotavirus diarrhea incidence, 4.50 (2.67-7.60) times greater risk of measles incidence, 5.99 (5.12-6.96) times greater risk of pneumococcal and Hib pneumonia incidence, and 1.48 (1.33-1.68) times greater risk of HPV incidence compared to Addis Ababa (the richest region with an entirely urban population). And the relative risk gradients, by both wealth quintile and region, are steeper when evaluating disease mortality. The relative VPD incidence and mortality gradients, in the baseline scenario in which vaccination coverage and treatment-seeking behavior are now included, indicate greater relative inequalities. Applying quintile-and region-specific vaccination coverage rates results in steeper relative incidence gradients (compared to the counterfactual scenario). This impact is compounded when examining relative mortality gradients, which also take into account quintile-and region-specific treatment-seeking behavior. Compared to the richest quintile, the poor have 3.34 (2.75-4.26) times the risk of mortality from rotavirus diarrhea, 4.23 (3.45-5.23) vaccination coverage scenario presents a scenario in which vaccination coverage in all population subgroups is increased to the level of vaccination coverage of the subgroup with the greatest vaccination coverage (e.g., all quintiles are modeled as having the same level of vaccination coverage as the richest quintile).
https://doi.org/10.1371/journal.pgph.0000819.g001 times the risk of mortality from measles, and 3.05 (2.84-3.32) times the risk of mortality from pneumococcal and Hib pneumonia. Regionally, compared to Addis Ababa, the population in Afar has 8.36 (5.83-13.38) times the risk of mortality from rotavirus diarrhea and 9.27 (7.87-11.07) times the risk of mortality from pneumococcal and Hib pneumonia.
In the maximized vaccination scenario, the relative risk gradients for VPD incidence return to the relative gradients observed in the counterfactual scenario, which demonstrates that increasing and equalizing vaccination coverage in Ethiopia would significantly improve health equity. As all quintiles and regions have the same vaccination coverage under this scenario, only the underlying inequalities in the distribution of risk and prognostic factors therefore contribute to the relative incidence gradients. Differences remain, however, for the relative mortality gradients, as treatment-seeking behavior was not modified in this scenario (i.e., treatment seeking was maintained similar across the baseline and maximized vaccination coverage scenarios).

Discussion
In this research we build and expand upon previous methods, which utilize the PAF and differences in the prevalence of risk and prognostic factors across population subgroups to estimate PLOS GLOBAL PUBLIC HEALTH Table 3. Estimated relative risk distribution of (A) rotavirus diarrhea incidence and mortality, (B) measles incidence and mortality, (C) pneumonia (pneumococcal and Hib) incidence and mortality, and (D) HPV incidence across three scenarios by wealth quintile and geographic region. The three scenarios include the distribution of risk and prognostic factors. The "counterfactual" scenario presents a scenario without vaccination or treatment-seeking behavior, the baseline scenario presents a scenario with current vaccination coverage and treatment-seeking behavior, and the maximized vaccination coverage scenario presents a scenario in which vaccination coverage in all population subgroups is increased to the level of vaccination coverage of the subgroup with the greatest vaccination coverage (e.g., all quintiles are modeled as having the same level of vaccination coverage as the richest quintile).

PLOS GLOBAL PUBLIC HEALTH
gradients are in the same general range as those found by Chang and colleagues [16], but the gradients are slightly greater. This makes sense at the individual level, since having multiple risk factors increases one's risk of morbidity and mortality compared to having a single risk factor, and at the population subgroup level, since the poor and those living in more rural regions are more likely to have multiple risk factors and often lack access to basic health care services. However, our modeling approach is limited by the way in which relative risks are measured. The standard methodology compares an individual's risk of morbidity or mortality with a given factor to a null scenario (i.e., no risk), whereas a more realistic evaluation would compare risk factors in various scenarios (e.g., comparing the risk related to having multiple factors with the risk related to having a single factor or a different set of factors). This model develops a more realistic depiction of how risk and prognostic factors can combine to produce morbidity and mortality without completely reworking the standard way in which relative risks are calculated. Additionally, this model does not include dynamic transmission of the infectious diseases examined within and across socioeconomic groups [35], and it therefore does not consider some meaningful factors that can contribute to disease incidence, such as household crowding, population density, and social contact among population subgroups. Lastly, quality of care can play an important role in mortality outcomes in addition to treatment-seeking behavior; however, data were not sufficiently available to include such a quality factor in our modeling.
Our results show that current unequal levels and distributions of vaccination and treatment-seeking patterns do not necessarily redress health inequalities across wealth quintiles and geographic regions. The poor and those in more rural regions are at a greater risk of incidence and mortality under the "counterfactual" scenario, because of the higher prevalence of risk and prognostic factors and lower health-seeking behavior among these population subgroups. In absolute terms, all populations subgroups will benefit from health interventions-in particular, the poorest will likely benefit the most in terms of absolute numbers of disease cases and deaths averted by vaccines. However, the relative distribution of health interventions across population subgroups may increase existing relative disparities (i.e., among the remaining cases and deaths not averted by vaccines). This is clearly not a reflection of the health interventions themselves but is rather due to unequal access to care (often pro-rich gradients) through vaccination coverage and treatment-seeking behavior.
This can be observed in the case of HPV: the poorest quintile is at an estimated 3.03 times greater risk of disease than the richest quintile in the baseline scenario (with vaccination and treatment-seeking behavior), compared to 1.13 times greater risk in the "counterfactual" scenario (no health interventions). For deaths, the greatest difference occurs for measles: the poorest quintile is at an estimated 4.20 times greater risk of death than the richest in the baseline scenario compared to 2.34 times greater risk in the "counterfactual" scenario. This is likely most apparent for these two VPDs, because of the corresponding vaccines' high efficacy (85% for MCV and 100% for HPV vaccine).
These findings highlight the urgent need to target both (1) the underlying risk factors (and corresponding social determinants) that contribute to VPD incidence and set the foundation for health inequalities in Ethiopia, and (2) low vaccination coverage and low treatment-seeking behavior, as lower rates among the poor and those in rural regions exacerbate relative inequalities. This analysis highlights that poor and rural children in Ethiopia stand to benefit the most from more equitable vaccination programs because they would have the most to gain. Even pro-rich vaccination coverage, in which the richest quintiles have higher coverage rates (as observed in Ethiopia), would likely result in greater health benefits in absolute terms (i.e., number of infections and deaths averted) given the more severe risk and prognostic profiles of the poor [9]. Ethiopia's vaccination program could become more pro-poor if coverage were more equal across wealth quintiles and the rural-urban divide, but might maintain relative health inequities if current pro-rich vaccine coverage rates remain.
The WHO's "Reach Every District" (RED) approach outlines key operational components specifically aimed at improving immunization coverage in every district [36]. Stronger implementation of this approach, which includes routine supervision of health facilities and monitoring for action, would likely improve disparities in Ethiopia's vaccination coverage. Indeed, in Ethiopia, DPT3 coverage has been found to be more than 4 times greater in zones (the administrative level below region) that conduct health facility supervision, provide written feedback to their health facilities, and have health facilities that monitor their immunization performance [37]. Beyond the health system, evidence in Ethiopia suggests that improving women's autonomy through education and gender-related initiatives can also lead to improved immunization coverage [38].
Our study presents a number of limitations. First is the somewhat outdated (from the DHS year 2016) and lacking evidence on risk and prognostic factors, especially for HPV, and the way in which these factors are estimated (i.e., compared to a null scenario, discussed above). Much has happened since the 2016 Ethiopian DHS (e.g., the COVID-19 pandemic, armed conflict); however, the 2016 Ethiopian DHS remains the most recent source of comprehensive disaggregated (by wealth group and geographic region) empirical data on risk and prognostic factors relevant to this analysis. Additionally, the epidemiological link between VPD cases and deaths is likely much more complicated than any of the three methods used by our modeling to estimate relative mortality risk gradients. The second relates to the model's optimization. The optimization was performed from a random starting point, and therefore there is no guarantee of stable results; and because the minimizer is not unique, there is no assurance that the scheme samples uniformly amongst the minimizers. Uniform sampling from high dimensional convex spaces can be computationally challenging, so we adopted a less sophisticated rearrangement method based on independent samples for each variable to generate starting points. Furthermore, the mathematical methods operate under the assumption that the proportional distribution across quintiles and regions are invariant under a scaling of the population from the survey sample. Rounding cutoffs and small perturbations of survey population were employed for computational reasons. Lastly and importantly, we lacked the empirical data (e.g., relative distributions in the incidence and mortality of those VPDs across socioeconomic groups in Ethiopia) to validate the model. Nevertheless, the methodology employed by this research takes simultaneously into account the joint nature of VPD risk and prognostic factors and uses the distribution of these factors to estimate relative gradients of VPD morbidity and mortality across population subgroups, wealth quintiles and geographic regions, in Ethiopia. While vaccination undoubtedly improves health outcomes in absolute terms, especially so among the poorest, the results of our modeling highlight the fact that unequal access to healthcare interventions can exacerbate pre-existing inequalities in the relative distribution of the risk and prognostic factors that give rise to VPD burden (among diseases cases and deaths not averted by vaccines). Current unequal levels and pro-rich gradients of vaccination and treatment-seeking patterns should be redressed so to significantly improve health equity across wealth quintiles and geographic regions in Ethiopia.