Spatio-temporal patterns of childhood pneumonia in Bhutan: a Bayesian analysis

Pneumonia is one of the top 10 diseases by morbidity in Bhutan. This study aimed to investigate the spatial and temporal trends and risk factors of childhood pneumonia in Bhutan. A multivariable Zero-inflated Poisson regression model using a Bayesian Markov chain Monte Carlo simulation was undertaken to quantify associations of age, sex, altitude, rainfall, maximum temperature and relative humidity with monthly pneumonia incidence and to identify the underlying spatial structure of the data. Overall childhood pneumonia incidence was 143.57 and 10.01 per 1000 persons over 108 months of observation in children aged < 5 years and 5–14 years, respectively. Children < 5 years or male sex were more likely to develop pneumonia than those 5–14 years and females. Each 1 °C increase in maximum temperature was associated with a 1.3% (95% (credible interval [CrI] 1.27%, 1.4%) increase in pneumonia cases. Each 10% increase in relative humidity was associated with a 1.2% (95% CrI 1.1%, 1.4%) reduction in the incidence of pneumonia. Pneumonia decreased by 0.3% (CrI 0.26%, 0.34%) every month. There was no statistical spatial clustering after accounting for the covariates. Seasonality and spatial heterogeneity can partly be explained by the association of pneumonia risk to climatic factors including maximum temperature and relative humidity.

. Map of Bhutan with districts and sub-districts with altitude. Map was created using ArcMap 10.5 software (ESRI, Redlands, CA).  14.59) higher in children aged < 5 years as compared to those 5-14 years. Females were 11.3% (95% CrI 10.0%, 12.5%) less likely to develop pneumonia compared to males. The burden of pneumonia decreased by 0.3% (95% CrI 0.26%, 0.34%) every month. Each 1 °C increase in the maximum temperature was associated with a 1.3% (95% CrI 1.27%, 1.4%) increase in pneumonia cases. Conversely, each 10% increase in   www.nature.com/scientificreports/ relative humidity was associated with a 1.2% (95% CrI 1.1%, 12.5%) decrease in the incidence of pneumonia (Table 2). There was no statistical evidence of spatial clustering after accounting for the covariates (Table 2 and Fig. 4). There was > 95% probability of a lower than the national average trend of pneumonia in 56/205 sub-districts, whereas 66/205 sub-districts had > 95% probability of a trend more than the national average. There was no clear spatial pattern, with sub-districts showing higher and lower average trends across all the 20 districts (Fig. 5).

Discussion
This is the first study to have explored the climatic and geographic effects on childhood pneumonia in Bhutan. Pneumonia was spatially and temporally heterogeneous across sub-districts of Bhutan during the study period.
There was a decreasing trend in disease burden from 2010 through 2018, in addition to a strong seasonal pattern. Pneumonia mainly affected males and children aged < 5 years. Maximum temperature were associated with an increased incidence of pneumonia, while relative humidity was associated with a decreased incidence of disease. In addition to climatic factors, spatial heterogeneity could be due to differences in the socio-demographic characteristics of sub-districts. Known risk factors responsible for exacerbation and spread of pneumonia in Bhutan include low birth weight, malnutrition, indoor air pollution, overcrowded households, lack of breastfeeding in infants, and poor personal and environmental hygiene 24,25 . This was evident from the two districts of Haa and Paro, which have the lowest poverty levels 26 , and also reported lowest SMR associated with pneumonia. Similar to the decreasing trend in the incidence of global childhood pneumonia 6 , the pneumonia trend in Bhutan decreased during the study period from 2015 onwards. This could be attributed to a decrease in exposure to key risk factors including poor housing conditions and overcrowding, incomplete immunisation and malnutrition 6 . www.nature.com/scientificreports/ Pneumonia is the single largest infectious cause of death in children worldwide. It accounts for 15% of all deaths of children < 5 years 27 . In this study, children < 5 years were at a much higher risk of pneumonia compared to those 5-14 years old. Infants (aged between 0-11 months) were reported to contribute up to 24.2% of pneumonia cases in another study from Bhutan 28 . The WHO and United Nations Children's Fund (UNICEF) initiated a global action plan for pneumonia and diarrhoea (GAPPD) to accelerate pneumonia control in children 29 . The GAPPD strategies include promoting: (1) exclusive breastfeeding and adequate complementary feeding to protect children from pneumonia; (2) vaccination; (3) hand hygiene; (4) reductions in household air pollution; (5) HIV prevention; (6) co-trimoxazole prophylaxis for HIV infected and exposed children; and (7) treatment of childhood pneumonia with antibiotics and oxygen. Strengthening GAPPD strategies should be considered in Bhutan, as is the case in other countries in South Asia (including Bangladesh and India). The introduction of pneumococcal conjugate vaccines (PCVs) in Bhutan in 2019 is timely in prevention of pneumonia, and the effect of introduction of PCV on the burden of childhood pneumonia in Bhutan must be evaluated 21,30 . Exclusive breastfeeding rates from birth until six months in Bhutan varies from 35.9 to 51.0% 31,32 . Increasing exclusive breastfeeding rates are likely to reduce pneumonia associated morbidity 33 .
Pneumonia incidence was highly seasonal, and was associated with climatic factors including temperature and relative humidity. The association of temperature with pneumonia has been reported in other studies 34,35 . A plausible explanation is the association of higher temperature with air pollution 36 which in itself is known risk factor and cause of pneumonia 7,37-39 . Most industries are located in the southern parts of Bhutan where air pollution is expected to be higher as compared to other districts. This was reflected by these sub-districts having higher SMR for pneumonia. Additionally, traditional methods of cooking in rural Bhutan using fire wood could also contribute to respiratory illness such as pneumonia 40,41 .
The incidence of pneumonia tends to be higher during the rainy season but not so in this study [42][43][44] . Rainfall may trigger socio-ecological behavioural changes such as increased contact between people and the distribution of pathogens. Further, heavy rainfall during the monsoon is likely to pollute drinking water, particularly surface water from streams, which is the main drinking water source for rural populations 45 . Unsafe drinking water and sanitation are important drivers of pneumonia 46 . Relative humidity was associated with a decrease in pneumonia incidence in this study which is in concordance with other studies 35,47 . Higher relative humidity decreases the survival of lipid-enveloped viruses such as influenza A, influenza b and respiratory syncytial virus 48,49 .
There are a number of limitations that need to be considered when interpreting the results of this study. First, the study used routine case reports to measure pneumonia incidence. Known issues exist surrounding completeness and representativeness of such data. Secondly, the causal organisms of pneumonia were not available and the association could be different based on the organisms. Thirdly, there was no reconciliation to accommodate different levels of aggregation of the climate variables (district) and the disease data (sub-district), and the climate conditions were assumed to be homogeneous within a district. Lastly, unaccounted risk modifiers were not included in the modelling due to a lack of available data. These important unmeasured factors, such as immunization coverage, air pollution level, living standards and socio-economic status, crowding, smoking, access to safe drinking water and latrine usage might have resulted in confounding, which was not able to be quantified 39,50,51 . www.nature.com/scientificreports/ Despite these limitations, the strengths of this study are the capacity to implement the spatial analysis at a relatively fine resolution, being the sub-district level, and over an extended time series (108 months). Traditionally, spatial patterns of infectious disease risk have been displayed at larger geographical units, such as a district, province, national, regional, and global scales 45,52,53 . Such low resolution can mask localized disease patterns due to averaging 54 .

Conclusion
Pneumonia is an important childhood disease and the introduction of PCVs to reduce the burden of this disease is timely. Pneumonia was highly seasonal and spatially heterogeneous across sub-districts. Seasonality can be explained by climatic factors including maximum temperature and relative humidity. The spatial and temporal variability of pneumonia should inform prevention and control efforts in Bhutan, through rational decision making and proper resources allocation.

Materials and methods
Study area. Bhutan 2% (452,178) of the population live in rural areas and practice subsistence farming. The altitude ranges from 75 m above sea level in the south to more than 7000 m in the Himalayas (Fig. 1).
Study design and data source. This is a retrospective study using secondary data on pneumonia from January 2010 to December 2018, stratified by sex and age (< 5 years and 5-14 years) at the sub-district level. The data were obtained from the National Acute Respiratory Infections surveillance system, hosted by the Bhutan Health Information and Management Systems (HIMS) under the Bhutan Ministry of Health. These data contain all pneumonia cases treated by health centres including hospitals and primary health care facilities and reported to the HIMS every month. Pneumonia is defined as "a patient with history of cough or reported breathing difficulty, and increased respiratory rate, (respiratory rate ≥ 50 breaths per minute in children aged two months or more and less than 12 months, or repiratory rate ≥ 40 breaths per minute in children aged 12 months or more and less than 60 months) or chest indrawing" 56 . Daily climatic variables (rainfall, relative humidity, minimum and maximum temperature) were obtained from the National Centre for Hydrology and Meteorology under the Ministry of Economic Affairs of Bhutan. Monthly average climatic variables were calculated at the district level because they were available only at the district level. Population estimates used in the study were obtained from the National Statistical Bureau, Bhutan 57 . Administrative boundary maps were downloaded from the DIVA-GIS website 58 .
Crude standardized morbidity ratios. An initial descriptive analysis of pneumonia incidence across the country was conducted. Crude SMR for each sub-district were calculated using the following formula: where Y is the overall SMR in sub-district i, O is the total number of observed pneumonia cases over the entire study period in the sub-district and E is the expected number of pneumonia cases in the sub-district across the study period. The expected number was calculated by multiplying the national incidence by the average population for each sub-district over the study period.
Exploration of seasonal patterns and inter-annual patterns. The time series of pneumonia incidence was decomposed using STL weighted regression to show: the seasonal pattern, inter-annual patterns and the residual variability. The STL model was structured as follows: where Y t represents numbers of local pneumonia cases with logarithmic transformation, S t is the additive seasonal component, T t is the trend, and R t is the "remainder component"; t is time in months 59,60 . Spatio-temporal model. A Bayesian statistical framework was deployed for spatial analysis. It provides a convenient framework for the simultaneous inclusion of covariates and spatial autocorrelation in a single model, while providing robust evaluation and expression of uncertainty. The posterior distributions can be used to quantify uncertainties in parameters of interest (e.g., covariate effects and spatial patterns of disease risk) 61 .
Initially, a preliminary bivariate Poisson regression of pneumonia cases was undertaken to select the covariates. The covariates were altitude and climatic variables including rainfall, minimum and maximum temperature and relative humidity without lag, and lagged up to 3 months. The covariates with a p-value of < 0.05 and the lowest Akaike's information criterion (AIC) were selected (Supplementary Table 1). The co-linearity of the selected climatic and environmental variables was tested using variance inflation factors (VIF) (Supplementary Table 2). In the final model, altitude, and rainfall, maximum temperature without lag and relative humidity lagged at 3 months were included.
Of the 83,315 observations stratified by sub-districts, age (< 5 and 5-14 years) and sex over an observation period of 108 months, there were 63,742 (72%) zero counts of pneumonia. Therefore, Zero-inflated Poisson (ZIP) regression was constructed in a Bayesian framework (Supplementary Table 3). The first model (Model I), www.nature.com/scientificreports/ assumed that spatial autocorrelation was not present in the relative risk of pneumonia. This model was developed with selected climatic factors (rainfall, maximum temperature and relative humidity), altitude, age (< 5 and 5-14 years) and gender as explanatory variables, and an unstructured random effect for sub-districts. The second model (Model II) contained a spatially structured random effect in addition to the aforementioned covariates. Model III, a convolution model, contained all of the components of the preceding two models. The model with the lowest DIC was selected as the final explanatory model. Model III assumed that the observed counts of pneumonia, Y, for the i th sub-district (i = 1 … 205) in the j th month (January 2010-December 2018) followed a Poisson distribution with mean (μ ij ), that is, where the expected number of cases in sub-district i, month j (acting as an offset to control for population size) was represented by E ij and θ ij is the mean log relative risk (RR). The intercept (α), and coefficients for age (5-14 years as reference), sex (male as reference), monthly trend, altitude, rainfall, relative humidity and maximum temperature are β 1 , β 2 , β 3 , β 4 , β 5 , β 6 and β 7 . The spatially unstructured and structured random effects are represented as u i and s i , respectively, with u i excluded from Model II and s i excluded from Model I. Spatiotemporal random effect with a mean of zero and variance of σ w 2 was denoted by w i as in other studies 62,63 . A conditional autoregressive (CAR) prior structure was used to model the spatially structured random effect. Spatial relationships between the sub-districts were based on a 'queen' contiguity matrix 64 . A weight of 1 was assigned to sub-districts sharing a border and 0 otherwise. A flat prior distribution was specified for the intercept, whereas non-informative normal prior distributions (mean = 0.0, precision = 0.00001) were used for the coefficients. The priors for the precision of unstructured and spatially structured random effects were specified using non-informative gamma distributions with shape and scale parameters of 0.5, 0.001.
Markov Chain Monte Carlo simulation was used to estimate model parameters 65 . The models were run for an initial 10,000 iterations for burn in, which were then discarded. Subsequently, visual inspection of posterior density and history plots were used to note convergence at intervals of 20,000 iterations. Convergence occurred at approximately 100,000 iterations for all models. Following convergence, posterior distributions from model parameters were stored for inference. Summaries of parameters were calculated, including posterior mean and 95% CrI. In the Bayesian analyses, RR with 95% CrI that excluded 1 were considered statistically significant.
Seasonality decomposition was carried out using the R statistical software, version 3.3.1 66 . The ZIP regression model was constructed using WinBUGS software, version 1.4.3 (MRC Biostatistics Unit 2008) 67 . ArcMap 10.5 software (ESRI, Redlands, CA) was used to generate maps of the posterior means of the unstructured and structured random effects and the spatiotemporal random effects.