Winter Is Coming: A Socio-Environmental Monitoring and Spatiotemporal Modelling Approach for Better Understanding a Respiratory Disease

Abstract

Chronic Obstructive Pulmonary Disease is a progressive lung disease affecting the respiratory function of every sixth New Zealander and over 300 million people worldwide. In this paper, we explored how the combination of social, demographical and environmental conditions (represented by increased winter air pollution) affected hospital admissions due to COPD in an urban area of Christchurch (NZ). We juxtaposed the hospitalisation data with dynamic air pollution data and census data to investigate the spatiotemporal patterns of hospital admissions. Spatial analysis identified high-risk health hot spots both overall and season specific, exhibiting higher rates in winter months not solely due to air pollution, but rather as a result of its combination with other factors that initiate deterioration of breathing, increasing impairments and lead to the hospitalisation of COPD patients. From this we found that socioeconomic deprivation and air pollution, followed by the age and ethnicity structure contribute the most to the increased winter hospital admissions. This research shows the continued importance of including both individual (composition) and area level (composition) factors when examining and analysing disease patterns.

1. Introduction

This research focuses on people suffering Chronic Obstructive Pulmonary Disease (COPD), a progressive lung disease affecting the respiratory function. In New Zealand (NZ) one in fifteen people suffer from COPD, creating a significant healthcare burden [1]. While cigarette smoking is usually the main factor causing COPD, long-term exposure to environmental irritants, such as air pollution, chemical fumes, or dust can also initiate and exacerbate the disease [2]. The effect of both chronic and acute air pollution exposure on human health has been broadly discussed in the scientific literature [3,4]. While the relationship between air pollution and COPD is not entirely understood, the majority of studies indicate air pollution as an important cause of the disease [5,6,7] or at least as a source of disease exacerbation (often in the combination with other causes) leading to increased hospitalisation rates and lung function breathing impairment [8,9]. Exposure of the population to air pollution is not the only causal environmental factor of COPD. Literature links increased COPD hospitalisation rates with colder temperatures and seasonality [10,11,12] as well as with urban environment [13]. Lifestyle, demographic, and socioeconomic factors, such as smoking, deprivation, and living conditions [14]; age, gender, and education [15]; ethnicity [9,16]; and the occupational exposure to lung irritants [15,17,18] have also been identified as causing or exacerbating COPD. Previous studies including in New Zealand have found that inequality in air pollution exposure often has a socioeconomic component [19,20,21,22]; while the housing quality is also a significant factor in poor respiratory health [23,24].

The study took place in the urban area of Christchurch (NZ) (43.53° S, 172.64° E) that is a regional capital of Canterbury. It is the largest city located in the New Zealand’s South Island and second most populous city in the country (381,500) [25]. The city is located on the east coast of the South Island, and although in the vicinity of Southern Alps, it is mostly flat except its southern part that is dominated by the Port Hills of Banks Peninsula. The city also lies near to active fault lines that make it vulnerable to earthquakes, e.g., the infamous series of earthquake events damaging the city in 2010 and 2011.

The city suffers a poor air quality during winter (June to August) with levels exceeding WHO guidelines for PM_2.5 (25 μg·m⁻³ 24-h average) [26]. The post-2010/2011 earthquakes and subsequent rebuild has arguably both exacerbated poor air quality, but also have provided a city with an opportunity to become smarter by incorporating environmental sensors into its rebuilt infrastructure. The air quality monitoring network that was established during the rebuild has provided high-precision, real-time environmental data allowing more detailed and fine-grained monitoring and environmental exposure modelling. This enhancement helps to improve the understanding of processes that lead to worse air quality and its subsequent health impacts.

This paper aims to investigate whether there are any existing associations between the number of COPD hospitalisations and air pollution conditions in the city of Christchurch in winter compared to the non-winter seasons, and identify how social and demographical factors might alter this relationship. We focused on the spatial, temporal dimension, and also on the place and people’s characteristics to provide overview picture of the situation. To do this we analysed aggregated health records, air quality monitoring data as well as socio-environmental and census data associated with COPD hospital admissions and identified populations most vulnerable to COPD, and COPD vulnerability ‘hot spots’ within the city. This procedure enables a more personalised approach to healthcare, advising patients, and vulnerable groups based on the latest environmental information; and longer term, helping to improve their quality of life and increase their ability to adjust behaviour to minimise COPD exacerbation.

2. Materials and Methods

2.1. Data

The primary dataset we analysed consisted of hospitalisation admissions records for the period August 2014 to October 2017 provided by the Canterbury District Health Board (CDHB). District Health Boards (DHBs) are government funded organisations responsible for ensuring the provision of health and disability services to populations within a defined geographical area, in this case Canterbury, within which Christchurch is located. Each record contained information about the ethnicity, age, gender, smoking, date of admission and also the address (if recorded) of individual patients hospitalised due to the acute COPD exacerbations. There were 3893 cases of hospitalisation documented during the study period of which 2350 patients had their permanent residence in the Christchurch study area, and so eligible for inclusion in the subsequent analysis. The majority of those hospitalised belonged to the European ethnic group (85.7%) or followed by Māori (the first indigenous people of New Zealand) (10.5%); females (55%) were hospitalised more often than males (45%), the only exception to this was observed in the Pacific Islanders ethnic group where males were hospitalised more frequently than females. The age structure of those hospitalised (Figure 1) shows that the increase in the number of hospitalisations is apparent for both males and females older than 50 years, as COPD is a disease affecting middle-aged and elderly persons [27,28]. June to August are designated winter months in the southern hemisphere, however, there was also an increase in the number of hospital admissions observed during the September as depicted (Figure 1b), therefore, all four months were considered winter for the purpose of the study. There were almost as many hospitalisations in those four months (1801) compared with the other eight months of the year (2092) that suggest the distinction of the winter season that requires a deeper understanding. There is also a deterioration in air quality in the winter months due to winter wood burning associated with home heating. The COPD hospital admissions data were aggregated into the Census Area Units (CAU) geography based on the residential addresses of patients. CAUs are census administrative areas that, in urban areas, contain 3000 to 5000 people. There were 133 CAUs in the study area in total (see Figure A1 in Appendix A for more details). The geographic distribution of season-specific hospital admissions due to COPD is depicted in Figure A2 located in the Appendix A.

The hospital admission data contains the information about the current smoking behaviour of the patient, which is important since smoking is usually considered the most important cause of COPD exacerbation. However, only 2.53% of people (2.26% of men/2.64% of women) hospitalised for COPD were coded at their admission as being current smokers, however, we believe this is likely to be an underestimate.

Looking at absolute numbers of hospital admissions between 2014 and 2016, the typical patient would be characterised as a female of European ethnicity over 60 years old and hospitalised during winter months.

During the study period, the 3893 hospitalisations due to COPD consisted of 2311 individual patients of which 1512 had one hospitalisation, therefore nearly one third of those in the dataset had multiple hospitalisations, with some individuals hospitalised even more than twenty times during the study period. Due to the focus on the hospitalisation rates and not only on disease incidence, multiple hospitalisations of one patient were treated as distinct events.

Demographic data used in the study came from the New Zealand Census in 2013, namely; age, gender and ethnic structure by CAU. The Census data were also used to calculate age and sex-adjusted incidence rates and Standardized Incidence Ratios (SIR) (based on the internally indirect standardized expected number of disease cases [29]) that served as the primary measure of the COPD hospital admissions in individual CAUs, allowing for the comparison among aggregated units overall, and between winter and non-winter months. The New Zealand Socioeconomic Deprivation Index NZDep2013 [30] was used as a proxy measure of the socio-economic status and quality of the neighbourhood represented by CAU.

A dense network of sensors was used to estimate average annual and season-specific PM_2.5 exposure in the study area. It consists of more than twenty high-quality sensors located in a regular grid (Figure A1 in Appendix A), supplemented and further validated by monitoring sites managed by local official agencies that provide data to the public. Air pollution data, specifically concentration of PM_2.5 (particulate matter smaller than 2.5 µm), came from two sources; the official monitoring network operated by Environment Canterbury (ECan) that is the regional council whose responsibilities include environmental management, and provided data from three sites—St Albans, Woolston and Kaiapoi (see the Figure A1 in Appendix A for the location of individual suburbs), and from the temporary experimental network of twenty additional monitoring sites aiming to densify the ECan network in order to improve the understanding of the air pollution and its behaviour within the city at finer spatial and temporal scales [31]. More details about the monitoring network are mentioned in [31,32]. While the network provided data at one-minute intervals, the data paired with either individual COPD cases or CAUs were aggregated to obtain annual averages, and winter and non-winter averages aggregated for CAUs (Figure 2). They were calculated as averages of 10-minute air pollution surfaces created using spatiotemporal kriging [33,34]. Specifically, a metric spatiotemporal variogram model was utilised (nugget = 0, partial sill = 1.5, range = 1200, space-time anisotropy = 15). The average PM_2.5 concentrations for CAUs are depicted in Figure 3. The differences between annual average (Figure 2a), winter season (Figure 2b) and other months (Figure 2c) are easily observable with little air pollution during the non-winter months and increased levels of winter PM_2.5 most noticeable in the eastern part of the study area.

2.2. Spatial and Spatiotemporal Patterns

We used global Getis-Ord G-statistic and local G_i* [35] to quantify global and local spatial autocorrelation in the data and to explore if any spatial patterns (hot and cold spots) of SIRs exist in the study area. The G and G_i* statistics are generally used as an indicator of local clustering that measures a concentration of spatial data [36]. The local G_i* allow for detecting and subsequent visualisation of the location of identified clusters of dependence that may not show up when using only a global statistics [35]. The method’s results distinguish between clusters concentrating high or low values, in the case of this study the clusters of high and low COPD hospitalisation rates (unhealthy and healthy clusters). The second order queen’s contiguity scheme of the neighbourhood was used to obtain spatial weights for individual CAUs, while the spdep package for R [37] was used for the processing of the data and computing of the measures.

The spatio-temporal scan statistics aimed to identify clusters of high and low rate areas over time and space and was processed using SaTScan 9.3 software [38] via the rsatscan package for R [39]. Input data consists of (1) hospitalisations aggregated into CAUs and grouped by week of admission, age, and gender of patients; (2) demographic structure of CAUs (age and gender); (3) coordinates of CAU centroids. The retrospective space-time analysis of high and low rate clusters was then based on age and gender-stratified data applying the Poisson probability model. SaTScan was set to find clusters within the circular window including up to 20% of the population, while the maximum time span of the cluster was set to four months. We opted to use the log-linear temporal trend adjustment to ensure the comparability of rates within various periods. Resulting indirectly standardised rates (expressed as the relative risk) for each identified geographic cluster were estimated, and only clusters significant at p-value lower than 0.05 remained in the results.

2.3. Correlations, Similarities, and Clustering

In the first stage, we carefully selected several factors likely to influence annual SIR or season specific SIR of events leading to hospital admissions of the patient. The correlations among all SIR rates and the deprivation of the area, proportion of elderly in the population in the area, PM_2.5 averages, and the ethnic structure of the population are visualised using a correlation matrix. We used a non-parametric Spearman’s rank correlation (Spearman’s ρ) to analyse associations between pairs of characteristics above. However, the correlation provides only a general overview of the associations that might differ locally. To compensate for this, the local correlations among the characteristics of neighbouring CAUs were computed. The local correlation coefficient between the pair of variables for individual spatial units provides local estimates of correlation based on the spatial neighbours of the units [40,41]. This method allows for the exploration and visualisation of the strength of associations that vary locally, showing the pattern of similar and dissimilar units. In this paper, we employed local Spearman’s ρ based on the first order queen’s contiguity of CAUs [42].

Cluster analysis could be described as a group of methods, whose purpose is to identify subsets of similar objects in the data leading to the formation of classes. While spatial clustering creates groups based primarily on spatial proximity (neighbourhood) or location, and one common variable, the methods of the multivariate clustering aim to categorise the set of objects with the emphasis on their quantitative and/or qualitative characteristics [43]. Nonetheless, several attempts at the combination of both approaches have appeared in recent years [44,45]. At the end of the clustering process, it is necessary to characterise and interpret the classes.

We clustered several characteristics of CAUs to get a picture of areas that show similar characteristics. Firstly, we applied Euclidean distance based metric measure of dissimilarity on the data as advised by simulations available in clustersim package for R [46]. These simulations provided the necessary sensitivity analysis of possible results of clustering that aimed to remove an excessive subjectivity from the selection of clustering parameters (number of clusters, type of dissimilarity measure, type of clustering). We then re-weighted the computed dissimilarity matrix spatially using a standardised distance matrix of CAU centroids, so the original measure of dissimilarity was spatially smoothed to ensure better contiguity and thus a local similarity of neighbouring spatial units. Lastly, the Ward algorithm [47] was utilised as the agglomeration method for the hierarchical clustering.

2.4. Geographically Weighted Poisson Regression

The most common method used for the prediction and modelling of relationships between dependent and independent variable(s) is a standard linear regression (OLS). However, when modelling count data, the Poisson generalised linear model (GLM) is a suitable modelling method. This is especially useful in case of rare events such as disease cases. Although OLS and GLM often fit non-spatial data well, when used with spatial data the methods can be problematic due to violations of several assumptions such as the normal distribution and homogeneity of residuals, and their independence (lack of autocorrelation). Researchers also assume that the relationships are constant over the space meaning the associations are spatially stationary [48]. However, spatial data often bring these behaviours as a result of their nature, in other words, the data can be (auto) correlated, introduce heterogeneity into linear models and the relations can be nonstationary in the study area. Also, as this study aims to explore local variations of the effect of independent variables on the outcome, the global modelling methods such as (generalised) linear regression models missing the necessary detail (local focus) since they focus mostly on the description of general global relations. The family of Geographically Weighted Regression (GWR) models [49] has been developed in order to capture the spatial non-stationarity and spatially varying associations, and it has been utilised in quantitative social and health sciences [48,50,51].

Motivated by the existing spatial pattern of the dependent variable (spatial autocorrelation existing in the data) and possible local variations of effects of exploratory variables on the target that this may pose, we opted for the usage of Geographically Weighted Poisson Regression (GWPR) model as implemented in self standing software GWR4 [52] and R packages spgwr [53] and GWmodel [54] that follow the original work described in [55]. The GWPR model enables the investigation of local spatial variations with regard to the number of COPD hospitalisation events in individual CAUs in the study areas. The GWPR model is defined by (Equation (1)) [48,55]:

$$y_i ~ Poisson\left(N_i exp\left({\displaystyle \sum }_k{\beta }_k\left(u_i,v_i\right)x_{k,i}\right)\right)$$

(1)

where y_i denotes a dependent variable (the total number of COPD hospitalisations in CAU), x_k,i k-th independent variable including the constant term, and N_i the offset variable corresponding to population size at risk at location I, and (u_i, v_i) is the geographic coordinate of the (centroid of) location i. The coefficients β_k (u_i, v_i) are assumed to be smoothly varying conditional on their location.

Following the technique of GWR, one can calibrate this model using a kernel regression methodology in which smoothed geographical variations of parameters are estimated with a spatial weighting kernel [55]. In this study, the adaptive bi-square kernel weighting function (Equation (2)) was used as it has a clear-cut range where kernel weights are non-zero, and it is more suitable for when one seeks a definitive local extent for model fitting [48].

$$W_{ij}=\{\begin{array}{c}{\left[1-{\left(\frac{d_{ij}}{d}\right)}^2\right]}^2\\ 0\end{array}\begin{array}{c}{d}_{ij}<d\\ otherwise\end{array}$$

(2)

where w_ij is the geographical weight assigned to the i-th observation at a regression based at point j in the study area (referred to as the j-th regression location), d_ij is the Euclidean distance between observation i and regression location j, and d is a non-linear parameter (bandwidth) size defined by a number of neighbours [48,51]. The choice of bandwidth size is of great importance since it controls the level of smoothing of the outputs as larger bandwidth size generally produces more smoothed estimates of coefficients [51]. A smaller bandwidth may also reduce the precision of local estimates since they are based on few observations, leading to large variance in the local estimates [48].

With a large dataset, the bandwidth selection can be made using a sample (%) of data points [56], or it can be selected subjectively based on the local knowledge or objectively based on the minimisation of some criterion as corrected Akaike Information Criterion (AICc) or cross-validation. The bandwidth of 52 neighbouring points for the adaptive bi-square kernel in this study was selected based on the AICc.

To identify the effects of social deprivation, air pollution, the proportion of elderly, and ethnicity on the number of COPD hospitalisations in Christchurch’s CAUs and to analyse the influence of the season, we created a pair of GWPR models. The first model focused on the assessment of the winter season while the second evaluated the rest of the year.

The common problem for explanatory variables in regression models is the presence of collinearity. As well as using (global and local) correlation, we also run the regression diagnostics in order to identify the presence of correlated independent variables. From these, we remove the variables with a high variance inflation factor (VIF < 10) that is recommended as an indicator of multicollinearity [57].

3. Results

3.1. Mapping

Look at the spatial distribution of hospital admissions due to COPD, there were several areas of a high number of hospital admissions located in a central city, southern centre and eastern suburbs, both season specific (Figure A2b,c in Appendix A) and annually (Figure 3a). However, the number of occurring events there might have been shaped by the population density. Thus, we calculated gender and age-adjusted incidence rate and standardised incidence ratio (SIR) that are more suitable for the direct comparison among areal units. Focusing on the differences between winter and non-winter months (Figure 3b,c), the incidence of COPD hospital admissions is generally higher in winter in nearly all CAUs with milder effects of winter in the west-central areas, and more obvious effects in the east and south-east suburbs of the city where the rates are also higher in non-winter months.

The SIR visualisation (Figure 4) shows higher winter hospitalisation rates even more clearly. While only the eastern part of the central city and Kaiapoi area in the north showed more cases than expected in non-winter months (Figure 4c), most of the city exhibited a higher number of cases than was expected in winter (Figure 4b). Generally, the city centre and (south) east of the city appears to be more vulnerable than surrounding areas, followed by the town of Kaiapoi in the north and some rather rural area units on the west. Nevertheless, the worsening is apparent for the winter months comparing overall situation, and the strength of the winter effect also skews the average SIR.

3.2. Spatial and Spatio-Temporal Patterns

The results of global Getis-Ord G-statistic detected statistically significant clustering processes existing in the study area (average SIR: G = 0.0091, expected G = 0.0076, p-value < 0.001; winter SIR: G = 0.0088, expected G = 0.0076, p-value < 0.001; non-winter SIR: G = 0.0096, expected G = 0.0076, p-value < 0.001). The results of local G_i* then shows all the significant clusters in Figure 5 showing Kaiapoi (north of the study area) and a small part of the east-central city as a high-value cluster in the non-winter season, while most of the south and south-east central city was indicated as high cluster value during the winter. The southernmost CAUs of the study area have been suggested as clusters of low values during the winter season and also in average without season specification. While central city and its eastern part, as well as Kaiapoi also appeared in the analysis without season specification, the situation in the densely populated CAUs of central city aggravated in winter. The non-winter hot spots (Figure 5c) are larger in overall area (especially in the north) when compared to winter situation (Figure 5b). However, an appearance of a hot spot in the central city during the winter makes the difference in the amount of population affected (47,955 people in winter vs. 16,960 in non-winter) as well as in the absolute number of cases of COPD admissions coming from these clusters (207 in winter vs. 81 in non-winter).

Up to this point, we have proceeded only with spatial analysis without the temporal dimension for the spatiotemporal evaluation of hot and cold spots (Figure 6). Using all the data provided and aggregated them into CAUs by individual months, four spatiotemporal clusters were identified. Three clusters of high values occurred during the winter seasons of 2016 and 2017 while the fourth—a cluster of low values—occurred in the January–March 2017 period (12 CAUs, relative risk (RR) = 0.11, the green cluster in Figure 6). Both the high values clusters located in the central to eastern suburbs of the city (21 CAUs, RR = 2.18, the purple cluster in Figure 6) and south-western suburbs (10 CAUs, RR = 2.08, the coral cluster in Figure 6) took place in the winter of 2017. In the winter of 2016, the smallest of significant high values clusters were located in Kaiapoi in the north part of the study area (5 CAUs, RR = 5.09, the blue cluster in Figure 6). The low-value cluster appeared just west of the central city adding to the initial findings coming from the map displays of SIRs (Figure 4).

3.3. Correlations and Similarities

It is not uncommon to find that human health is affected by a mix of environmental and social conditions rather than just a single factor possibly triggering the outbreak of a certain disease or health impairment. The situation is no different in Christchurch, and although our main aim was to explore the effect of winter air pollution on the COPD hospital admissions, we also considered other, primarily socioeconomic conditions. Firstly, we assessed if PM_2.5, socioeconomic deprivation, the proportion of elderly in the population, and ethnicity were associated with the season-specific standardised incidence ratios of COPD hospital admissions in CAUs. The correlation matrix (Figure 7) shows the non-spatial Spearman’s correlation coefficients for the aforementioned factors. In general, the socioeconomic deprivation correlates positively with the SIR of COPD hospitalisations regardless of the season. Other significantly positive correlations can be found for the proportion of Māori population in the CAU and the case of the average level of PM_2.5, although this association is even weaker. The proportion of European population appeared to be negatively associated with the winter SIRs. Finally, the proportion of elderly in the population of CAU correlates negatively with SIR, although the correlation is rather weak.

The Spearman’s correlation coefficients explain only a general level of association between investigated area level characteristics, and thus it is not suitable to describe spatial variations of these associations. The local correlations as presented by [41] introduces the improvement of the method by calculating the local associations in the spatial neighbourhood using the contiguity scheme. This allows a spatial assessment and geovisualisation of the variation of the association between the pairs of characteristics. A map matrix representing the local correlations of season-specific SIRs and selected characteristics is shown in Figure 8. The figure compares a set of geographically weighted correlation pairs that were later used in the spatial regression and clustering. Positive local correlations are shown in shades of red, while negative correlations in shades of green colour.

The overall results of local correlations correspond with the results of Spearman’s correlation highlighting socioeconomic deprivation as an important correlate of SIR regardless of the season. However, the differences between winter and non-winter season were identified in the central and southern parts of the study area. Analysis found that the deprivation acts as important direct factor in the winter, while there are more CAUs showing negative local correlations in the non-winter period. Also, the local correlations of SIRs and PM_2.5 exposure indicate different patterns in winter and non-winter season. While there is a distinctive area of the negative local correlation in the city centre in the winter, this is not evident in the rest of the year. The proportion of elderly people in the populations is positively associated with the COPD SIR mainly in the northern part of the study area (specifically Kaiapoi) whereas it is more of a negative local correlation in the rest of the study area. The role of the ethnicity does not differ much when comparing general and season-specific SIRs. The proportion of European and Māori population in CAU performs inversely. The higher proportion of Māori (or lower the proportion of Europeans) in the population the higher is SIR in the area (and vice versa).

3.4. Clustering

In order to further investigate the potential associations between COPD hospitalisations and social and environmental characteristics, we assessed possible similarities between area units that would allow for forming groups (clusters) and the meaningful description of that could distinguish between these groups.

Here, the analysis focused on the possible grouping of area units (CAUs) into which we aggregated individual level data and calculated rates. It is based on the multivariate characteristics of CAU represented by the attributes and characteristics extracted by the combination of hospital admission records, census records, socioeconomic deprivation index, and averaged season specific levels of PM_2.5. We identified five clusters of CAUs in the study area. A summary of these CAU clusters is shown in

Table A1. Basic statistical characteristics of clusters of CAUs. Underlined numbers indicate the highest value of the characteristics amongst groups; bold numbers denote the lowest values.

	Group 1 (n = 17)			Group 2 (n = 34)			Group 3 (n = 16)			Group 4 (n = 48)			Group 5 (n = 17)
	Mean	SD	Median	Mean	SD	Median	Mean	SD	Median	Mean	SD	Median	Mean	SD	Median
Deprivation	2.53	1.94	1.00	7.41	1.02	7.00	7.75	1.24	8.00	2.77	1.48	2.50	3.71	1.96	3.00
PM2.5 winter	9.21	0.86	9.24	8.92	0.81	8.98	10.54	0.41	10.52	9.46	0.78	9.47	8.63	0.35	8.65
PM2.5 non-winter	5.54	0.33	5.44	5.44	0.23	5.43	5.98	0.32	5.92	5.56	0.32	5.57	5.07	0.27	5.08
Elderly proportion	16.02	4.10	15.38	11.74	3.96	11.48	16.18	3.46	16.66	16.91	4.21	16.46	16.08	3.15	16.26
SIR winter	0.71	2.13	0.00	2.57	1.72	2.24	2.07	0.89	2.06	1.11	0.66	0.96	0.98	0.88	0.62
SIR non-winter	1.65	3.31	0.58	1.15	0.83	0.94	1.47	0.91	1.25	0.55	0.33	0.47	0.46	0.29	0.42
European	0.91	0.04	0.89	0.74	0.08	0.75	0.78	0.07	0.78	0.88	0.03	0.88	0.73	0.05	0.75
Māori	0.05	0.03	0.05	0.11	0.03	0.11	0.13	0.03	0.13	0.06	0.02	0.07	0.05	0.01	0.05
Asian	0.04	0.03	0.04	0.10	0.06	0.09	0.06	0.04	0.06	0.05	0.03	0.04	0.19	0.04	0.19
Winter admissions	0.12	0.14	0.00	0.57	0.17	0.54	0.42	0.14	0.43	0.50	0.11	0.50	0.49	0.15	0.53

located in Appendix A where numbers in bold mark the lowest mean value of individual characteristics amongst all groups while the grey cells with underlined numbers denote the highest mean value of characteristics amongst all groups. Figure 9 then presents the graphical visualisation of groups represented by their mean values. The location of the clusters of CAUs in the study area is then depicted in the map in Figure 10. The colours of groups in the table match with colours of groups in the boxplot and map.

CAUs grouped within the first identified cluster (yellow areas in Figure 10) are located mainly around the outskirts of the study area creating a spatially contiguous region only in the southern part of the study area. This cluster is characterised by the highest proportion of the population of European ethnicity and the lowest deprivation as well as the highest non-winter hospitalisation rate, but the winter SIR and the ratio of winter months’ hospital admissions are the lowest. The fourth group (green areas in Figure 10) has very similar characteristics as the first group, however, the major differences between these two clusters are the higher proportion of elderly in the population and an inverse regime of the hospitalisations with a low non-winter SIR and significantly higher winter SIR. CAUs within this cluster are also located around the outskirts of the central city creating two regions, one to the south and another in the northern suburbs of the city. The second cluster covers most of the city centre and is distinctive to the other group mainly due to exhibiting the highest winter SIR and ratio of winter hospital admissions; these locations are also deprived, but with the lowest proportion of elderly in the population. The most deprived and most affected by the air pollution is the third cluster (red areas in Figure 10) with the highest proportion of Māori population. The winter and non-winter SIRs for cluster three are the second highest. CAUs belonging to this cluster are spatially concentrated in the eastern suburbs of the central city or in Kaiapoi on the northern boundary of the study area. The fifth cluster (purple areas in Figure 10) is the most spatially homogenous covering a contiguous area in the western part of the city. This cluster is the least affected by PM_2.5 also having the lowest non-winter SIR amongst all groups, and there is also the greatest proportions of the Asian population.

3.5. Geographically Weighted Poisson Regression

The ordinary Poisson regression models have been commonly utilised for the analysis of disease events when investigating counts and rates. However, as a global method, it does not allow for the exploration of local variations of relationships between dependent and explanatory variables. For this reason, the GWPR was used in this study. Nevertheless, the pair of season-specific global Poisson regressions was constructed for two reasons, (1) to find whether there is any strong collinearity present in the data, (2) to find out whether the GWPR model brings any improvements in explanatory part when compared to the global regression model.

There were statistically significant correlations on global as well as local level identified in the data (Figure 9 and Figure 10) possibly introducing multicollinearity into the regression model. The VIF (variance inflation factor) was used as indicator of multicollinearity in the fully specified ordinary Poisson regression model, where number of COPD hospitalisations acts as dependent outcome while socioeconomic deprivation, PM_2.5, Proportion of elderly in the area, Proportion of Europeans in the area, Proportion of Māori population in the area, Proportion of Asians in the area and the interaction term between socioeconomic deprivation and PM_2.5 are the independent variables Since we were dealing with counts of COPD hospitalisation, and COPD is a disease with an impact on the older ages of the population, we used the (a natural logarithm of) number of people aged 65 or more living in CAU as an offset for the model. Observed VIF values in the models ranging from suggested 1.4 to 296 suggested that strong collinearity was existed in the models. As a result of highest VIF values, the interaction terms between socioeconomic deprivation and PM_2.5, and Proportion of Asian population were removed from the explanatory variables causing that the VIF values of explanatory variables in reduced models were lower than 5 (1.3–4.0) suggesting no remaining collinearity concerns.

The explanatory variables from the reduced ordinary Poisson regression model served as a basis for GWPR in which the local variations of relationships regarding the COPD hospitalisations were investigated. The bandwidth of 52 neighbouring points for the adaptive bi-square kernel was selected based on the AICc of the winter specific model. Although the non-winter model bandwidth selection process evaluated 75 neighbouring points as optimal, we have chosen the threshold of 52 neighbours to ensure the comparability of both models.

The summaries of season-specific global Poisson regression models, as well as summarised results of GWPRs, are shown in

Table 1. Results of global Poisson regression model and Geographically Weighted Poisson Regression (GWPR) model parameters for winter and non-winter COPD hospitalisations in Christchurch. Bold numbers represent significant relationship at p < 0.05 level. GLM: generalised linear model; AICc: corrected Akaike Information Criterion; GWR: geographically weighted regression.

Independent Variable	Global Model GLM Poisson Winter				Global Model GLM Poisson Non-winter				GWR Poisson (n = 52) Winter					GWR Poisson (n = 52) Non-winter
	Est.	Std. Error	z-Value	Pr (>|t|)	Est.	Std. Error	z-Value	Pr (>|t|)	Min.	Q1	Median	Q3	Max.	Min.	Q1	Median	Q3	Max.
(Intercept)	−3.92	0.46	−8.47	0.00	−6.57	0.68	−9.68	0.00	−12.83	−5.59	−3.06	−1.69	−0.77	−17.29	−8.05	−6.35	−3.18	2.34
Deprivation	0.14	0.02	6.32	0.00	0.09	0.03	2.83	0.00	−0.02	0.07	0.12	0.14	0.30	−0.04	0.02	0.05	0.11	0.25
PM2.5	0.03	0.05	0.59	0.56	0.56	0.12	4.61	0.00	−0.17	−0.04	0.04	0.21	0.74	−0.56	0.14	0.57	0.82	1.13
Elderly ratio	−0.03	0.01	−3.72	0.00	−0.04	0.01	−3.24	0.01	−0.09	−0.05	−0.03	−0.01	0.04	−0.07	−0.03	−0.02	−0.01	0.03
European	−0.61	0.54	−1.12	0.26	−1.22	0.72	−1.69	0.09	−4.58	−3.01	−1.74	−0.26	5.72	−6.12	−3.89	−2.22	−0.66	12.25
Maori	1.05	1.38	0.76	0.44	2.39	1.56	1.53	0.08	−13.49	−2.95	0.62	4.67	14.13	−3.97	−0.24	2.21	4.79	10.69
(pseudo) R2	0.44				0.48				0.62					0.65
AICc	365.17				196.40				300.15					186.65

. The local parameters are described by summary presenting values of minimum, lower quartile (Q1), median, upper quartile (Q3), and maximum. The spatial distribution of local parameters of explanatory variables is shown in Figure 11 and Figure 12 for the winter season and non-winter season respectively. All parameters in the map were classified to sextiles, hence the colours and individual classes are comparable between the two pictures. It is also worthwhile to mention that GWPR estimates are log-odds (due to the definition of Poisson regression model), so negative values of parameter denote negative relationship while positive values indicate positive relationships between number COPD hospitalisations and explanatory variables.

According to the results of a global Poisson regression model, the socioeconomic deprivation increases the number of COPD hospitalisations in both winter (log-odds 0.14) and non-winter season (log-odds 0.09). However, when comparing the local relationship between socioeconomic deprivation and number of COPD hospitalisations, it is apparent that the association is much stronger during the winter in most of the study area (Figure 11a and Figure 12a) and especially in the eastern part of the study area, where most of CAUs are also denoted as statistically significant based on pseudo t-value (t > 2.032). The local intensification [11] of the association in winter is noticeable, particularly in the southern and eastern part of the study area, with an apparent borderline of change of the association leading through the city centre. The opposite process of weakening of the association in the winter is apparent in the northernmost part of the study area (Kaiapoi).

One of the main goals of this study was to increase the understanding of the influence of winter air pollution on the number of COPD hospitalisations in the study area. Previous studies [58,59] found a relationship between peaks of air pollution and delayed increased hospital admissions. We further elaborate on this finding, and using the GWPR we demonstrated the local and season-specific variations of this relationship. Average levels of PM_2.5 proved to be globally significant only in the non-winter season as well as the positive direct effect of air pollution is locally stronger in non-winter months. Locally, this is particularly visible in central to south-western CAUs and also in the eastern CAUs (Figure 11b and Figure 12b). Although weaker in winter, the air pollution is able to aggravate the situation in the areas dominated by deprivation as well as in other CAUs located in the southwest of the studied region. The dominating effect of socioeconomic deprivation [14] could explain these findings since the bad quality housing in higher deprived areas can potentially overshadow the effect of increased levels of air pollution in the winter season, the effect also mentioned by other studies [1,11,60]. Although the relationship between COPD hospitalisations and air pollution might seem to be weaker in the winter, there is the remaining progressive effect of air pollution caused by increased winter levels of air pollution (particulate matter) when compared to the non-winter season. Nonetheless, the air pollution still has a multiplicative effect on the number of COPD hospitalisations regardless of the season as mentioned regularly in the literature [16,59,61].

The link between the proportion of elderly people [3,10] (representing an increased average age of CAU) living in the area and number of COPD hospitalisations is globally statistically significant in both seasons. Its local importance (Figure 11c and Figure 12c), although rather weak, is noticeable in the west of the study area and its spatial pattern is stable in both seasons following the direction from south-east to north-west of the study area. The proportion of people of NZ European ethnicity living in the CAU showed globally negative association with COPD hospitalisations insignificant in both seasons. This “protective“ effect is visible in most of the study area except the northern outskirts (light green areas Figure 11d and Figure 12d), and it is especially strong in the east. Its pattern in both seasons is rather stable. The last explanatory variable is the proportion of Māori population. The global effect of the variable was not evaluated as significant in global models. However, according to the results of GWPR, its relationship to the total number of COPD hospitalisations in both seasons can be locally significant. Positive effect can be seen particularly in the western central part of the study area, while there is a negative effect in eastern suburbs both in winter and non-winter. The pattern of the values of the parameters is similar (Figure 11e and Figure 12e), but in overall the effect of Māori is stronger (positive and negative) in the winter season. In winter, CAUs denoted as significantly related to COPD hospitalisations extended the area where there was a weaker influence of deprivation, i.e., if there is more Māori living in the area the number of hospital admissions grew, while in summer the “Māori effect” multiplies the effect of air pollution. Although rather subtle, especially in winter, this relationship should not be lost as there is a burden of COPD for Māori and Pacific peoples in New Zealand [62,63,64]. GWPR explains COPD hospitalisation better in southern and eastern CAUs of the study area in both winter and non-winter (Figure 11f and Figure 12f).

4. Discussion

Based on the individual and aggregated geographical data and area-level characteristics, we detected high-risk health hot spots related to increased risk of COPD hospital admission both overall and season-specific with patterns of increased SIR being stronger during the winter season. Spatiotemporal scan statistics supported this by finding a statistically significant hot spot only in winter 2016 and 2017. The central city was identified as one of the most vulnerable areas in most of the spatial analyses, and it was also classified as an individual group (the fourth cluster) by spatiotemporal scan statistics.

We show the continued importance of including both individual and area level factors when examining and analysing disease patterns. However, as in other health and geohealth research, it becomes salient that researchers must find a thin threshold in the study design and the study results presentation, that on the one hand allows the most efficient data manipulation, analysis and overall usability of the study, while ensuring the data confidentiality and privacy of individuals on the other. For this reason, we decided to use reasonably large census area units (CAU) to aggregate the data. CAUs also provides a practical trade-off between the detail (spatial granularity) of a study and its usability. The use of individual records would provide further insight into the mechanism of social, geographical, and environmental interactions affecting people who suffer COPD. However, in this study, we decided to focus rather on the general effect of the place on the people with COPD rather than focus on individual patients.

Analyses also differentiated only between winter and non-winter season, not considering the immediate temperature changes or monthly temperature averages. This was done for several reasons, (1) we did not operate with enough monitoring sites recording temperature within the city as not all the air pollution monitoring sites measured the temperature, so the coverage of the temperature grid would not be sufficient when compared to the other datasets in the study; (2) the temperatures differences in Christchurch do not have to be that significant when comparing individual months, i.e., according to long-term monthly averages May is colder than September, however there is a “winter effect” [11,12] in September when the symptoms and medical conditions of patients cumulate over the winter; (3) we were interested in how the studied associations of the winter season differ from non-winter.

It is also of value to draw the attention to the clustering process itself. In both approaches presented in this paper, the evaluation of spatial (and spatiotemporal) pattern and multivariate clustering, necessitated a degree of subjectivity was needed to analyse and interpret the results. Decisions had to be made about the spatial and attribute weights, methods of similarity measures and/or a clustering type. That is why we conducted a sensitivity analysis on every step of the analyses to minimise the risk of misleading results. Since the spatial autocorrelation exists in the data, we employed the generalised Poisson GWR that allows accounting for the heterogeneity introduced by spatial autocorrelation and also allows for the exploration of local variations of parameters and associations. However, although commonly used in scientific studies, there is a critique existing on the GWR, since it sometimes can identify spatial patterns even in the situations where there are not any [65,66,67]. We tried to precede this by looking at the spatial patterns first, so we knew there were hot and cold spots of outcomes existing in the study area.

We carefully selected several explanatory variables based on the previous research and literature considering environmental conditions and lung irritants (particulate matters), socioeconomic factors and quality of life (deprivation), and demographic structure (age and ethnicity structure) of CAUs. Yet, we are sure that all possible conditions affecting COPD patients were not included, so the constructed model and clustering has to be taken as an approximation of reality. For reasons mentioned earlier, we did not include personal information (including other health, ethnicity and social data) and physical geographical factors as temperature and humidity. Some of the interactions and explanatory variables were also removed due to strong correlations causing multicollinearity in final global and local models. Geographically weighted Poisson regression model was selected as a suitable tool for exploration of spatial variations of the relationship between COPD hospitalisations and selected explanatory variables. The bandwidth and kernel type used in the study were chosen based on the AICc optimisation and recommendations provided by the software, however, if one decides for another criterion, then the resulting optimal bandwidth could be different, and because of that the local parameters estimates and possibly even their interpretation could change as well. We are also cognisant of a post-earthquake context linked to census results in the Christchurch area that reflects changes in the city’s urban and social structure. Nevertheless, the data are the most complete available to researchers. Earthquakes were also a reason why the winter air pollution situation in Christchurch has been improving. The primary source of the winter air pollution in Christchurch is smoke from solid fuel-burning domestic heating systems [68,69,70]. Firstly, many chimneys were destroyed in the earthquakes and people made the decision not to replace them, and secondly wood burners are now more tightly regulated by the changes in local policy environment relating to the type of burners being installed in new and renovated houses. Although smoking is often mentioned as the leading trigger of COPD, we did not include it in the analyses, as it is likely to be underreported in the hospital records (only about 2.5% of cases) due to self-reporting and often unknown smoking history of patients.

5. Conclusions

The primary focus of this study was to examine the impact of air pollution on patterns of COPD hospital admissions, and, our results support the finding mentioned in previously conducted studies [5,6] that highlight the important role of air pollution as a risk factor for COPD. This study therefore contributes further evidence that mitigating air pollution would potentially lead to lower rates of hospitalisations for those with COPD. We also argue that the winter season is closely associated with the increase of hospitalisation especially in the areas with high socioeconomic deprivation due to colder temperatures and cumulative effect of winter [12]. However, the air pollution proved to be interconnected with COPD hospitalisations even in non-winter season showing the importance that may be boosted even more in case of extreme air pollution events or people moving in or through an environment with high concentrations of lungs irritants.

We found a complex picture and in it what is significant is not only a single factor, but rather a group of closely related socio-environmental conditions that includes measures of deprivation, influence of ethnicity and proportion of elderly in the area unit or neighbourhood all of which affecting hospitalisation rates locally. We have also demonstrated that an understanding of the spatial distribution of factors that affect COPD is vital in order to gain a more complete understanding of the precise mechanisms that may be affecting COPD hospital admissions. Our contribution in this instance is to emphasise that a geographic perspective aids understanding of patterns and processes of respiratory disease and it is a valuable addition to more traditional non-spatial studies.