# Two-stage interrupted time series analysis with machine learning: evaluating the health effects of the 2018 wildfire smoke event in San Francisco County as a case study

Arnab Dey, Yiqun Ma, Gabriel Carrasco-Escobar, Changwoo Han, François Rerolle, Tarik Benmarhnia

► To cite this version:

Arnab Dey, Yiqun Ma, Gabriel Carrasco-Escobar, Changwoo Han, François Rerolle, et al. Two-stage interrupted time series analysis with machine learning: evaluating the health effects of the 2018 wildfire smoke event in San Francisco County as a case study. American Journal of Epidemiology, 2025, 194 (10), pp.2936-2944. 10.1093/aje/kwaf147. hal-05200007

HAL Id: hal-05200007

https://hal.science/hal-05200007v1

Submitted on 9 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
OXFORD

<figure><img src="image_2.png" /></figure>

JOHNS HOPKINS

BLOOMBERG SCHOOL
of PUBLIC HEALTH

# Two-stage interrupted time series analysis with machine learning: evaluating the health effects of the 2018 wildfire smoke event in San Francisco County as a case study

Arnab K. Dey<fnref n="1" />,<fnref n="2" /> ID, Yiqun Ma<fnref n="1" /> ID, Gabriel Carrasco-Escobar<fnref n="3" />, Changwoo Han<fnref n="1" />,<fnref n="4" /> ID, François Rerolle<fnref n="1" />,<fnref n="5" />, Tarik Benmarhnia*<fnref n="1" />,<fnref n="6" /> ID

## Abstract

Randomized controlled trials (RCTs) are considered a key identification strategy for establishing causal relationships between exposures and outcomes. When evaluating the health impacts of extreme weather events, however, RCTs are generally infeasible due to ethical issues, costs, and the lack of a suitable control group. Quasi-experimental designs capitalizing on the timing of natural experiments, such as Interrupted Time Series (ITS), offer a valuable alternative to estimate causal effects when control groups are not available. This paper explores the application of a 2-stage ITS framework that compares traditional autoregressive integrated moving average (ARIMA) models and 2 machine learning algorithms: Neural Network Autoregressive (NNETAR) and Prophet-Extreme Gradient Boosting (XGBoost). As a case study, we assess the impacts of the 2018 wildfire smoke event on respiratory hospitalizations in San Francisco County, CA. We split the data into pre- and postevent periods to train and evaluate the models, perform cross-validation for hyperparameter tuning, and predict hospitalizations under the counterfactual scenario. Data and R code are provided for reproducibility. In the case study, the Prophet-XGBoost shows the best model performance and was used to generate the counterfactual trends. We estimate that the 2018 smoke event resulted in a total of 92 (95% empirical CI, 24-125) excess respiratory hospitalizations (12.5% of the observed hospitalization count during the event period). Our proposed approach offers a powerful tool for assessing the effects of extreme weather events and can be broadly applied to other epidemiological contexts, such as public health policy evaluation.

**Key words:** causal inference; interrupted time series analysis; wildfire smoke; respiratory hospitalization; California.

## Introduction

Randomized controlled trials (RCTs) are widely regarded as a key identification strategy to infer causal effects. However, RCTs are not always desirable or feasible due to ethical concerns, costs, or simply because the intervention of interest has already been implemented. This is particularly true when evaluating the health impacts of extreme weather events, which often occur on a large scale and cannot be randomly and ethnically assigned to specific populations. In this context, alternative quasi-experimental designs have been advanced, which capitalize on the timing of natural experiments to mimic randomization, enabling the estimation of causal effects using observational data to evaluate the impacts of extreme weather events or public health interventions.<fnref n="1" />⁻<fnref n="3" />

Most quasi-experimental methods compare an outcome of interest in a population exposed to an event, such as an extreme

weather event, with a control group that is not exposed to the event.<fnref n="2" />,<fnref n="3" /> However, control groups are not always available. This can be due to limited data availability or because all available geographical units were exposed at the same time. In this context, there are some approaches that compare the effects of an exposure within a single exposed population.<fnref n="1" /> Such designs are valuable in contexts where interventions are implemented arbitrarily (ie, without randomization) or where exposures occur in a whole population (ie, without a comparison group).<fnref n="1" />,<fnref n="4" /> One such design is the Interrupted Time Series (ITS) analysis, which is well suited to evaluate the effects of exposures occurred at a clearly defined time point.<fnref n="5" />

The ITS approach assumes that an event will change either the magnitude or direction of an outcome following its occurrence and leverages the randomness in the timing of the event to detect these changes. This is achieved by collecting time series data on

<footnote n="1">¹Scripps Institution of Oceanography, University of California, San Diego, CA, United States</footnote>

<footnote n="2">²Herbert Wertheim School of Public Health and Human Longevity Science, University of California, San Diego, CA, United States</footnote>

<footnote n="3">³Health Innovation Laboratory, Institute of Tropical Medicine 'Alexander von Humboldt', Universidad Peruana Cayetano Heredia, Lima, Peru</footnote>

<footnote n="4">⁴Department of Preventive Medicine, Chungnam National University College of Medicine, Daejeon, South Korea</footnote>

<footnote n="5">⁵Francis I. Proctor Foundation, University of California, San Francisco, CA, United States</footnote>

<footnote n="6">⁶Irset Institut de Recherche en Santé, Environnement et Travail, UMR-S 1085, Inserm, University of Rennes, EHESP, Rennes, France</footnote>

<footnote>*Corresponding author: Tarik Benmarhnia, Scripps Institution of Oceanography, University of California, San Diego, 9500 Gilman Dr, La Jolla, 92093 CA (tbenmarhnia@ucsd.edu)</footnote>

<footnote>**Author contributions:** Arnab K. Dey and Yiqun Ma are co-first authors, with names sequenced alphabetically to reflect equivalent contribution.</footnote>
an outcome of interest at multiple time points before and after an exposure event. This approach then extrapolates the time series data from the pre-event phase to generate a counterfactual trend that represents what would have been observed (under specific assumptions) in the absence of the exposure event. The effect of an event is then estimated by comparing the actual time series in the postevent period with the counterfactual time series.⁴ By studying a single population over a period of time, the ITS design can avoid issues related to differences between groups, like selection bias or unmeasured confounders.⁴

Generating counterfactual trends for ITS typically involves methods such as autoregressive integrated moving average (ARIMA),⁶ which assumes that the time trends in the pre- and postevent periods can be modeled as a linear combination of parameters. Such traditional ITS designs have been shown to generate biased estimates when certain underlying assumptions are not met.⁷ In addition, most of the ITS literature (both empirical examples and tutorials) focused on a single-stage model where all data units (pre- and post-data points) are included, and models are optimized based on a stepwise approach (eg, using the Akaike criterion). In parallel, some researchers proposed a 2-stage ITS framework,⁸,⁹ where the identification of causal effects is decomposed into 2 distinct steps: (1) using the pre-event data only, optimize a predictive model that will fit the observed data as best as possible based on calendar/seasonal and time-varying variables; and (2) once the predictive model in optimized in step 1, applying it to postevent data to compare with observed outcomes and infer the potential effects related to the exposure over time. A 2-stage approach offers a few advantages over traditional 1-stage ITS approaches. First, it allows for a careful examination of outcome trends in the pre-event period using multiple modeling approaches. This initial stage enables researchers to check whether the models and available data can adequately capture the underlying temporal patterns of the outcome before the event, which is essential for generating a valid counterfactual trend in the postevent period within an ITS design. In cases where pre-event trends cannot be adequately captured, researchers may opt for alternative causal inference methods that utilize different strategies for constructing counterfactuals. Second, this approach does not rely on the parametrization of changes in intercepts and slopes, as treatment effects can be estimated directly and for each time unit after the occurrence of the event, which offers more modeling flexibility.

Many recent machine learning (ML) approaches have demonstrated high performance in predicting trends through data-driven processes. These approaches allow nonlinear interactions between covariates and complex functional forms, and they have been shown to outperform predictions generated by linear approaches.¹⁰ While some ML approaches have been successfully used in prior studies in conjunction with ITS,¹¹,¹² they have not been used or compared in the context of 2-stage ITS analyses. In this paper, we aim to describe how ML algorithms can be used in the context of 2-stage ITS and compare the performance of multiple algorithms. Specifically, we will focus on Neural Network Autoregressive (NNETAR) and Prophet-Extreme Gradient Boosting (XGBoost) hybrid models and compare their performances to the traditional ARIMA model. We demonstrate the use of these methods, estimating the impacts of the 2018 wildfire smoke event on respiratory hospitalizations in San Francisco County, CA, as an illustrative case study.

The 2018 Camp Fire in Butte County, Northern California, which began on November 8, 2018, and burned for more than 2 weeks, was one of the most destructive wildfires in the state's history.

It directly caused 85 fatalities, displaced over 50 000 people, and destroyed more than 18 000 structures, resulting in an estimated $16.5 billion in damage.¹³ Smoke from the Camp Fire traveled south and west into the San Francisco County, driving the county's 24-hour average fine particulate matter (PM₂.₅) concentrations to exceed 35 µg/m³—above the US Environmental Protection Agency (EPA)'s primary standard—for 12 consecutive days, with a peak of 105 µg/m³. Wildfire PM₂.₅ has been shown to affect respiratory health.¹⁴,¹⁵ However, the impact of this 2018 wildfire smoke event on respiratory hospitalizations in San Francisco County has not been quantified yet.

This paper introduces a hybrid methodology that integrates multiple ML algorithms within a 2-stage ITS analysis, aiming at enhancing the robustness of the ITS approach for causal inference. We demonstrated the application of this hybrid approach through a case study that assessed the impacts of the 2018 wildfire smoke event on respiratory hospitalizations in San Francisco County, CA.

# Methods

## Study overview

### Methodological framework

This study introduced an observational 2-stage ITS design coupled with multiple algorithms to generate counterfactual trends. As a case study, we conducted a longitudinal analysis using daily data on respiratory hospitalization and environmental variables from 2009 to 2018 in San Francisco County, CA.

First, we used data before the 2018 wildfire smoke event to train and test 3 models: a traditional ARIMA model, a neural network designed for time series autoregression (NNETAR), and a hybrid model that combines Prophet algorithm with XGBoost (Prophet-XGBoost). Then, the model with the best performance was selected as the primary model to predict expected respiratory hospitalizations during and after the smoke event, representing the counterfactual scenario in which the smoke event did not occur. Finally, the difference between the observed and expected respiratory hospitalizations during the event period was calculated to determine the excess hospitalizations attributable to the smoke event. Figure 1 presents an overview of the methodology of the study.

### Wildfire smoke event period

In the main analysis, we defined the 2018 smoke event period as the 12 consecutive days (November 9-20, 2018) during which the daily mean PM₂.₅ concentration in San Francisco County exceeded the US EPA's primary standard of 35 µg/m³ (Figure 2). Considering that the Camp Fire began in its source location on November 8, 2018, we selected an alternative event period of November 8-20, 2018, for secondary analysis. Additionally, to account for the potentially delayed health effects of smoke exposure, we estimated the excess hospitalizations in the week following the smoke event as well (November 21-27, 2018).

# Data sources

### Respiratory hospitalization data

Our study analyzed respiratory hospitalization data from 2 complementary sources maintained by the California Department of Health Care Access and Information, collectively referred to hereafter as unscheduled hospitalizations: the Patient Discharge Data, which captures inpatient admissions, and the Emergency Department Data, which records emergency department visits.¹⁶ Using both databases enabled us to comprehensively capture the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Methodology framework for estimating excess respiratory hospitalizations due to the smoke event in San Francisco County.</figcaption></figure>

full spectrum of unscheduled respiratory healthcare encounters, including cases treated solely in emergency settings and those requiring inpatient admission. We calculated daily unscheduled respiratory hospitalization counts at the county level by aggregating ZIP code-level data based on patients' residential addresses. To identify respiratory hospitalizations, we utilized primary diagnosis codes for Diseases of the Respiratory System from the primary International Classification of Diseases (ICD)-9 and ICD-10 coding systems. The specific diagnosis codes are detailed in Appendix S1. This study was approved by the California Health and Human Services Agency's Committee for the Protection of Human Subjects (project number: 2021-116).

## Environmental data

Dail ZIP Code Tabulation Area (ZCTA)-level total and wildfire-specific PM<sub>2.5</sub> concentrations at population centroids were estimated using an ensemble-based statistical approach in a previous study and were averaged to the county level.<sup>17</sup> Nonsmoke PM<sub>2.5</sub> concentration was calculated as the difference between total PM<sub>2.5</sub> and wildfire-specific PM<sub>2.5</sub>. Daily maximum temperature, minimum temperature, dew point temperature, total precipita-

tion, and nonsmoke PM<sub>2.5</sub> concentration were used as predictors for our models. Maximum temperature, minimum temperature, and dew point temperature data at 4 × 4 km<sup>2</sup> were obtained from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) Climate Group.<sup>18</sup> Precipitation data at 4 × 4 km<sup>2</sup> were obtained from the Gridded Surface Meteorological (gridMET) reanalysis product.<sup>19</sup> To be consistent with the wildfire-specific PM<sub>2.5</sub> data, we first extracted values of these meteorological variables at the population-weighted centroids of each ZCTA and then calculated the average across all ZCTAs within San Francisco County.

## Data analysis

We implemented 3 modeling approaches: ARIMA, NNETAR, and Prophet-XGBoost. Our full dataset comprised daily counts of respiratory hospitalizations and associated predictor variables in San Francisco County from January 1, 2009, through December 31, 2018. We divided the dataset into 3 segments: a training period (January 1, 2009, to November 7, 2016), a testing period (November 8, 2016, to November 7, 2018), and a postevent period (during and after the smoke event; November 9, 2018, to December 31, 2018)
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Daily mean PM<sub>2.5</sub> concentrations in San Francisco County, before, during, and after the 2018 wildfire smoke event. This figure shows the time series of the daily mean total PM<sub>2.5</sub>, wildfire smoke PM<sub>2.5</sub>, and nonsmoke PM<sub>2.5</sub> concentrations (µg/m<sup>3</sup>) in San Francisco County, before, during, and after the 2018 wildfire smoke event. The horizontal dashed line represents the US Environmental Protection Agency (EPA)'s 24-hour average standard for PM<sub>2.5</sub> (35 µg/m<sup>3</sup>). The shaded rectangle marks the period of the smoke event (November 9-20, 2018), when the daily mean total PM<sub>2.5</sub> concentrations exceeded the EPA's standard. The lines of total PM<sub>2.5</sub> and nonsmoke PM<sub>2.5</sub> overlap on days without smoke PM<sub>2.5</sub>.</figcaption></figure>

(Figure 1). All statistical analyses were conducted with R software (version 4.4.0) using the tidymodels and modeltime packages. The datasets and R scripts used in this study are publicly available at https://github.com/benmarhnia-lab/Two_Stage_ML_ITS.

## Modeling approaches

**ARIMA.** We first used a traditional ARIMA model to capture the main components of the time series like the trend and seasonality. The ARIMA approach assumes that the future value of a variable is a linear function of several past observations, past forecast errors, and a random error term.<sup>20</sup> The ARIMA model can be mathematically defined as follows:

$$y_t = c + \sum_{i=1}^{p} \phi_i y_{t-i} + \sum_{j=1}^{q} \theta_j \epsilon_{t-j} + \epsilon_t$$

where c is a constant, $\phi_i$ are the coefficients of the autoregressive terms for $i = 1, 2, \dots, p$ (the order of the autoregressive part), $\theta_j$ are the coefficients of the moving average terms for $j = 1, 2, \dots, q$ (the order of the moving average part), and $\epsilon_t$ represents the error term related to random noise at time t.<sup>20</sup>

The ARIMA model tuning explored both nonseasonal and seasonal components. For nonseasonal parameters, the autoregressive (AR) order and moving average (MA) order were both tuned between 1 and 3, while differencing was tuned between 0 and 2. For seasonal components, both AR and MA orders were tuned between 0 and 3, with seasonal differencing limited to 0 and 1. The tuning process used a space-filling design with 100 combinations. The model also incorporated Fourier terms for yearly seasonality with 3 components to handle periodic patterns.

**Neural network autoregressive model.** Neural networks are flexible frameworks for modeling nonlinear problems.<sup>21</sup> Specifically, we proposed using a simple feed-forward neural

network. This model consists of a network of n processing units (neurons) organized in layers and connected by weighted links. The output of the network can be expressed by the following equation for a single hidden layer network:

$$y_t = \alpha_0 + \sum_{j=1}^{s} (\alpha_j g(\beta_{0j} + \sum_{i=1}^{r} \beta_{ij} x_{t-i})) + \epsilon_t$$

Where $x_{t-i}$ represents the input at time t - i, $\alpha_j$ (for $j = 0, 1, 2, \dots, s$) are the output layer weights, and $\beta_{ij}$ (for $i = 1, 2, \dots, r$; $j = 1, 2, \dots, s$) are the weights of the connections from the input layer to the j-th neuron in the hidden layer. The term r denotes the number of input nodes, and s is the number of hidden nodes. The function g(x), often chosen as the logistic sigmoid function, is the activation function for the hidden layer and is defined as:

$$g(x) = \frac{1}{1 + \exp(-x)}$$

This function introduces nonlinearity into the model, enabling it to capture complex patterns in the data beyond what is possible with linear models. The error term $\epsilon_t$ represents random noise.

The NNETAR tuning focused on both network architecture and time series components. The network architecture parameters included hidden units (8-20 nodes), number of networks (40-100), and penalty term (0.01-0.1) for regularization. For time series aspects, it tuned seasonal AR (1-4) and nonseasonal AR (1-6) components. The training process was configured with 50-200 epochs. The tuning process used a space-filling design with 75 combinations, allowing for efficient exploration of the hyperparameter space. The model also incorporated additional temporal features implemented via time series signatures.
**Prophet-XGBoost.** Finally, we used a hybrid approach combining the Prophet and XGBoost algorithms. According to Zhang,²¹ hybrid methodology suggests that by combining different models, the different underlying patterns of the time series of respiratory hospitalizations Y_t can be captured. This approach suggests that every model can be expressed in the following way:

$$Y_t = L_t + N_t$$

where L_t denotes the linear component and N_t denotes the nonlinear component. This allows great flexibility to create coupled models whose components can be optimally chosen. Thus, initially, from the adjustment of the linear model, we would obtain residuals that capture the nonlinear patterns. Let e_t denote the residual at time t from the linear model, then

$$e_t = Y_t - \hat{L}_t$$

A suitable model can be used to analyze these residuals to identify and represent the nonlinear patterns as follows:

$$\hat{N}_t = f(e_{t-j}) + \epsilon_t \quad \text{for } j = 1, \dots, n$$

where f is a nonlinear function defined by the nonlinear model and ε_t is the random error. Given all the above, the values predicted by the hybrid model would be as follows:

$$\hat{Y}_t = \hat{L}_t + \hat{N}_t$$

For this study, Prophet and XGBoost models were combined under the hybrid methodology framework to improve the prediction performance in the pre-event periods.

Prophet is an additive regression model developed by Facebook for forecasting time series data. It decomposes time series into 3 main components: trend, seasonality, and holidays. The trend component is modeled using a logistic growth curve to accommodate saturated growth patterns, while seasonal effects are captured through Fourier series to model cyclical behavior. The holiday component allows for the incorporation of specific events or occasions that can influence the time series data.²² The model can be expressed as:

$$y(t) = g(t) + s(t) + h(t) + \epsilon_t$$

Where g(t) represents the trend function, s(t) is the seasonal component that captures periodic changes, h(t) accounts for the effects of holidays, and ε_t is the error term.

XGBoost represents an advanced implementation of gradient boosting, where predictive models are built as an ensemble of weak learners, typically decision trees. It operates by consecutively adding predictors to an ensemble, with each new predictor correcting its predecessor's errors through gradient descent optimization. XGBoost minimizes a regularized objective function that balances prediction error and model complexity, ensuring both accuracy and simplicity in the final model.²³ The predictive model in XGBoost is formalized as follows:

$$\hat{y}_i = \sum_{k=1}^{K} f_k(x_i), \quad f_k \in \Phi$$

where ŷ_i is the prediction for the i-th instance, x_i represents the features of the i-th instance, K is the number of trees in the model, f_k represents the individual decision trees, and Φ is the space of all possible regression trees.

The hybrid Prophet-XGBoost model combines these approaches into a robust forecasting framework. Prophet decomposes the time series into trend, seasonality, and holiday components, providing a structured representation of temporal patterns. Then, XGBoost augments this framework by modeling the residuals from the Prophet model through an ensemble of decision trees, minimizing a regularized objective function that balances prediction error and model complexity. The hybrid implementation leverages Prophet's strength in capturing explicit temporal patterns while using XGBoost's ability to identify complex relationships in the residuals and additional features.

In our analysis, this hybrid approach was configured with tunable parameters including Prophet's changepoint_range (0.4-0.8), prior_scale_changepoints (0.01-0.5), and prior_scale seasonality (0.1-2.0), with fixed linear growth and yearly seasonality. The XGBoost component was configured for regression with gradient boosting (implied by the engine choice), tuning key parameters including number of variables per split (mtry: 4-25), minimum node size (min_n: 1-12), tree depth (8-20), learning rate (0.001-0.1), loss reduction (10⁻¹² to 10²), and early stopping iterations (10-30). The tuning process used a space-filling design with 100 combinations. The model incorporated additional features through time series signatures and 7-day lags of nonsmoke PM₂.₅ and considered Fourier terms for yearly seasonality with 3 components to handle periodic patterns.

## Model training and performance assessment

To maintain temporal dependencies, we implement a structured expanding-window time series cross-validation framework on the training data (January 1, 2009 to November 7, 2016). It started with a 2.9-year initial training window followed by a 12-month validation window. We generated 5 sequential folds with nonoverlapping validation window positioned at different times during the training period. The training data in each subsequent fold incorporates all previous training data plus additional 12-month observations (Figure S1).

We evaluate our 3 modeling approaches. Each methodology underwent extensive parameter optimization, exploring 100 parameter configurations for ARIMA, 75 for NNETAR, and 100 for Prophet-XGBoost. The optimal parameter set for each methodology is determined by minimizing the Root Mean Square Error (RMSE) across the cross-validation folds.

Then, for each modeling approach, we evaluated each optimized model using the testing data (November 8, 2016 to November 7, 2018). Model performance was comprehensively evaluated based on coefficient of determination (R²), mean absolute error (MAE), RMSE, mean absolute percentage error (MAPE), and symmetric mean absolute percentage error (SMAPE).

## Excess hospitalization calculation and uncertainty measurement

The best performing modeling approach in the testing period was then used to predict respiratory hospitalizations for the whole study period. To quantify prediction uncertainty, we employed a Moving Block Bootstrap (MBB) approach,²⁴,²⁵ performing 1000 model iterations. In brief, the observed series of daily hospitalizations during the pre-event period were decomposed into predicted values by the model and the residuals. The residuals of pre-event data were divided into n - L + 1 overlapping blocks,
<table><caption>**Table 1.** Performance metrics of the ARIMA, NNETAR, and Prophet-XGBoost model.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Training (January 1, 2009 to November 7, 2016)</th><th colspan="3">Testing (November 8, 2016 to November 7, 2018)</th></tr><tr><th>ARIMA</th><th>NNETAR</th><th>Prophet-XGBoost</th><th>ARIMA</th><th>NNETAR</th><th>Prophet-XGBoost</th></tr></thead><tbody><tr><td>R<sup>2</sup></td><td>0.71</td><td>0.84</td><td>0.76</td><td>0.65</td><td>0.58</td><td>0.69</td></tr><tr><td>MAE</td><td>7.08</td><td>5.36</td><td>6.49</td><td>8.65</td><td>9.44</td><td>8.64</td></tr><tr><td>RMSE</td><td>9.19</td><td>6.88</td><td>8.39</td><td>11.82</td><td>13.03</td><td>11.15</td></tr><tr><td>MAPE</td><td>0.14</td><td>0.11</td><td>0.13</td><td>0.16</td><td>0.16</td><td>0.17</td></tr><tr><td>SMAPE</td><td>0.13</td><td>0.10</td><td>0.12</td><td>0.16</td><td>0.17</td><td>0.16</td></tr></tbody></table>

Abbreviations: MAE, mean absolute error; MAPE, mean absolute percentage error; R², coefficient of determination; RMSE, root mean square error; SMAPE, symmetric mean absolute percentage error.

where *n* is the length of the pre-event data (3958 days) and *L* is the size of each block. In each iteration, we sampled *n*/L blocks with replacement and concatenated them to form a bootstrapped residual series. These bootstrapped residuals were added back to the model predictions to yield a new pseudo-time series to which the pretrained model was reapplied to generate predictions after the event occurrence. The 95% empirical CIs (eCIs) were derived from the empirical distribution of these 1000 predictions, taking the 2.5th and 97.5th percentiles as the lower and upper bounds, respectively. Note that this approach does not enforce symmetry around the original point prediction, so the central estimate may not lie exactly at the midpoint of the intervals, but its simplicity and minimal assumptions have made it a common choice in time series applications.²⁶⁻²⁸ We selected *L* = 14 based on a nonparametric plug-in rule proposed by Lahiri et al.²⁹ This choice was further confirmed when we ran the MBB procedure for a range of block size (1-100), inspected how the width of the 95% CIs change, and chose the smallest *L* that yields stable interval widths (Figure S2). By sampling blocks of contiguous observations rather than individual observations, the MBB approach allows the bootstrap samples to retain the correlation structure within blocks, making it more suitable for time series data.²⁴,²⁵

The predicted hospitalizations during and after the smoke event represents the expected hospitalization counts under the counterfactual scenario in which the smoke event did not occur. The excess hospitalizations resulting from the smoke event were then calculated by the difference between the observed and expected hospitalizations. The 95% eCIs of the estimated excess hospitalizations were calculated based on the 95% eCIs of the model predictions. We calculated both daily excess hospitalizations and the total excess hospitalizations during the whole period of the smoke event.

## Results

During the 2018 wildfire smoke event (November 9-20, 2018), a total of 738 respiratory hospitalizations were recorded in San Francisco County, CA. The average daily total PM₂.₅ concentration was 75.6 µg/m³, of which approximately 88% was from wildfire smoke.

Table 1 presents the performance metrics for the 3 models evaluated on both training and testing datasets. In the training dataset, the NNETAR model showed an overall top performance, with an R² of 0.84 and consistently the lowest error measures, while the ARIMA model showed the poorest performance. In the testing dataset, Prophet-XGBoost showed the best performance, with the highest R² (0.69) and the lowest MAE, RMSE, and SMAPE. Overall, although both ARIMA and NNETAR models showed reasonable performance, they were outperformed by the Prophet-XGBoost model in the testing dataset. Therefore,

Prophet-XGBoost was selected to predict the counterfactual scenario in the postevent period.

The hybrid Prophet-XGBoost model predicted a total of 646 (95% eCI, 613-714) expected respiratory hospitalizations in San Francisco County from November 9 to 20, 2018, under the counterfactual scenario in which the smoke event did not occur. Therefore, by comparing the observed and expected hospitalizations, we estimated a total of 92 (95% eCI, 24-125) excess respiratory hospitalizations resulting from the 2018 wildfire smoke event in San Francisco County (approximately 8 excess hospitalizations per day), accounting for approximately 12.5% of the total observed hospitalization count during the event period (Figure 3; Table 2). The estimated daily excess respiratory hospitalizations during the 2018 wildfire smoke event were also listed in Table 2. The highest attributable fraction of respiratory hospitalizations occurred on November 16 (approximately 23.5%), followed by the first day of the smoke event (November 9; approximately 22.7%) and November 12 (approximately 20.1%).

In a secondary analysis, we included November 8, 2018 (the day the Camp Fire began in its source location) as the beginning of the wildfire smoke event period. The estimated excess respiratory hospitalizations during this extended period were 88 (95% eCI, 15-125), which is similar to the estimate obtained using the primary definition of the smoke event period. In the week following the smoke event (November 21-27, 2018), the estimated difference between the total observed and expected hospitalizations was null (−1 [95% eCI, −42 to 19]).

## Discussion

In this paper, we proposed a hybrid approach combining ML and a 2-stage ITS analysis. We applied this methodology to assess the health impacts of the 2018 wildfire smoke event in San Francisco County, CA, serving as a case study. Using the hybrid Prophet-XGBoost model to generate the hospitalization trend under counterfactual scenario, we estimated that this 12-day smoke event resulted in a total of 92 (95% eCI, 24-125) excess respiratory hospitalizations in this county, which accounted for about 12.5% of the observed hospitalizations. This innovative method augments traditional ITS analysis by integrating ML algorithms to generate robust counterfactual trends, allowing for more accurate causal inference of the impacts of extreme weather events. Additionally, using ML models to generate counterfactuals allows for more flexible assumptions as compared to the stringent assumptions of linearity and seasonality under the traditional ITS framework.

Our findings on the effects of wildfire smoke events on respiratory health are overall consistent with the broader literature.¹⁴,¹⁵ Specifically, the magnitude of the effect estimate in our study aligns with previous research. For example, a previous study investigating the health effects of the same 2018 wildfire
<table><caption>**Table 2.** Estimated excess respiratory hospitalizations during the 2018 wildfire smoke event in San Francisco County.</caption><thead><tr><th>Date/period</th><th>Observed hospitalizations (count)</th><th>Expected hospitalizations (count, 95% eCI)</th><th>Excess hospitalizations (count, 95% eCI)</th><th>Attributable fraction<sup>a</sup> (%, 95% eCI)</th></tr></thead><tbody><tr><td colspan="5"><strong>Daily estimates during the event</strong></td></tr><tr><td>November 9, 2018</td><td>63</td><td>49 (46-53)</td><td>14 (10-17)</td><td>22.7 (15.2-27.7)</td></tr><tr><td>November 10, 2018</td><td>58</td><td>54 (51-59)</td><td>4 (-1 to 7)</td><td>7.6 (-1.0 to 12.7)</td></tr><tr><td>November 11, 2018</td><td>68</td><td>56 (54-62)</td><td>12 (6-14)</td><td>17.9 (9.0-20.8)</td></tr><tr><td>November 12, 2018</td><td>73</td><td>58 (54-63)</td><td>15 (10-19)</td><td>20.1 (13.9-25.4)</td></tr><tr><td>November 13, 2018</td><td>54</td><td>51 (49-57)</td><td>3 (-3 to 5)</td><td>4.8 (-6.2 to 9.4)</td></tr><tr><td>November 14, 2018</td><td>63</td><td>50 (48-57)</td><td>13 (6-15)</td><td>20.1 (10.2-23.9)</td></tr><tr><td>November 15, 2018</td><td>44</td><td>49 (47-56)</td><td>-5 (-12 to -3)</td><td>-12.3 (-27.0 to -6.8)</td></tr><tr><td>November 16, 2018</td><td>67</td><td>51 (48-57)</td><td>16 (10-19)</td><td>23.5 (15.5-28.3)</td></tr><tr><td>November 17, 2018</td><td>59</td><td>56 (53-61)</td><td>3 (-2 to 6)</td><td>5.4 (-4.1 to 10.6)</td></tr><tr><td>November 18, 2018</td><td>62</td><td>59 (56-64)</td><td>3 (-2 to 6)</td><td>4.4 (-3.9 to 9.6)</td></tr><tr><td>November 19, 2018</td><td>61</td><td>59 (57-65)</td><td>2 (-4 to 4)</td><td>3.3 (-6.9 to 7.3)</td></tr><tr><td>November 20, 2018</td><td>66</td><td>53 (51-60)</td><td>13 (6-15)</td><td>19.7 (9.2-22.8)</td></tr><tr><td colspan="5"><strong>Total estimates during the event</strong></td></tr><tr><td>November 9-20, 2018</td><td>738</td><td>646 (613-714)</td><td>92 (24-125)</td><td>12.5 (3.2-17.0)</td></tr></tbody></table>

The expected hospitalizations were predicted using the Prophet-XGBoost model.

event in California using a synthetic control approach reported approximately 98 increased respiratory hospitalizations in San Francisco County during the event.³⁰ Another study utilized a 2-stage ITS design to estimate the healthcare utilization effects of the December 2017 Lilac Fire in San Diego County among pediatric patients, and found that the Lilac Fire was associated with approximately 16 excess respiratory visits per day (underlying population size: approximately 690 000 children aged 0-19),⁹ which is higher than our estimate of about 8 excess hospitalizations per day (underlying population size: about 800 000 people of all age groups). However, this direct comparison of estimated excess hospital visits across studies focusing on different smoke events should be interpreted with caution due to multiple factors, such as

the difference in wildfire characteristics, potential heterogeneity in wildfire effects across age groups, whether the region was directly burned or only exposed to smoke, and possible differences in the chemical composition of the smoke. The overall consistency between results from existing methods and the newly proposed hybrid approach provides robust support for the acute respiratory impacts of wildfire smoke. Our hybrid methodology, which combines ML and a 2-stage ITS analysis, can be utilized alongside other causal inference designs to generate more robust estimates of the impacts of extreme weather events or the effectiveness of interventions.

Our study presents a novel application of ML techniques within an ITS framework to evaluate the health impacts of extreme

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Daily observed and expected respiratory hospitalizations in San Francisco County, before, during, and after the 2018 wildfire smoke event. This figure illustrates the time series of the daily observed (orange line) and expected (the counterfactual; purple line) respiratory hospitalizations in San Francisco County, before, during, and after the 2018 wildfire smoke event. The expected hospitalizations were predicted using the Prophet-XGBoost model. The shaded rectangle marks the period of the smoke event (November 9-20, 2018). The shaded area around the expected hospitalizations represents the 95% empirical CI for the expected hospitalizations.</figcaption></figure>

<footnote>ªAbbreviations: Attributable fraction, fraction of observed hospitalizations attributable to the smoke event; eCI, empirical CI.</footnote>
weather events. We showcase the 2-stage hybrid approach and
demonstrate its flexibility and robustness of causal inference.
This method can be widely applied to other contexts, such as
evaluating the effectiveness of policies and interventions. The 2-
stage approach can also be easily adapted to staggered interven-
tions³¹ (when no control groups are available), which is a very
active area of research in the context of Difference-in-Differences
(DID) methods. It can be extended to assess heterogeneity in effect
estimates across time or treated units (if multiple treated units
receive the treatment in a sequential way).³² We also believe this
approach can be readily applied to evaluate interventions across
diverse health issues and geographical contexts, highlighting its
potential as an important tool in public health research.

While our approach marks a significant methodological advancement, it diverges from other analytical methods, such as DID or synthetic control methods. One caveat in applying this approach is the potential overfitting of the ML models. Although overfitting was not an issue in our case study, it may arise in other contexts. However, since the primary purpose of training the ML model is to generate a counterfactual scenario, the potential close alignment of pre-event trends with actual data is of less concern. In addition, unlike standard single-stage ITS models, although the nonparametric and complex modeling strategies in our proposed ML-integrated 2-stage ITS models make them more flexible to unexpected patterns in the data, it makes it more difficult for researchers to explicitly state and test their assumptions about the impact model and the temporal patterns of effects in the postevent period. Further work to address this limitation is warranted in future studies.

Conclusion

There is an urgent need for more robust and flexible methods to evaluate the impacts of extreme weather events, particularly as climate change exacerbates a range of existing and emerging public health threats. Our approach, which integrates ML techniques with traditional quasi-experimental methods, provides researchers with an enhanced tool for generating counterfactual scenarios. This methodological framework enhances our capacity to assess the impact of events or interventions in real-world situations where obtaining a direct counterfactual may not be possible. The approach presented in this paper offers an innovative method for assessing the effects of extreme weather events and public health interventions, which can provide robust scientific evidence to better inform public health strategies amid evolving environmental and epidemiological landscapes.

Supplementary material

Supplementary material is available at the *American Journal of Epidemiology* online.

Funding

This work was supported by Climate Action 2023 UC Innovation and Entrepreneurship Awards R02CE6859 and by the National Institute on Aging (RF1AG080948).

Conflict of interest

The authors declare they have no actual or potential competing financial interests.

Data availability

The datasets and R scripts used in this study are publicly available
at https://github.com/benmarhnia-lab/Two_Stage_ML_ITS.

References

1. Bonell CP, Hargreaves J, Cousens S, et al. Alternatives to randomisation in the evaluation of public-health interventions: design challenges and solutions. *J Epidemiol Commun Health.* 2011;65(7):582-587. https://doi.org/10.1136/jech.2008.082602

2. Nianogo RA, Benmarhnia T, O'Neill S. A comparison of quasi-experimental methods with data before and after an intervention: an introduction for epidemiologists and a simulation study. *Int J Epidemiol.* 2023;52(5):1522-1533. https://doi.org/10.1093/ije/dyad032

3. Sheridan P, McElroy S, Casey J, et al. Using the generalized synthetic control method to estimate the impact of extreme weather events on population health. *Epidemiology.* 2022;33(6):788-796. https://doi.org/10.1097/EDE.0000000000001539

4. Bernal JL, Cummins S, Gasparini A. Interrupted time series regression for the evaluation of public health interventions: a tutorial. *Int J Epidemiol.* 2017;46(1):348-355. https://doi.org/10.1093/ije/dyw098

5. Craig P, Katikireddi SV, Leyland A, et al. Natural experiments: an overview of methods, approaches, and contributions to public health intervention research. *Annu Rev Public Health.* 2017;38(1):39-56. https://doi.org/10.1146/annurev-publhealth-031816-044327

6. Hartmann DP, Gottman JM, Jones RR, et al. Interrupted time-series analysis and its application to behavioral data. *J of App Behav Analysis.* 1980;13(4):543-559. https://doi.org/10.1901/jaba.1980.13-543

7. Baicker K, Svoronos T. Testing the Validity of the Single Interrupted Time Series Design. National Bureau of Economic Research. 2019.

8. Feder KA, Mojtabai R, Stuart EA, et al. Florida's opioid crackdown and mortality from drug overdose, motor vehicle crashes, and suicide: a Bayesian interrupted time-series analysis. *Am J Epidemiol.* 2020;189(9):885-893. https://doi.org/10.1093/aje/kwaa015

9. Leibel S, Nguyen M, Brick W, et al. Increase in pediatric respiratory visits associated with Santa Ana wind-driven wildfire smoke and PM₂.₅ levels in San Diego County. *Ann Am Thorac Soc.* 2020;17(3):313-320. https://doi.org/10.1513/AnnalsATS.201902-150OC

10. Makridakis S, Spiliotis E, Assimakopoulos V, et al. Statistical, machine learning and deep learning forecasting methods: comparisons and ways forward. *J Oper Res Soc.* 2023;74(3):840-859. https://doi.org/10.1080/01605682.2022.1128629

11. Linden A, Yarnold PR. Using machine learning to identify structural breaks in single-group interrupted time series designs. *J Eval Clin Pract.* 2016;22(6):855-859. https://doi.org/10.1111/jep.12544

12. Linden A, Yarnold PR. Using machine learning to evaluate treatment effects in multiple-group interrupted time series analysis. *J Eval Clin Pract.* 2018;24(4):740-744. https://doi.org/10.1111/jep.12966

13. Zimmerman, J. Service Assessment: November 2018 Camp Fire. National Weather Service: Salt Lake City, UT, USA. 2020. https://www.weather.gov/media/publications/assessments/sa1162SignedReport.pdf
14. Reid CE, Brauer M, Johnston FH, et al. Critical review of health impacts of wildfire smoke exposure. *Environ Health Perspect.* 2016;124(9):1334-1343. https://doi.org/10.1289/ehp.1409277

15. Liu JC, Pereira G, Uhl SA, et al. A systematic review of the physical health impacts from non-occupational exposure to wildfire smoke. *Environ Res.* 2015;136:120-132. https://doi.org/10.1016/j.envres.2014.10.015

16. California Department of Health Care Access and Information (HCAI). Data Documentation [Internet]. HCAI; 2024 [cited 2024 March 14]. Available from: https://hcai.ca.gov/data-and-reports/request-data/data-documentation/

17. Aguilera R, Luo N, Basu R, et al. A novel ensemble-based statistical approach to estimate daily wildfire-specific PM<sub>2.5</sub> in California (2006–2020). *Environ Int.* 2023;171:107719. https://doi.org/10.1016/j.envint.2022.107719

18. PRISM Group, Oregon State University, https://prism.oregonstate.edu, data created 4 Feb 2014, accessed 21 Dec 2023.

19. Abatzoglou JT. Development of gridded surface meteorological data for ecological applications and modelling. *Int J Climatol.* 2013;33(1):121-131. https://doi.org/10.1002/joc.3413

20. Box GE, Jenkins GM, Reinsel GC, et al. *Time Series Analysis: Forecasting and Control.* John Wiley & Sons; 2015.

21. Zhang GP. Time series forecasting using a hybrid ARIMA and neural network model. *Neurocomputing.* 2003;50:159-175. https://doi.org/10.1016/S0925-2312(01)00702-0

22. Taylor SJ, Letham B. Forecasting at scale. *Am Stat.* 2018;72(1):37-45. https://doi.org/10.1080/00031305.2017.1380080

23. Chen T, Guestrin C. Xgboost: a scalable tree boosting system. In: *Proceedings of the 22nd ACM sigkdd International Conference on Knowledge Discovery and Data Mining* 785-794. Association for Computing Machinery. 2016.

24. Künsch HR. The jackknife and the bootstrap for general stationary observations. *Ann Stat.* 1989;17(3):1217-1241.

25. Liu RY, Singh K. Moving blocks jackknife and bootstrap capture weak dependence. In: Lepage R, Billard L, eds. *Exploring the Limits of Bootstrap*. Vol. 225. John Wiley; 1992:248.

26. Kioumourtzoglou M-A, Coull BA, Dominici F, et al. The impact of source contribution uncertainty on the effects of source-specific PM<sub>2.5</sub> on hospital admissions: a case study in Boston, MA. *J Expo Sci Environ Epidemiol.* 2014;24(4):365-371. https://doi.org/10.1038/jes.2014.7

27. Zhang F, Wagner AK, Soumerai SB, et al. Methods for estimating confidence intervals in interrupted time series analyses of health interventions. *J Clin Epidemiol.* 2009;62(2):143-148. https://doi.org/10.1016/j.jclinepi.2008.08.007

28. Mudelsee, M. & Mudelsee, M. Bootstrap confidence intervals. In: *Climate Time Series Analysis: Classical Statistical and Bootstrap Methods*. 2014: 61-104. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-04450-7_3

29. Lahiri S, Furukawa K, Lee Y-D. A nonparametric plug-in rule for selecting optimal block lengths for block bootstrap methods. *Stat Methodol.* 2007;4(3):292-321. https://doi.org/10.1016/j.stamet.2006.08.002

30. Letellier N, Hale M, Salim KU, et al. Applying a two-stage generalized synthetic control approach to quantify the heterogeneous health effects of extreme weather events: a 2018 large wildfire in California event as a case study. *Environ Epidemiol.* 2025;9(1):e362. https://doi.org/10.1097/EE9.0000000000000362

31. Batomen B, Benmarhnia T. Staggered interventions with no control groups. *Int J Epidemiol.* 2024;53(6):dyae137. https://doi.org/10.1093/ije/dyae137

32. Wing C, Yozwiak M, Hollingsworth A, et al. Designing difference-in-difference studies with staggered treatment adoption: key concepts and practical guidelines. *Annu Rev Public Health.* 2024;45(1):485-505. https://doi.org/10.1146/annurev-publhealth-061022-050825