## Abstract

**Background** Numerous epidemiological studies have documented the adverse health impact of long-term exposure to fine particulate matter (PM_{2.5}) on mortality even at relatively low levels. However, methodological challenges remain to consider potential regulatory intervention’s complexity and provide actionable evidence on the predicted benefits of interventions. We propose the parametric g-computation as an alternative analytical approach to such challenges.

**Method** We applied the parametric g-computation to estimate the cumulative risks of non-accidental death under different hypothetical intervention strategies targeting long-term exposure to PM_{2.5} in the Canadian Community Health Survey cohort from 2005 to 2015. On both relative and absolute scales, we explored benefits of hypothetical intervention strategies compared to the natural course that 1) set the simulated exposure value at each follow-up year to a threshold value if exposure was above the threshold (8.8 µg/m^{3}, 7.04 µg/m^{3}, 5 µg/m^{3}, and 4 µg/m^{3}); and 2) reduce the simulated exposure value by a percentage (5% and 10%) at each follow-up year. We used the three-year average PM_{2.5} concentration with one-year lag at the postal code of respondents’ annual mailing addresses as their long-term exposure to PM_{2.5}. We considered baseline and time-varying confounders including demographics, behavior characteristics, income level, and neighborhood socioeconomic status. We also included the R syntax for reproducibility and replication.

**Results** All hypothetical intervention strategies explored led to lower 11-year cumulative mortality risks than the estimated value under natural course without intervention, with the smallest reduction of 0.20 per 1000 participants (95% CI: 0.06 to 0.34) under the threshold of 8.8 µg/m^{3}, and the largest reduction of 3.40 per 1000 participants (95% CI: -0.23 to 7.03) under the relative reduction of 10% per interval. The reductions in cumulative risk, or numbers of deaths that would have been prevented if the intervention was employed instead of maintaining status quo, increased over time but flattened towards the end of follow-up. Estimates among those ≥65 years were greater with a similar pattern. Our estimates were robust to different model specifications.

**Discussion** We found evidence that any intervention further reducing the long-term exposure to PM_{2.5} would reduce the cumulative mortality risk, with greater benefits in the older population, even in a population already exposed to low levels of ambient PM_{2.5}. The parametric g-computation used in this study provides flexibilities in simulating real world interventions, accommodates time-varying exposure and confounders, and estimates adjusted survival curves with clearer interpretation and more information than a single hazard ratio, making it a valuable analytical alternative in air pollution epidemiological research.

## Introduction

As collective efforts in previous decades have successfully reduced the level of fine particulate matter (PM_{2.5}) globally, quantifying the effectiveness of policies that further reduce ambient PM_{2.5} is becoming particularly important in supporting evidence-based policymaking. Indeed, previous studies found consistent evidence of deleterious associations between long-term exposure to low levels of PM_{2.5} (e.g., below the current health-based standards or guidelines) and risk of mortality^{1–6} and morbidity,^{7–9} suggesting potential reductions in health burden if the PM_{2.5} level were to be further reduced. While the evaluation of exposure-response functions in existing studies provides important information in understanding the potential effectiveness of policies, further methodological considerations are required to estimate the potential benefits of realistic interventions.

First, evidence suggested that the risk associated with the changes in acute exposure to PM_{2.5} could vary with time,^{10–13} potentially due to changes in chemical compositions of PM_{2.5} with different toxicity and population susceptibility towards PM_{2.5}.^{14,15} Similar disparity in toxicity across long-term exposure to PM_{2.5} components was also observed,^{16,17} suggesting that such temporal changes could exist in risk associated with long-term exposure to PM_{2.5}. In other words, it is important to use analytical methods flexible enough to incorporate such temporal changes in estimation of related health burdens. However, existing studies of health impacts of long-term exposure to PM_{2.5} generally considered time-fixed exposure and confounders (see Table S1 for a narrative review of recent studies on health impact of long-term exposure to PM_{2.5} and their methodological considerations). Furthermore, the most widely used estimate for exposure-response function in this field is a single hazard ratio (HR) for the follow-up period estimated with a standard Cox-proportional hazard model (Table S1), which is assumed to be constant over time and precludes consideration of temporal changes. Although extension of a Cox-proportional hazard model could provide period-specific HRs that incorporate temporal changes,^{18} recent developments in causal inference literature raise concern about the ambiguity in the causal interpretation of HR and period-specific HRs.^{19} Specifically, period-specific HRs have a built-in selection bias because susceptible people exposed to higher PM_{2.5} are more likely to die early if PM_{2.5} truly increases risk of mortality, and are removed from the susceptible population in later time.^{20} This differential depletion of susceptible subjects over time can lead to artificially diminished or even reversed period-specific HR later in study even when the cumulative survival is still lower among those exposed to higher PM_{2.5}, violating the proportional assumption and hindering interpretation.^{21}

Second, calculation of health burden related to long-term exposure to PM_{2.5} commonly employed exposure-response function previously estimated with the static intervention strategy, where a fixed change of exposure value was assigned to the entire population.^{22(chap19)} However, the more flexible and realistic dynamic intervention strategy, where the exposure value was assigned based on individuals’ history of covariates including exposure, is hard to apply when existing exposure-response functions are used.^{22(chap19)} Methods capable of incorporating dynamic intervention strategy to imitate complexities in actual regulatory interventions are needed to provide direct evidence on effectiveness of air pollution control policies.^{23} To fill this gap in knowledge translation, we propose the parametric g-computation as an analytical alternative in air pollution epidemiological research, which could better predict the effectiveness of hypothetical policies while being more flexible in resembling real world interventions.

G-computation (also known as g-formula) is a generalization of non-parametric standardization developed under the potential outcome framework for causal inference,^{24} and parametric g-computation is a variation that employs parametric modeling. Under the consistency (i.e., the exposure is defined unambiguously, and all exposed individuals receive the same version of treatment),^{22(chap3),25} exchangeability (i.e., no unmeasured confounding or informative censoring),^{25} and positivity (i.e., probability of receiving every exposure conditioning on confounders is greater than zero) assumptions,^{22(chap3)} and a time-to-event outcome setting, g-computation can provide marginal causal risk estimates at each follow-up time point under hypothetical intervention strategies (i.e. adjusted survival curves), while allowing other population characteristics to be altered according to the intervention.^{26} Particularly, parametric g-computation excels in estimating adjusted survival curves under dynamic intervention strategies. In other words, g-computation can directly answer causal questions such as: how many lives could we save if we promulgate a policy that further reduces air pollution to levels lower than the current standard among those whose exposure were above the current standard, compared to maintaining the status quo? Although parametric g-computation has been widely applied in other fields of epidemiology,^{27–30} application in air pollution studies remains limited. Previous applications in this field either focused on a small cohort in occupational settings,^{31–33} or modelled simple air pollution changes on asthmatic outcomes among children (i.e. not considering time-varying confounding nor changes in effect estimates over time).^{34,35}

In this study, we aim to demonstrate the use of parametric g-computation to evaluate the effectiveness of hypothetical intervention strategies targeting long-term exposure to PM_{2.5} on reducing mortality using a Canadian cohort experiencing low PM_{2.5} exposure from 2005 to 2015. This analytical alternative can account for previously unaddressed complexities, refine the effect estimates with less restrictive identification conditions and provide estimates more intuitive to policy makers.

## Methods

### Study population

We created a retrospective cohort with respondents to the Canadian Community Health Survey (CCHS) from three enrolling cycles in the years of 2000/2001, 2003 and 2005, respectively.^{36–38} CCHS is a national cross-sectional survey collecting health status, heath care utilization and health determinants information of the Canadian population, covering the population 12 years and over in the ten provinces and the three territories. The survey excluded individuals living on reserves and other Aboriginal settlements, full-time members of the Canadian Forces, the institutionalized population, and residents of certain remote regions.

Among CCHS respondents who gave permission to share and link their information with other administrative datasets, we obtained their mailing address history and death records through December 31^{st}, 2015 via Statistics Canada’s Social Data Linkage Environment, using probabilistic methods based on common identifiers.^{2,39} We focused on non-accidental death as outcome (International Classification of Diseases ninth revision codes 001 to 799 and International Classification of Diseases tenth revision codes A to R) in this study. To facilitate pooling of results across cycles, we restricted the cohort to participants who were alive on January 1^{st}, 2005 and used this date as the start of follow-up for all cycles. We also restricted our cohort to individuals older than 25 years and younger than 80 years in 2005 thus all cohort participants were adults and followed for 11 years or till death. Besides, we dropped respondents without data for covariates including exposure in 2005. This study was approved by the Health Canada-Public Health Agency of Canada Research Ethics Board.

### Exposure assessment

To estimate respondents’ long-term exposure to PM_{2.5}, we utilized the ground-level PM_{2.5} concentrations from V4.NA.02.MAPLE of the Atmospheric Composition Analysis Group of Washington University,^{40} which covers all of North America below 70°N. The 0.01° × 0.01° (roughly equivalent to 1×1km^{2} at the latitudes where most Canadians live) annual estimates of PM_{2.5} from 2001 to 2015 were derived using satellite retrievals of aerosol optical depth and chemical transport model simulations, and calibrated with ground-based observations using geographically weighted regression.^{41} The annual estimates of PM_{2.5} closely agree with long-term cross-validated ground measurements at fixed-site monitors (n=2,312) across North America (R^{2}=0.70).^{41} Using the ground-level PM_{2.5} concentration surfaces described above, we first assigned the annual PM_{2.5} concentration of the grid cell into which the postal code centroid falls as the postal code specific annual PM_{2.5} concentrations. Then we calculated respondents’ annual long-term exposure to PM_{2.5} as three-year average postal code specific PM_{2.5} concentrations with one-year lag based on their mailing address history (e.g., a respondent’s long-term exposure to PM_{2.5} in 2013 is the average of their postal code specific PM_{2.5} concentrations in 2010, 2011 and 2012).^{2} We utilized three-year average with one-year lag to represent long-term exposure of PM_{2.5} so that the exposure always precedes the outcome and the timeframe is consistent with the regulatory review of Canadian Ambient Air Quality Standards for annual PM_{2.5}.^{42} This metric of long-term exposure to PM_{2.5} was widely used in previous Canada based studies of long-term health impacts of PM_{2.5}.^{2,43,44}

### Covariates other than exposure

In this section we summarized the data sources and meaning of covariates in this study while the covariate selection to control for in our model will be discussed in the statistical analysis section. We used covariates to describe the collection of exposure, time-fixed confounders, and time-varying confounders in this study. Baseline characteristics of respondents were ascertained at the time of enrollment into CCHS via self-report and were processed using the same method as previous studies,^{2,43} including sex, age (converted to value in 2005), body mass index, marital status, immigrant status, visible minority, indigenous status, smoking status, alcohol consumption, consumption of fruits and vegetables, leisure physical activity, working status, and educational attainment (details of variable categorization in Table 1). By using 2005 as the start of follow-up time for all individuals, we assumed that all baseline characteristics other than age ascertained at the time of enrollment would remain the same through the entire follow-up period.

We also obtained characteristics of the respondents and their neighborhoods through linkage with administrative datasets using similar methods as previous studies.^{2,43} Specifically, we obtained annual income quintile of respondents through linkage with tax data based on common identifiers.^{43} For person-years with missing annual family income, we imputed them with the nearest prior values and the proportions of missing are 5.21%, 4.97% and 4.69% for Cycle 2000/2001, 2003 and 2005, respectively. We also obtained annual characteristics of neighborhoods through linkage with census data from the nearest census year based on respondents’ mailing address postal codes, including community size at census metropolitan area level and four Canadian Marginalization Index at census dissemination area level. The Canadian Marginalization Index summarizes dissemination area-level socioeconomic status into four dimensions using principal component analysis to reduce the dimensionality of census data: the immigration and visible minority index combines information on proportion of recent immigrants and proportion of people self-identify as visible minority; the households and dwellings index combines information on types and density of residential accommodations and family structure characteristics; the material resources index combines information on access to and attainment of basic material needs; and the age and labour force index combines information on participation in labour force and proportion of seniors.^{45} Last, we obtained airshed (six distinct regions of Canada that cut cross jurisdictional boundaries and showed similar air quality characteristics and air movement patterns within each region) to capture large scale spatial variation,^{46} and urban form information of respondents’ neighborhoods in 2005 to capture urbanicity of participants’ residence, through linkage with census data.^{2}

### Hypothetical intervention strategies

In this study, we explored three types of intervention strategies: 1) applying the simulated value of time-varying covariates without any intervention (natural course); 2) setting the simulated long-term exposure to PM_{2.5} value at each follow-up year to a threshold value if PM_{2.5} was higher than the threshold (threshold intervention); and 3) reducing the simulated PM_{2.5} value by a fixed percentage at each interval (i.e., follow-up year) (relative reduction intervention).

Threshold values explored are the current Canadian Ambient Air Quality Standards for PM_{2.5} of 8.8 µg/m^{3}, 80% of the current Canadian Ambient Air Quality Standards for PM_{2.5} (or 7.04 µg/m^{3}), the new World Health Organization (WHO) air quality guideline of 5 μg/m^{3}, and a PM_{2.5} level that was further below the WHO guideline (4 µg/m^{3}). Interval-specific relative reduction values explored are 10% and 5% per interval. To avoid extensive model extrapolation, we restricted the relative reduction intervention so that subjects with exposure below 1.8 µg/m^{3}, the background PM_{2.5} level in Canada,^{47} will not be further reduced. The first type of intervention strategy represents the predicted covariates based on the observed data without intervening and serves as the reference for other strategies. The second and the third are dynamic intervention strategies that incorporate the impact of intervention on covariates during earlier time points while simulating covariates in later time points.

### Statistical analysis

We applied parametric g-computation with different hypothetical intervention strategies targeting long-term exposure to PM_{2.5} to understand the benefits of intervention strategies on cumulative risk of non-accidental death. We conducted g-computation analysis for each enrollment cycle separately and pooled the results across cycles using meta-regression. Briefly, we estimated the cumulative mortality risk at each follow-up year standardized to the distribution of the confounders and long-term exposure to PM_{2.5} in the study population, with all time-varying covariates (confounders and PM_{2.5}) conditioned on covariates history, with and without intervention on PM_{2.5} (i.e., adjusted survival curves). Next, we calculated the differences in cumulative mortality risks between the natural course and other intervention strategies on both absolute and relative scales to provide estimates for the benefits of hypothetical intervention strategies compared to maintaining status quo. We pooled results with fixed-effect meta-regression, which calculates a weighted average of cycle specific estimates with weights equal to the inverse of the variance using the “meta” package.^{48}

The proof of parametric g-computation are described extensively elsewhere,^{22(chap21),29} and detailed description of how to implement such an approach in a setting similar to our study was previously published,^{28} with available R package for easy implementation.^{49} However, since the application of parametric g-computation is limited in air pollution studies, we include a diagram (Figure 1) to summarize the four steps that carry out the g-computation in a time-to-event setting with time-varying exposure and confounders, and describe the steps in details below.

### Steps to implement parametric g-computation

In Step 1, we fitted a pooled logistic model (i.e., discrete-time hazard model) and adjusted for baseline characteristics, time-varying characteristics, quadratic function of year and interaction between long-term exposure to PM_{2.5} and categorical year. The pooled logistic model estimated the probability of death during the year conditioning on survival till the start of the year given all covariates (including PM_{2.5}), which allowed the conditional probability of death and its association with PM_{2.5} to vary over year. We chose confounders to control for in the outcome model based on substantive knowledge of the relationship between long-term PM_{2.5} and mortality as summarized in the simplified directed acyclic graph (Figure S1). We included a full list of covariates in Table 1 with specific forms of covariates in Table 2. We included both individual socioeconomic status indicators (e.g., education and family income) and community socioeconomic status indicators (e.g., Canadian Marginalization Index for dissemination area) to fully capture the variation in socioeconomic status among cohort participants, which is a major source of residual confounding. We also included individual behavior indicators like dietary and exercise patterns, which are strong risk factors for mortality, precede the exposure, and might share common unmeasured causes with the exposure, even though they might not directly cause the exposure.^{50} To note, in the setting when only time-fixed covariates were used, we could estimate marginal adjusted survival curves directly using outputs from this pooled logistic model by predicting the probability of death standardized to the distributions of covariates under the intervention of interest (e.g., setting the baseline level of exposure to a specific value while keeping all baseline covariates the same as observed for all participants).^{19,51} However, in our study setting of time-varying covariates and time-to-event outcome, we also need to model time-varying covariates (including PM_{2.5}) so that we could simulate time-varying covariates at all follow-up years for all participants, especially for periods after participants’ death.^{28,29}

In Step 2, we modeled the time-varying covariates (including PM_{2.5}) using linear regressions while including variables such as previous-year value of the covariate of interest, baseline characteristics, same-year values of time-varying covariates set to occur before the covariate of interest, and quadratic function of time. The choice of independent variables in covariate models are based on substantive knowledge as summarized in the simplified directed acyclic graph (Figure S1). We summarized the list of all covariates in Table 1 and the specific forms of covariates included in covariate models in Table 2. We set the sequence of time-varying covariates as community size, income, immigration and visible minority, material resources, households and dwellings, age and labour force and long-term exposure to PM_{2.5}. Since previous studies using different cycles of CCHS found a supra-linear PM_{2.5}-mortality association,^{2,43,52,53} we used natural logarithm transformed long-term exposure to PM_{2.5} as the independent variable in both the outcome and covariate model in the main analysis.

In Step 3, we simulated new datasets based on the intervention strategies. For each intervention, we randomly sampled 10,000 subjects from the cohort with replacement (i.e., Monte Carlo sampling) and created an empty dataset of all sampled subjects for all follow-up years till the end of period of interest (normally the end of follow-up as in this study but extrapolation is possible with extra assumptions). We only simulated new datasets for 10,000 subjects instead of the number of participants in study cohorts (~60,000 participants in each cohort) to save computation time and similar practice was conducted before with smaller cohort.^{29} Next, we assigned the baseline values of all covariates (values of baseline covariates and values of time-varying covariates at start of follow-up) in each simulated dataset to the same as its original dataset, then altered the relevant covariate values based on the intervention strategy (e.g., setting the baseline long-term exposure to PM_{2.5} to 5 μg/m^{3} if it is higher than 5 μg/m^{3} in the threshold intervention of 5 μg/m^{3} but could include other covariates if needed). Last, we simulated time-varying covariates at each year after baseline based on their history with covariate models estimated in the second step and altered the covariates based on the intervention strategy.

In Step 4, with the simulated datasets and outcome model from the first step, we calculated for each subject the probability of dying during each year conditioning on surviving to the beginning of this year, standardized to the distribution of the confounders and long-term exposure to PM_{2.5} under the intervention strategies, regardless of their observed outcomes. Next, we calculated for each subject the cumulative mortality risk at each year as the cumulative sum of the abovementioned conditional probability of mortality times the probability of surviving till the beginning of the time interval. The estimated cumulative morality risk at each year is the average of estimates from all subjects for all hypothetical interventions. We also calculated the absolute difference in cumulative morality risk and percentage change in cumulative morality risk with estimated cumulative morality risk from natural course as reference.

Besides, we calculated the 95% confidence intervals (CI) for all estimates using standard errors from 200 bootstrap iterations.^{54} In each iteration, we randomly sampled the same number of participants as in the original cohort with replacement and ran the four steps described above to calculate cumulative mortality risks under intervention strategies. We chose this number of iteration because we were constrained by available computational resources (>1h of computational time for each bootstrap iteration), and the original application of parametric g-computation in time-varying covariates and time-to-event setting used the same number.^{29} Future studies with more computational resources might consider larger number of bootstrap iteration.

### Sensitivity analyses

To test the robustness of our results to model misspecification, we considered a number of different model specifications for both outcome and covariate models including: 1) reordering the sequence of time-varying covariates in covariate models by moving age and labour force to before the other Canadian Marginalization Index, moving income to after Canadian Marginalization Index, and moving PM_{2.5} to the first place among all covariates; 2) extending the extent of history modeled by including previous year and two-year previous values of all the time-varying covariates in the covariate models other than just the previous-year value of the covariate of interest; and 3) including time-varying covariates other than long-term PM_{2.5} as categorical in the outcome model and using the multinomial logistic model for them in covariate model instead of modeling them as continuous with bounds using linear model (see Table 2 for details of model specifications for each time-varying covariates in main analysis). We also visually evaluated the simulated and observed adjusted survival curves and histories of covariates under no intervention in the main analysis as a heuristic check of model misspecification.^{27}

Next, to facilitate comparison with previous studies, we used long-term PM_{2.5} in original scale in all models as sensitivity analysis, which assumed the same log-linear PM_{2.5}-mortality association used in other cohorts^{4,7} instead of the supra-linear one supported by previous studies of different cycles of the CCHS cohort.^{2,43,52,53} Besides, we also ran a Cox-proportional hazard model with the same specification as the outcome model in our main analysis except that we included no time variable and used long-term PM_{2.5} in original scale, which assumed a log-linear PM_{2.5}-mortality association.

Last, since most deaths occurred among older individuals and age could modify the health impact of long-term exposure to PM_{2.5}, we conducted a subset analysis restricted to cohort participants aged ≥65 years at the time of enrollment. Since it took up to 24 hours to run one round of sensitivity analysis without bootstrapping, we were unable to perform bootstrapping to calculate CIs for sensitivity analyses due to computational constraints and did not present variances for our estimates. We pooled cycle-specific estimates from sensitivity analyses by averaging the estimates in each cycle. All analyses were done in R version 4.0.5^{55} with the “gfoRmula” package.^{49} The R code to replicate these analyses and a simulated dataset are available at the following link: https://github.com/suthlam/cchs_g_computation.git.

## Results

We observed 6,475 (10.4%), 6,525 (10.5%), and 6,135 (9.2%) non-accidental deaths in the 11 years of follow-up starting from 2005 among the three cycles of CCHS cohorts of 62,365, 62,380, and 66,385 participants, respectively (Table 1). The three cycle cohorts were comparable in all descriptive statistics (Table 1). Without any hypothetical intervention, the observed average long-term exposure to PM_{2.5} in three cycles of our cohort decreased slightly from 6.4 ± 2.2 µg/m^{3}, 6.5 ± 2.3 µg/m^{3}, 6.5 ± 2.3 µg/m^{3} in 2005 to 5.8 ± 2.0 µg/m^{3}, 6.0 ± 2.0 µg/m^{3}, and 6.0 ± 2.0 µg/m^{3} in 2015, respectively (Table 1).

All hypothetical intervention strategies explored in this study led to lower 11-year cumulative mortality risks than the estimated value under a natural course without intervention, 102.5 per 1000 participants (95% confidence interval (CI): 100.3 to 104.8 per 1000 participants). The reductions in 11-year cumulative mortality risks from the natural course were 0.20 per 1000 participants (95% CI: 0.06 to 0.34 per 1000 participants) under the threshold of 8.8 µg/m^{3}, 0.63 per 1000 participants (95% CI: 0.18 to 1.07 per 1000 participants) under the threshold of 7.04 µg/m^{3}, 1.87 per 1000 participants (95% CI: 0.53 to 3.21 per 1000 participants) under the threshold of 5 µg/m^{3}, 3.08 per 1000 participants (95% CI: 0.85 to 5.31 per 1000 participants) under the threshold of 4 µg/m^{3}, 1.68 per 1000 participants (95% CI: -0.15 to 3.51 per 1000 participants) under the relative reduction of 5% per interval, and 3.40 per 1000 participants (95% CI: -0.23 to 7.03 per 1000 participants) under the relative reduction of 10% per interval. To note, the reduction in 11-year cumulative mortality risks could also be interpreted as the number of deaths that would have been prevented if the intervention was employed instead of maintaining status quo. Changes in relative scale showed similar pattern (Table 3). To fulfill the four threshold intervention strategies, averages of 18.7%, 38.3%, 72.0% and 91.4% of subjects experienced change in their natural course exposure every year, respectively, while 100% had their exposure changed under the two relative reduction intervention strategies (Table 3). The corresponding reductions in average simulated PM_{2.5} from the start of follow-up to the end of year 11 ranged from 0.13 to 1.87 µg/m^{3} for threshold intervention strategies, and 1.25 to 2.18 µg/m^{3} for relative reduction intervention strategies (Table 3). Across all strategies, we observed steady expansions in differences of yearly cumulative mortality risks between the natural course and other strategies until the 7^{th} year of follow-up, after which the differences remain constant and shrink during the last year of follow-up (Figure 2 with numeric results in Table S2). In the main analysis, we pooled estimates of yearly cumulative mortality risks across cycles using random-effect meta-regression and pooled estimates of differences (absolute and relative scale) in cumulative mortality risks using fixed-effect meta-regression. Cycle-specific results with corresponding I^{2} values are summarized in Figure S2 with numeric results in Table S3.

Heuristic checks of model fitting in the main analysis support the robustness of our estimates: 1) the cumulative mortality risk estimated by parametric g-computation under the natural course closely tracked the value observed (Figure S3); and 2) the observed mean values of all time-varying covariates were similar to those simulated under the natural course over time (Figure S3). To note, since participants had no time-varying covariates recorded after their death while we simulated participants’ time-varying covariates for all years, differences between observed and simulated covariates are expected, especially later in the study period. Furthermore, sensitivity analyses with different model specifications (different sequence of time-varying covariate, extent of history modeled, and parametrization of time-varying confounders) resulted in similar estimates as the main analysis (Figure 3, numeric results in Table S4).

When assuming a log-linear PM_{2.5}-mortality association in the sensitivity analysis (compared to the supra-linear association assumed in main analysis by log-transforming the exposure), reductions in 11-year cumulative mortality risks comparing other intervention strategies to the natural course ranged from 0.01 per 1000 participants to 1.65 per 1000 participants, slightly smaller than main analysis assuming a supra-linear PM_{2.5}-mortality association (log-transformed PM_{2.5} used as exposure in modeling) (Table S4). The Cox model assuming a log-linear PM_{2.5}-mortality association found 15.6% (95% CI: 4.0 to 28.5%) increase in hazard of death per 10 µg/m^{3} increase in PM_{2.5}. When focusing on cohort participants ≥65 years at the start of follow-up, reductions in 11-year cumulative mortality risks comparing other intervention strategies to the natural course ranged from 0.49 per 1000 participants to 7.07 per 1000 participants (Table S4), which is larger than the main analysis using the general population ≥35 years.

## Discussion

In this study, we applied the parametric g-computation as an analytical alternative that is particularly valuable for air pollution epidemiological research, especially for evaluating specific intervention strategies. With application in a large Canadian cohort, we demonstrated how to incorporate consideration of complex time structure in the data and how to calculate causally interpretable cumulative risk estimates over the follow-up time (i.e., adjusted survival curves) with parametric g-computation. We described that any intervention further reducing the long-term exposure to PM_{2.5} would reduce the cumulative mortality risk, even in a region with relatively low levels of ambient PM_{2.5}. Such reduction in cumulative risk increased over time and flattened towards the end of follow-up on both relative and absolute scales. The older population also experienced greater benefits from the explored hypothetical intervention strategies than the general population.

Numerous observational studies found positive associations between long-term exposure to PM_{2.5} and non-accidental mortality. A meta-analysis reported a pooled effect estimate of 6% (95% CI: 4 to 8%) increase in hazard of death per 10 µg/m^{3} increase in PM_{2.5} (HR-1).^{5} A recent study in a similar Canadian cohort from 2000 to 2012 found 11% (95% CI: 4 to 18%) increase in hazard of non-accidental death per 10 µg/m^{3} increase in chronic exposure to PM_{2.5} with a Cox proportional hazard model.^{2} Our sensitivity analysis using Cox model without time-varying coefficients found similar numeric results [15.6% (95% CI: 4.0 to 28.5%)]. Although we can’t directly compare our estimates from the main analysis to previous results given the difference in interventions explored, the consistent reductions in cumulative mortality risk over follow-up time across intervention strategies when compared to natural course in this study lend further support to previous findings that PM_{2.5} is detrimental to health even at levels below current standards.

For example, we identified a 0.19% (95% CI: 0.05 to 0.33%) decrease in 11-year cumulative mortality risk comparing the hypothetical intervention strategy with threshold of 8.8 µg/m^{3} to natural course, which provided evidence of health benefits from policies that further reduce the air pollution level to below current Canada standard of 8.8 µg/m^{3}, which is lower than the 12 µg/m^{3} standard of U.S. explored by previous studies.^{4,9} To facilitate comparison with previous studies assuming a log-linear PM_{2.5}-mortality association, we included sensitivity analysis using PM_{2.5} on the normal scale and found reduced cumulative mortality risks in all hypothetical interventions compared to maintaining status quo but the numeric values are smaller than those in the main analysis. The observed difference in the numeric values of analysis assuming log-linear association and analysis assuming supra-linear association is a combination of difference in how the exposure-response relationship is modeled and how the exposure model performs. However, given the existing evidence in Canadian cohorts and similarity between observed survival curve and estimated survival curve using parametric g-computation under no intervention in the main analysis,^{2,43,52,53} we have confidence in the validity of results assuming a supra-linear association.

More importantly, we demonstrated how to incorporate more flexibilities in simulating real world interventions with g-computation in this study and provide intuitive estimates for benefits of such interventions. Taking the hypothetical intervention strategy with threshold of the current Canadian Ambient Air Quality Standards as an example, the average long-term exposure to PM_{2.5} in 2005 was around 6.5 µg/m^{3}, below the standard of 8.8 µg/m^{3}. However, some cohort participants were exposed to PM_{2.5} levels higher than 8.8 µg/m^{3} during some years of follow-up and our hypothetical intervention only affected these subject-years by reducing their exposure to 8.8 µg/m^{3}, which represented a three-cycle average of 18.7% of subjects across all years. Since the observed PM_{2.5} levels decreased without any intervention in our study, fewer subjects were directly intervened on in later years under threshold intervention strategies, which explained the flattened differences in cumulative risks between intervention strategies in later years. However, all time-varying covariates after the intervention on PM_{2.5} would change accordingly due to the intervention on PM_{2.5}, thus influencing future outcomes as well. Such dynamic intervention strategy incorporated considerations of people who could be intervened on and are more realistic than the static intervention strategy commonly employed in health burden estimation with traditional exposure-response function, which sets change in exposure at a fixed value for all individuals throughout the period of interest. Besides, although we only provided differences in cumulative risk as compared to the natural course, it is easy to estimate differences between any two hypothetical intervention strategies.

Furthermore, the estimated cumulative risks over the follow-up time by g-computation (i.e., adjusted survival curves) and corresponding comparisons of values between different hypothetical interventions provided clearer causal interpretation and more information than a single HR or period-specific HRs, which is generally used in air pollution studies (Table S1). In the context of health impacts from chronic exposure to PM_{2.5}, HR can change over time since the toxicity of PM_{2.5} (e.g., chemical composition of PM_{2.5}) and susceptibility of population to PM_{2.5} could change over time, while standard Cox model assumed constant HR and period-specific HR from extensions of Cox had built-in bias that led to ambiguity in causal interpretation.^{56} On the other hand, the cumulative mortality risks estimated in this study avoided such ambiguity in interpretation while also demonstrating the change of intervention effect over time.^{19} Also, obtaining the casually interpretable absolute differences in cumulative risks between hypothetical intervention strategies over time could be particularly helpful for comparing different scenarios regarding public health benefits.^{57} Besides, if policies affecting air pollutants such as PM_{2.5} could further affect prognostic covariates influencing both future PM_{2.5} levels and health outcomes (commonly referred to as exposure-confounder feedback), traditional regression methods based on adjustment in a multivariable model would fail and lead to biased estimates for the effect while g-computation is designed to particularly solve this problem.^{24,26,58} An example of such exposure-confounder feedback is that people might move due to high level of PM_{2.5} in their current community and subsequently change the characteristics of their community of residence, while the characteristics of their current community also affect the level of PM_{2.5} and probability of death. Controlling for such community characteristics is necessary for confounding control but doing so with traditional methods will remove the indirect effect mediated through community characteristics and introduce collider-stratification bias^{59} if any unmeasured confounder of the community characteristics and death exists.^{58} However, making moving decision based on community level of PM_{2.5} is unlikely in countries with relatively low PM_{2.5} like Canada and exposure-confounder feedback is expected to be negligible in our study, but it is possible in countries with higher PM_{2.5}.

This study has a few limitations that need to be acknowledged. First, parametric g-computation can only account for measured confounders and lack of conditional exchangeability (i.e., residual confounding) might exist due to unmeasured confounders or measurement errors of existing confounders, regardless of our extensive list of confounders considered based on substantive knowledge on risk factors of PM_{2.5} exposure and death (Figure S1). For example, we assumed many individual behavior, demographic, and socioeconomic variables to be time-invariant (e.g., employment status and body mass index) due to data availability (these variables were only reported once at the time of enrollment) while they likely changed over the study period.

However, we also included time-varying individual income and community demographic and socioeconomic variables in our models, which mitigated the concern of residual confounding from these baseline variables. Besides, like other cohort studies of air pollution, we utilized postal-code level PM_{2.5} levels as surrogates for individual exposure to PM_{2.5}, which might introduce exposure misclassification.^{60} Recent studies showed that such bias may either not bias effect estimates ^{61} or bias these estimates towards the null,^{62} making our estimates more conservative.

Second, although we explored different model specifications and found similar results in sensitivity analyses, we cannot rule out the possibility of model misspecification, especially given the fact that parametric g-computation requires correct model specification of both the outcome and covariate models to achieve unbiased results. Notably, McGrath et al. demonstrated that the “g-null paradox”, a form of model misspecification that was traditionally believed to cause false rejection of null hypothesis under sharp null effect,^{63} is unavoidable in parametric g-computation even when the sharp null hypothesis does not hold, and recommended more flexible models in analysis.^{64} However, the magnitude of bias depends on the extent of exposure-confounder feedback and time-varying confounding. In the context of this study, we would expect relatively small exposure-confounder feedback thus less concern over g-null paradox.

Also, consistent results from sensitivity analysis using more flexible models supported the robustness of our results.

Third, being an active research field, the existing R package for parametric g-computation does not support features like incorporation of spline functions of time-varying covariates in the model, direct estimation of randomized interventional strategy,^{65} model fit checking with significance tests, or bias analysis. However, the current package provided enough flexibility for our study to employ flexible models that mitigated the possibility of violating the positivity assumption via model extrapolation. For example, we were able to incorporate flexible supra-linear PM_{2.5}-mortality association and temporal changes in the conditional probability of mortality in the estimation as supported by previous studies, incorporate restricted cubic spline function of baseline age in all models, and conduct sensitivity analysis with categorical time-varying confounders. Besides, although not relevant to our cohort since we had the all-cause mortality as the outcome and no loss to follow-up, right censoring and informative loss to follow-up could be handled by parametric g-computation and the existing R package by simulating data on participants as though they had not been censored.^{66}It is worth mentioning that other methods could also handle the methodological considerations that g-computation addresses—consideration of complex time structure and reporting of adjusted survival curves— and have been applied in air pollution epidemiological research, including Inverse Probability of Treatment Weighting (IPTW).^{6} Furthermore, some recent approaches such as the targeted maximum likelihood estimation can also be used to directly evaluate individualized dynamic intervention strategies of continuous exposures and provide doubly robust estimates that are less vulnerable to model misspecification with valid statistical inference when data-adaptive/machine-learning methods are incorporated.^{67,68}

Finally, PM_{2.5} is a mixture of varying chemical components and toxicity and is generated from different sources of emissions (e.g. traffic, industries, and wildfires). In this paper, we focused on PM_{2.5} without distinguishing the PM_{2.5} composition nor the sources of emissions. This potentially violated the consistency assumption (i.e. no-multiple-versions-of-treatment and all exposed individuals received the same version of treatment). If there is any unmeasured confounder of the “version of treatment” and outcome relationship, the effect estimates could be biased according to a recent simulation study, with magnitude and direction of such bias depending on the strength of confounding.^{69} In future studies, it would be important to consider the possible differential toxicity of PM_{2.5} components and define hypothetical interventions targeting different sources of PM_{2.5} emissions separately.

## Conclusion

This study demonstrated the benefits of using parametric g-computation as an analytical alternative for air pollution epidemiological research, especially for evaluating the potential effects of realistic dynamic intervention strategies in the time-to-event setting with time-varying exposure and confounders. With a large Canadian cohort, we calculated causally interpretable cumulative risk estimates over the follow-up time and corresponding benefits compared to maintaining status quo. We also found that any intervention further reducing the long-term exposure to PM_{2.5} would reduce the cumulative mortality risk from maintaining the status quo, even in a population already exposed to relatively low levels of ambient PM_{2.5}.

## Data Availability

The original data could not be shared to protect the privacy of study participants. All data produced in the present work are contained in the manuscript.

## Footnotes

↵* Co-senior authors

**Conflicts of interest:**The authors declare they have nothing to disclose.**Funding:**This study was funded by Health Canada (#810630). The funder supported the study design, data collection and analysis in this study.