Inequalities in exposure to the air pollutants PM2.5 and NO2 in Australia

Vulnerable subpopulations may be exposed to higher levels of outdoor air pollution than the rest of the population. Due to the potential for this to exacerbate their existing health burden, concerns about disparities in subpopulations’ air pollution exposure have motivated international public health researchers to examine this topic. In Australia, such research is lacking to date, despite heterogeneity in air pollution at multiple spatial scales across the continent. This study aimed to investigate disparities in exposure to two health-relevant outdoor air pollutants: particulate matter <2.5 μm (PM2.5) and nitrogen dioxide (NO2). We used national land-use regression models to estimate annual average concentrations of PM2.5 and NO2, and area-level census data on ethnicity, age and socio-economic status (SES) to calculate the bivariate associations between each census-derived variable with the concentration of air pollutants. We also used multivariable models including specific measures of SES as covariates to assess to what extent associations were explained by SES. Associations were calculated separately for rural and urban areas using generalised additive models which accounted for spatial autocorrelation. Bivariate results showed significant nonlinear associations (p < 0.001) between vulnerable subpopulations and pollutant concentration. These associations suggested that areas with greater socio-economic disadvantage, a higher proportion of ethnic minorities, and elderly people are exposed to higher concentrations of PM2.5 and NO2, although differences in the magnitude of exposure were small overall. Our multivariable models showed that the associations between ethnic minorities and pollutant concentration appear to be substantially affected by area-level SES. Our results suggested that these vulnerable subpopulations are inequitably exposed to PM2.5 and NO2. While the magnitude of differences in exposure were generally small, the predicted differences in exposure among vulnerable subpopulations could contribute to a potentially avertable health burden at a population-level.


Introduction
Outdoor air pollution is recognised as one of the greatest global threats to human health, with 4.2 million deaths per year in 2016 attributed to outdoor air pollution (WHO 2018). Of the wide range of air pollutants, PM 2.5 and NO 2 are among the most ubiquitous and extensively studied, and they are known to be associated with adverse health effects (Weinmayr et al 2010, Rückerl et al 2011, Hoek et al (Chakraborty et al 2011, Taylor 2011, Hajat et al 2015, Holifield et al 2017. Despite substantial international research on vulnerable subpopulations' exposure to air pollution, there is a paucity of comparable research in Australia (Chakraborty andGreen 2014, Masterman-Smith et al 2016). This is notable because, although Australia has relatively clean air by global standards, levels of PM 2.5 and NO 2 observed in previous studies are associated with negative health effects (Johnston et al 2007, 2011a, Hanigan et al 2008, O'Loingsigh et al 2017, Hanigan et al 2019. Disparities in exposure are often caused by socio-economic and political factors which have been found-in a range of environmental contexts-to cluster vulnerable communities into heavily polluted areas (Schlosberg 2007, Mohai et al 2009, Walker 2009, Holifield et al 2017. To understand if these patterns might be similar in Australia, our study explored whether this exposure disparity to PM 2.5 and NO 2 exists for vulnerable subpopulations located in rural and urban areas across the continent.

Background
The air pollutant PM 2.5 consists of all airborne particles below 2.5 μm in diameter. It has a range of anthropogenic and natural sources and consequently is both concentrated in major urban centres and widely, and patchily, distributed across the continent (Kelly and Fussell 2016, Knibbs et al 2018a, Park et al 2018. The air pollutant NO 2 is primarily formed from traffic emissions and as such, has a skewed urban distribution that contrasts markedly with the distribution of PM 2.5 over the continent, but is also similarly concentrated in major urban centres (WHO 2000, Weinmayr et al 2010, Cowie et al 2016, SoE 2016. Air pollution in rural and remote areas is primarily caused by natural sources, such as bush fires and dust storms, or from industrial activity (Johnston et al 2011b, SoE 2016, Knibbs et al 2018a. This air pollution tends to be seasonal and often consists of short-term acute peaks, unlike the longer-term chronic air pollution found in cities (CSIRO 2001, Johnston et al 2007, Hanigan et al 2008, Williamson et al 2013.
The pollutants PM 2.5 and NO 2 are associated with a range of health effects including respiratory and cardiovascular diseases, and in extreme cases, mortality (Brugge et al 2007, Rückerl et al 2011, Hamra et al 2015, Kelly and Fussell 2016, Khreis et al 2017, Schraufnagel et al 2019a, 2019b. The associations between negative health effects and the concentration of pollutants are significant for acute and chronic exposure, although the exact associations vary between exposures (Rohr and Wyzga 2012, Hoek et al 2013, Atkinson et al 2014, Faustini et al 2014, Wyzga and Rohr 2015, Li et al 2016. There is evidence that PM 2.5 and NO 2 have a steeper exposure-response curve at low concentrations, indicating that increases in low concentrations intensifies negative health effects more than equivalent increases at higher concentrations (Burnett et al 2018, Hanigan et al 2019.
World Health Organisation (WHO) guidelines state that average PM 2.5 concentrations should not exceed 10 μg m −3 annually and 25 μg m −3 over 24 h; and, that average NO 2 concentrations should not exceed 40 μg m −3 annually and 200 μg m −3 h −1 (WHO 2006). Air pollution concentration in Australia is generally low relative to these standards, as the Australian population-weighted average PM 2.5 concentration for 2013 was estimated to be 5.93 μg m −3 (Brauer et al 2016), and no exceedances of hourly or annual air quality standards for NO 2 have occurred since the 1990s (SoE 2016). However, these guidelines do not reflect current research which has found no safe concentration levels of these pollutants (Kelly and Fussell 2016, Hanigan et al 2019, Schraufnagel et al 2019b. Despite Australia's seemingly low pollution levels, ambient and short-term concentrations of PM 2.5 and NO 2 have been correlated with negative health outcomes in a number of Australian studies . This is problematic because of the few studies that have, Indigenous people have been found to be more susceptible to respiratory and ischaemic heart disease from, for example, exposure to PM 2.5 (Johnston et al 2007, Hanigan et al 2008. Sources of air pollution in rural areas, particularly natural sources such as bush fires and dust storms, may inequitably affect Indigenous communities, compounding their health burden from air pollution exposure (Knibbs and Sly 2014, Masterman-Smith et al 2016, Melody et al 2016. Our changing climate may serve as an additional factor in exacerbating these health burdens on Indigenous Australians, as well as intensifying emissions from natural sources in rural areas (Webb et al 2014, Green et al 2015. In major urban areas, emissions from transport and industry may result in vulnerable subpopulations experiencing inequitable chronic exposure to numerous air pollutants. A nation-wide study of industrial polluting sources (Chakraborty and Green 2014), and a study of annual average NO 2 concentration in urban areas (Knibbs and Barnett 2015), found that socio-economically disadvantaged and Indigenous people were inequitably exposed to air pollution compared to the rest of the population.
A significant body of international quantitative Environmental Justice research has found that vulnerable subpopulations-particularly low-income communities and ethnic minorities-are exposed to higher concentrations of air pollution compared to the rest of the population (Chakraborty et al 2011, Taylor 2011, Hajat et al 2015, Holifield et al 2017. A small number of studies have examined age as well, with some evidence of disparities in exposure to air pollution between age groups, although these findings are not consistent across studies (Clark et al 2014, Gurram et al 2015, Clark et al 2017, Gurram et al 2019. Other studies have examined disparities in exposure separately for rural and urban areas, allowing their analyses to account for varied exposure, either acute or chronic, from different sources of air pollution and different factors impacting residential patterns by urbanicity relative to the distribution of pollution (Clark et al 2014, 2017, Choi et al 2016, Rosofsky et al 2018. Disparities in exposure are typically caused by socio-economic and political factors including housing prices, historical land-zoning patterns and poor political representation (Mohai et al 2009, Cooper et al 2018, although studies suggest historical and institutional discrimination can also contribute to them (Pulido 2000, Taylor 2011, Holifield et al 2017. These factors lead to polluting sources being sited near disadvantaged communities, or the migration of disadvantaged people to co-locate near polluting sources, resulting in these subpopulations spatially clustering in heavily polluted areas (Mohai et al 2009, Mohai andSaha 2015).
The findings of this body of research have major implications for public health, as many inequitably exposed subpopulations are also potentially vulnerable to exposure to air pollution. Due to their physiology, young and elderly people are particularly susceptible to air pollution (Rückerl et al 2011, Sacks et al 2011, Schraufnagel et al 2019a. Socio-economically disadvantaged people have also been found to be vulnerable to exposure, possibly due to poorer quality housing, which can increase their exposure to air pollution-or compound its effects (Evans and Kantrowitz 2002, Rauh et al 2008, Rückerl et al 2011, Wang et al 2017, Nieuwenhuijsen et al 2018. Similarly, ethnic minorities have been shown to be more vulnerable to exposure from air pollution (Sacks et al 2011, Bell et al 2014, Wang et al 2017. This may be partly explained by their socio-economic status (SES), cultural factors and institutional discrimination, which can make it harder for ethnic minorities to access resources to reduce their exposure or address health issues (Johnstone and Kanitsaki 2008, Durey et al 2012, Sherwood 2013, Kelaher et al 2014.
As inequitable exposure to air pollution has considerable public health ramifications, it is important to investigate large scale disparities in exposure, which can reveal which subpopulations are most exposed to various air pollutants and their sources. It is also important to analyse possible causes of inequitable exposure, particularly SES, to inform policy responses addressing environmental inequalities. To inform appropriate policy responses, this study investigates continental-level disparities in exposure to air pollution in Australia by answering the following research questions: • What are the continental level associations between vulnerable subpopulations and exposure to PM 2.5 and NO 2 in urban and rural areas of Australia?
• What is the influence of SES on the associations between vulnerable subpopulations by age and ethnicity and exposure to both pollutants?

Data-air pollution
We obtained estimates of annual average concentrations of PM 2.5 and NO 2 from satellite-based land-use regression (LUR) models (Knibbs et al , 2018a (Knibbs et al , 2018a. As we did not have more recent values, we used only the data from the year closest to the most recent Australia Census conducted in 2016, as the most current measure of vulnerable subpopulations' exposure to both pollutant, 2015 for PM 2.5 and 2011 for NO 2 . Previous studies showed that predicted and measured pollutant levels varied little year-to-year, and that spatial patterns of both pollutants were consistent over decadal and longer times scales in Australia (Knibbs et al 2018a(Knibbs et al , 2018b.
Consequently, it is reasonable to assume that the data for both pollutants are representative of current exposure.

Data-population variables
We acquired data from the most recent 2016 Australian Census, to determine the distribution of vulnerable subpopulations defined by ethnicity, age and SES (ABS 2016a). Subpopulations defined by these factors have been frequently found to be inequitably exposed to air pollution and/or vulnerable to air pollution compared to the rest of the studied population (Rückerl et al 2011, Sacks et al 2011, Bell et al 2014, Holifield et al 2017, Wang et al 2017. Data for each measure of population vulnerability were available for each SA1 unit. SA1s are the smallest census unit containing detailed demographic information, with ∼57 000 SA1s across Australia, varying in size from 0 to 328 000 km 2 and hosting an average of 407 people (ABS 2016a). In contrast, mesh blocks, which were the same census area unit that our datasets of PM 2.5 and NO 2 pollutant estimates were available for, only contain information on the number of people living in each area (figure 1).
We assessed potentially vulnerable subpopulations by age using census data of the percentage of people aged under 15 years old (percent under 15) and the percentage aged over 65 years old (percent over 65) in each SA1. For ethnicity, we measured the percentage of people in each SA1 who were Indigenous (percent Indigenous-defined as individuals that have declared their Indigenous status as Aboriginal and/or Torres Strait Islander). We also calculated the percentage of people of minority ethnicity (percent minority ethnicity-this is a measure of the number of non-Anglo European and non-Indigenous people that are within a culturally distinct ethnic minority, based on demographic research on waves of migration to Australia post-WWII) ( (table 1).
For each SA1, we quantified SES using the Socio-Economic Indexes for Areas (SEIFA) variables, which are four indices representing different aspects of SES, including economic resources and education. The four SEIFA variables are: the Index of Relative Socio-economic Advantage and Disadvantage (IRSAD), the Index of Relative Socio-economic Disadvantage (IRSD), the Index of Economic Resources (IER) and the Index of Education and Occupation (IEO). Each index is constructed using a complex weighted combination of values of different variables including measures of unemployment, income, education, occupation, house prices and internet quality (ABS 2013). These variables are selected and weighted based on the aspect of SES quantified by each index, with the number of variables used to construct each index varying from 9 to 25 (ABS 2013). For example, IEO is calculated from measures of education and occupation including the percentage of people that are over 15 and are unemployed, have completed high school and have a minimum diploma level educational qualification, each with a different weighting to calculate the index value (ABS 2013). To examine certain aspects of SES in greater detail, we examined specific measures of SES for education, employment, income and housing alongside the SEIFA indices. Many of the variables we examined had been used to construct the SEIFA indices, such as the percentage of people over 15 years with a higher school certificate, and the percentage of people in the labour force who were unemployed (table 2).
To minimise issues with low population counts, SA1s with fewer than 10 people, or which had missing values for any variable assessed were excluded from the final dataset (Fraser and Wooten 2005). This excluded between 2600 and 2850 SA1s overall depending on the pollutant variable (table 3). Index measuring financial aspects of socio-economic advantage and disadvantage, particularly income and wealth-calculated by ABS (2013) IEO Index measuring educational and occupational aspects of socio-economic advantage and disadvantage-calculated by ABS (2013) Education High School % Percentage of residents above 15 years old with a high school certificate qualification and/or equivalent Diploma % Percentage of residents above 15 years old with at least a diploma level tertiary education qualification Employment Unemployed % Percentage of residents above 15 years old that are unemployed Labour % Percentage of residents above 15 years old that are in the labour force Blue-collar % Percentage of residents above 15 years old with a blue-collar job-typically lower-paying jobs associated with trades-as defined by the census (ABS 2011) Blue-collar2% Percentage of employed residents with a blue-collar job-typically lower-paying jobs associated with trades-as defined by the census (ABS 2011) White-collar % Percentage of residents above 15 years old with a white-collar job-typically higher-paying managerial jobs and similar occupations-as defined by the census (ABS 2011) Income House-low inc. % Percentage of households with a low income as defined by the census-less than $26 000 a year (ABS 2013) House-high inc. % Percentage of households with a high income as defined by the census-over $78 000 a year (ABS 2013) Indiv-low inc. % Percentage of individuals above 15 with a low income as defined by the census-less than $15 000 a year (ABS 2013) Indiv-high inc. Associations between pollutant concentrations and vulnerable subpopulations were estimated using Generalised Additive Models (GAMs) with smoother splines. GAMs can account for nonlinear and non-monotonic relationships (Ruppert et al 2003, Wood 2017, which we anticipated would be present in the associations, based on our previous work (Knibbs and Barnett 2015).
Smoother splines were applied to each population variable using natural cubic basis functions, as well as Robust Estimation of Maximum Likelihood to calculate the smoothing parameter and prevent overfitting.
We first estimated the bivariate associations between each population variable and each pollutant. We calculated separate GAMs for rural and urban areas using the ABS Greater Capital City Statistical Areas categorisation to define census areas because we were interested in potential rural/urban differences that may reflect underlying differences in pollutant sources, subpopulation characteristics, or both (ABS 2016b). Based on this categorisation we designated ∼19 700 SA1s as rural areas and ∼35 500 SA1s as urban areas, with an average area of 360 and 1.5 km 2 (Q 1 : 0.23 and 0.11 km 2 ; Q 3 : 4.3 and 0.27 km 2 ) and an average population of 387 and 442 respectively. The number of large rural SA1s with relatively small populations meant that a more fine-scale analysis would likely provide limited insight and consequently was not pursued.
Prior to analysis, NO 2 values were log-transformed to satisfy statistical assumptions regarding residuals. For the SAT-W (PM 2.5 ) GAMs, a gamma distribution was used to account for the non-normal distribution of SAT-W values. We accounted for spatial autocorrelation in the residuals of the GAMs by including interaction terms for the geographic coordinates of each data value in each model, as we found evidence of spatial autocorrelation in the original GAMs (p-value <0.001, chapter S3 section 3.3). Further details on how we accounted for spatial autocorrelation are also included in chapter S3 section 3.3.
We then fitted multivariable models with covariates for SES to determine if the bivariate associations estimated between ethnicity and age with exposure to PM 2.5 and NO 2 were still present when the influence of SES was accounted for. This analysis was carried out because previous studies have shown that disparities in exposure experienced by vulnerable subpopulations may be caused by the status of those subpopulations (Chakraborty et al 2011, Taylor 2011, Holifield et al 2017. We used specific measures of SES as covariates rather than the SEIFA indices as there was a high degree of multicollinearity between indices. In addition, each index value did not explain as much variance in pollutant values as the specific variables did, suggesting that each index on its own did not account for the full influence of SES on the associations between other vulnerable subpopulations and exposure to either pollutant. We categorised the specific SES variables into 5 groups based on the aspect of SES they measure as listed in table 3: education, employment, income, housing and population. For each group, we measured the variance in pollutant concentration explained by combinations of SES variables and the multicollinearity between them. We pooled together variables from each group that explained the most variation in pollution concentration with the least multicollinearity. We then tested the multicollinearity and variance explained between different combinations of variables within this larger group. Based on these results, we included at least one variable from each SES group as covariates for the multivariate models to account for the influence of each aspect of SES on pollutant exposure by age and ethnicity (table 4). The only exceptions were SES variables for population, as these had significant problems with multicollinearity in most models.
We calculated standard statistical measures, including p-values, for all models (details of which can be found in section 2 of S1 is available online at stacks.iop.org/ ERL/14/115005/mmedia). As we had a large sample size which would favour statistically significant findings of unknown practical relevance, and as we were more interested in the shape of the associations-particularly non-monotonicity-we instead focused mainly on the plots of the predicted relationship between pollution concentrations and each population variable. The predicted concentrations were calculated by summing the intercept of the model with the fitted smoother values generated in the model for the vulnerable subpopulation predictor. We generated 95% confidence intervals by summing or subtracting 1.96 standard deviations, calculated by the model for the vulnerable subpopulation predictor, from the predicted concentrations. To make the results more intuitive, we truncated our plots at the first and last 25 values of the curve to reduce problems associated with small numbers leading to excessive uncertainty in the estimated curves.
All analyses were conducted in R version 3.4.0 using the spdep (Bivand and Wong 2018) and mgcv (Wood 2017) packages.

Results
Bivariate models-urban Urban locations: The bivariate associations between PM 2.5 and NO 2 concentration with each vulnerable subpopulation are shown in figure 2. All results for specific measures of SES are in chapter S3 section S3.5 (PM 2.5 ) and section 3.6 (NO 2 ). Due to the similarity of our results for both modelled estimates of PM 2.5 , we included the results for one PM 2.5 dataset (SAT-W) in chapter S3 section 3.7, as although the results were very similar with SAT, the overall model fit for GAMs using this dataset were not as well-fitted.
We found positive curvilinear associations between the percentage of people of Indigenous status with the concentrations of both PM 2.5 and NO 2 until the 95th percentile value of percent Indigenous, where the association flattened. We also found positive curvilinear associations between percent minority ethnicity with the concentrations of both PM 2.5 and NO 2. As values of percent Indigenous and percent minority ethnicity increased between the 5th and 95th percentiles, the predicted concentration of PM 2.5 increased by approximately 0.2 and 0.6 μg m −3 and the predicted NO 2 concentration increased 0.5 and 1.5 ppb respectively.
For age, we observed negative curvilinear associations between percent under 15 years and the concentration of both pollutants. Predicted PM 2.5 and NO 2 concentration decreased by 0.25 μg m −3 and 1 ppb, respectively, between the 5th and 95th percentile values of percent under 15 years. Despite the associations being statistically significant (tables S1-6 and S1-7), we only calculated a small u-shaped association between percent over 65 years and the concentration of PM 2.5 , and found no clear overall positive or negative association between percent over 65 years and NO 2 concentration based on our plots.
For SES, between the 5th and 95th percentiles of each SEIFA index, we found negative associations between all SEIFA indices and the concentration of PM 2.5 and NO 2 . The highest predicted concentrations of PM 2.5 and NO 2 were in areas with a low SES, especially around the 5th percentile of SEIFA values. Between the 5th and 95th percentile of SEIFA values, the predicted concentration of PM 2.5 and NO 2 decreased by 0.4-0.6 μg m −3 and 1.5-2 ppb respectively depending on the index.
These results were consistent for specific measures of SES, as we found similar associations between the 5th and 95th percentile values of each measure of SES with the concentrations of both PM 2.5 and NO 2 , (these can be seen in figures S3.2 and S3.4).
Bivariate models-rural Rural locations: the bivariate associations between PM 2.5 and NO 2 concentration with each vulnerable subpopulation are shown in figure 3. We found positive associations between percent Indigenous with the concentration of PM 2.5 and NO 2 between the 5th and 95th percentile values of percent Indigenous, and negative associations after the 95th percentile. Between the 5th and 95th percentile values for percent Indigenous, the concentration of predicted PM 2.5 and NO 2 increased approximately 0.3 μg m −3 and 1 ppb respectively. We calculated positive curvilinear associations between percent minority ethnicity with the concentration of PM 2.5 and NO 2 , with the predicted concentration of PM 2.5 and NO 2 increasing 0.6 μg m −3 and 1.8 ppb respectively as values of percent minority ethnicity rose between the 5th and 95th percentiles.
We found negative curvilinear associations between percent under 15 years and the concentrations of PM 2.5 and NO 2 . As values of percent under 15 years increased between the 5th and 95th percentiles, predicted PM 2.5 concentrations decreased by 0.2 μg m −3 and predicted NO 2 concentrations decreased by 0.8 ppb. We calculated positive curvilinear associations between percent over 65 years and the concentration of both pollutants, with the predicted concentrations of PM 2.5 and NO 2 increasing For SES, we calculated predominately negative associations between all SEIFA indices with PM 2.5 and NO 2 concentration between the 5th and 95th percentile values of each SEIFA index. Low SES areas around the 5th percentile of SEIFA values were exposed to the highest predicted concentrations of both pollutants. High SES areas and the very lowest SES areas-those below the 5th percentile value for each SEIFA index-were exposed to the lowest concentrations. The difference between the lowest and highest predicted concentrations was 0.2-0.4 μg m −3 for PM 2.5 concentration and 1-1.5 ppb for NO 2 depending on the SEIFA index.
These results were consistent for specific measures of SES, as we found similar associations between most specific measures of SES and the concentrations of both pollutants, particularly between the 5th and 95th percentiles of each population variable (figures S3.3 and S3.5).

Multivariable models-urban
Urban locations: figure 4 shows the multivariable associations between PM 2.5 and NO 2 concentration with each vulnerable subpopulation when specific SES variables were included as covariates. Positive curvilinear associations between percent Indigenous and percent minority ethnicity with pollutant concentration were observed up to the 95th percentile of population values, similar to the bivariate models shown in figure 2. However, predicted increases in the concentration of PM 2.5 and NO 2 were considerably smaller than for the bivariate models, and decreases in predicted pollutant concentrations after the 95th percentile of either population variable were much larger. Between the 5th and 95th percentile values of percent Indigenous predicted PM 2.5 and NO 2 concentrations increased only 0.02 μg m −3 and 0.1 ppb respectively, and areas above the 95th percentile value of percent Indigenous experienced the lowest concentrations of both pollutants. As percent minority ethnicity increased, predicted concentrations of PM 2.5 and NO 2 increased 0.35 μg m −3 and 0.8 ppb respectively between the 5th and 95th percentile values of percent minority ethnicity. The concentrations of both pollutants were highest in areas with the highest percentage of people from a non-Indigenous ethnic minority.
For age, we found negative curvilinear associations between percent under 15 years with the concentration of PM 2.5 and NO 2 which were similar to the bivariate models. As values of percent under 15 years increased between the 5th and 95th percentiles, the predicted concentration of PM 2.5 and NO 2 decreased approximately 0.5 μg m −3 and 2 ppb respectively. We found flat associations between percent over 65 years and the concentration of PM 2.5 and NO 2 up to the 95th percentile value for percent over 65 years, which were similar to the bivariate models. However, unlike the bivariate models, after the 95th percentile value for percent over 65 years we found negative curvilinear associations between percent over 65 years and the concentration of PM 2.5 and NO 2 .
Multivariable models-rural Rural locations: figure 5 plots the multivariable associations between PM 2.5 and NO 2 concentration with each  vulnerable subpopulation when specific SES variables were included as covariates. We found that the association between percent Indigenous with PM 2.5 concentration was not statistically significant (table S1-12). We also found that changes in predicted NO 2 concentration with increased percent Indigenous were particularly small, even though between the 5th and 95th percentile values of percent Indigenous we found a positive association between percent Indigenous and the concentration of NO 2 , similar to the bivariate model. Positive curvilinear associations were calculated between percent minority ethnicity with the concentration of PM 2.5 and NO 2 which were similar to the bivariate models, but with smaller increases in predicted pollutant concentrations. As values of percent minority ethnicity increased between the 5th and 95th percentiles, predicted PM 2.5 concentration increased 0.5 μg m −3 and NO 2 concentration increased 1 ppb.
We found a negative curvilinear association between percent under 15 years with NO 2 concentration, which was very similar to the bivariate model. We calculated a flat association between percent under 15 years and PM 2.5 concentration between the 5th and 95th percentile values for percent under 15, which was considerably different from the negative association calculated by the bivariate model. We calculated curvilinear positive associations between percent over 65 years and the concentration of both pollutants, which was similar to the bivariate models except for substantial increases in predicted pollutant concentration observed after the 95th percentile value of percent over 65 years in the multivariable models. Overall, between the minimum and maximum values of percent over 65 years, the predicted concentration of PM 2.5 and NO 2 increased 0.9 μg m −3 and 1.75 ppb respectively.

Discussion
Our study analysed the associations between the concentration of PM 2.5 and NO 2 with the distribution of vulnerable subpopulations in rural and urban areas of Australia. Our analysis found that there were significant associations between some vulnerable subpopulations and the annual average concentration of both pollutants. These associations were complex, with many nonlinear and non-monotonic relationships estimated between population vulnerability and pollutant concentration. Nonetheless, the results of our study suggest that in urban and rural areas, the most vulnerable subpopulations may be exposed to higher concentrations of both pollutants compared with less vulnerable subpopulations.
The results of our analysis showed that people of a non-Indigenous ethnic minority are likely to be exposed to higher concentrations of PM 2.5 and NO 2 in urban and rural areas, compared to the rest of the population. In regard to age, similar results were observed for people over 65 years old living in rural areas, for whom we found evidence of inequitable exposure to PM 2.5 and NO 2 , although we did not find evidence of similar inequalities in urban areas. Further, we found that people under 15 were the least likely to be exposed to high concentrations of either pollutant in urban and rural areas. Areas with greater socio-economic disadvantage were generally found to be exposed to inequitably high concentrations of PM 2.5 and NO 2 , particularly in urban areas. Our results also indicated that Indigenous people were inequitably exposed to PM 2.5 and NO 2 in rural and urban areas, although the predicted disparities in exposure to both pollutants were notably smaller compared to other potentially vulnerable subpopulations. Overall, associations between vulnerable subpopulations and exposure to PM 2.5 and NO 2 were consistent between rural and urban areas, particularly for inequalities in exposure by ethnicity and SES. The results of our multivariable models demonstrated that SES significantly influenced some bivariate associations between censusderived indicators of other vulnerable subpopulations and PM 2.5 and NO 2 concentration. In particular, SES had a significant influence on the disparities in exposure experienced by Indigenous and non-Indigenous ethnic minorities to both pollutants, although non-Indigenous ethnic minorities still experienced relatively large inequalities in exposure when the influence of SES was accounted for. Concerning age, SES significantly affected associations between people under 15 and exposure to PM 2.5 and NO 2 in rural areas, although the initial disparities in exposure were weak and showed that people under 15 were less likely to be exposed to both pollutants compared to the rest of the population. Furthermore, SES had little influence on disparities in exposure experienced by people over 65 years. The results of these models indicate that although SES does not explain all inequitable exposure to air pollution by age and ethnicity, it explains a significant degree of the disparities in exposure experienced by ethnic minorities compared to the rest of the population.
The results of our study are generally consistent with international and (albeit limited) Australian studies on vulnerable subpopulations' exposure to air pollution. Similar to this analysis, existing studies have consistently found that socio-economically disadvantaged people and ethnic minorities, particularly Indigenous communities in Australia, are inequitably Our results showed that inequalities in exposure to air pollution by age may exist in Australia, particularly for elderly people, but that these inequalities are not as pervasive as disparities in exposure by SES and ethnicity. Moreover, existing Australian studies have not investigated and compared inequalities in air pollution exposure between rural and urban areas, with most studies purely examining exposure in urban areas (Knibbs andBarnett 2015, Cowie et al 2016). International studies that have compared inequalities in air pollution exposure between rural and urban areas have typically found that associations between vulnerable subpopulations and air pollutants vary by urbanicity. Similar to age, our results for inequitable exposure to PM 2.5 and NO 2 in urban and rural areas are consistent with the findings of these international studies. In accordance with existing studies, our analysis found some variation in the associations between vulnerable subpopulations and exposure to PM 2.5 and NO 2 by urbanicity, especially for people over 65 who only experienced inequitable exposure in rural areas. At the same time, there was consistent evidence of inequitable exposure to both air pollutants by SES and ethnicity in both urban and rural areas, suggesting that these inequalities are prevalent across Australia. Finally, existing Australian studies that have documented inequitable exposure to air pollution have typically used limited measures of air pollution, either examining exposure to a single pollutant (Knibbs andBarnett 2015, Cooper et al 2018), or a particular source of pollution, specifically anthropogenic sources (Chakraborty and Green 2014, Cowie et al 2016). Our analysis expands on these studies by finding disparities for both NO 2 and PM 2.5 across rural and urban areas. As NO 2 is mostly emitted from traffic while PM 2.5 is emitted from various anthropogenic sources in urban areas and natural sources in rural areas, our results demonstrate that both anthropogenic and natural sources of PM 2.5 may contribute to disparities in exposure experienced by vulnerable subpopulations. Overall, the results of this study are consistent with international research on environmental justice, while expanding on existing Australian research and demonstrating that inequities in air pollution exposure are more pervasive than previously understood.
While our results provide considerable evidence of inequitable exposure to air pollutants, it is not clear whether the disparities in exposure computed in our analysis will translate to disparities in health outcomes. Changes in the predicted concentration of annual average PM 2.5 and NO 2 were small, as all predicted increases in PM 2.5 and NO 2 concentration were less than 1 μg m −3 and 2 ppb respectively for areas with the most vulnerable subpopulations. These differences in predicted pollutant concentration are small compared to the prediction error of the pollution datasets used in our study, as RMSE values for the PM 2.5 and NO 2 LUR models were ∼1.2 μg m −3 and ∼2 ppb respectively when evaluated independently (Knibbs et al , 2018a, suggesting that the results of our analysis may be influenced by error in the initial PM 2.5 and NO 2 data. Nonetheless, these changes in exposure may still be of concern, as studies indicate that negative health outcomes can occur from chronic exposure at concentrations of PM 2.5 and NO 2 that are lower than Australian (and WHO) standards, and suggest that there is no lower threshold concentration of PM 2.5 or NO 2 where no health outcomes occur. In addition, studies suggest that increases in low concentrations of PM 2.5 and NO 2 may intensify negative health outcomes more than equivalent increases at higher concentrations (Burnett et al 2018, Hanigan et al 2019). The health effects from exposure to PM 2.5 and NO 2 may also be confounded or, conversely, compounded by the presence of other pollutants such as heavy metals in PM 2.5 or ozone, as well as the presence of each other. Previous studies have found that the overlap between PM 2.5 and NO 2 in Australia is quite weak overall, with previous LUR models estimating weak to fair correlations between estimates of both pollutants nationally (r: 0.09-0.5) (Knibbs et al 2018a). However, it is possible that there is some overlap of both pollutants in major urban centres for this study, as these are the areas most likely to experience the highest concentrations of both PM 2.5 and NO 2 (Knibbs et al , 2018a, thereby exacerbating the potential for disparate health impacts from the inequalities measured in urban areas for this study. Vulnerable subpopulations are more susceptible to experiencing negative health outcomes from exposure to air pollution, making small differences in pollution concentration more likely to cause disparities in negative health outcomes between these subpopulations and the rest of the population as well. Furthermore, it is possible that disparities in exposure to air pollutants measured in this study could be greater at smaller scales, as local variations in exposure between subpopulations may be attenuated when analysed at large scales (Baden et al 2007). Locations with relatively high concentrations of air pollution, such as around major industrial areas and motorways in urban areas or within mining cities in rural areas, have the most potential for inequalities in air pollution exposure to be obscured by analyses at a national level. For these reasons, our results, as well as previous studies in Australia, indicate that there is a need for further research to quantify the health effects from disparate exposure to air pollution. The results of this and previous studies also suggest that public health policy could consider vulnerable subpopulations' exposure to air pollution independently of overall health impacts of air pollution on the entire population.
To inform policies addressing inequalities in exposure to air pollution, studies have traditionally focused on the influence of societal factors on exposure, particularly SES, as they have been found to be a frequent cause of inequitable exposure experienced by vulnerable subpopulations. The results of our multivariable analysis were consistent with these studies, indicating that public health policy could focus on SES to address disparities in air pollution exposure, due to its influence on ethnic minorities' exposure. However, how policy should be directed to address this depends on the sources of pollution causing disparities in exposure, especially as the major sources of PM 2.5 and NO 2 are a mix of anthropogenic and natural sources which are likely to require significantly different policy responses. Anthropogenic sources, such as traffic and industry which are primarily concentrated in urban areas, can be regulated directly through environmental and planning policy. Examples of relevant policy measures include setting limits on acceptable levels of emissions in the licensing agreements of industrial sources (Dobbie and Green 2015) and legislating fuel efficiency standards (Smit 2019). In contrast, emissions from natural sources of PM 2.5 , such as bush fires and dust storms in rural areas, cannot be regulated directly via policy due to the large number of environmental and meteorological factors that affect their occurrence, location and intensity. Emissions from these natural sources are projected to potentially increase due to our changes in Australia's climate over time (Bradstock 2010, Goudie 2014. Consequently, all else being equal, disparities in exposure to PM 2.5 in rural areas (and some urban areas) may worsen, further exacerbating existing disparities in health outcomes from pollution exposure. The results of our study indicate the need for further public health and environmental research to assess and quantify disparities in vulnerable subpopulations' exposure to PM 2.5 and NO 2 , their health impacts and the causes of inequalities in exposure. Moreover, the results of our study indicate that policy responses may be necessary to address vulnerable subpopulations' exposure to air pollution and that these policy responses would need to be nuanced and multi-faceted to address the different sources of pollution. Our results show that particular attention should be given to the influence of SES in causing disparities in pollution exposure, but that this must be considered in the context of the nature and sources of emissions. In particular, when addressing natural sources of emissions, the influence of climate change on these sources should be accounted for to address long-term disparities in exposure.

Limitations
A possible limitation of our analysis is that we did not compare the results of models that accounted for spatial autocorrelation differently. We specifically presented the models that were the most accurate for our analysis based on their model fit and reduction in spatial autocorrelation, as shown in section 3 of S1. Nonetheless, associations between census-derived vulnerable subpopulation variables and pollutant concentration can vary significantly-including changes in statistical significance, predicted differences in exposure, and even changes the direction of associations-depending on how comprehensive the smoothing spline is for geographic coordinates in the model. This is shown by the models in S5, which used smoother splines with far less knots (k=5) for geographic location compared to the models in our analysis (k=25-50) and that accounted for spatial autocorrelation less extensively. This ultimately meant that the results of the analysis depend heavily on the number of knots selected for the geographic location smooth spline, which was guided by the authors' prior experience.
Another limitation is that the predicted differences in annual average concentration for PM 2.5 and NO 2 in our analysis were smaller than the independent validation prediction error-measured as absolute RMSEof the LUR models used in our study (Knibbs et al , 2018a. Additionally, the accuracy of pollutant estimates in rural areas are more difficult to validate due to the lack of monitoring data, potentially leading to greater error than was measured or accounted for in the models. Consequently, the predicted increases in PM 2.5 and NO 2 concentration calculated in our analysis, which are smaller than the prediction error, need to be interpreted cautiously.
The use of annual average concentration estimates from LURs to quantify exposure may not represent the effects of acute exposure from important sources like fires and dust storms in Australia. Values for annual average concentration were the most temporally refined pollutant estimates available for this analysis, and are an adequate measure of overall exposure between subpopulations, particularly as chronic exposure to PM 2.5 and NO 2 has well-documented health effects. Nonetheless, the use of more temporally refined data would allow the models to more accurately quantify disparities in exposure to acute emissions.
The large scale of our analysis meant that only disparities in exposure experienced by broad measures of vulnerable subpopulations and the causes of inequitable exposure, specifically SES, could be investigated in this study. Analysis of other possible factors explaining disparities in exposure by age or ethnicity-such as cultural or political disadvantage-require more refined data on vulnerable subpopulations and, therefore, were beyond the scope of this study. Future studies at finer spatial scales could use more refined data on vulnerable subpopulations to investigate the influence of these factors on inequalities exposure to air pollution in Australia.
Finally, due to the ecological fallacy, which is possible when applying this type of analysis, the use of aggregated pollutant concentration data and census data may not be representative of individual exposure to outdoor air pollution. Personal exposure to air pollutants may be heavily affected by the mobility patterns of individuals and the quality of housing they live in, which were not investigated in this analysis. Previous studies have found that SES has a significant influence on housing quality (Evans and Kantrowitz 2002, Cooper et al 2018) and patterns of behaviour (Myers 2009). These studies provide some evidence indicating that the influence of these factors may accentuate socio-economically disadvantaged people's exposure to air pollution beyond what is measured with geospatial measures of pollutant concentration. It is possible that these patterns may extend to ethnic minorities, who are often more likely to experience socio-economic disadvantage compared to the rest of the population. Nonetheless, as there is limited research of similar patterns by age or ethnicity in Australia, it is unknown if accounting for housing quality and mobility patterns would significantly compound geospatial measures of exposure to PM 2.5 and NO 2 experienced by these subpopulations. As a result, we have attempted to minimise the limitations of a geospatial analysis by using estimates of PM 2.5 and NO 2 concentration at the mesh block scale, the smallest census area units available, before populationweighted averaging the values to the next largest census area unit to represent average exposure for our analysis.

Conclusion
Our study found that vulnerable subpopulations are likely to be exposed to higher concentrations of the hazardous air pollutants PM 2.5 and NO 2 in urban and rural areas of Australia, compared to the rest of the population. Specifically, ethnic minorities, socio-economically disadvantaged people living in rural and urban areas, and elderly people (only in rural areas), experienced the greatest disparities in exposure to these pollutants. SES was found to be a major cause of these inequalities in exposure. These results are consistent with previous research on vulnerable subpopulations' exposure to air pollution in Australia. Whether this disparity in exposure for these vulnerable subpopulations would translate to similar disparities in negative health outcomes is harder to identify. However, these findings suggest that this issue should be investigated further by researchers to properly understand the public health burden of inequitable exposure, and should be considered by policy makers when assessing the impacts of air pollution. Further, the consistency of our results across pollutants and across urban and rural areas indicate that anthropogenic and natural sources of pollution are worthy of policy attention. As climate change directly and indirectly impacts air pollution (through increasing fire weather and increasing the need for hazard reduction burns), these issues may become more pressing in the near future.