Spatio-temporal patterns of pneumonia in Bhutan: A Bayesian analysis

26 Pneumonia is one of the top 10 diseases by morbity in Bhutan. This study aimed to investigate 27 the spatial and temporal trends and risk factors of pneumonia in Bhutan. A multivariable Zero- 28 inflated Poisson regression using a Bayesian Markov chain Monte Carlo simulation was 29 undertaken to quantify associations of age, sex, rainfall, maximum temperature and relative 30 humidity with monthly pneumonia incidence and identify underlying spatial structure of the 31 data. Overall pneumonia incidence was 96.5 and 4.57 per 1,000 populations over nine years in 32 people aged <5 years and ≥5 years, respectively. Children <5 years or being a female are more 33 like to get pneumonia than ≥5 years a nd males. A 10mm increase in rainfall and 1°C increase 34 in maximum temperature was associated with a 7.2% (95% (credible interval [CrI] 0.7%, 35 14.0%) and 28.6% (95% CrI 27.2%, 30.1%) increase in pneumonia cases. A 1% increase in 36 relative humidity was associated with a decrease in the incidence of pneumonia by 8.6% (95% 37 CrI 7.5%, 9.7%). There was no evidence of spatial clustering after accounting for the 38 covariates. Seasonality and spatial heterogeneity can partly be explained by the association of 39 pneumonia risk to climatic factors including rainfall, maximum temperature and relative 40 humidity. 41


43
Pneumonia is a major cause of morbidity and mortality worldwide 1 . Each year, pneumonia 44 accounts for over 12 million hospital admissions and 1.3 million deaths in children aged less 45 than 5 years worldwide 2,3 . In 2017, pneumonia was the fourth-leading cause of death and it is 46 estimated that it will be the third-leading cause of death by 2040 4 . The World Health 47 Organization (WHO) estimates that respiratory infections account for 6% of the total global 48 burden of disease. This accounts for a higher percentage compared to the burden of diarrheal 49 disease, cancer, human immunodeficiency virus (HIV) infection, ischemic heart disease or 50 malaria 5 . 51 Pneumonia is a potentially life-threatening illness with a particularly high burden in South Asia 52 and sub-Saharan Africa 3,6,7 . It is not only a major cause of morbidity and mortality but is also 53 associated with a substantial economic burden on healthcare systems 8,9 and household 54 income 10 . Pneumonia often has a complex aetiology involving multiple pathogens, including 55 many that are transmitted person-to-person. Past time-series analyses have identified various 56 pneumonia and influenza outcomes to be temporally seasonal, demonstrating highly consistent 57 peaks in winter months and troughs in summer months 11,12 . Other studies have found that 58 pneumonia admissions were highly spatially clustered 13 , driven by contact with infected people 59 during indoor activities 14 . 60 Pneumonia continues to be an important communicable disease in Bhutan-locate in the Eastern 61 Himalayas 15-17 (Fig. 1). In 2019, pneumonia was one of the top-ten ranked diseases in terms of 62 morbidity and accounted for 19% of the overall disease burden 18 . Every year the Bhutanese 63 government spends a huge amount on the treatment and management of pneumonia. In the 64 financial year 2017-2018, 7.1% of current health expenditure was spent on treating infectious 65 respiratory diseases 19,20 . Despite the importance of pneumonia, and the infectious nature of the 66 disease, there have been no previous studies to understand the underlying ecological drivers of 67 a 7.2 % (95% CrI 0.7%, 14.0%) increase in incidence of pneumonia. Similarly, a maximum 92 temperature increase of a 1°C was associated with a 28.6% (95% CrI 27.2%, 30.1%) increase 93 in pneumonia cases. However, a 1% increase in relative humidity was associated with a 94 decrease in the incidence of pneumonia by 8.6% (95% CrI 7.5%, 9.7%) ( Table 2). 95 There was no evidence of spatial clustering after accounting for the covariates ( Table 2 and 96 Fig. 4). There was >95% probability of a higher than the national average trend of pneumonia 97 in 56/205 sub-districts, whereas 67/205 sub-districts had >95% probability of a trend less than 98 the national average. There was no clear spatial pattern, with sub-districts showing higher and 99 lower average trends across all the 20 districts (Fig. 5). 100 101

102
Pneumonia was spatially and temporally heterogeneous across sub-districts of Bhutan during 103 the study period. There was a decreasing trend, in addition to a strong seasonal pattern during 104 the study period. Pneumonia mainly affected children aged <5 years and females. Rainfall and 105 maximum temperature were associated with an increased incidence of pneumonia while 106 relative humidity was associated with a decrease incidence. 107 In addition to climatic factors, spatial heterogeneity could be due to differences in the socio-108 demographic characteristics of sub-districts. The risk factors responsible for exacerbation and 109 spread of pneumonia in Bhutan were low birth weights, malnutrition, smoky and overcrowding 110 households, bottle-feeding of infants and poor personal and environmental hygiene 23,24 . This 111 was evident from the two districts of Haa and Paro which has the lowest poverty 25 , and also 112 reported lowest SMR of pneumonia. Similar to the decreasing trend in the incidence of global 113 childhood pneumonia 26 , the national pneumonia trend decreased during the study period. This 114 could be attributed to a decrease in exposure to key risk factors including poor housing 115 conditions and overcrowding, incomplete immunisation and malnutrition 26 . 116 Pneumonia is the single largest infectious cause of death in children worldwide. It accounts for 117 15% of all deaths of children <5 years 27 . In this study, children <5 years were at a much higher 118 risk of pneumonia compared to those ≥5 years. Infants (aged between 0-11 months) was 119 reported to contribute up to 24.2% of cases in another study 28 . The WHO and United Nations 120 Increasing exclusive breastfeeding rates are likely to reduce pneumonia associated morbidity 33 . 131 Pneumonia was highly seasonal and was associated with climatic factors including 132 temperature, rainfall and relative humidity. The association of temperature with pneumonia has 133 been reported in other studies 34,35 . A plausible explanation is the association of higher 134 temperature with air pollution 36 which in itself is known risk factor and cause of pneumonia 37-135 40 . Most industries are located in the southern parts of Bhutan where air pollution is expected 136 to be higher as compared to other districts. This was reflected by these sub-districts having 137 higher SMR for pneumonia. Additionally, traditional methods of cooking in rural Bhutan using 138 fire wood could also contribute to respiratory illness such as pneumonia 41,42 . 139 The incidence of Pneumonia tends to be higher during the rainy season 43-45 . Rainfall may 140 trigger socio-ecological behavioural changes such as increased contact between people and the 141 distribution of pathogens. Further, heavy rainfall during the monsoon is likely to pollute 142 drinking water, particularly the surface water from streams, which is the main drinking water 143 source for rural populations 46 . Unsafe drinking water and sanitation are important drivers of 144 pneumonia 47 . Relative humidity was associated with a decrease in pneumonia incidence in this 145 study which is in concordance with other studies 35,48 . Higher relative humidity decreases the 146 survival of lipid-enveloped viruses such as influenza A, influenza b and Respiratory Syncytial 147 There are a number of limitations that need to be considered when interpreting the results of 149 this study. First, the study used routine case reports to measure incidence of pnuemonia. Known 150 issues exist surrounding completeness and representativeness of such data. Secondly, the causal 151 organisms of pneumonia were not available and the association could be different based on the 152 organisms. Thirdly, there was no reconciliation to accommodate different levels of aggregation 153 of the climate variables (district) and the disease data (sub-district), and the climate conditions 154 were assumed to be homogeneous within a district. Lastly, unaccounted risk modifiers were 155 not included in the modelling due to a lack of available data. These important unmeasured 156 factors, such as immunization coverage, air pollution level, living standards and socio-157 economic status, crowding, smoking, access to safe drinking water and latrine usage might 158 have resulted in confounding, which was not able to be quantified 39,51,52 . 159 Despite these limitations, the strengths of this study are the capacity to implement the spatial 160 analysis at a relatively fine resolution, being the sub-district level, and over a long time series 161 (108 months). Traditionally, spatial patterns of infectious disease risk have been displayed at 162 larger geographical units, such as a district, province, national, regional, and global 163 scales 46,53,54 . Such low resolution can mask localized disease patterns due to averaging 55 . 164

165
Pneumonia is an important childhood disease and the introduction of pneumococcal conjugate 166 vaccines to reduce the burden of this disease is timely. Pneumonia was highly seasonal and 167 spatially heterogeneous across sub-districts. Seasonality can be explained by climatic factors 168 including temperature, rainfall and relative humidity. The spatial and temporal variability of 169 pneumonia should inform in better targeting of its prevention and control in the country through 170 rational decision making and proper resources allocation. 171 live in rural areas and practice subsistence farming. The altitude ranges from 75m above sea 177 level in the south to more than 7000m in the Himalayas (Fig. 1). 178

Study design and data source 179
This is a retrospective study using secondary data on pneumonia from January 2010 to 180 An initial descriptive analysis of pneumonia incidence across the country was conducted. 196 Crude SMR for each sub-district were calculated using the following formula: 197 Where Y is the overall SMR in sub-district i, O is the total number of observed pneumonia 199 cases over the entire study period in the sub-district and E is the expected number of pneumonia 200 cases in the sub-district across the study period. The expected number was calculated by 201 multiplying the national incidence by the average population for each sub-district over the 202 study period. 203

Exploration of seasonal patterns and inter-annual patterns 204
The time series of pneumonia incidence was decomposed using STL weighted regression to where expected number of cases in sub-district i, month j (acting as an offset to control for 241 population size) was represented by Eij and θij is the mean log relative risk (RR). The intercept 242 (α), and coefficients for age (≥5 as reference), sex (male as reference), monthly trend, rainfall, 243 relative humidity and maximum temperature are β1, β2, β3, β4, β5 and β6. The spatially 244 unstructured and structured random effects are represented as ui and si, respectively, with ui 245 excluded from Model II and si excluded from Model I. Spatiotemporal random effect with a 246 mean of zero and variance of σw 2 was denoted by wi as in other studies 63,64 . 247 A conditional autoregressive (CAR) prior structure was used to model the spatially structured 248 random effect. Spatial relationships between the sub-districts were based on a 'queen' 249 contiguity matrix. A weight of 1 was assigned to sub-districts sharing a border and 0 otherwise. 250 A flat prior distribution was specified for the intercept, whereas a non-informative normal prior 251 distribution was used for the coefficients. The priors for the precision of unstructured and 252 spatially structured random effects were specified using non-informative gamma distributions 253 with shape and scale parameters equal to 0.01. 254 The model was run for an initial 10,000 iterations, which were then discarded. Subsequently, 255 visual inspection of posterior density and history plots were used to note convergence at 256 intervals of 20,000 iterations. Convergence occurred at approximately 100,000 iterations for 257 all models. Following convergence, posterior distributions from model parameters were stored 258 for inference. Markov Chain Monte Carlo simulation was used to estimate model parameters 259 65 . Summaries of parameters were calculated, including posterior mean and 95% credible CrI. 260 In all analyses, an α-level of 0.05 was adopted to indicate statistical significance (as indicated 261 by 95% CrI for relative risks (RR) that excluded 1). 262 Seasonality decomposition was carried out using the R statistical package, release 3.     Figure 1 Map of Bhutan with districts and sub-districts with altitude. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. ARISupplementary.pdf