Bayesian spatial modelling for quasi-experimental designs: An interrupted time series study of the opening of Municipal Waste Incinerators in relation to infant mortality and sex ratio

Background: There is limited evidence on potential health risks from Municipal Waste Incinerators (MWIs), and previous studies on birth outcomes show inconsistent results. Here, we evaluate whether the opening of MWIs is associated with infant mortality and sex ratio in the surrounding areas, extending the Interrupted Time Series (ITS) methodological approach to account for spatial dependencies at the small area level. Methods: We specified a Bayesian hierarchical model to investigate the annual risks of infant mortality and sexratio (female relative to male) within 10 km of eight MWIs in England and Wales, during the period 1996–2012. We included comparative areas matched one-to-one of similar size and area characteristics. Results: During the study period, infant mortality rates decreased overall by 2.5% per year in England. The opening of an incinerator in the MWI area was associated with −8 deaths per 100,000 infants (95% CI −62, 40) and with a difference in sex ratio of −0.004 (95% CI −0.02, 0.01), comparing the period after opening with that before, corrected for before-after trends in the comparator areas. Conclusion: Our method is suitable for the analysis of quasi-experimental time series studies in the presence of spatial structure and when there are global time trends in the outcome variable. Based on our approach, we do not find evidence of an association of MWI opening with changes in risks of infant mortality or sex ratio in comparison with control areas.


Introduction
In recent years, incineration of municipal waste in Europe has increased in response to legislation to divert waste from landfills, leading to the construction of new municipal waste incinerator (MWI) plants. Emissions from MWIs include particulates, SO 2 , nitrogen oxides (NO X ), heavy metals, polycyclic aromatic hydrocarbon (PAH) and polychlorinated biphenyls (PCBs), as well as polychlorinated dibenzo-pdioxins and furans (PCDD/Fs) which are reported to have a range of biological effects including enzyme induction, immunotoxicity, developmental toxicity and with reproductive and endocrine effects (Birnbaum, 1994). Exposures to particulates as a proxy for MWI emissions are very low from modern well-regulated incinerators (Douglas et al., 2017).
Despite high public concern about potential effects of waste management processes on birth and infant mortality outcomes, there are relatively few epidemiological studies available, particularly of modern incinerators. Some previous studies have found associations with infant mortality (Tango et al., 2004), miscarriages (Candela et al., 2015), sex ratio (Williams et al., 1992), preterm delivery (Candela et al., 2013), and congenital anomalies (Dummer et al., 2003;Cordier et al., 2010), with congenital anomalies and pre-term birth being the leading causes of infant mortality. The established reproductive toxicity of dioxins and heavy metals, present in emissions from incinerators, provides plausible biological mechanisms for these observed associations, although doses received by individuals at ground level are likely to be very low (Ashworth et al., 2014). In addition, sex ratio may behave as a sentinel for the endocrine disruption due to chemicals in MWI emissions including dioxins and furans (Candela et al., 2013;Van Larebeke et al., 2008).
A systematic review (Ashworth et al., 2014) found that the evidence for adverse birth and neonatal outcomes in relation to waste incineration was limited and inconclusive. Problems identified were the small numbers of available studies and low statistical power in each study, potential exposure misclassification, and residual confounding. Additionally, the current evidence-base was found to be difficult to synthesise due to the heterogeneity between studies in their design, location, time period, number of incinerators and the birth outcomes explored. Thus, the authors recommended further epidemiological investigation, particularly multi-site studies to increase statistical power, with better control for confounding, more accurate exposure assessment and robust methodology (Ashworth et al., 2014).
To help address these recommendations, we consider the opening of multiple new MWIs across England as a natural experiment to study annual trends in birth outcomes, with a focus on infant mortality and sex ratio, using an interrupted time series (ITS) approach (Zhang et al., 2014).
In an ITS study, a time series is used to establish an underlying trend, which is 'interrupted' by an intervention at a known point in time. Such methods have been applied in: epidemiology (Gasparrini et al., 2009;James Lopez, 2017), econometrics (Shadish et al., 2002), drug use research (Jandoc et al., 2015), evaluation of health care quality (Penfold and Zhang, 2013) and to assess impact of changes in public health policies (Kontopantelis et al., 2015) e.g. the effect of a smoking ban on birth outcomes (Bakolis et al., 2016). The ITS approach has also become popular in the road safety literature, where authors have focused on deriving indices to describe the effect of safety interventions both in an empirical Bayesian and fully Bayesian framework (Park et al., 2010;El-Basyouny and Sayed, 2011).
In its original formulation, ITS evaluates only the changes in the time series of the outcome of interest before and after the intervention, failing to account for time trends independent of it; therefore inference may be biased if the outcome exhibits a time trend that might confound the intervention effect (Linden and Adams John, 2010). Li et al. (2013) proposed a two-stage Bayesian hierarchical model to evaluate the effect of an intervention to reduce crime, allowing for data sparsity and including control areas to account for global time trends. In the first stage, the control areas (matched to the exposed areas) were used to assess the pre-intervention trend and then the estimates were fed to a second stage, where the effect of the intervention was evaluated in the exposed areas. However, the two-stage formulation prevented the propagation of uncertainty from exposed to control areas; also, no structured local spatial dependency was included.
Here, we develop the approach used by Li et al. (2013) and propose a framework for applying the ITS approach to small area data, defining the opening of a new MWI as the intervention, and infant mortality and sex ratio in the surrounding areas as the outcome.
We consider important confounders including area-level measures of ethnicity and deprivation, both known to be associated with early life outcomes. For example, a systematic review of UK studies found that infant mortality is associated with greater area-deprivation and lower social class (Weightman et al., 2012). Infants of South Asian and Black ethnicity are reported to have higher infant mortality, partially explained by the increased prevalence of congenital anomalies in these populations (Li et al., 2018). In addition, sex ratio (males over total) declines in association with socio economic status, possibly because of higher stress levels with lower socio-economic status (Grech, 2017), and differences have been reported in newborn sex ratio according to paternal ethnicity (Grunebaum and Chervenak, 2016).
Our modelling framework extends the previous statistical work in the field in the following ways: (i) we include control areas matched on size and area characteristics to allow for global trends; (ii) we specify the spatial dependency between areas; (iii) we perform joint Bayesian inference of exposed and control areas which allows uncertainty to be propagated across the model; and (iv) we define an index to quantify the before-after difference in infant mortality risks comparing the MWI and control areas.
This work is part of a larger national investigation of 22 MWIs in relation to reproductive outcomes in Great Britain during 2003-10, where we estimated exposure to particulate matter from MWI emissions in Great Britain (Douglas et al., 2017) and further used both modelled emissions concentrations and distance from incinerators to examine associations with fetal outcomes (Ghosh et al., 2019).

Study design
We included all eight new MWIs coming into operation in England and Wales between 2003 and 2010 operating under the EU WID. Infant mortality and sex ratio data were available for the period 1996-2012, at Middle layer Super Output Area (MSOA) Census geography (average population 7500).

Incinerator (exposed) areas
We defined the exposed area for each incinerator as the set of MSOAs whose centroids lie within 10 km of a MWI (Fig. 1). We chose 10 km as the cut-off for consistency with screening criteria used for implementing the Habitats Regulations (1992/43/EEC). These regulations state that incineration plants that are within 10 km of a European Site require an assessment of their impact for short-range air emissions (Ashworth et al., 2013). For each MWI, we considered five years before and after the opening date (Table 1), except for Isle of Wight (3 years) and Grundon (Lakeside) (2 years) where opening dates were respectively 2009 and 2010, because as noted outcomes were available only to 2012.

Comparator areas
For each MWI we selected one comparator area (composed by MSOAs), as follows: (i) For each MWI, we identified the potential set of comparator MSOAs, as all MSOAs within the same Government Office Region as that of the MWI. For Grundon more than half of the MSOAs were in London (the rest in South East) therefore we used the London Region as reference (Table 1). (ii) For each MSOA in the reference region, we then constructed a 10 km circular buffer around each MSOA centroid. We excluded MSOAs within 10 km of an existing or new MWI or sharing a boundary with such an MSOA. (iii) For each of these potential comparator areas, we selected the same opening years as the MWI, and then constructed a vector that included: • annual number of live singleton births for the pre-opening period, • population density (Census 2011), • percentage of population of non-white ethnicity (Census 2011) • Carstairs deprivation index input variables: proportion of male unemployment, proportion of households without a car, proportion of overcrowded households, proportion of persons in private households with an economically active head of household in social class IV or V (partially skilled and unskilled occupations) (Carstairs and Morris, 1990). (iv) We computed the Euclidean multivariate distance between the MWI and comparator vectors and selected the comparator area with the smallest distance.
For the five incinerators located in the South East, we computed the comparator area for each of them independently; 28% of MSOAs were shared between control areas (see Fig. 1). We depict the steps for the choice of the comparator area (Supplementary material Fig. 1).

Outcome and confounder data
Data on infant mortality and sex ratio were obtained from the Office for National Statistics; counts of annual infant mortality were defined as the sum of deaths occurring between 0 and 365 days of life. Sex ratio is defined as the sum of the female births divided by sum of male births at area level.
We accounted for the following potential MSOA-level confounding variables: • Area-level socio-economic deprivation defined using the Carstairs 2011 score in quintiles (Carstairs and Morris, 1990) • Population density from the 2011 Census • Ethnicity (2011 Census) as the percentage of the MSOA non-white female population ages 16-49 years old, split in three categories (Black, Asian and Other ethnicity), with cutoffs defined using the England and Wales average (12.1%) and twice the average (24.2%)  • NO 2 concentrations at 200 m × 200 m resolution, postcode headcount weighted to MSOA level .
In sensitivity analyses, we investigated the use of MSOA centroids within 4 km of the MWI, for consistency with previous published studies on incinerators (Candela et al., 2015;Candela et al., 2013). We also restricted the time window to include only those MWIs with data available for 5 years before/after the opening. Further, we allowed area-specific intercepts for all the MWI and comparator areas to account for possible heterogeneity of the MWIs.

Bayesian interrupted time series
We used a Bayesian hierarchical model to jointly fit the data for MWI (exposed) and comparator areas. We defined y ikt as the counts of infant deaths, distributed as:  where i = 1, … , N MSOAs, t is time (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and k = 1, … , K MWI and E ikt is the number of live singleton births. Then on the link function log(λ ikt ) we specified a linear regression: where D i indicates if the MSOA centroid lies within 10 km of an MWI (D = 1) or a comparator area (D = 0); t 0k t t (1 ) k 0 is the opening year for the k − th MWI and Z i is the vector representing potential confounders at the MSOA level.
The regression coefficients are normally distributed centred on zero, with a variance equal to 10 4 . The area-specific random effects b i are here defined as exchangeable, i.e. b i~N (0, σ b 2 ), with the hyperparameter 1/σ b 2~G amma(1,0.00005) (Martins et al., 2013;Rue et al., 2009). We extended the specification below to allow for local spatial dependency.
In the Eq. (1) α indicates the log baseline risk for infant deaths; β 1 represents the difference in risk between MWI and comparator areas before the incinerator opened; β 2 indicates the general time trend, which we assumed was linear in this setting, but could be extended to be non-linear if appropriate; β 3 identifies a potential change in the timetrend in the comparator areas after the incinerator opening and β 4 is the coefficient describing the change in linear trend between MWI and comparator areas before the incinerator opened. Finally, β 5 indicates the change in time trend between MWI and comparator areas beforeafter the MWI opened. We present the model for count data using a Poisson distribution; for the log transformed sex ratio the distribution is Gaussian but the interpretation of the parameters remains the same. We extend the modelling framework in Eq. (1) to allow for spatially structured random effects (Freni-Sterrantino et al., 2018;Riebler et al., 2016;Simpson et al., 2017), appropriate when there is localised spatial dependence (Supplementary materials Appendix A).

Calculation of summary index
We constructed a summary index to describe the difference in before-after effects between MWI and comparator areas, for infant mortality: where μ CB |y, μ CA |y, μ EB |y, μ EA |y (E = Exposed, C = Comparator, A = After, B = Before) are the posterior distributions obtained from Eq.
(1) (see Supplementary material Table 1) and n EA is the total number of deaths occurred in MWI areas after the opening. The µ EB µ µ CA CB gives an estimate of the effect of MWI opening adjusted for the comparator area trends. We report γ as the difference in deaths per 100,000 infants and as ratio difference for sex ratio.

Results
Overall from 1996 to 2012 there were 526,642 single live births and 2462 infant deaths recorded across the regions included in the study, across 304 and 276 MSOAs in MWI and comparator areas, respectively.
Descriptive statistics are shown in Table 2. Comparing periods before and after the opening of an MWI, infant mortality rates decreased from 4.8 to 3.9 (per 1000) in the MWI areas, 5.0 to 4.6 (per 1000) in the comparator areas. The inter-quartile ranges (IQRs) of infant mortality were higher before opening in the MWI and comparator areas (9.62 and 9.43 per 1000 respectively), compared with after opening (IQR 7.93 and to 9.01 per 1000 respectively). We observed a similar trend for England and Wales over the same period; rates decreased from 6.0 to 4.1 (per 1000), rates decreased overall by 2.5% per year (Office for National Statistics).
For sex ratio, there was no difference between before and after MWI opening; in the exposed areas, average sex ratio was 0.96 (IQR 0.28) and 0.97 (1QR 0.26) respectively and it was 0.96 (IQR 0.25) and 0.97 (1QR 0.26) respectively in the comparator areas. Table 3 shows the results of the regression model. The summary index γ indicates there were 8 fewer infant deaths per 100,000 live, single births (95% CI −62, 40) when comparing the before and after periods, adjusted for the global trend in the comparator areas. 27% of the variability in infant mortality was accounted for by the spatially structured random effect b (the posterior mean of the mixing parameter ϕ), and this was robust to selection of priors. The posterior estimates show that infant mortality was higher by 48% (95% CI 24, 76) for the more deprived areas and by 27% (95% CI 10, 47) for areas where the proportion of non-white ethnicity was above the national average.
For sex ratio, the summary index γ as ratio difference is −0.004 (95%CI −0.02, 0.01), and 46% of variability observed was accounted for by a spatially structured random effect. Fig. 2 depicts the posterior mean of the linear predictors and credible intervals for the MWI and comparator areas. The lines are almost parallel, and the credible intervals overlap, indicating that the decrease in infant mortality risk is similar for both areas as is sex ratio.
By restricting the MWI boundary to 4 km in sensitivity analyses, there were 569 deaths among 113,940 live births across 59 MSOAs and 63 comparator MSOAs. We observed decreasing rates in both areas, from 5.4 to 4.9 (per 1000) before/after MWI opening for MWI areas and 5.0 to 4.5 (per 1000) in the comparator areas (Supplementary material  Tables 3, 4, 5). Model estimates show the same spatial variability as observed for the 10 km analyses and a summary index γ of 6 deaths per 100,000 infants (95% CI −126, 126) in the MWI areas corrected for trends in the comparator areas. Mean sex ratio between MWI and comparator areas was the same as for the 10 km buffer and the summary index γ is 0.009 (95% CI −0.03, 0.07); spatial structured variability accounted for 21% of the variance.
Sensitivity analysis for 4 and 10 km buffers restricted to incinerators with data five years before and after the opening, did not alter the findings (Supplementary material Tables 6 and 7).

Discussion
We present a statistical framework that extends the ITS model to account for spatial dependence applied to the opening of MWIs. We applied the developed framework to evaluate the effect of the opening of a MWI on infant mortality and sex ratio at birth, and found no difference when comparing MWI areas with comparator areas. For infant mortality we accounted for the strong decreasing global trends seen across England and Wales over the study period.
A strength of our study is the use of all available data on MWIs and infant mortality and sex ratio from 1996 to 2012, for England and Wales, to investigate possible effects of MWI opening on these outcomes, hence avoiding selection effects (bias). Few studies to date have published results on infant mortality in relation to MWIs. Tango et al. (2004) reported no significant excess of infant deaths near 63 MWIs in Japan, but found a small statistically significant peak decline in risk with distance up to 10 km from incinerators. In addition, Dummer et al. investigated stillbirths and neonatal deaths and found no associations with incineration (MWI and crematorium) (Dummer et al., 2003).
In Ghosh et al. (2019) we previously reported neonatal and postneonatal mortality for 2003-10 within 10 km of MWIs and found no association related to modelled emissions or distance. Here we use an alternative approach to investigate the association of infant mortality with incinerators, taking advantage of time series to investigate potential 'before/after' effects.
For sex ratio, Williams et al. (1992) reported a significant excess in female births in one of three "at risk" districts compared to the Scottish average, at aggregated data level. No significant results were found with modelled incinerator generated dioxins since 1992 for one incinerator located in Taipei (Taiwan) (Lin et al., 2006); also no significant excess was found near 63 MWIs in Japan (Tango et al., 2004). In Italy Candela et al. (2013) evaluated the effect of 8 MWIs in the Emila Romagna Region using dispersion models to estimate particulate matter (PM10) exposure, and found no significant results for sex ratio (Candela et al., 2013). Similarly, in UK, we previously reported no significant association with sex ratio (Ghosh et al., 2019), but using a dispersion model and distance in relation to 22 MWI, and not the ITS approach.
The ITS framework is well suited for policy evaluation of time-series studies in a quasi-experimental setting. Most previous studies using the ITS approach have not allowed for the possibility of analysing global trends, which is critical to interpretation of any before-after effects as these are non-randomised interventions. It offers advantages over difference in difference methods because it not only controls implicitly for baseline mean and trends but provides a model to deal with autocorrelation present in the data.
We aimed to maximize the similarity between exposed and comparator areas by including time-varying variables (i.e. live births) in the period before opening as well as socio-economic deprivation indicators. We investigated the inclusion of the MSOA rural/urban classification and NO 2 (a proxy for road traffic) as additional matching variables, but we found no association on the selected comparator areas. The fact that there was only a small posterior mean difference (and credible interval including zero) for β 5 (the parameter that indicates a potential   difference in the trend before-after the incinerator opening, between MWI and comparator areas), suggests not only that there is no effect from the incinerator but also that the matching has been done appropriately. Nonetheless, despite having matched areas for socio-economic deprivation, we included spatial structured random effects to account for any potentially unmeasured confounders at small area level and an area-specific intercept to account for differences in MWI characteristics. The inclusion of a structured spatial effect resulted in a smaller Deviance Information Criterion (model selection index) compared to a spatial unstructured random effect (data not shown) and should improve inference in the ITS approach (Ning et al., 2019). Area-specific intercept did not improve the DIC, as the intercepts had similar values. The difference before-after opening is represented as a step or shift in the intercept and slope in the regression models. We found no evidence of temporal autocorrelation or seasonality in the residuals of our model, suggesting that a linear time effect was appropriate in our data, although the proposed framework is extremely flexible and can deal with non-linearity via semi or non-parametric functions (e.g. splines or Gaussian processes) if needed.
Our approach has limitations. The official opening date for each MWI was provided by the Environment Agency but may not have coincided with the actual operating date. Douglas et al. (2017), showed that MWIs were not all operating continuously over this time period, due to repairs, upgrades, and closures. This might have resulted in some misclassification of the before-after opening period.
Other limitations include lack of gestational age for the time period 1996-2006 and lack of knowledge of maternal address before birth, as 24.4% women move during pregnancy (Hodgson et al., 2015), which meant that we could not analyse possible exposures during pregnancy, as an alternative exposure window. Furthermore, the confounders used in our analyses were obtained from the Census 2011, and we assigned them equal weight before/after MWI opening in the matching process as they were not time varying. Although the population characteristics are likely to remain similar over this relatively short period, we cannot exclude possible residual confounding.
Distance is a relatively crude measure of exposure to an MWI, but it may capture exposures not reflected in the stack emissions modelling (e.g. transport of waste to the MWI) as well as non-airborne exposure pathways. We selected a 10 km boundary around each MWI, as has been used in previous studies in the UK (Reeve et al., 2013;Elliott et al., 1996;Gouveia and Prado, 2010;Ghosh et al., 2019), also including the dispersion modelling study of MWI emissions carried out by Douglas et al. (2017) for all 22 British MWI between 2003 and 2010. In a sensitivity analysis, we also selected a 4 km boundary based on previous studies conducted in Italy (Candela et al., 2013;Candela et al., 2015;Ranzi et al., 2011), with results similar to those we observed for 10 km. Smaller buffer size of 1-3 km were not considered because of limited statistical power.
In conclusion, we did not find an association between the opening of a new MWI and changes in infant mortality trends or sex ratio at birth for 10 and 4 km buffers, using distance as proxy of exposure, after taking into account temporal trends in comparator areas and potential confounding factors.
The framework we propose, applying ITS to small areas, can readily be adapted to take account of different types of data, structures, spatial and temporal dependencies. The method may be relevant to applications that investigate a potential causal effect of an intervention or policy, taking advantage of the quasi-experimental design.

Conflict of interest
Anna Hansell declares a Greenpeace membership but has not received any money from the organisation nor has been involved in campaign, nor other relationships or activities that could appear to have influenced the submitted work. All other authors declare no conflict of interest.

Source of funding
The study was funded by a grant from Public Health England (PHE), by a grant from the Scottish Government, and funding from the MRC-PHE Centre for Environment and Health and from the National Institute for Health Research Health Protection Research Unit (NIHR HPRU) in Health Impact of Environmental Hazards at King's College London and Imperial College London in partnership with Public Health England (PHE). The work of the UK Small Area Health Statistics Unit is funded by Public Health England as part of the MRC-PHE Centre for Environment and Health, also funded by the UK Medical Research Council (MR/L01341X/1). PE is Director of the MRC-PHE Centre for Environment and Health and acknowledges support from the NIHR Imperial Biomedical Research Centre. This work used the computing resources of the UK MEDical BIOinformatics partnership -aggregation, integration, visualisation and analysis of large, complex data (UK MED-BIO) which is supported by the Medical Research Council (MR/ L01632X/1).
The funders had no role in study design, analyses, interpretation of the data, or decision to submit results.

Contributors
AFS is the guarantor of this paper. All authors had full access to all of the data in the study and can take responsibility for the integrity of the data and the accuracy of the data analysis.

Data and computing code
Births and deaths data were obtained from the Office for National Statistics (ONS) National Mortality, Births and Stillbirth registers for England and Wales.
2011 Census aggregate data came from the Office for National Statistics (for England and Wales): UK Data Service (Edition: June 2016).
Health data are available from the data providers on application with appropriate ethics and governance permissions, but we do not hold data provider, ethics, or governance permissions to share these datasets with third parties. The computing code is available upon request to the corresponding author.