Mapping subnational HIV mortality in six Latin American countries with incomplete vital registration systems

Background Human immunodeficiency virus (HIV) remains a public health priority in Latin America. While the burden of HIV is historically concentrated in urban areas and high-risk groups, subnational estimates that cover multiple countries and years are missing. This paucity is partially due to incomplete vital registration (VR) systems and statistical challenges related to estimating mortality rates in areas with low numbers of HIV deaths. In this analysis, we address this gap and provide novel estimates of the HIV mortality rate and the number of HIV deaths by age group, sex, and municipality in Brazil, Colombia, Costa Rica, Ecuador, Guatemala, and Mexico. Methods We performed an ecological study using VR data ranging from 2000 to 2017, dependent on individual country data availability. We modeled HIV mortality using a Bayesian spatially explicit mixed-effects regression model that incorporates prior information on VR completeness. We calibrated our results to the Global Burden of Disease Study 2017. Results All countries displayed over a 40-fold difference in HIV mortality between municipalities with the highest and lowest age-standardized HIV mortality rate in the last year of study for men, and over a 20-fold difference for women. Despite decreases in national HIV mortality in all countries—apart from Ecuador—across the period of study, we found broad variation in relative changes in HIV mortality at the municipality level and increasing relative inequality over time in all countries. In all six countries included in this analysis, 50% or more HIV deaths were concentrated in fewer than 10% of municipalities in the latest year of study. In addition, national age patterns reflected shifts in mortality to older age groups—the median age group among decedents ranged from 30 to 45 years of age at the municipality level in Brazil, Colombia, and Mexico in 2017. Conclusions Our subnational estimates of HIV mortality revealed significant spatial variation and diverging local trends in HIV mortality over time and by age. This analysis provides a framework for incorporating data and uncertainty from incomplete VR systems and can help guide more geographically precise public health intervention to support HIV-related care and reduce HIV-related deaths. Supplementary Information The online version contains supplementary material available at 10.1186/s12916-020-01876-4.


Background
Human immunodeficiency virus (HIV) continues to be a large contributor to morbidity and mortality across the globe [1,2]. While the burden of the HIV epidemic is most concentrated in sub-Saharan Africa, HIV remains a public health priority in Latin America. In 2017, the Global Burden of Disease Study (GBD) estimated over 30,000 deaths from HIV/AIDS-related causes in the region [1,2]. To combat the epidemic, the Joint United Nations Program on HIV/AIDS (UNAIDS) Fast-Track strategy emphasizes the need to reduce HIV-related deaths to less than 500,000 deaths worldwide by 2020-a 75% reduction from deaths in 2010 [3]. UNAIDS further targets a 90% reduction in HIV-related deaths by 2030 [3]. Despite increased access to antiretroviral therapy in many countries in Latin America [3,4], few countries show substantial reduction in HIV mortality since 2000 [1,2]. Continued efforts are needed to track progress towards meeting the UNAIDS Fast-Track goals with respect to HIV mortality.
Country-level estimates of HIV mortality in Latin America are available from a variety of sources [1,5] and estimates of mortality exist at the state level in select countries such as Brazil and Mexico [1]. Beyond this, however, detailed subnational estimates at the second administrative level in many countries in Latin America are absent. This lack of subnational estimates is alarming given the localized burden of HIV in urban areas and among high-risk subgroups such as people who inject drugs, sex workers, and men who have sex with men (MSM) [6][7][8][9][10]. Additionally, inequalities in HIV burden across local geographies likely occur given that many underlying drivers of HIV infection and death-such as poverty, incarceration, undernutrition, distribution of health practitioners, and access to health services-vary across geographic areas and through time [11][12][13][14]. Previous studies have confirmed substantial within-country variation in mortality rates but are limited to select countries, areas, and years [15][16][17][18].
The paucity of evidence on subnational HIV mortality is likely due to several methodological challenges associated with a granular spatial modeling of HIV mortality in Latin American countries. Deaths attributable to HIV are inherently small in number in areas with small populations, adding stochastic noise to direct estimates [19]. Past approaches have used Bayesian models that apply small area methods that borrow strength across age, time, and space to produce stable estimates of mortality rates in areas with a small number of deaths [15,[20][21][22]]. An additional complication in HIV mortality estimation is that-due to stigma or misdiagnosis-HIV deaths may be misclassified and coded to other underlying causes of death such as tuberculosis, endocrine disorders, meningitis, or encephalitis [23][24][25]. Moreover, vital registration (VR) systems in many countries in Latin America are incomplete and not all deaths are recorded in official statistics [23,[25][26][27]. While a variety of methods have been proposed to estimate the completeness of death registration at the country level [28,29], standard methods rely on stable population and sex pattern assumptions that often do not hold in small subnational areas [29,30]. To resolve the difficulties associated with small numbers of deaths and VR completeness, Schmertmann and Gonzaga proposed a Bayesian model framework for small area life expectancy estimation in countries with incomplete VR systems [31]. This method incorporates a novel functional form for mortality that is informed by prior distributions for VR completeness coverage based on empirical evidence.
In this analysis, we address these challenges by utilizing comprehensive cause of death assignment and applying a small area estimation framework that incorporates prior information on VR completeness to produce estimates of HIV mortality and deaths due to HIV by age and sex at the municipality level in six countries in Latin America: Brazil, Mexico, Guatemala, Costa Rica, Colombia, and Ecuador. Our modeling approach leverages information from neighboring areas across space and time to produce estimates across all years of available data. These six countries were selected based on availability of VR data, but not all contained the same range of available years:

Overview
This analysis complied with the Guidelines for Accurate and Transparent Health Estimates Reporting (GATHER) [32]. Our ecological study estimated the HIV mortality rate and the number of HIV deaths by age group and sex at the municipality level for all years of available VR (Additional file 1: Figure S1). All analyses were carried out at the second administrative unit level, which we refer to as municipalities for convenience unless referencing country-specific results, where we use the appropriate national nomenclature for the administrative subdivision (municipality for Brazil, Colombia, Guatemala, Mexico; canton for Costa Rica and Ecuador). Municipalities were combined as needed to create stable units of analysis over the study period, reducing the total number of areas analyzed in select countries (Table 1; Additional file 1: Table S1). In the results, all presented rates are age-standardized for comparison between countries, unless otherwise stated. We used standard age weights produced by GBD 2017 for age standardization [1].

Vital registration data
Vital registration (VR) mortality data consisted of anonymized individual-level records from all deaths reported in each country's VR system occurring between the years of study (Additional file 1: Table S2). These records were tabulated by municipality of residence, age group (0-4, 5-9, 10-14, …, 75-79, and ≥ 80 years), sex, and the underlying cause of death according to the tenth revision of the International Classification of Diseases (ICD-10) [33]. We standardized VR data using methods developed for the GBD [1]. This process requires all deaths to be attributed to a single underlying cause of death following ICD guidelines and fits within a hierarchy of mutually exclusive and collectively exhaustive causes. Deaths that were coded with ICD-10 codes that could not be an underlying cause of death, as well as deaths that were coded to non-specific causes of death, were redistributed to most detailed causes of death by age, sex, municipality, and year according to a framework developed by Naghavi et al. [34] and updated for GBD 2017 [1]. This includes an HIV correction step that rectifies deaths assigned to comorbidities such as tuberculosis, endocrine disorders, meningitis, or encephalitis that diverge from locations without HIV epidemics [1].

Location hierarchy
We created country-specific location hierarchies that list all subnational administrative units for each year in the specified time period and match each corresponding death from the VR system to the municipality level (Additional file 1: Figure S2). For each country, municipalities were geo-matched to shapefiles provided by the Global Administrative Unit Layers (GAUL) [35] (Brazil, Costa Rica, Guatemala) or the Humanitarian Data Exchange [36] (Colombia, Ecuador, Mexico). In all selected countries, municipality boundaries changed over time, reflecting new boundary designations across the years of study ( Table 1). Municipalities that underwent a boundary change during the period of the analysis were merged to create a stable unit across the period of observation, and the municipality-level shapefiles were manually edited to match the split hierarchy using ArcMap version 10.6 [37]. Merged units that included multiple municipalities were modeled as one area, and in the results share the same estimates of HIV mortality rate. Details of these shifts are provided in Additional file 1: Table S1.

Covariates and population
We included several available covariates to help inform estimates of HIV mortality: population density [38], night-time light brightness [39], urbanicity [40], and travel time to the nearest settlement of more than 50,000 inhabitants [41] (Additional file 1: Table  S3). These covariates were selected because they are factors or proxies for factors previously identified in the literature as associated (not necessarily causally) with HIV mortality. Specifically, these four variables were included as measures or proxies for connectedness and urbanicity as HIV historically spread among high-risk groups in urban areas [6,42] and is typically found to be higher in more urban compared to more rural locations. Each covariate was obtained in a raster format at a 5 × 5-km resolution and required aggregation to the modeled municipality level for inclusion in our modeling framework. This aggregation was done fractionally: raster cells that crossed municipality borders were fractionally allocated to municipalities in proportion to the covered area.
We created age-and sex-specific populations for each municipality unit by aggregating the WorldPop [38] raster to the modified shapefile, utilizing the same fractional aggregation process. The age-and sex-specific populations for municipalities were then scaled to the national population estimates derived from the GBD [1]. To do so, for each country, sex, age group, and year, we defined a population raking factor as the ratio of the GBD population estimate for that same sex, age group, and year to the sum of the WorldPop population for all municipalities within the country, and then multiplied the WorldPop population estimates for each municipality within the country by this raking factor. This resulted in age-and sex-specific population estimates for each municipality which aligned with the GBD national population sizes and structures.

Statistical model Vital registration completeness
Expanding on previous literature [31], we used a Bayesian framework that bypasses a lack of identifiability between the mortality rate and completeness estimate by incorporating an informed prior on VR completeness. In this analysis, we incorporated information from GBD [1] on subnational (for Brazil and Mexico) and national VR completeness (for remaining countries) as well as geographic patterns in under-5 VR completeness from past analyses [43] to generate priors on municipality-level VR coverage by two age groups (< 15 years and 15+ years) and year (Additional file 1: Figure S3). We selected these two age groups based on the available national VR completeness estimated in GBD and established literature and expert opinion [31,44]. The supplemental methods outlined in Additional file 1 summarize our process for generating informed priors on VR completeness in greater detail. We did not model VR completeness for adults if national GBD completeness estimates exceeded 95% in all years of available VR data (Costa Rica and Colombia). Similarly, we did not model under-15 VR completeness if GBD estimates of child completeness were greater than 90% in all years of available VR data (Costa Rica, Guatemala, Mexico). We therefore modeled adult completeness in Ecuador, Guatemala, Mexico, and Brazil, and child completeness in Ecuador, Colombia, and Brazil.

Modeling framework
We estimated HIV mortality separately by sex using a small area estimation framework built upon a model developed in prior modeling studies [15,45]. This Bayesian hierarchical generalized linear model used a Poisson data likelihood to model the number of HIV deaths in a municipality, year, and age group (supplemental methods in Additional file 1). The Poisson distribution was characterized by a parameter that multiplied the mortality rate, and population by municipality and VR completeness by municipality (Colombia, Ecuador, and Guatemala) or state (Brazil and Mexico). Completeness priors added probabilistic information about VR coverage that allowed estimation of mortality rates given counts of registered deaths. We modeled the log of the mortality rate as a linear combination of terms including random effects with conditional autoregressive distributions to smooth over age, year, and municipality. We also included covariates as fixed effects (see supplemental methods in Additional file 1).
Models were estimated separately for each country and sex and fit using the TMB package [46]. One thousand draws were sampled from the approximated posterior distributions of each modeled parameter and used to construct 1000 draws of HIV mortality (m j, t, a ) for each municipality j, year t, and age group a. We calculated point estimates from the mean of these draws, and the lower and upper bounds of the 95% uncertainty interval from the 2.5th and 97.5th percentiles, respectively, for each age, sex, year, and municipality. Municipality-level estimates for each age, year, and sex were aggregated to the state and national level using a population-weighted average. In Brazil and Mexico, estimates were calibrated to GBD at the state level and estimates were calibrated to national HIV mortality estimates for the remaining four countries. To accomplish this, we calculated the ratio of the national-or state-level estimate from GBD to the mean national estimate derived from population-weighting m j, t, a , and multiplied all draws of m j, t, a by this ratio. We generated the number of HIV deaths for each age-sex-year-municipality by multiplying the mean, lower, and upper bounds of our mortality estimates by the corresponding WorldPop population estimate. We quantified the relative inequality as the mortality rate ratios for municipalities in the 90th percentile versus those in the 10th percentile of mortality rate by year. We calculated the absolute inequality as the difference in mortality between municipalities within a country in the 90th percentile and those in the 10th percentile in terms of mortality rate by year. Throughout our analysis, we qualify statements as statistically significant if the posterior probability of that statement exceeds 95%. We completed our analysis using R version 3.6.3 [47].

Model assessment
To assess if including VR completeness in our statistical framework improved model estimates, for the five countries where we applied completeness priors to either children under-15 or adults we also fit a model where the completeness term in the statistical model, π k;t;a Ã , was removed. We then compared the ratio of annual national HIV mortality in children under 15 and adults from GBD to 1000 draws of national estimates of annual HIV mortality from the standard and completeness models. This ratio, known as the raking factor, is plotted in Additional file 1: Figure S4-S9. A raking factor closer to 1 indicates better alignment with GBD, and inclusion of completeness priors generally resulted in closer alignment with national GBD mortality estimates.

Concentrated deaths due to HIV
We estimated that a large proportion of HIV deaths were concentrated in a small number of geographical areas with large populations. In all countries, over half of the HIV deaths were located in less than 10% of municipalities in the latest year of study (Fig. 7)

Absolute and relative inequality over time
Relative inequality-the mortality rate ratio for municipalities in the 90th percentile versus those in the 10th percentile-varied from 5.0 (4.9-5.2) and 4.8 (4.6-5.0) among men and women in Brazil in 2017, to 63.4 (31.6-113.0) and 167.6 (54.6-403.5) among men and women in Costa Rica in 2016 (Fig. 8). The estimated relative geographic inequality in HIV mortality increased in all countries from the first to the last year of study, and this increase was statistically significant in all countries barring Guatemala and Costa Rica. The largest percent increase in relative inequality for each sex over the study period was 49.0% in Colombian men (from 6.5 Absolute inequality-the difference between the mortality rate in the 90th percentile and the 10th percentile-showed less temporal variation in Brazil, Mexico, Fig. 8 Relative and absolute inequality among municipalities in HIV mortality. a Relative inequality, defined as the ratio of estimated HIV mortality for municipalities in the 90th percentile versus 10th percentile, by year with 95% uncertainty intervals. Costa Rica is omitted from this panel because its estimated relative inequality was > 50. b Absolute inequality, defined as the difference in HIV mortality rates for municipalities in the 90th versus 10th percentile, by year with 95% uncertainty intervals. Selected countries are differentiated by color and line type and Costa Rica: the difference in estimated absolute inequality between the first and last year of study was less than 1.5 deaths per 100,000 among both men and women. Male absolute inequality increased in Colombia by 40.9% (from 6.9 [6.4-7.4]  Local disparities in median age group among those who died from HIV The estimated median age group among men who died varied substantially at the municipality level in the latest year of study: by 15 years in Brazil, Ecuador, and Mexico, and by 10 years in Colombia, Costa Rica, and Guatemala (Fig. 9). Among women, estimated median age group among those who died in the latest year of study varied at the municipality level by 15 years in Brazil, Guatemala, and Mexico, by 10 years in Colombia, and by only 5 years in Costa Rica and Ecuador. Differences in median age group among those who died also shifted over time. In Brazil, the estimated median age group among male decedents rose in 99.6% of municipalities from 2000 to 2017, while in Guatemala only 21.5% of municipalities saw an increase in estimated median age group among male decedents from 2009 to 2017. An increase in estimated median age was also observed among women: in Mexico, Ecuador, Colombia, and Brazil, the median age group among female decedents rose in > 97% of all municipalities in each country. Additional file 1: Figure S16-S21 show estimated HIV mortality by age group and sex for each country in the last year of study.

Discussion
There are few past analyses that assess HIV mortality in Latin America using VR data at a subnational level, largely due to the statistical challenges of incorporating incomplete VR systems and estimating mortality in areas with small populations and small numbers of deaths. The few analyses that integrate estimates of VR completeness into their modeling framework are often done at the national or state level or are limited to a single country and year [15][16][17][18]. In this analysis, we expand on previously described methods [31] that include prior estimates of VR completeness and demonstrate the utility of estimates that combine uncertainty from both incomplete registration systems and from statistical methods designed to leverage information across space, time, and age to inform mortality rates in areas with small numbers of HIV deaths.
Our estimates revealed large-scale spatial heterogeneity in HIV mortality across the six Latin American countries considered in our analysis. We also reveal divergent national trends in HIV mortality in the six countries across the study period, and variable relative change within countries at the municipality level. From the first to the last year of study, HIV mortality decreased among men and women in all countries, with the exception of women in Colombia from 2000 to 2017 and both men and women in Ecuador from 2004 to 2014. Despite the progress in reducing HIV mortality among both sexes at the national level in Brazil, Guatemala, Costa Rica, and Mexico, inequalities in municipality-level HIV mortality persist and relative inequality increased over time in all countries. This analysis highlights uneven progress towards reducing HIV mortality and reaching UNAIDS Fast-Track goals, and emphasizes an alarming trend in Ecuador, where over 95% of cantons experienced increases in estimated mortality among both sexes from 2004 to 2014. Nonetheless, it also underlines stories of success: all countries contained municipalities with an estimated decrease in HIV mortality. Further evaluating municipalities with the greatest decreases in HIV mortality within a country may help decision-makers recognize successful strategies that could be implemented in municipalities experiencing increases or slower declines.
There are likely a multitude of factors that contribute to the spatial and temporal patterns in HIV mortality observed in our analysis. Consistent with past analyses, we found higher rates of HIV mortality among men compared to women and slightly different spatial patterns by sex [15,48,49]. The consistently elevated levels of mortality among men is likely partially because men who have sex with men (MSM) continue to be one of the populations with the highest prevalence throughout Latin America, and a group that suffers a higher level of stigma and discrimination [48,[50][51][52]. These trends may also reflect prevalent gender norms [49], unequal access to timely diagnosis and treatment [48], or differences in disease burden from comorbidities between genders [1]. The spatial distribution of high-risk groups possibly also contributes to the HIV epidemic remaining concentrated in large urban centres [6]. For example, one important driver of spatial differences in HIV mortality is prison populations, where HIV transmission is high due to overcrowding, violence, and lack of information on the risk of HIV acquisition [53]. In Brazil, there is evidence of low adherence to antiretroviral therapy (ART) and a higher proportion of primary and secondary resistance among prison populations, which are predominantly male [54]. Municipalities with large prison populations-such as several in the state of Sao Paulo, Brazil-show higher levels of HIV mortality rates and number of deaths due to HIV. Another potential driver of spatial heterogeneity is population migration. Political conflicts and economic hardships across the region, notably in Venezuela and Central America, have fostered waves of migration that can affect HIV prevention, treatment, and care program [55]. Furthermore, difficulties in acquiring HIV treatment and ART shortages spurs regional migration that can differentially impact HIV care and control programs in bordering countries [55,56].
A key driver of temporal trends in HIV mortality is the implementation of ART treatment programs, which have been incorporated to varying degrees in all Latin American countries and generally led to increases in ART coverage [55]. ART treatment is often a central focus of national HIV programs, but differences in priority and ability to commit resources have likely impacted progress in reducing HIV mortality. Access to HIV treatment for people living with HIV is country-dependent and has shifted over time. Brazil was the first middleincome country to offer free ART treatment to people living with HIV (PLHIV) in 1996 [57], with Costa Rica following soon after in 1998. Within the ensuing decade, Mexico [58] and Colombia [59] adopted similar policies of universal treatment, which may have contributed to the observed reduction in national HIV mortality as these programs matured. In recent years, Guatemala has also expanded ART treatment options through joint support from the national government and The Global Fund to Fight AIDS, Tuberculosis, and Malaria [60]. Given the concentration of HIV in urban high-risk groups, many national programs-such as those in Costa Rica and Guatemala-have emphasized a combination prevention strategy that focus on HIV testing, STI diagnosis, and linkage to care in vulnerable populations [61]. In Ecuador, the only country in our analysis where we estimated increases in national HIV mortality among men and women over the study period, there has historically been a paucity of information available on research findings concerning HIV/AIDS burden [62]. Notably, our analysis in Ecuador only extends to 2014 and patterns of HIV mortality may have changed with the recent emphasis on community testing and treatment in metropolitan areas [63]. While all countries have expanded access to treatment and documented increases in ART coverage, albeit at different time periods, persistent disparities in access to quality health services and adherence to ART remain [6,64,65] and may contribute to differences in HIV mortality declines observed in this analysis. Further, in all countries in our analysis, communities with socioeconomic and health disadvantages-such as indigenous communities, sex workers, and transgender populations-often have unequal access to treatment and are an emerging or establish public health priority [66,67].
Our analysis also revealed substantial variation in median age group among those who died from HIV at the municipality level, and an increase over time in median age group among those who died of HIV. These results agree with past research that demonstrates an accelerating growth in the number of people living with HIV that are above 50 years of age [68,69]. Several countries, notably Brazil and Mexico, contained municipalities with a 15-year difference in median age groups among men and women who died from HIV in the latest year of study. While attributing these trends to specific drivers is outside the scope of this analysis, there are several factors that could influence an increase in median age of death, including changes in HIV incidence in specific age groups [70], increases in life expectancy among people living with HIV [70], access to ART [71], migration, and the age distribution of the population across the study period.
This analysis provides novel subnational estimates of HIV mortality that convey important information to policymakers and could inform future action. Knowledge of local differences in HIV mortality can help guide scale-up of ART where mortality might reflect suboptimal coverage. Our estimates highlight how deaths due to HIV are concentrated in a low proportion of municipalities. In the longer term, HIV mortality measures could be used to highlight areas that might benefit from programmatic interventions that target HIV prevention such as pre-exposure prophylaxis (PrEP). As countries in Latin America consider expanding access to PrEP, studies have demonstrated that prioritization of PrEP to those at highest risk could save money and lives [72,73]. Cost-effective interventions are especially important in Guatemala, Ecuador, Colombia, and Costa Rica, where HIV programs depend on donor funding [74]. Furthermore, subnational differences in HIV burden have already been used to develop localized strategies for HIV prevention and elimination in sub-Saharan Africa [75,76].

Limitations
This analysis is subject to a number of limitations. First, the VR data that inform our estimates are subject to misclassification biases. HIV is generally under-reported as a cause of death [23]. While the HIV correction methodology from the GBD employed in this analysis corrects for biases in HIV deaths classified to other underlying causes of death, there may be additional country-specific biases not addressed by our methodological approach. While we matched decedents to their municipality of residence as provided by the VR mortality databases, these data could overrepresent urban areas with larger health care facilities where individuals may have died after seeking treatment. It is possible that in some cases these individuals may have been mistakenly recorded as residing in that municipality if their municipality of residence was not available. Second, our method for correcting for incomplete VR across space and time makes several crucial assumptions. Our analysis incorporates estimates of VR completeness for children under 15 and adults 15 and over based on previous analyses, but VR completeness may vary within these age ranges. Further, we use the geographic variation in completeness identified by comparing reported under-5 all-cause deaths to previous estimates, and the variation in under-5 mortality may not be comparable to patterns in adult VR incompleteness. Additionally, variation in all-cause completeness may differ from patterns in HIVspecific VR completeness. Third, we use population estimates from WorldPop in this analysis that are subject to error, especially in sparsely populated areas. While WorldPop estimates include census data as inputs [77], depending on timing and data accessibility, estimates may differ from the underlying census measures and may not utilize the most recent census or the most detailed tabulations. Fourth, population migration in response to political conflict or economic instability in the region, including Central America and Venezuela, may not be properly captured in WorldPop estimates or recorded in vital registration systems. Fifth, our analysis is subject to large uncertainty. This reflects uncertainty both due to the small number of HIV deaths at the municipality level and the need to estimate completeness. While we believe this method better captures major sources of uncertainty, care must be taken when interpreting results. Sixth, we use custom shapefiles that are matched to country-level administrative subdivisions, and differences in administrative divisions between GAUL [35] or the Humanitarian Data Exchange [36] and an individual country's designation of administrative areas may affect the accuracy of results, especially in our estimates of the number of HIV deaths by municipality. Seventh, our small area estimation models smooth over space and time by making assumptions about the temporal and spatial structure of HIV mortality that may not always hold. Eighth, VR data availability varied across the countries selected in our analysis, and comparison between temporal trends in HIV mortality may be difficult to assess for countries with different years of data availability. Finally, it is difficult to directly assess violations of our modeling assumptions or quality issues in the underlying data sources given that VR completeness cannot be verified. Nonetheless, comparisons to GBD national estimates (Additional file 1: Figure S4-S9) provide reassurance in overall country trends.

Future directions
There is considerable opportunity to expand this analysis. First, access to VR data over more years of study, or in neighboring countries in Latin America, could provide valuable benchmarks for more direct comparisons and allow additional information across space and time to potentially improve our models. Our current study uses four available covariates that serve as proxies for urbanization and development, but in the future, availability of other drivers of HIV mortality at the municipality level such as socioeconomic status, healthcare infrastructure, high-risk group concentration, and ART treatment availability could improve our estimates. Further, the technique we used to include uncertainty and information on subnational VR completeness could be extended to other countries where VR systems are not complete. Finally, this small area estimation framework could be used to estimate all-cause and cause-specific mortality due to other causes at local levels in the six modeled Latin American countries.

Conclusion
Our analysis finds large-scale variation in HIV mortality among municipalities in six Latin American countries, both in the latest year of study as well as over the entire study period. Our estimates demonstrate the need to assess HIV burden at a granular geographic scale in Latin America, given that the HIV epidemic is concentrated in high-risk groups and select urban areas. The methods developed in this analysis provide a framework for incorporating prior information on VR completeness into subnational estimates of HIV burden. This analysis could be used to identify areas that have successfully reduced HIV mortality and areas of high HIV burden, as well as to inform the rollout of preventive interventions that are required to help countries progress towards achieving UNAIDS targets and advance health equity.
Additional file 1:. Supplemental methods, GATHER checklist, Supplemental Figure S1-S21, and Supplemental Tables S1-S4. Figure S1. Analytical process overview. Figure S2. Analytical process for VR data.  Figure S16. Estimated HIV mortality in Brazil by age group, 2017. Figure S17. Estimated HIV mortality in Colombia by age group, 2017. Figure S18. Estimated HIV mortality in Costa Rica by age group, 2016. Figure S19. Estimated HIV mortality in Ecuador by age group, 2014. Figure S20. Estimated HIV mortality in Guatemala by age group, 2017. Figure S21. Estimated HIV mortality in Mexico by age group, 2017. Table S1: Merged municipalities by country to form stable geographical units. Table S2: Vital Registration data. Table  S3: Covariate data sources. Table S4: National HIV mortality rates among men and womenfwil.

Funding
This work was primarily supported by grant OPP1132415 from the Bill & Melinda Gates Foundation. The funder of the study had no role in study design, data collection, data analysis, data interpretation, writing of the report, or decision to publish. The corresponding authors had full access to all the data in the study and had final responsibility for the decision to submit for publication.

Availability of data and materials
Our study follows the Guidelines for Accurate and Transparent Health Estimates Reporting (GATHER). Estimates can be further explored at national, state, and municipality level by age group, sex, and year through our online visualization tools (https://vizhub.healthdata.org/lbd/hiv-mort-la). The source code used to generate estimates, as well as the outputs of the study (including full sets of estimates at the state and municipality levels), are publicly available online via the Global Health Data Exchange (http://ghdx. healthdata.org/record/ihme-data/latin-america-hiv-mortality-estimates-2000-2 017). All maps presented in this study were generated by the authors and no permissions are required to publish them.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Competing interests
Dr. Singh reports personal fees from Crealta/Horizon, Medisys, Fidia, UBM LLC, Trio health, Medscape, WebMD, Clinical Care options, Clearview healthcare partners, Putnam associates, Focus forward, Navigant consulting, Spherix, Practice Point communications, the National Institutes of Health and the American College of Rheumatology, personal fees from Simply Speaking, other from Amarin, Viking, Moderana and Vaxart pharmaceuticals, nonfinancial support from FDA Arthritis Advisory Committee, non-financial support from Steering committee of OMERACT, an international organization that develops measures for clinical trials and receives arm's length funding from 12 pharmaceutical companies, non-financial support from Veterans Affairs Rheumatology Field Advisory Committee, non-financial support from Editor and the Director of the UAB Cochrane Musculoskeletal Group Satellite Center on Network Meta-analysis, outside the submitted work. Dr. Krishan reports grants from UGC Centre of Advanced Study, CAS II, awarded to the Department of Anthropology, Panjab University, Chandigarh, India, outside the submitted work. Prof. Postma reports grants and personal fees from various pharmaceutical industries, all outside the submitted work. Prof Postma holds stocks in Ingress Health and Pharmacoeconomics Advice Groningen (PAG Ltd) and is advisor to Asc Academics, all pharmacoeconomic consultancy companies. Dr. Ancuceanu reports he received consultancy and speakers' fees from various pharmaceutical companies. Walter Mendoza is a Program Analyst in Population and Development at the United Nations Population Fund-UNFPA Country Office in Peru, an institution which does not necessarily endorse this study. Dr. Pandi-Perumal reports a non-financial relationship with Somnogen Canada Inc. and occasional royalities from publishing houses, outside the submitted work.