Cost-effectiveness of rotavirus vaccination in children under five years of age in 195 countries: A meta-regression analysis

Background Rotavirus caused an estimated 151,714 deaths from diarrhea among children under 5 in 2019. To reduce mortality, countries are considering adding rotavirus vaccination to their routine immunization program. Cost-effectiveness analyses (CEAs) to inform these decisions are not available in every setting, and where they are, results are sensitive to modeling assumptions, especially about vaccine efficacy. We used advances in meta-regression methods and estimates of vaccine efficacy by location to estimate incremental cost-effectiveness ratios (ICERs) for rotavirus vaccination in 195 countries. Methods Beginning with Tufts University CEA and Global Health CEA registries we used 515 ICERs from 68 articles published through 2017, extracted 938 additional one-way sensitivity analyses, and excluded 33 ICERs for a sample of 1,418. We used a five-stage, mixed-effects, Bayesian metaregression framework to predict ICERs, and logistic regression model to predict the probability that the vaccine was cost-saving. For both models, covariates were vaccine characteristics including efficacy, study methods, and country-specific rotavirus disability-adjusted life-years (DALYs) and gross domestic product (GDP) per capita. All results are reported in 2017 United States dollars. Results Vaccine efficacy, vaccine cost, GDP per capita and rotavirus DALYs were important drivers of variability in ICERs. Globally, the median ICER was $2,289 (95% uncertainty interval (UI): $147–$38,993) and ranged from $85 per DALY averted (95% UI: $13–$302) in Central African Republic to $70,599 per DALY averted (95% UI: $11,030–$263,858) in the United States. Among countries eligible for support from Gavi, The Vaccine Alliance, the mean ICER was $255 per DALY averted (95% UI: $39–$918), and among countries eligible for the PAHO revolving fund, the mean ICER was $2,464 per DALY averted (95% UI: $382–$3,118). Conclusion Our findings incorporate recent evidence that vaccine efficacy differs across locations, and support expansion of rotavirus vaccination programs, particularly in countries eligible for support from Gavi, The Vaccine Alliance.


Introduction
Mortality due to diarrheal diseases has declined by 60% since the year 2000, but diarrheal diseases remain one of the leading causes of death for children under 5 every year, and were responsible for an estimated 497,434 childhood deaths in 2019 [1]. Rotavirus is a leading etiology of diarrheal disease mortality globally, responsible for 30% percent of diarrheal deaths among children under 5 [2].
Some of the recent decline in mortality has been attributed to rotavirus vaccination, which averted an estimated 28,000 deaths in 2016 [2]. By the end of 2018, Gavi, The Vaccine Alliance, supported rotavirus vaccination programs in 45 of 57 eligible countries, covering 39% of children in those countries with a full course of a rotavirus vaccine [3]. By April 2020, a total of 107 countries introduced rotavirus vaccines to their routine immunization schedules [4].
As more countries consider introducing rotavirus vaccination into their routine immunization programs, cost effectiveness of the vaccine can inform their decisions. However, cost effectiveness analyses (CEAs) are not available in 53 countries. Among the 142 countries with existing CEA, 49 have estimates that are based on a single study, meaning that the results are sensitive to the modeling assumptions of the study authors and do not consider the cumulative evidence. For countries with more than one CEA, results can vary greatly. For example, the 20 published incremental cost-effectiveness ratios (ICERs) for Peru vary from $132 to $2,438 per DALY averted. Similarly, in Bangladesh, the 19 published ICERs vary from $23 to $1,543 per DALY averted. Because of this variability, efforts to synthesize CEA results and transfer them to other settings are growing. Two systematic reviews synthesized published evidence of CEA on rotavirus vaccines [5,6]. The first compared methodological approaches between high-income (HIC) and lowincome (LIC) countries [5]. The second summarized study characteristics and findings, identifying vaccine efficacy as an important source of variability [6].
More recent cost-effectiveness research relied on meta-analysis to synthesize evidence from studies conducted in multiple countries. Haider et al. used a meta-regression approach to estimate the cost-effectiveness of the rotavirus vaccine across 29 LIC and lower-middle income (LMIC) countries [7]. They produced pooled estimates for the LICs and LMICs under consideration, rather than providing country-specific estimates. Further, their analysis incorporated only three covariates (Gross Domestic Product [GDP], vaccine coverage, and literacy rate) [7]. Jit et al. investigated costeffectiveness across five European countries using a single model and incorporating country-specific rotavirus burden. They found that results varied across settings due to differences in rotavirus burden and vaccine price, among other factors [8]. Rosettie et al. developed and applied a meta-regression approach to CEA to the human papilloma virus (HPV) vaccine [9]. They exploited oneway sensitivity analyses extracted from published studies and used covariates for study methods, vaccine characteristics (coverage, cost, type, booster, and target group), and country-specific variables for GDP per capita and HPV burden to predict incremental cost effectiveness ratios (ICERs) across 195 countries.
In this paper we apply meta-regression methods developed by Rosettie et al. to estimate the cost-effectiveness of rotavirus vaccination in 195 countries. We identify and quantify sources of heterogeneity in published CEAs from study methods, and vaccine characteristics, well as economic and epidemiologic conditions in each country (i.e. GDP per capita, rotavirus burden). We focus specifically on vaccine efficacy, because it differs by location [1,10] and has been identified as a driver of heterogeneity [6], but has not yet been used in any meta-analyses of a vaccinepreventable disease.

Cost-effectiveness data
We used cost-effectiveness data through 2017 from two comprehensive registries of all published CEA maintained by Tufts University's Center for Evaluation of Risk and Value in Health: (1) CEA registry with cost per quality-adjusted life year (QALY) results [11], and (2) Global Health CEA registry with cost per disabilityadjusted life year (DALY) results [12]. As previously reported, they searched PubMed for English-language articles using keywords ''QALYs", ''quality-adjusted" [11] and ''cost-utility analysis" for the CEA registry, and ''disability-adjusted" or ''DALY" for the Global Health CEA registry [12]. Abstracts were screened to identify original estimates, and assigned an International Classification of Dis-eases, 10th Revision (ICD-10) code to each article, [11,12] which we subsequently mapped to 2017 Global Burden of Disease, Injury and Risk Factor Study (GBD) causes and etiologies. For rotavirus,  the CEA registry contains 30 studies measuring the cost per QALY,  while the Global Health CEA registry contains 41 studies measuring the cost per DALY (Supplementary Table S1). Combined, these two registries include 519 incremental cost-effectiveness ratios (ICERs) across 142 countries.
This study complies with the Guidelines for Accurate and Transparent Health Estimates Reporting (GATHER) (Supplementary Table S2).

Data standardization, mapping, and extraction
Building on the registry data, we performed several tasks to create covariates and extend the dataset, as described in the Supplementary Sections 3 and 4, and previously reported [9]. Briefly, we extracted data on 5 study methods variables with missing values in the registry data whenever possible: cost discount rate, QALY/DALY discount rate, health outcome measure, perspective, and time horizon. We collapsed perspectives into two categories: healthcare payer and health sector versus all other categories. To account for economic and epidemiological differences across countries, each ICER was mapped to GBD categories for at least one age group, sex, location, and cause, which were used to pull GDP per capita and rotavirus burden measured as DALYs from diarrhea attributable to rotavirus per 100,000 population for children under 5 for the study year [2]. DALYs combine the effects of mortality and morbidity where a death from rotavirus in a child aged 2 to 4 years is measured as 81.81 years of life lost, and an episode of rotavirus is measured as roughly 0.002 years lived with disability, i.e. 5/365 years duration times an average disability weight of 0.12. We extracted data on four vaccine characteristics: vaccine type (monovalent, pentavalent, or both), coverage, cost, and efficacy. All comparisons were relative to a null comparator defined as no intervention (n = 1,245, 92.6% of final sample), standard of care (n = 98, 7.3%) or placebo (n = 2, 0.1%). When registry entries did not correspond to these categories, we extracted data to recalculate them. Finally, we extracted ICERs and covariates from oneway sensitivity analyses, identifying the covariate that differed, and the reference analysis. ICER values, vaccine costs, and GDP per capita were adjusted to 2017 US dollars ($US).

Inclusion and exclusion criteria
Of the 519 initial cost-effectiveness results from the Tufts database, we excluded two articles that compared modeling methods (two ratios), or did not report vaccine type (two ratios) (Fig. 1). We extracted an additional 938 cost-effectiveness results from one-way sensitivity analyses for seven variables. We also excluded 35 ratios due to missing data, reporting errors, or because the ICER was 0 (Fig. 1).
We conducted a logistic regression analysis (see below) to estimate the probability that a cost-effectiveness result was costsaving using 1,418 cost-effectiveness results, of which 73 were cost-savings. We conducted the meta-regression analysis with a total of 1,345 ICERs (excluding the cost-savings results) from 142 countries and 68 studies.

Statistical analysis
To estimate ICERs across countries, we adopted a Bayesian mixed-effects meta-regression framework that employs regularization and outlier detection to select covariates and estimate their associations with the observed ICER. We used those estimated associations to predictions ICERs for 195 countries. This framework consists of five stages described in Supplementary Section 5 and previously reported [9,13].
The first stage estimated prior distributions for seven preselected covariates in which one-way sensitivity analyses were done in the published literature: cost discount rate, DALY/QALY discount rate, perspective, vaccine type, vaccine coverage, vaccine cost, and vaccine efficacy. The aim of the sensitivity analyses was to reduce omitted variable bias and identify prior distributions for variables that were highly correlated, such as log-vaccine cost and log-GDP per capita. The one-way sensitivity analyses differed by one covariate and no unmeasured covariates from their corresponding reference analysis. We therefore matched each sensitivity analysis with its corresponding reference analysis and fit separate linear models to estimate the effect of the difference between covariate values in reference and sensitivity analyses on the difference in the corresponding log-ICERs. The estimated regression coefficients and standard errors from these models were then used as Gaussian prior distributions for these covariates in all subsequent models.
The second stage estimated the association between log-ICER and log-GDP per capita, using a spline ensemble that can estimate a non-linear response curve. The estimates included log-rotavirus burden as a covariate. This stage had a robust statistical approach to outlier detection, with 10% of observed log-ICERs being classified by the model as outliers and trimmed from the dataset ahead of further modeling steps. The resulting spline-transformed GDP per capita variable was used in subsequent stages.
In the third stage, we used a generalized Lasso approach to select additional covariates to include in the final model, thereby balancing covariate selection from a priori knowledge versus statistical approaches. The seven covariates from the first stage, and the spline-transformed log-GDP per capita variables were preselected. Candidate covariates for the model were: whether or not the study used a lifetime time horizon, whether the outcome mea- sure was QALYs or DALYs, and log-rotavirus burden, all of which were selected for the final meta-regression model. In the fourth stage, we created Gaussian priors for the spline transformed GDP per capita variable, and the two variables selected in stage 3 to safeguard against overfitting in the final estimates. Given that we had coefficient estimates for these variables from Stage 3, we standardized the variables, and used a grid-search and 10-fold cross validation to select the standard deviation for Gaussian priors. The selected prior standard deviation for the standardized covariates minimized the mean squared error for predicting data in the hold-out set.
In the fifth stage, we fit a linear mixed model that accounted for between-study variability and the fact that studies often reported multiple ICERs as part of sensitivity analyses. The model has a random intercept for each article, and covariates selected in the third stage of the analysis using priors estimated from the first and fourth stages, respectively. The parameter estimates for the covariates were used to predict ICERs for 195 countries.
Using this modeling framework, we explore three different parameterizations of efficacy to determine whether or not its inclusion improves model fit. First, we include efficacy as a main effect. Second, we define effective coverage as the product of coverage and efficacy (scaled to be between 0 and 100%), where covered population is the primary beneficiary of the intervention. Finally, we estimate a model without efficacy. All model comparisons are based on R-squared and root mean square error (RMSE) fit statistics, with the best fitting model used for final results. The model with efficacy as the main effect had the best fit (Rsquare = 0.962, RMSE = 0.607), and was used for final results.
We estimated the probability that rotavirus vaccination was cost-saving using a mixed effects logistic regression as described in Supplementary Section 6. The model was fit subsequent to the meta-regression and included the same final set of covariates, but used the full data set with 1,418 results instead of the subset of the data selected in the trimming stage of the analysis. As with the meta-regression, we estimate a random intercept for each article.
We use the estimates from the meta-regression and logistic regression to predict both the ICER and probability the vaccine was cost-savings for 195 countries as a function of country-level covariates. The predicted ICER with adjustment is the product of the predicted ICER and (one minus the cost-saving probability). The predicted cost-saving probabilities were small, as were the adjustments to the predicted ICERs.
Our predictions assume 90% vaccine coverage (the GAVI target), lifetime time horizon, healthcare payer perspective, 3% burden and cost discount rates, monovalent vaccine type, and DALYs as the outcome measure. Our estimate of country-specific vaccine efficacy comes from GBD [1]. Because the rotavirus vaccine does not have a single market price, we used the cost of all required doses based on the 2017 cost per dose reported in the WHO's Market Information for Access (MIA4) and aggregated by Linksbridge [14,15] (see Supplementary Section 7). We used the listed UNICEF price for the 57 Gavi-eligible countries [16], and Pan American Health Organization (PAHO) for 22 countries eligible for PAHO Revolving Fund.
We predict both the mean ICER and 95% UI, which are then adjusted by the probability that the intervention is cost-saving. We report two ratios: 1) the adjusted mean predicted ICER to GDP per capita, and the upper bound of the UI of the adjusted predicted ICER to GDP per capita. While GDP per capita has been criticized as a cost-effectiveness threshold as outdated [17,18], and not an accurate measure of a population's preference for spending on health improvements [19], the ratios may provide a useful screening tool [19] and place the results in the context of the country's economy.
In meta-regression analysis, vaccine efficacy and log-DALYs per capita were negatively associated with the log-ICER, and logvaccine cost and log-GDP per capita were positively associated with it (Supplementary Table S5). A 10 % increase in burden would reduce the ICER by 4.2%, and a 10% increase in vaccine price would increase the ICER by 6.4%. Including vaccine efficacy in the analysis slightly improved model fit; R 2 increased from 0.960 to 0.961, while Root mean-squared error fell from 0.613 to 0.607. The random effects variance for the model including efficacy was also lower (0.576 including efficacy versus 0.593). We did not observe similar improvements for models of effective coverage (Supplementary Table S6). We report estimates for the model including efficacy in the main text, and estimates for the model without efficacy in Supplementary Table S9.
The ranges of observed ICERs within countries were broad due to differences in methods and intervention characteristics. The ranges of observed ICERs within super-regions were broad for the same reasons as well as differences in rotavirus burden and GDP among countries. Despite these differences, the final metaregression model fit the data well, as shown in the comparison of the observed and fitted ICERs globally and by super-region (Fig. 2).
Among the seven GBD super-regions (Table 2), Sub-Saharan Africa and South Asia had the lowest population-weighted mean ICERs among GBD super-regions. Among the 46 countries making up the sub-Saharan African region, the mean predicted ICER was $251 per DALY averted (90% UI 38-903), while for the 5 countries making up the South Asia region, the mean predicted ICER was $294 per DALY averted (95% UI 45-1 062). The 34 countries making up the high-income region had the highest ICER ($40 914 per DALY averted; 95% UI 6 382-153 130).
Predicted ICERs varied across countries (Table 3; Fig. 3A). Among all 195 countries, the lowest mean adjusted ICERs were observed in Central African Republic (2017 US$ 85 per DALY averted; 95% UI 13-302), and Chad ($120 per DALY averted; 95% UI 18-428), which have the two highest burdens worldwide. The highest mean predicted ICERs were observed in the US ($70,599 per DALY averted; 95% UI 11,030-263,858), and Luxembourg ($46,158 per DALY averted; 95% UI 7,256-169,808), which are in the bottom 8% of rotavirus burden worldwide. Predicted ICERs may be outside the range of observed ICERs, due to differences in parameters in the published and predicted results, such as parameters such as vaccine cost and efficacy.
Considering these results in the context of each country's economy, 21 countries had adjusted ICERs with upper UIs less than onehalf times GDP per capita (Fig. 3c). Among them, 15 were in the sub-Saharan Africa region, 4 in Southeast Asia, East Asia, and Oceana, one in South Asia (India), and one in North Africa and the Middle East (Sudan). An additional 43 countries had upper UIs that were between 0.5 and one times the GDP per capita, including 17 in sub-Saharan Africa, 13 in Latin America and the Caribbean, six in Southeast Asia, East Asia, and Oceana, three in North Africa and the Middle East, two in South Asia, and one each in high-income (Norway) and Central Europe, Eastern Europe, and Central Asia (Tajikistan).

Discussion
To our knowledge, this is the most comprehensive metaregression analysis of CEA of rotavirus vaccination conducted, building on previous work that focused on 29 LICs and LIMCs [7], and previous efforts to leverage a meta-regression approach to synthesize all available evidence on HPV and provide costeffectiveness estimates for 195 countries worldwide [9]. An asterisk (*) denotes that the total number of articles may exceed 68 since some articles examined multiple regions, vaccine characteristics, and cost-effectiveness analyses characteristics. QALY = quality-adjusted life-year, DALY = disability-adjusted life-year, IQR = interquartile range.      Unlike previous approaches, our approach facilitates the transferability of published estimates across settings, while accounting for differences between settings. Currently, a Ministry of Health in a country where there are no cost effectiveness estimates available may look to estimates from neighboring countries. However, these estimates may vary considerably within a country, and there would be no quantitative way to adjust them to account for differences in economic activity, burden of disease, vaccine efficacy or cost. For example, there are no published cost-effectiveness estimates in Equatorial Guinea, while neighboring Cameroon has estimates that range from $10 to $87 per DALY averted (Fig. 2). Relying on estimates from Cameroon would likely underestimate the costeffectiveness of the intervention because of differences in GDP per capita ($15,803 in Equatorial Guinea vs $1,671 in Cameroon) and corresponding vaccine costs ($22.53 vs $4.61, respectively). This is consistent with the findings from Jit et al, who noted that different vaccine prices contributed to variability in CEA results [8]. Using a meta-regression approach and accounting for the drivers of variability in cost-effectiveness, our estimates are $1,714 per DALY averted (95% UI: $267 -$6,306) in Equatorial Guinea.
Our approach also facilitates decision-making for policy-makers in a country equipped with multiple estimates. For example, Bangladesh's 19 published estimates ranged from $23 to $1,543 per DALY averted owing to differences in modeling assumptions that cannot be easily synthesized in a policy-making setting. Our approach quantifies the effects of important drivers of variability in ICERs, producing an estimated ICER and simple linear equation that allows policy-makers to leverage the estimated effects and   Our results provide encouraging evidence in favor of introducing rotavirus vaccination worldwide. All countries had predicted mean ICERs below three times GDP per capita, with the exception of the United States and Somalia. Taking uncertainty into account, the upper bound of the UI of the predicted ICER was below one times GDP per capita in 64 countries, but exceeded three time GDP per capita in five countries (the US, UK, Portugal, Greece, and Somalia).
Nevertheless, much progress needs to be made in order to meet the Global Vaccine Action Plan's target of 90% coverage by 2030. Of the 195 countries considered here, only 24 had reached 90% coverage in 2017 [20]. Twelve countries (Zambia, Eritrea, Senegal, The Gambia, Uzbekistan, Sao Tome and Principe, Ghana, Burkina Faso, Mozambique, Burundi, Rwanda, and Nicaragua) are LICs, two (Palestinian Territories and Morocco) are LMICs, five (Armenia, Fiji, Guyana, Mauritius, Paraguay) are UMICs, and five (Luxembourg, Norway, Saudi Arabia, Bahrain, and Qatar) are HICs. These countries all have different economies but a shared commitment to rotavirus vaccination, suggesting that the large scale-up and delivery of vaccinations across the rest of the world is possible.
Our work here has a number of limitations. First, we are unable to account for all differences between the studies used in our metaanalysis. Doing so would require more detailed descriptions of methodologies used in each article. Second, although we used the full sample of results to estimate the probability that the result was cost-saving, and adjusted all predicted ICERs by the predicted cost-saving probability, we did not propagate the correlation between the logistic regression and the meta-regression models. The probability that the rotavirus vaccine was cost-saving at current vaccine prices was low however, meaning that the correlation would not have substantially changed our results. Third, we applied Lasso in Stage 3 of the analysis framework to identify essential features from a subset instead of all covariates. Future research may improve the estimates by extracting sensitivity analyses for a more covariates, applying selection criteria to the variables in Stage 1, or improving methods for variable selection in the presence of collinearity in Stage 3. Fourth, we did not assess validity of our predictions in terms of model calibration and discrimination. CEA researchers working with decision analytic models have initiated research on predictive validity, and this field will have implications for the quality of the published results that are the data for in our meta-regression analysis. Fifth, uncertainty of reported ICERs is largely based on sensitivity analyses of different assumed values for input parameters.
Despite these limitations, our work is the most comprehensive set of country-specific estimates of rotavirus vaccination. In producing these estimates, our meta-regression approach synthesizes all available evidence, quantifies uncertainty to facilitate decisionmaking, and incorporates country-specific estimates of rotavirus vaccine efficacy. Our results support introducing or expanding rotavirus vaccination, which remains a critical step in reducing rotavirus burden.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.