Trends in Untreated Tuberculosis in Large Municipalities, Brazil, 2008–2017

We adapted a mathematical modeling approach to estimate tuberculosis (TB) incidence and fraction treated for 101 municipalities of Brazil during 2008–2017. We found the average TB incidence rate decreased annually (0.95%), and fraction treated increased (0.30%). We estimated that 9% of persons with TB did not receive treatment in 2017.

Finally, we obtained municipality-level sociodemographic covariates describing wealth level and healthcare access, 2 factors thought to be associated with TB incidence and the fraction receiving treatment. We used GDP per capita to describe wealth level (6). As an indicator of healthcare access we used municipal-level Family Health Strategy coverage (7); the Family Health Strategy is Brazil's method of delivering primary care (8). Finally, we obtained municipal-level population estimates for each study year from the Brazilian Institute of Geography and Statistics (Instituto Brasileiro de Geografia e Estatistica, IBGE) (6).
All data were deidentified and extracted from publicly available sources.

Model Description
The model estimates incidence as the number of persons exiting a state of undetected and untreated TB. Persons can exit this state by treatment initiation, death (before treatment initiation), or recovery (before treatment initiation). We estimated TB incidence, the fraction of cases receiving treatment, and the number of untreated cases (Appendix Table 1) as well as the average annual rates of change in these measures over the study period (Appendix Table 2).
We specified Poisson likelihood functions for SINAN case notification data and SIM mortality data.
For municipality i in year j, where γij represents population size, αij represents the modeled TB incidence rate, βij represents the modeled fraction of treated cases, δij represents the probability of death after treatment initiation, represents the probability of recovery without treatment, πh represents the estimated coverage of SIM (calculated at the state-level, denoted h), and ρij is an adjustment for misreporting of TB deaths in the SIM database, described below.
We specified exponential and inverse logit functions for incidence ( ) and fraction treated ( ), respectively: For municipality i in year j, where φ0 and ω0 are intercepts, Xij is a vector of the standardized municipal-level covariates (primary care access and log GDP per capita), and φ and ω are the associated vectors of regression coefficients for incidence and fraction treated respectively. The inclusion of these variables allows for partial pooling among municipalities with similar sociodemographic characteristics. Additionally, and are municipality-time random effects for incidence and fraction treated respectively. For each municipality these random effects are assumed to follow a random walk: For municipality i in year j, where ψ0 and ϕ0 are intercepts, ψ1i and ϕ1i are demeaned municipal-level random effects; ψ2ij and ϕ2ij are demeaned autoregressive municipality-year effects, set equal to zero at j = 1; and 1 , 2 , 1 , and 2 are standard deviation terms.
We estimated the probability of death among persons who initiated treatment ( ): Where bij is the probability that the treatment outcome is "death," cij is the probability that the treatment outcome is "loss to follow up," and τ is the probability that an individual dies given that they were lost to follow up. Values for bij and cij were estimated via logistic regression functions fitted to data for persons with a treatment outcome recorded (97.3% of all treated persons): with their associated regression coefficients. These regression estimates were used in preference to raw values to reduce the stochastic variation in the reported rates of these measures. We conducted a sensitivity analysis in which we coded all persons without a recorded treatment outcome as having abandoned treatment. While this led to an increase in the probability of abandonment, it did not meaningfully change estimates of incidence or fraction treated.
Finally, we estimated the systematic underreporting of TB as a cause of death ( ): For municipality i in year j, θ0 is the intercept, ηi is a demeaned municipal-level random effect, is the random-effects variance, θ1 is a linear time trend, xij is the percentage deaths in SIM attributed to a poorly-defined cause, and θ2 is the associated regression coefficient.
There is substantial uncertainty around true values for several model parameters. We used a Bayesian approach to represent and propagate this uncertainty through the analysis. We used prior probability distributions to summarize existing evidence on all model parameters (Appendix Table 3). The prior probabilities of μ, and θ3 were elicited through an expert opinion survey, described below. The prior distribution for τ was assumed to be Beta Stan code for the model is available at github.com/mel-hc/bayesian_subnat_est.

Model Priors from Expert Opinion
Prior distributions for the probability of recovery without treatment and underreporting of TB as a cause of death were created based on an expert opinion survey described by Chitwood et al . To review, the prior distribution for recovery without treatment was created based on median values for lowest, highest, and best-guess estimates of respondents. Incorporating expert opinion into the death adjustment was slightly more complex. Respondents were asked to estimate the rate of TB death misclassification among HIV negative persons in SIM in 2017 as it relates to the quality of cause of death reporting in a state or municipality, where the rate of ill-defined cases of death is an indicator for cause of death reporting quality. Two scenarios were presented: Scenario A, in which ≈1% of deaths an ill-defined cause, and Scenario B, where ≈15% of deaths had an ill-defined cause. Estimates were summarized as Beta distribution parameters. We incorporated these estimates into the death adjustment (ρij) estimate as follows: With the prior distributions: The parameters θ0, ηi, and θ2 are used to estimate the overall death adjustment, described above.

Model Performance
The total model run time was ≈37 minutes. There were no divergent transitions or iterations that saturated the maximum tree depth of 12. Key parameter prior distributions and posterior means, confidence intervals, effective sample sizes, and R-hat values are presented in Appendix Table 3.

Sensitivity Analysis
We tested the sensitivity of results to parameter prior distributions Altering individual model priors had little impact on model diagnostics or run time; while we observed differences in posterior distributions between the main run and sensitivity runs (Appendix Table 4), we did not observe meaningful changes to the distribution of our outcomes of interest (Appendix Figure   3). In addition to testing weaker priors, we ran the model with probability of survival without treatment fixed at 30% or 70%. While point estimates for several outcomes changed, they fell within the confidence intervals reported in the main analysis, and the relative burden of disease among modeled cities did not change.

Notes on Outliers
In 2 cities, Rio de Janeiro and São Vicente, model outcomes were indicative of biased case notification or TB-related death data. In Rio de Janeiro, we observed a decrease in TB-