Arima

Including a control series in ITS analysis improves causal inference, as ITS cannot exclude the possibility that any observed change was due to the intervention of interest, or another co-intervention or event. A control series is one that is not impacted by the intervention; selection of an appropriate control is described elsewhere [3 (link)]. As with ITS in segmented regression, including a control series involves running an ARIMA model for the series of interest, and separately for the control series [17 (link)]. If a change is observed in the intervention series but not the control series, this provides evidence that the impact was specific to the intervention.

Chinese HFRS incidence data from 1975 to 2008 was obtained from the Chinese Center for Disease Control and Prevention. All HFRS cases were initially diagnosed by clinical symptoms. Patient blood samples were also collected and sent to local Centers for Disease Control and Prevention (CDC) laboratories for serological confirmation. Finally, data were collected by case number according to the sampling results. There might be admission rate bias in the disease report, but this has been reduced as much as possible. In China, HFRS is a nationally notifiable disease and hospital physicians must report every case of HFRS to the local health authority within 12 hours. Local health authorities later report monthly HFRS case totals to higher the national level CDC for surveillance purposes. Due to mandatory reporting, it is believed that the degree of compliance in disease notification over the study period was consistent.
We used the Box-Jenkins approach to ARIMA (p, d, q) modeling of time series [20 ]. This model-building process is designed to take advantage of associations in the sequentially lagged relationships that usually exist in periodically collected data [21 (link)]. The following were the parameters selected when fitting the ARIMA model: p, the order of autoregression; d, the degree of difference; q, the order of moving average.
The annual data used in this study did not show seasonal pattern, so the series was differenced at the non-seasonal level to induce stationarity. Autocorrelation function (ACF) graph and Partial autocorrelation function (PACF) graph were used to identify the order of moving average (MA) and autoregressive (AR) terms included in the ARIMA model. Estimates of the model's parameters were obtained by the conditional least squares method. Diagnostic checking including residual analysis and the Akaike Information Criterion (AIC) was used to compare the goodness-of-fit among ARIMA models. The Ljung-Box test was used to measure the ACF of the residuals. In addition, we used the mean absolute percentage error (MAPE) and fitting effect diagram to assess forecast accuracy.
, where x_t and denote observed and fitted values at time point t. The MAPE value was calculated based on observed values and fitted values from 1978 to 2008. A lower MAPE value indicates a better fit of the data. Finally, the fitted ARIMA model was used for short-term forecasting of HFRS incidence for years 2009 to 2011. All analyses were performed using SAS9.1 with a significant level of p < 0.05.

We used the mortality model previously published by Foreman and colleagues,²² (link) and extended it to 2100 with slight modifications. Briefly, the cause-specific model included three components: the underlying mortality, modelled as a function of the Socio-demographic Index (SDI), time, and additional cause-specific covariates where appropriate; a risk factor scalar that captured the combined risk factor effects for specific causes, based on the GBD 2017 cause-risk hierarchy and accounting for risk factor mediation;²⁴ (link) and an autoregressive integrated moving average (ARIMA) model²⁵ that accounted for unexplained residual mortality.
To accommodate long-range forecasts, we removed the spline on SDI and used a random walk with attenuated drift for the ARIMA model. Foreman and colleagues found that our mortality model had better out-of-sample predictive validity than the most widely used demographic forecasting model.²² (link) The method used to develop reference scenario values for each of the independent drivers in the mortality model was not modified from Foreman and colleagues.²² (link)

All data were stored on a secure drive accessible only to the study team. Statistical analyses were performed in Rstudio version 4.0.1 (Rstudio Team 2022 ). Using the NEISS-Work dataset, national ED-treated occupational injury count estimates were produced using the R packages “survey” and “srvyr” (Ellis et al. 2021 ; Lumley 2021 ) using the aforementioned NEISS-Work survey weights. ED-treated occupational injury count estimates were generated for all injuries and by injury event type, a categorical variable denoting the way an injury was incurred and is based on the aforementioned BLS OIICS v 2.01 classification system (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021b ); all analyses were conducted both for total injury rate estimates and stratified by injury event type. ED-treated occupational injury rates were calculated per 10,000 FTE using Current Population Survey (CPS) estimates which were generated using NIOSH’s Employed Labor Force (ELF) query system; as NEISS-Work includes all work-related ED-treated injuries, FTE estimates were generated for all jobs (as opposed to “primary” or “secondary” jobs only) (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021c ). Standard errors (SE) for FTE estimates were generated using generalized variance functions provided by BLS; standard errors were used to calculate monthly FTE variances by multiplying the square of the SE by corresponding ELF-generated monthly FTE estimates (i.e., the corresponding monthly sample size) (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021c ). Variances of both numerator (injury count estimates) and denominator (FTE) data were used to calculate 95% confidence intervals (CI) for ED-treated occupational injury rate estimates based on Taylor series expansion (National Institute for Occupational Safety and Health (NIOSH) Division of Safety Research 2021d ) and were reported as injury rate estimates ± margin of error.
Seasonality of injury rate estimates was assessed by calculating seasonality indices per month. Seasonality indices were calculated by dividing the mean rate for each month by the mean monthly occupational injury rate for the entire dataset; seasonality indices of greater and less than one indicate higher than and lower than expected injury rates for a given month, respectively (Zhang et al. 2014 ).
To assess linear trends in injury rates over time, we fit a linear regression model to monthly injury rate estimates and adjusted for autocorrelation and serially correlated error terms using autoregressive integrated moving average (ARIMA) modeling. This analysis was conducted using both monthly total injury rate estimates and monthly estimates stratified by injury event type. In data violating the linear regression assumption of no autocorrelation, ARIMA models are used to control for serial correlation (e.g., seasonality) by including lagged dependent variable values and errors, including in studies of injury data (Box et al. 2016 ; Zhu et al. 2015 (link)). An ARIMA model takes the form ARIMA(p,d,q)(P,D,Q)_m, where p is the order of autocorrelation, d is the number of differences applied to the data, q is the order of moving average terms, P, D, and Q are the seasonal versions of these terms, and m is the order of seasonality (e.g., 12 for annually seasonality in monthly data) (Hyndman and Athanasopoulos 2018a ). ARIMA models were fit to monthly injury rates by examining autocorrelation and partial autocorrelation plots. A lagged regression estimate was included if it showed statistical significance (p < 0.05) and was necessary to control for serial correlation. Finally, significance of each model’s Ljung-Box Q statistic was observed to ensure proper model fit, with a non-significant value considered a properly fit model (Ljung and Box 1978 ). The conditional sum of squares method was used to estimate all models. To assess temporal trends, a trend regressor with slope of one was included in each ARIMA model as a covariate and reported with 95% CIs (Hyndman and Athanasopoulos 2018b ). A total percent decrease in injury rates throughout the study period was estimated by multiplying this term by 96 (i.e., the total number of months in the study period) and calculating the percent difference from the model’s intercept; an analogous calculation using each trend parameter’s 95% CI was performed to determine each percent decrease’s 95% CI.

We calculated incidence rates as number of cases (COVID-19 or admissions) divided by each country's population. We also calculated 7-day smoothed rolling average rates to reduce the effects of lower reporting on weekends.
We performed the first analysis on the 7-day smoothed data using NBSR. We also considered ARIMA models for autocorrelation.
To further strengthen the results and given that England did not implement the COVID-19 certification when it was effective in the other three countries, we used its data as a counterfactual for Difference-in-Differences (DiD) models. To help visualize this method, a plot of the difference and cumulative difference of the incidence rates for cases and hospitalizations of all countries is provided in Supplementary Figure 5. The numbers shown are essentially what constitute the basis of the DiD model. The differences have been calculated extracting England's incidence rates from the other countries' incidence rates. A decreasing trend in the difference's plots (on the left) is associated with a protective effect of the intervention date on the outcome.
We performed all the analyses in R v4.3 and used the packages epiR, tidyverse, forecast, ggplot2, MASS and lmtest. Code is available in https://github.com/KimLopezGuell/Covid-passport.

We performed Negative binomial segmented regression (NB) and Autoregressive integrated moving average (ARIMA) models as a preliminary and sensitivity analysis of the main model, Difference-in-Differences. Detailed methods results and output for NB and ARIMA can be found in the Supplementary material.
Difference-in-Differences (DiD) methods compare the mean of the variable of interest for an exposed and control group, before and after a certain interruption point, providing insight on the changes of the variable for the exposed countries relative to the change in the negative outcome group (15 (link)). We cannot draw causal conclusions by simply observing before-and-after changes in outcomes, because other factors might influence the outcome over time. DiD methods overcome this by introducing a comparison between two similar groups exposed to different conditions. First, DiD takes the difference of the variable of interest of both groups before and after the intervention. Then it subtracts the difference of the control group to the difference of the exposed one to control for time varying factors, therefore giving a result which constitutes a difference of the differences. This approximates the clean impact of the intervention. In essence, the DiD estimating equation is the following,
where Y_gt is the outcome for an individual in group g and treated unit t, P_t is a binary time variable indicating whether the observation belongs to the period before or after the intervention and T_g is a binary variable indicating whether the observation belongs to the exposed or the controlled group. In this setting, the treatment effect is estimated with the coefficient β₃ from the regression.
For this method to be rightly used, all the typical OLS assumptions must be met. The parallel trends assumption, which requires both groups to present similar trends before the intervention time point (16 (link)), must also be satisfied. We tested all these assumptions, and the latter can be visually inspected in Figures 1–3.
DiD models produce estimates which consider a counterfactual group, therefore adjusting for unmeasured confounding. This cannot be done by neither of the two previous models. Its biggest limitation is that, in the end, the measured effect can only be attributed to the timepoint chosen. If that is due to the intervention placed then, or to other underlying reasons around the same time, cannot be known by design.