Variation in smoking attributable all-cause mortality across municipalities in Belgium, 2018: application of a Bayesian approach for small area estimations

Abstract

Background

Smoking is one of the leading causes of preventable mortality and morbidity worldwide, with the European Region having the highest prevalence of tobacco smoking among adults compared to other WHO regions. The Belgian Health Interview Survey (BHIS) provides a reliable source of national and regional estimates of smoking prevalence; however, currently there are no estimates at a smaller geographical resolution such as the municipality scale in Belgium. This hinders the estimation of the spatial distribution of smoking attributable mortality at small geographical scale (i.e., number of deaths that can be attributed to tobacco). The objective of this study was to obtain estimates of smoking prevalence in each Belgian municipality using BHIS and calculate smoking attributable mortality at municipality level.

Methods

Data of participants aged 15 + on smoking behavior, age, gender, educational level and municipality of residence were obtained from the BHIS 2018. A Bayesian hierarchical Besag-York-Mollie (BYM) model was used to model the logit transformation of the design-based Horvitz-Thompson direct prevalence estimates. Municipality-level variables obtained from Statbel, the Belgian statistical office, were used as auxiliary variables in the model. Model parameters were estimated using Integrated Nested Laplace Approximation (INLA). Deviance Information Criterion (DIC) and Conditional Predictive Ordinate (CPO) were computed to assess model fit.

Population attributable fractions (PAF) were computed using the estimated prevalence of smoking in each of the 589 Belgian municipalities and relative risks obtained from published meta-analyses. Smoking attributable mortality was calculated by multiplying PAF with age-gender standardized and stratified number of deaths in each municipality.

Results

BHIS 2018 data included 7,829 respondents from 154 municipalities. Smoothed estimates for current smoking ranged between 11% [Credible Interval 3;23] and 27% [21;34] per municipality, and for former smoking between 4% [0;14] and 34% [21;47]. Estimates of smoking attributable mortality constituted between 10% [7;15] and 47% [34;59] of total number of deaths per municipality.

Conclusions

Within-country variation in smoking and smoking attributable mortality was observed. Computed estimates should inform local public health prevention campaigns as well as contribute to explaining the regional differences in mortality.

Introduction

Smoking is one of the main causes of preventable mortality and morbidity worldwide. According to the World Health Organization (WHO), globally 12% of all deaths among adults aged 30 years and over were attributed to tobacco [1]. Tobacco use has a proven causal relationship with a number of chronic conditions, like cancer, cardiovascular disease, chronic respiratory disease and others [2]. In the COVID-19 pandemic, current smokers who contracted coronavirus were shown twice as likely to be hospitalized and tended to report more symptoms than non-smokers [3], although some controversy remains whether smokers are less likely to get infected with this virus [3].

Despite notable successful efforts by WHO and national governments worldwide to limit tobacco use [4], smoking prevalence in European countries remains a public health concern. Among the WHO regions, Europe has the highest prevalence of tobacco smoking among adults (28% in overall population, 38% in males and 19% in females) and one of the highest prevalence of tobacco use by adolescents [5]. It is thus not surprising that, compared to the rest of the world, the WHO European Region has one of the highest proportions of deaths attributable to tobacco use. WHO has estimated that tobacco use is currently responsible for 16% of all deaths in adults over 30 in European Region, with many of these deaths occurring prematurely [5].

Smoking prevalence in Belgium is slightly lower than the EU-15 average and was estimated at 15% in 2018. This proportion is higher in men (18%) than in women (12%). On the positive side, the prevalence of daily smoking has decreased by 40% between 1997 and 2018 [6]. Substantial regional variation in smoking rates was reported with higher prevalence in Wallonia (18%), followed by Brussels (16%) and Flanders (13%) and socio-economic disparities in smoking are persistent in Belgium similar to other countries [7,8,9]. Differences at lower geographical level (e.g. municipalities) may thus also exist while these have never been studied. These differences, if substantial, may be highly relevant for public health policies and prevention campaigns. Absence of municipality level statistics on smoking prevalence further hinders estimates of spatial distribution of smoking attributable mortality (i.e., the number of deaths that can be attributed to tobacco) in support of local health policy-making and understanding regional variations in mortality.

Since 1997, the Belgian Health Interview Surveys (BHIS) has been regularly collecting data on tobacco use in Belgium. The sampling design allows reliable estimates of many health indicators including smoking on national and subregional level, however, as only approximately one fourth of nearly 600 municipalities has been effectively sampled each survey round, direct estimates of smoking at municipality level are precluded. Recently, novel statistical methods have been proposed to address the growing demand for estimating small-area level indicators in situations when few or no data in some area precludes direct estimations [10]. In such cases, indirect estimates may be computed. Such estimates borrow strength from related areas through linking models based on auxiliary data such as recent census and administrative records [11]. Accounting for the complex survey design in computing these indirect estimates poses additional methodological challenges.

The aim of this study was to obtain estimates of smoking prevalence in each Belgian municipality using the BHIS 2018 data and to calculate smoking attributable mortality at municipality level using all-cause mortality statistics.

Methods

The Belgian Health Interview Survey (BHIS)

The BHIS is a state-funded cross-sectional population survey carried out every several years since 1997. In 2018, the sixth BHIS has been conducted. Stratified multistage clustered sampling is used to compose the sample: a limited number of municipalities are selected in which the survey is conducted. Within each municipality, a sample of households is drawn so that groups of 50 individuals are interviewed in total. Clustering also takes place at the household level since members of the same household are more alike than persons not belonging to the same household. To ensure that the final sample reflects the composition of the population, post stratification weights are calculated, considering population data from the National Register on age, gender, residence and household size.

Study outcome and independent variables

For this study, data on smoking status (daily smoker, occasional smoker, former smoker and never smoked) and amount of cigarettes per day (categorized as “light smoking” in case of < 10 cigarettes per day, “moderate smoking” in case of 10 to 20 cigarettes per day, and “heavy smoking” for ≥ 20 cigarettes per day), age (in years), gender and educational level (no education/primary education, lower secondary education, higher secondary education, and highest education) of individuals aged 15 + were used. Statistics on age, gender and education level distribution for all Belgian municipalities were obtained from the 2011 census. The all-cause mortality data at municipality level in 2018 by age and gender were requested and obtained from Statbel.

Statistical analyses

Direct and smoothed design based estimators for municipality level smoking prevalence

Our study interest is to estimate smoking prevalence per municipality using the BHIS data. A design-based direct and smoothed estimators are discussed in this section. This study deals with a binary outcome (smoking yes vs. no) denoted as an indicator variable ${\mathrm{y}}_{\mathrm{ik}}$ for the event of interest on the $\mathrm{k}$-th individual ($\mathrm{k}=1,\dots ,{\mathrm{N}}_{\mathrm{i}}$) in the $\mathrm{i}$-th municipality ($\mathrm{i}=1,\dots ,589$). The survey is conducted with sampling probabilities for person $\mathrm{k}$ from municipality $\mathrm{i}$ being ${\uppi }_{\mathrm{ik}}$. We denote by ${\mathrm{s}}_{\mathrm{i}}$ the set of individuals who are sampled from municipality$\mathrm{i}$. The design weights are an inverse of the sampling probabilities, i.e.${\mathrm{w}}_{\mathrm{ik}}= \frac{1}{{\uppi }_{\mathrm{ik}}}$. Horvitz-Thompson proposed the following formula for an estimator of prevalence ${\mathrm{P}}^{\mathrm{HT}}$[12]:

$\mathrm P_{\mathrm i}^{\mathrm{HT}}={\sum\limits_{\mathrm{k{\in}s}_{\mathrm{i}}}\frac{{\mathrm y}_{\mathrm{ik}}}{{\mathrm\pi}_{\mathrm{ik}}}}$ 

When the number of samples in each area is large and all areas are sampled, the design-based variance is usually small and Horvitz-Thompson estimates work well. However, in the BHIS, area sample sizes are small and samples are not available from each area, therefore smoothing over the areas may be considered. Bayesian hierarchical models offer a flexible approach to fitting such models and are used in this setting [12], in particular, a so-called Besag-York-Mollie (BYM) model that assigns a Conditional Autoregressive (CAR) distribution to the random effect to account for proximity of municipalities that share a common boundary [13,14,15]. A penalizing Complexity (PC) prior, using a scaled spatially structured component and an unstructured component [15] was used for the random effects. Compared to the default Gamma prior which is commonly used in BYM model, PC priors, as a weak information prior, have been suggested as useful, understandable and conservative [15]. A sensitivity analysis with alternative parameters of the PC prior assigning more prior mass on smaller variance of the random effects were conducted to assess the robustness of results to the choice of the prior. Further details can be found in Supplementary file 1 (Technical details on design-based estimators).

To take into account the survey design, the logit transformation of the design-based direct estimates of the prevalence was modelled as $\widehat{\mathrm{p}}$_i ~ N(${\uptheta }_{\mathrm{i}}$_,$\widehat{{\mathrm{V}}_{\mathrm{i}}}$), where ${\uptheta }_{\mathrm{i}}$ is as specified above and $\widehat{{\mathrm{V}}_{\mathrm{i}}}$ is the asymptotic variance on the logit scale. A series of models including combinations of area-level covariates (age, gender and education distribution of the population at municipality level) were fitted. We estimated prevalence for two indicators, namely, current smoking vs non-current smoking (former or never), and ever smoking (current and former) vs never smoking. Prevalence and 95% credibility interval of former smokers were obtained by bootstrapping as difference between the prevalence of current and former smokers and prevalence of former smokers. Similarly, estimates of heavy, moderate, light, former and never smoking were obtained. Differences computed during bootstrapping were bounded to be positive (i.e., negative differences were set to zero and respective prevalence estimates were rescaled so that sum of categories remained equal to 100%). For the purpose of these analyses, individuals with missing data (19%, see Supplementary file 1 (Table S1) for details) on smoking status were excluded given the methods have not yet been developed to be applied to multiple imputed data.

BYM belongs to a family of Bayesian hierarchical models. Bayesian inference is useful to model the spatial and spatio-temporal data as it allows for complete flexibility in how estimates borrow strength across space and time and allow to account for similarities based on sharing borders or on the distance [13, 16]. Integrated Nested Laplace Approximation (INLA) was used to obtain estimates of model parameters. INLA has been shown as a more efficient and less time-consuming alternative to Markov chain Monte Carlo (MCMC) methods.

Fitted models were assessed and compared using the deviance information criterion (DIC) and the conditional predictive ordinates (CPO). The DIC is a criterion based on the trade-off between the fit of the data to the model and the corresponding complexity of the model. The CPO is the probability density of an observed response based on the model fit to the rest of the data. Large values indicate a better fit of the model to the data, while small values of the CPO represent unexpected response values. CPOs were plotted and inspected visually.

Smoking attributable all-cause mortality

Relative risks for all-cause mortality from smoking

Ovid Medline was searched (02/02/2021) for published systematic reviews of original studies of relative risk (RR) estimates for all-cause mortality related to smoking using filters for publication type (“systematic review”) and a combination of search terms “smoking” and “all-cause mortality”. The following information was extracted from eligible studies: author and date of publications, region/country of included studies, age and gender characteristics of study population, and RR risks for all relevant exposures.

Population attributable fraction and smoking attributable mortality

The population attributable fraction (PAF) is the proportion of cases of an outcome of interest that can be attributed to a given risk factor among the entire population [17]. When several categories of exposure to the risk factor exist, the formula takes the form:

${\mathrm{PAF}}_{\mathrm i}=\frac{\sum_{\mathrm j=1}^{\mathrm n}{\mathrm P}_{\mathrm i}\left(\mathrm{RR}-1\right)}{\sum_{\mathrm j=1}^{\mathrm n}{\mathrm P}_{\mathrm i}\mathrm{RR}},$ 

where $\mathrm{i}=1,\dots ,\mathrm{k}$ is the number of municipalities, and $\mathrm{j}=\mathrm{1,2},\dots ,\mathrm{n}$ is the number of exposure categories.

The PAF is usually expressed as the percentage of disease cases/deaths attributable to exposure. Smoking attributable mortality (SAM) can be obtained as.

${\mathrm{SAM}}_{\mathrm{i}}$ = ${\mathrm{D}}_{\mathrm{i}}$ * ${\mathrm{PAF}}_{\mathrm{i}}$,

where ${\mathrm{D}}_{\mathrm{i}}$ is the number of observed deaths in municipality $\mathrm{i}$.

Smoking attributable mortality was calculated using two scenarios. In a first scenario, we accounted for the differences in population demographic structure between municipalities by performing calculations using the number of deaths and smoking prevalence standardized to the age and gender distribution of total Belgian population. In a second scenario, given that prevalence of smoking behavior varies between age and gender groups, the calculations were performed in six strata by age and gender (age 15–39, 40–64, 65 + , males and females) and subsequently summed up to a total number of deaths per municipality that can be attributed to smoking.

In each scenario, mean estimates and 95% credibility intervals estimates were obtained by performing 1000 computations based on 1000 samples from the posterior distribution of smoking prevalence and 1000 samples from gamma distribution assumed for the relative risks. Number of computations was determined empirically, i.e. increasing number of computations from 500 to 1000 resulted in minor changes in computed estimates therefore 1000 computations were considered sufficient. Gamma distribution was chosen given its properties are compatible with the expected shape of the underlying distribution of RR: it is continuous, non-negative and puts relatively less weight on the tails of the distribution compared to e.g. log-normal or uniform distribution.

All analyses were conducted in R (R version 4.0.4, R studio version 1.4.1103). Smoothed estimates were obtained with SUMMER package (version 1.1.0) [18, 19].

Results

Descriptive statistics

The BHIS 2018 data included 9,753 respondents from 154 municipalities with 7,829 included in the analyses. The number of observations per municipality varied from 2 to 249 (Fig. 1). The mean age of participants (accounting for survey weights) was 48 years and 48% were males. Half of the respondents (50%) had high education and one third (32%) had secondary school diploma. Fifteen percent were daily smokers, 4% smoked occasionally, 23% reported to have quit smoking and 58% have never smoked. Heavy smokers constituted 2% of the total population, with another 10% being moderate smokers and 7% light smokers (Table 1). Proportion of missing data in smoking status was higher among lower educated individuals, with no differences by age or gender (23% vs 16%, see Supplementary file 1 (Table S1) for details).

Fig. 1

Belgian Health interview survey: in- and off-sample municipalities

Table 1 Survey variables –weighted estimates, based on 7,829 respondents

Design-based and smoothed municipality level prevalence of smoking

The design-based Horvitz-Thomson estimator of smoking prevalence for in-sample municipalities varied between 0 and 100% (Fig. 2 A). This large variation can be attributed in particular to small sample sizes in many municipalities. Smoothed prevalence obtained from the hierarchical model without adjustment for area-level demographic and socio-economic covariates revealed less variability in smoking prevalence between municipalities (Fig. 2 B). The DIC of this model (254.08) however revealed worse fit compared to models that accounted for socio-demographics (DIC 248.23 for a model with gender and education). Adjusting for municipality level covariates (education, gender and age) further calibrated the estimates of smoking prevalence (Fig. 2 C & D). Smoothed estimates for current smoking ranged between 11–27% per municipality, and for former smoking between 4 and 34% (Table 2). Sensitivity analyses with alternative parameters for the PC prior showed that the model was robust to the choice of priors (Supplementary file 1 (Fig. S1).

Fig. 2

Prevalence of smoking behavior (current smoker) per municipality. A Horvitz–Thompson estimator design-based estimates, B Smoothed estimates without adjustment for covariates, C Smoothed estimates adjusted for education and gender at municipality level, D Smoothed estimates adjusted for education, gender and age at municipality level

Table 2 Smoothed prevalence of current (heavy, moderate and light smokers), former smokers and persons who never smoked with 95% credible interval from model adjusted for municipality level age, gender and education

Comparing design-based estimators with smoothed prevalence estimates (for the 154 in-sample municipalities) revealed a great level of smoothing, in particular for municipalities with few observations (Supplementary file 1 Fig. S2). CPO diagnostics of fully adjusted model revealed few municipalities with very low values (Supplementary file 1 Fig. S3). Data check revealed that these were the municipalities with very high (e.g. 46% in municipality of Mechelen) or very low (e.g. 7% in municipality of Maasmechelen) design-based smoking prevalence which may be attributed to random variability due to small number of observations. These estimates were namely smoothed to 19% and 18% for Mechelen and Maasmechelen, respectively.

Smoking attributable all-cause mortality

The literature search yielded 74 papers, of which four meta-analyses provided pooled estimates of relative risk for all-cause mortality from smoking (Table 3). Smoking attributable mortality calculated according to the two outlined scenarios is presented in Table 4 and Fig. 3. Different assumptions made across the scenarios resulted in differences in number of deaths attributable to smoking indicating the uncertainty inherent to such estimates. In terms of fraction of total number of deaths, smoking was responsible for between 10 and 47% of total mortality in municipalities.

Table 3 Relative risks for all-cause mortality from smoking

Table 4 Smoking attributable all-cause mortality per municipality, in absolute number of deaths

Fig. 3

Smoking attributable mortality (as absolute number of deaths and proportion of total deaths) according to the two scenarios

Discussion

The objective of this study was to apply the Bayesian hierarchical model to obtain indirect small area estimations of smoking prevalence in Belgium and compute municipality-level smoking attributable mortality. The final estimates of current/former smoking prevalence ranged between 11 [Credible Interval 3;23] and 27% [21;34] between the municipalities. Across Belgian municipalities, smoking was assessed to be the cause of 10 [7;15] to 47% [34;59] of all deaths. This demonstrates that smoking poses a serious public health burden and, given that smoking-related deaths are avoidable, it also presents a prevention opportunity to improve quality and duration of life of many people.

To our knowledge, this is the first study to compute municipality level estimates for smoking and smoking attributable mortality in Belgium. Information about spatial distribution of these health indicators should support local policy-making as well as contribute to explaining the regional differences in mortality. The methodology used in the paper can be broadly applicable in the context of national health surveys in other countries.

While a number of small-area estimation methods have been proposed in the literature, it was out of scope of this work to provide a comprehensive comparative analysis of these methods and we rather focused on practical application of the method proposed by Mercer et al. [12]. This approach had the advantage in a sense that it is suitable for complex survey design with large amounts of off-sample areas. A crucial question is the validity of the estimates in the off-sample areas as well as in areas with very few respondents.

With respect to off-sample areas and in an absence of the true estimates in each municipality it is hard to make strong claims about the validity of the produced estimates. Existing methodological literature including simulation studies suggests this method yields valid estimations [12], however, it is reasonable to expect for its accuracy to decrease with increasing number of off-sample areas. In the current setting, data were available from 154 of 589 municipalities, which implied a lot of smoothing and extrapolation. To our knowledge, no simulation studies are available to assess the performance of this method in context where almost three quarters of the areas were off-sample (Watjou et al. have omitted up to 56% of the areas in their simulation study [10]). To improve the accuracy of estimations, auxiliary data from three municipality-level socio-economic covariates (age, gender and education) were used. These factors have been shown to be associated with smoking in the BHIS data but also from previous research [24], and thus, accounting for distribution of municipality population across age, gender and educational groups was expected to further fine-tune the estimates. In fact, only minor improvements in the fitted models were observed according to DIC which may be explained by ecological fallacy – area-level covariates may fail to pick up the existing associations with the outcome on individual level. Incorporating individual level covariates in the model would likely improve the predictions; however, more methodological research is needed to accommodate challenges related to complex survey data (e.g. how to incorporate the individual weights when part of areas is off-sample) and missing data (e.g. missing outcome or covariate data in survey data and in off-sample areas) in this context. In current analyses, provided estimates of smoking prevalence are likely conservative given lower educated respondents (how are known to be more frequently smokers) were more likely to leave the question on smoking blank. For missing survey data multiple imputation (MI) could provide a remedy under assumption of ‘missing at random’ (MAR) [25]. Smoothed estimates method should be further developed to accommodate the theoretical and practical implementation of MI in this context. Last but not least, given the increasing policy relevance of small area estimates, future rounds of the BHIS survey may reconsider the tradeoff between increasing the number of in-sample areas, which improve possibilities of small area analyses, and the cost implications of such design changes.

With respect to areas with very few respondents, model diagnostics revealed substantial smoothing towards the mean prevalence. Where direct smoking prevalence estimates (based on 2 to 249 respondents) ranged between – likely unplausible – 0 and 100%, the smoothed indirect estimates had a range of 11 to 27%. This is subject to large uncertainty if the sample size is very small. On the other hand, omitting important covariates can result in over smoothing, so obtained estimates should be interpreted with caution and validated in future studies.

Another important methodological step in this work was calculating the population attributable fractions that are subsequently used to obtain smoking attributable mortality. These widely used in public health measures are subject to strong limitations. First, these calculations rely on assessment of relative risk of smoking for all-cause mortality from existing studies across the world or indirect methods such as Peto-Lopez method that bases relative risks on the excess rates of lung cancer mortality [26]. For this study, we have used the estimates from four meta-analyses of existing studies [20,21,22,23] as meta-analyses provide more robust and stable estimates of relative risks compared to single studies. However, the obtained estimates did not cover in sufficient detail all age and gender groups of general adult population. Analyses under the two scenarios (standardizing the population for age and gender vs summing up estimates from six strata's by age and gender) resulted in a wide range of estimates with wide confidence intervals reflecting substantial uncertainty (Table 4). Second, while relative risks were obtained after adjusting for several important confounders, the calculations of population attributable fraction may still fail to account for a number of possible exposures such as, for example, air quality or work-related exposure. Third, important considerations were made in respect of standardization of number of deaths per municipality to age and gender structure of general Belgium population. On the one hand, estimates of relative risks may have already captured the age and gender distribution in the study cohort (which may or may not be comparable to Belgian population), on the other hand, obtaining PAF estimates on standardized or stratified number of death per age and gender group was supposed to explicitly eliminate differences in number of deaths between relatively ‘green’ (with large proportion of young population) and ‘grey’ (with large proportion of older population) municipalities. To further improve the accuracy of PAF calculations, data on age and gender specific relative risks of all-cause mortality from smoking in Belgian population are needed. Last but not least, the estimates of smoking prevalence and death statistics from 2018 were used, while there is a known large lag between the onset of smoking, length of exposure to smoking and death. Given that smoking rates have been declining substantially over the last 20 years [6], provided estimates of smoking attributable mortality are likely to be conservative as deaths occurring in 2018 are related to past smoking.

In conclusion, substantial within-country variation in smoking and subsequently smoking attributable mortality was observed. This study provided an illustration of the methodological approach to obtain small area estimates from national health survey that can be used for a broad range of health indicators.