Estimating Under Reporting of Leprosy in Brazil using a Bayesian Approach

Data source

We used data from the Brazilian National Notifiable Diseases Information System (SINAN). We collected the total and new number of leprosy cases for each of the 558 mainland microregions of Brazil, from 2007 up to 2015. New cases are defined, by the Ministry of Health, as cases with no previous treatment (Brasil 2016). Confirmed cases of leprosy are all cases that presented one of the following conditions: (a) lesion(s) and/or area(s) of the skin with changes in thermal and/or painful and/or tactile sensibility; (b) peripheral nerve thickening associated with sensory and / or motor and / or autonomic changes; (c) presence of Mycobacterium leprae bacilli, confirmed by intradermal smear microscopy or skin biopsy. New cases of leprosy can be reported through an active detection, by performing epidemiological investigation of contacts and collective exams, such as surveys and campaigns (rapid serological tests), or by passive detection, when individuals spontaneously search health care systems.

The detection of the disease depends on the manifestation of symptoms in infected individuals. However, differences in the resistance of their immunological systems may contribute to the appearance of symptoms only after a long period after incubation of the bacterium, ranging from 2 up to 7 years (Brasil 2001). Additionally to such individual conditions, endemic levels and unfavorable socioeconomic factors, as well as access to public health services and a high rate of housing occupation, influence the pattern of the disease incidence in the country. Based on these factors, socioeconomic indicators and living conditions of the population in each microregion are considered to account for space variability of the disease incidence. Therefore, to conduct our analysis we collected the following information: (1) the proportion of examined household contacts with relation to the total registered household contacts (x₁), obtainded from SINAN; (2) coverage of the cash transfer program Bolsa Familia in terms of the proportion of the population at risk that is assisted (x₂); (3) coverage of the Family Health Program in terms of the proportion of registered population (x₃); (4) average number of people per household (x₄), obtained from the Atlas of Human Development in Brazil (2010); and (5) proportion of people living in urban areas (x₅), obtained from the Brazilian Institute of Geography and Statistics - IBGE (2010 demographic census).

The evolution of symptoms will be reflected into disabilities, ranging from low, when patients have no problems in eyes, hands or feet, to high consequences, as visible deformations in hands, feet or eyes or even severe visual impairment. Therefore, new cases of patients presenting high levels of disabilities is an indirect indicator of diagnosis delay (Nobre et al. 2017), which leads to underreporting. Thus, in this work, the proportion of diagnosed new leprosy cases with Grade 2 of physical disability (w), collected from SINAN, was considered in the characterization of the underreporting mechanism throghout the Brazilian microregions.

Data Curation

To peform our analises, we negleted missing informations. Complementary to this, from the total 311,970 new cases, with indentification of microregion, we excluded Fernando de Noronha from the analysis, once it is an island that stays more than 500 km from the mainland Brazil and might generate problems for the spatial analysis. Therefore, for the spatial analysis, we considered 557 mainland Brazilian microregions. After that, it remained 311,968 new cases from the period 2007 a 2015.

To obtain the numerator to calculate the proportion of diagnosed new leprosy cases with GD2, we exclude from the analysis the null and ”not evaluated” values from the variable type of diagnosis. They represented 10.35% of the total number registered. Therefore, 279,719 had information about the type of diagnosis.

The proportion of household examined contacts was obtained as the ratio between the number of examined contacts and the number of identified contacts in each microregion. We excluded the null values and the cases with less than 43 contacts in those variables. We obtained a total of 267,160 with information.

Statistical model for underreported leprosy data in Brazil

In disease surveillance and epidemiological studies, two main statistical approaches have been considered to deal with underreport counts. Inference under the first one is based on a censored Poisson likelihood function, allowing the estimation of both the disease rate and the probability of the observed count being underreported in each area (Bailey et al. 2005; Oliveira et al., 2017). The second and potentially more flexible approach relies on the specification of a hierarchical Poisson model from which is possible to estimate both the disease rates and the proportion of reported cases in each area (Stoner et al. 2019; Dvorzak and Wagner 2015; Shaweno et al. 2017; Stamey et al. 2006; Whittemore and Gong 1991; Papadopoulos and Silva 2012). Based on the sample information we have available, we consider the Bayesian framework proposed by Stoner et al. 2019 to fit the Brazilian leprosy data.

Under such a framework, in each microregion i (for i = 1,…, 557), the reported (observed) count Y_i is modeled as a Binomial random variable, where the number of trials is an unobserved Poisson variable T_i corresponding to the true number of cases that has been incompletely recorded. The true count generating process is modeled through the mean of the Poisson variable, denoted by μ_i, and the reporting mechanism is modeled through the Binomial probability, denoted by ϵ_i. Then, the basic structure for the hierarchical model is: for i = 1,…, 557. Since T_i is not observed, statistical inference is based on the marginal distribution of Y_i obtained from the joint model given in equation (1), which is

We assume that the mean expected new cases of leprosy, μ_i, depends on the five socioeconomic covariates x₁ to x₅ previously described, such that

We also include the logarithm of the total exposure population, log(P_i), as an offset, so that parameter μ_i can be interpreted as the leprosy incidence rate at area i. Additionally, in order to capture any residual variation in the leprosy, we include a spatially structured random effect ϕ_i and a local unstructured random effect δ_i in the log-linear predictor of μ_i.

To model ϵ_i, the probability of reporting a new leprosy case at each area i, we made use of the covariate w: the percentage of diagnosed new leprosy cases with GD2 of physical disability, such that

As previously discussed, covariate w acts as a proxy of the appropriate variable that accounts for notification efficiency of leprosy new cases in each microregion. An unstructured random effect γ_i is included in the logistic regression model for ϵ_i in order to account for potential effects of unobserved covariates that may influence the detection of leprosy.

In the statistical literature, the joint model defined by equations (2), (3) and (4) is called by the Poisson-Logistic (or Pogit) model. It worth noticing that the Pogit model provides a straightforward predictive analysis for the proportion of leprosy new cases that were not observed, denoted by Z_i = T_i − Y_i. Given and ϵ_i, Z_i can be predicted from the distribution for all i.

Despite its appealing features, it is well-known in the literature that the Pogit model suffers with the lack of identifiability. This occurs because only the product η_i = μ_iϵ_i is identified from the observed data since any other parameter combination, say and , such that yields the same likelihood function. In practice, that concept means that additional information must be introduced in the model in order to distinguish between parameters μ_i and ϵ_i, thus conducting to a meaningful inference about the true count T_i. Such extra information can be provided by validation datasets, active search surveys and experts’ opinion. The source of information to be used will depend on which one is available for the specific practical situation one is dealing with.

To overcome the identifiability issue when fitting the Pogit model to the Brazilian leprosy count data, we follow the approach of Stoner et al. 2019. Whenever the regression models of μ_i and ϵ_i do not share any common covariate, the lack of identifiability of the Pogit model relies on the confounding between the two intercepts α₀ and β₀. Then, special attention must be given to their prior distributions. As discussed in Stoner et al. 2019, by taking x₁ to x₅ and w as centered covariates it provides that β₀ and α₀ are, respectively, interpreted as the mean reported number of leprosy new cases (on the log scale) and the mean reporting rate (on the logistic scale) when the covariates are at their centering value. In this context, the appealing interpretation of α₀ and β₀ can be appropriately used to elicit an informative prior distribution for one of them, thus providing an identifiable Pogit model.

We found that there is no study seeking to estimate the detection rate of leprosy over the country between the period of 2007 to 2015. If such a study was available, we could center the covariate w with respect to its observed mean and then elicit an appropriate informative prior for parameter α₀, which, in this case, could be interpreted as the overall mean reporting rate. In the work by Cunha et al. 2015, the authors present results of an active search survey, performed in 2012, in some municipalities of the state of Amazon, Northern of Brazil. Based on their findings, the microregion Rio Preto da Eva (RPE), composed by the municipalities Rio Preto da Eva and Presidente Figueiredo, showed reporting rates of approximately 91% when comparing the total cases registered by the local health centers. In particular, the microregion of RPE is located in the metropolitan region of Manaus, capital of the State of Amazonas, and it is one of the most developed regions in northwest of Brazil.

As far as we know, this is the most trustful available study about the level of leprosy underreporting in a Brazilian microregion. Therefore, we rely on the information provided for this microregion to define an informative prior distribution for parameter α₀. To do so, we center the covariate w with respect to its observed value for Itacoatiara so that the parameter α₀, on the logarithmic scale, represents the average level of notification of leprosy cases in such microregion. Moreover, we assume that the information for the level of notification for the period from 2007 to 2015 can be represented by findings of Cunha et al. 2015 for year 2012. We then specified a N(2.5, 0.3) as the prior distribution for the parameter, which provides an a priori average level of reporting of approximately 91% for this microregion, on the logistic scale. The choice for the variance term was guided by the sensitivity studies presented by Stoner et al. 2019. They suggest that the model is robust in terms of quantifying uncertainty, as long as the practitioner specifies a prior for α₀ that is informative without the need of being a degenerated one (a prior with null variance). Through an a priori predictive analysis in our application, the variation expressed by the prior distribution N(2.5,0.3) provides values for ϵ in the microregion RPE with high concentration (about 90%) between 88.2% and 95.1%, reflecting well our belief regarding the average level of reporting used as a reference for the period.

We now provide the prior distribution for parameter β₀. As an informative prior was already elicited for the other model intercept, α₀, there is no need for doing the same in relation to β₀. Nevertheless, to avoid the generation of unrealistic values for the leprosy incidence rate (especially quite elevated values), we also follow the approach proposed in Stoner et al. 2019. Each covariate x₁ to x₅ was centered with relation to its observed mean such that is interpreted as the mean reported number of leprosy cases, on the log scale. Then, we assume a prior N(−8,1) for β₀ to represent our belief that is not plausible a very high value (such as over 1 million) for the total number of new leprosy cases. For each of the remaining regression coeficients α₁ and β₁, …, β₅ we elicit a relatively noninformative Normal prior distribution N(0,10²).

The model specification is completed with the following prior distributions for the random effects and their precision terms. As it is usual in the literature, the prior distribution of ϕ_i, with precision parameter ν, is represented by an intrinsic conditional autoregressive (iCAR) model (Besag et al. 1991). Here, a neighbor of an area i was defined as any i′ ≠ i sharing a geographical boundary with i. Both the unstructured effect γ_i and the local effect δ_i are assumed to have a Normal prior distribution , and , respectively. For each of the precision parameters v, and we elicited a Gamma distribution G(1,1).

The model was implemented using package NIMBLE from software R (R Core Team, 2015). For the MCMC scheme two chains were used, each ran for a total of 600K iterations. After discarding 200K iterations as burn-in and using a lag of 200, the potential scale reduction factor (PSRF) proposed by Brooks and Gelman (1998) was computed as less than 1.05 for all regression coefficients and precision parameters, indicating convergence to the target posterior distribution.

To summarise, we applied the Bayesian model to correct the underreporting of leprosy new cases, from which we estimate the posterior mean (Mean), the posterior standard deviation (SD) and the 95% highest posterior density (HPD) intervals (95%-HPD) for all regression parameters defined on Equations (3) and (4) (Table 1). We also provide the incidence rate ratio (IRR) for Poisson parameters and the odds ratio (OR) for logistic parameters with their respective 95% HPD intervals (95%-HPD).

For each model coefficient, we provide the IRR (or OR) and its 95% HPD interval. When such interval contains the value 1 it can be considered that the covariate effect is not significant. The same if the 95% HPD interval for the coefficient contains the value zero. The IRR indicates the effect in the mean incidence when a unity is increased in the explanatory variable.