Distributed lag interrupted time series model for unclear intervention timing: effect of a statement of emergency during COVID-19 pandemic

Abstract

Background

Interrupted time series (ITS) analysis has become a popular design to evaluate the effects of health interventions. However, the most common formulation for ITS, the linear segmented regression, is not always adequate, especially when the timing of the intervention is unclear. In this study, we propose a new model to overcome this limitation.

Methods

We propose a new ITS model, ARIMAITS-DL, that combines (1) the Autoregressive Integrated Moving Average (ARIMA) model and (2) distributed lag functional terms. The ARIMA technique allows us to model autocorrelation, which is frequently observed in time series data, and the decaying cumulative effect of the intervention. By contrast, the distributed lag functional terms represent the idea that the intervention effect does not start at a fixed time point but is distributed over a certain interval (thus, the intervention timing seems unclear). We discuss how to select the distribution of the effect, the model construction process, diagnosing the model fitting, and interpreting the results. Further, our model is implemented as an example of a statement of emergency (SoE) during the coronavirus disease 2019 pandemic in Japan.

Results

We illustrate the ARIMAITS-DL model with some practical distributed lag terms to examine the effect of the SoE on human mobility in Japan. We confirm that the SoE was successful in reducing the movement of people (15.0–16.0% reduction in Tokyo), at least between February 20 and May 19, 2020. We also provide the R code for other researchers to easily replicate our method.

Conclusions

Our model, ARIMAITS-DL, is a useful tool as it can account for the unclear intervention timing and distributed lag effect with autocorrelation and allows for flexible modeling of different types of impacts such as uniformly or normally distributed impact over time.

Introduction

Interrupted time series (ITS) has become a popular study design to evaluate the effects of health interventions in the field of public health and epidemiology [1, 2]. The effect of the intervention is estimated by comparing it with the “counterfactual" that is estimated by the expected trend assuming the absence of the intervention. ITS has been regarded as one of the best designs to estimate the causality of the intervention, especially when no control population is available, or randomized controlled trials (RCTs) are not feasible [3, 4]. Fretheim et al. (2015) demonstrated that ITS could generate effect estimates similar to those of RCTs [5]. Previous studies extensively examined the strengths and limitations of ITS and provided a practical guideline for its application [3, 6,7,8].

More formally, let Y_t be an outcome of interest at time t, and T₀ be the time of the intervention of interest. For statistical inference, we assume that the sampled time series data $\{Y_{t}\}_{t=1}^{T}$ is available. Standard ITS can be expressed as a simple segmented linear regression formulation, which is given by:

$$Y_{t} = \beta_{0} + \beta_{1}t + \beta_{2} D_{t} + \beta_{3} {tD}_{t} + \varepsilon_{t},$$

(1)

where ε_t is Gaussian white noise with constant variance $\sigma ^{2}_{\varepsilon }$ (i.e., independent of t) and D_t is a dummy variable indicating the post-intervention interval, coded as 0 in the pre-interruption period and 1 in the post-interruption period. β₀ represents the baseline level at t=0, β₁ denotes the change in outcome associated with a one-time increase and is regarded as the underlying pre-intervention trend, and β₂ and β₃ represent the changes in the intercept and slope of the trend after the intervention, respectively. Note that although we tentatively use the linear regression formulation expressed in Eq. (1), it can be easily extended to the generalized linear regression formulation.

One issue that has not yet been covered in detail in the prior literature is how the unclear timing of the intervention should be modeled. For example, as we discuss in “Application data analysis” section, consider the case where the statement of emergency (SoE) during the coronavirus disease 2019 (COVID-19) pandemic is the intervention for controlling the spread of infection. Although the date of the SoE declaration is fixed (i.e., nominal timing of the intervention is fixed and clear), the actual effect of the SoE begins long before the date of declaration because the media broadcast the declaration in advance, leading to a change in people’s behavior accordingly. In such a case, determining D_t in Eq. (1) becomes difficult as the actual timing of the intervention is unclear. In other words, the current ITS approach can only be employed when the timing of the intervention is clear and the nominal timing is the same as the actual timing of the intervention. To overcome this limitation, in this study, we aim to develop a new model that combines the ITS model in Eq. (1) with the 1) autoregressive integrated moving average (ARIMA) model and 2) distributed lag functional terms. The motivation to include the ARIMA technique is 1-A) to model autocorrelation, which is frequently observed in time series data but cannot be treated well in Eq. (1), and 1-B) to model the decaying cumulative effect of the intervention. By contrast, the motivation to include the distributed lag functional terms is 2-A) to represent the idea that the effect of the intervention does not start at a fixed time point but that the intervention timing is distributed over a certain interval (thus, the intervention timing is unclear), and 2-B) to model the distributed effect of the intervention over a certain period that includes the timing of the intervention.

The remainder of this study is structured as follows. In “Methods” section, we introduce the idea of the ARIMA ITS model with distributed lag functional terms to model the unclear intervention timing and its distributed effect and then describe the inference procedure. In “Application data analysis” section, we outline practical data analysis procedures of our method for the users by applying several human mobility datasets during the COVID-19 pandemic in Japan. In “Discussion” section, we discuss and propose possible further developments.

Methods

ARIMAX iTS with distributed lag functional terms

In this section, we first explain an ARIMA model with exogenous variables (ARIMAX) [9]. Then, by extending the exogenous variables to our new functional variables (distributed lag functional terms), we propose our model, ARIMAITS-DL, to simultaneously model both the unclear intervention timing and the distributed effect of the intervention.

Define the sampled time series data as $\{Y_{t}, \boldsymbol {X}_{t}\}_{t=1}^{T}$. The general class of the ARIMAX (p,d,q) model without a constant takes the form of; [9, 10]

$$\delta(B)\Delta^{d} Y_{t} = \boldsymbol{X}_{t}\boldsymbol{\beta} + \theta(B)\varepsilon_{t},$$

(2)

or equivalently,

$$\Delta^{d}Y_{t} = \delta(B)^{-1}\boldsymbol{X}_{t}\boldsymbol{\beta} + \frac{\theta(B)}{\delta(B)}\varepsilon_{t},$$

(3)

where B is a lag operator such that ${BY}_{t} = Y_{t-1},\ \delta (B) = 1 - \delta _ {1} B -\dots -\delta_{p}B^{p}=1-\sum_{i=1}^{p}\delta_{i}B^{i},\ \delta_{i}<1,\ \Delta = (1-B),\ \theta (B) = 1 - \theta_ {1} B -\dots -\theta_{q}B^{q}=1-\sum_{i=1}^{q}\theta_{i}B^{i},\ \delta_{i}$ is the parameter for the autoregressive (AR) part, θ_i is the parameter for the moving average part (MA), and p and q are the orders of the AR and MA parts, respectively. For simplicity, we set d=0, and thus, ARIMAX(p,d,q) reduces to ARMAX (p,q). Even when the model where d>0 is needed, the following procedure is still valid after differencing Y_t before fitting the model. The example of an analytical procedure in this case is explained in “Application data analysis” section. When the covariate vector X includes a dummy variable indicating the post-intervention intervals, which corresponds to D_t in Eq. (1), we can extend the simple ITS model to the ARIMA model [11].

To model both the unclear intervention timing and the distributed effect of the intervention over time, we propose the following ARIMAX ITS model with distributed lag functional terms with the covariates X_t, denoted by ARIMAITS-DL (p,q,l₁,l₂). We define T₀ as the nominal timing of the intervention.

$$(1-\sum\limits_{i=1}^{p}\delta_{i}B^{i})Y_{t} = \sum\limits_{k=0}^{l_{1}+l_{2}}w_{k}F_{t-k}+\boldsymbol{X_{t}\beta}+(1-\sum\limits_{i=1}^{q}\theta_{i}B^{i})\varepsilon_{t},$$

(4)

where $\{\delta _{i}\}_{i=1}^{p}, \{w_{k}\}_{k=0}^{l_{1}+l_{2}}, \boldsymbol {\beta }^{T}$, and $\{\theta _{i}\}_{i=1}^{q}$ are unknown parameters,

$$F_{t-k} = \left\{\begin{array}{ll} f(t-k-T_{0}+l_{1}) \ \ & \text{if}\ t\ \in [T_{0}-l_{1}+k, T_{0}+l_{2}] \\ 0 \ \ & \text{Otherwise}, \end{array}\right.$$

(5)

l₁ controls the duration before the effect of the intervention appears (i.e., the start timing of the effect), and l₂ controls the duration of the effect (i.e., the end timing of the effect). Thus, the effect of the intervention is assumed to last from T₀−l₁ to T₀+l₂. Given this formalization, we can model the unclear intervention timing during t∈[T₀−l₁,T₀+l₂]. The distributed lag functional term F_t−k represents how the effect of the intervention is distributed over the time during t∈[T₀−l₁+k,T₀+l₂]. The function f() is a probability density function (pdf) of time t that represents the proportion of the distributed effect of the intervention. Note that Eq. (4) can be considered one special case of the transfer function model popularized by Box and Jenkins [10].

Figure 1 illustrates examples of f(): ex. if the effect is assumed to be distributed uniformly and symmetrically with a peak at T₀ or asymmetrically around T₀, f() might be formulated as a uniform distribution (green line in Fig. 1), (truncated) normal distribution (blue or red lines in Fig. 1), or (truncated) log-normal distribution (purple in Fig. 1), respectively. The choice of f() is discussed in the following section with practical examples.

Fig. 1

Examples of f() corresponding to the proportion of the effect of intervention at time of T₀ between t∈[T₀−l₁,T₀+l₂]

Detailed explanation on distributed lag functional terms

For simplicity, we set p=1 and q=0 and drop X_t from Eq. (4). Thus, our model reduces to the ARX(1,0) model with distributed lag functional terms. Then, Eq. (4) is the reduced simple form that can be expressed as:

$$Y_{t} = \delta Y_{t-1} + \sum\limits_{k=0}^{l_{1}+l_{2}}w_{k}F_{t-k}+\varepsilon_{t},$$

(6)

where |δ|<1 is assumed to make Y_t a stationary process. Now, t is assumed to indicate a day. F_t−k indicates the proportion of the distributed effect of the intervention and w_k estimates the (lagged) effect. Thus, the term w_kF_t−k indicates the distributed lag effect on Y_t. Another advantage of the formulation of Eq. (6) is that because of the AR parameter, δ, we can model the cumulative effect of the intervention from the previous time point. For a clearer understanding, we explain Eq. (6) with examples as follows.

From the definition of F_t−k, when t<T₀−l₁, Eq. (6) becomes

$$Y_{t} = \delta Y_{t-1}+\varepsilon_{t}.$$

(7)

Consider the case where t∈[T₀−l₁,T₀+l₂] (i.e., during the intervention period). When t=T₀−l₁+m (m∈[0,l₁+l₂]), the iterative use of Eq. (6) results in the general form as

$$Y_{t} = \delta^{m-1}Y_{T_{0}-l_{1}-1} + \sum\limits_{j=0}^{m}\sum\limits_{i=0}^{m-j}\delta^{m-j-i}w_{i}F_{T_{0}-l_{1}+j}+\sum\limits_{i=0}^{m}\delta^{m-i}\varepsilon_{T_{0}-l_{1}+i}.$$

(8)

This formulation represents the cumulative effect of the intervention since a day before the intervention becomes effective (i.e., the day of T₀−l₁−1). The first term is the cumulative effect of the outcome $Y_{T_{0}-l_{1}-1}$ with decay δ^k−1, the second term is the delayed and cumulative effect of the intervention, which is explained below, and the third term is an error term.

Finally, we consider the post-intervention period (i.e., t>T₀+l₂). When t>T₀+l₂+n (n∈[0,l₁+l₂]), Eq. (6) becomes

$$Y_{t} = \delta Y_{t-1} + \sum\limits_{k=n}^{l_{1}+l_{2}}w_{k}F_{t-k}+\varepsilon_{t},$$

(9)

and when t>T₀+l₁+2l₂, Eq. (6) becomes the following simple form again as:

$$Y_{t} = \delta Y_{t-1} + \varepsilon_{t}.$$

(10)

Now by using Eq. (6) iteratively, we obtain the general form Eq. (8). For a clearer understanding, we introduce the following example calculations, (T0-T2), which correspond to the first three days after the actual intervention becomes effective (i.e., the actual timing of the intervention while the nominal intervention timing is still T₀). T0. When t=T₀−l₁, the outcome $Y_{T_{0}-l_{1}}$ is affected by the white noise $\varepsilon _{T_{0}-l_{1}}$, the decayed outcome at one day ago $\delta Y_{T_{0}-l_{1}-1}$, where |δ|<1 represents the decay rate and the proportion of the effect of intervention at the first date $F_{T_{0}-l_{1}}=f(0)$, and w₀, which represents the effect size as follows:

$$Y_{T_{0}-l_{1}} = \delta Y_{T_{0}-l_{1}-1} + \underbrace{w_{0}F_{T_{0}-l_{1}}}_{\substack{\text{effect of 0 day} \\ w_{0}:\ \text{0-day-delayed effect} \\ F_{T_{0}-l_{1}}: \%\text{ of effect at 0 day}}} + \varepsilon_{T_{0}-l_{1}}.$$

T1. When t=T₀−l₁+1, the outcome $Y_{T_{0}-l_{1}+1}$ is affected by the white noise $\varepsilon _{T_{0}-l_{1}+1}$, the decayed outcome at $\delta Y_{T_{0}-l_{1}}$, the proportion of the effect of intervention on the first and second dates $F_{T_{0}-l_{1}+1} =f(1), F_{T_{0}-l_{1}}=f(0)$, and w₀ and w₁, which represent the 0- or 1-day delayed effect), respectively, as follows:

$$\begin{array}{*{20}l} Y_{T_{0}-l_{1}+1} &= \delta Y_{T_{0}-l_{1}} + \underbrace{w_{0}F_{T_{0}-l_{1}+1}}_{\substack{\text{effect of}\ +1 \text{day} \\ w_{0}:\ \text{0-day-delayed effect} \\ F_{T_{0}-l_{1}+1}: \%\text{ of effect at \(+1\) day}}}+\underbrace{w_{1}F_{T_{0}-l_{1}}}_{\substack{\text{effect of 0 day}\\ w_{1}:\ \text{1-day-delayed effect}\\ F_{T_{0}-l_{1}}: \%\text{ of effect at 0 day}}} + \varepsilon_{T_{0}-l_{1}+1}. \end{array}$$

By plugging $Y_{T_{0}-l_{1}}$ in T0, we can obtain the following equation:

$$\begin{array}{*{20}l} Y_{T_{0}-l_{1}+1}&=\delta (\delta Y_{T_{0}-l_{1}-1} + w_{0}F_{T_{0}-l_{1}}+\varepsilon_{T_{0}-l_{1}})+w_{0}F_{T_{0}-l_{1}+1}\\ & +w_{1}F_{T_{0}-l_{1}}+\varepsilon_{T_{0}-l_{1}+1}\\ &= \delta^{2}Y_{T_{0}-l_{1}-1} + \underbrace{w_{0}F_{T_{0}-l_{1}+1}}_{\text{effect of}\ +1 \text{day}}+ \underbrace{\underbrace{(\delta w_{0}+w_{1})}_{\substack{\delta w_{0}:\ \text{0-day-delayed effect with decay}\ \delta \\ w_{1}:\ \text{1-day-delayed effect}}}}_{\text{effect of 0 day}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! F_{T_{0}-l_{1}}\\ &+\varepsilon_{T_{0}-l_{1}+1}+\delta\varepsilon_{T_{0}-l_{1}} \end{array}$$

T2. When t=T₀−l₁+2, we obtain $F_{T_{0}-l_{1}+2} =f(2), F_{T_{0}-l_{1}+1}=f(1), F_{T_{0}-l_{1}}=f(0)$. Using Eq. (6) iteratively, the model is reduced to:

$$\begin{array}{*{20}l} Y_{T_{0}-l_{1}+2} &= \delta Y_{T_{0}-l_{1}+1} + \underbrace{w_{0}F_{T_{0}-l_{1}+2}}_{\substack{\text{effect of}\ +2\ \text{day} \\ w_{0}:\ \text{0-day-delayed effect} \\ F_{T_{0}-l_{1}+2}:\ \%\text{ of effect at \(+2\) day}}}+ \underbrace{w_{1}F_{T_{0}-l_{1}+1}}_{\substack{\text{effect of}\ +1 \text{day} \\ w_{1}:\ \text{1-day-delayed effect} \\ F_{T_{0}-l_{1}+1}:\ \%\text{ of effect at}\ +1 \text{day}}} \\ &+\underbrace{w_{2}F_{T_{0}-l_{1}}}_{\substack{\text{effect of 0 day} \\ w_{2}:\ \text{2-day-delayed effect} \\ F_{T_{0}-l_{1}}:\ \%\text{ of effect at 0 day}}}+\varepsilon_{T_{0}-l_{1}+1} \\ &=\delta^{3} Y_{T_{0}-l_{1}-1}+\underbrace{w_{0}F_{T_{0}-l_{1}+2}}_{\text{effect of}\ +2 \text{day}}+\underbrace{\underbrace{(\delta w_{0}+w_{1})}_{\substack{\delta w_{0}:\ \text{0-day-delayed effect with decay}\ \delta \\ w_{1}:\ \text{1-day-delayed effect}}}F_{T_{0}-l_{1}+1}}_{\text{effect of}\ +1 \text{day}} \\ &+ \underbrace{\underbrace{(\delta^{2}w_{0}+\delta w_{1}+w_{2})}_{\substack{\delta^{2} w_{0}:\ \text{0-day-delayed effect with decay}\ \delta^{2} \\ \delta w_{1}:\ \text{1-day-delayed effect with decay}\ \delta \\ w_{2}:\ \text{2-day-delayed effect}}}F_{T_{0}-l_{1}+1}}_{\text{effect of 0 day}} \\ &+\varepsilon_{T_{0}-l_{1}+2}+\delta\varepsilon_{T_{0}-l_{1}+1}+\delta^{2}\varepsilon_{T_{0}-l_{1}}. \end{array}$$

(11)

For example, Eq. (11) provides an intuitive explanation of the cumulative and delayed effect of the intervention. The first term is the decayed outcome at T₀−l₁−1; the second term is the effect on the same day (i.e., T₀−l₁+2), which is decomposed into the 0-day-delayed effect w₀ and the proportion of effect on the same day $F_{T_{0}-l_{1}+2} =f(2)$; the third term is the cumulative and delayed effect from the previous day, which is decomposed into the 1-day-delayed effect w₁, the decayed 0-day-delayed effect δw₀, and the proportion of effect of the previous day $F_{T_{0}-l_{1}+1} =f(1)$; and the fourth term is the cumulative and delayed effect from 2-day ago, which is decomposed into the 2-day-delayed effect w₂, the decayed 0- and 1-day-delayed effects δ²w₀,δw₁, and the proportion of effect from 2-day ago $F_{T_{0}-l_{1}} =f(0)$. This example explains that we can model the proportion of the intervention on a given day and the cumulative and delayed effects of the intervention from the previous days.

Restrictions for estimation

The simplest procedure for estimating the parameters in Eq. (4) is using the maximum likelihood method as with the ordinary ARIMA modeling. Unfortunately, the precision of the estimates of $\{w_{k}\}_{k=0}^{l_{1}+l_{2}}$ is known to sometimes become poor because of the high correlation among $\{F_{t-k}\}_{k=0}^{l_{1}+l_{2}}$, resulting in multicollinearity in the model [12, 13]. To obtain stable estimates of w_ks, we impose the following constraints: We redefine $\sum_{k=0}^{l_{1}+l_{2}}w_{k}F_{t-k}$ in a matrix form as:

$$\begin{array}{*{20}l} \sum\limits_{k=0}^{l_{1}+l_{2}}w_{k}F_{t-k} = \boldsymbol{F}^{T}\boldsymbol{C\eta}, \end{array}$$

(12)

where $\boldsymbol {F}=(F_{t},\dots,F_{t-l_{1}-l_{2}})^{T}\in \mathbb {R}^{l_{1}+l_{2}+1},\ \boldsymbol {C\eta }=(w_{0},\dots,w_{l_{1}+l_{2}})^{T}\in \mathbb {R}^{l_{1}+l_{2}+1},\ \boldsymbol {C}\in \mathbb {R}^{(l_{1}+l_{2}+1)\times h}$ includes the basis variables derived from the specific constraint on w_ks, and $\boldsymbol {\eta }\in \mathbb {R}^{h}$ is a vector of unknown parameters. For example, a constant decline during the lag interval is modeled by:

$$\begin{array}{*{20}l} \boldsymbol{C\eta} = \left(\begin{array}{c} \frac{l_{1}+l_{2}}{l_{1}+l_{2}+1} \\ \vdots \\ \frac{1}{l_{1}+l_{2}+1} \end{array}\right) \eta_{0}, \end{array}$$

(13)

a moving average in the previous L period is modeled by (in the case of L=2):

$$\begin{array}{*{20}l} \boldsymbol{C\eta} = \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \\ 1/2 & 1/2 \\ 1/4 & 3/4\\ \vdots & \vdots \\ \cfrac{2^{l_{1}+l_{2}-1}-(-1)^{l_{1}+l_{2}-1}}{3\times 2^{l_{1}+l_{2}-1}} & 1-\cfrac{2^{l_{1}+l_{2}-1}-(-1)^{l_{1}+l_{2}-1}}{3\times 2^{l_{1}+l_{2}-1}} \end{array}\right) \left(\begin{array}{cc} \eta_{0} \\ \eta_{1} \end{array}\right), \end{array}$$

(14)

and a polynomial smoothing proposed by Schwartz (2000) [11] and Rondeau et al. (2004) [14] is modeled by:

$$\begin{array}{*{20}l} \boldsymbol{C\eta} = \left(\begin{array}{cccccc} 1 & 0 & 0 & \dots & 0 & 0 \\ 1 & 1 & 1 & \dots & 1 & 1 \\ 1 & 2 & 2^{2} & \dots & 2^{s-1} & 2^{s} \\ & & \ddots \\ 1 & (l_{1}+l_{2}) & (l_{1}+l_{2})^{2} & \dots & (l_{1}+l_{2})^{s-1} & (l_{1}+l_{2})^{s} \end{array}\right) \left(\begin{array}{c} \eta_{0} \\ \eta_{1} \\ \vdots \\ \eta_{s} \end{array}\right), \end{array}$$

(15)

where s denotes the degree of the polynomial function. Other examples with non-linear functions, such as splines, can be found in Gasparrini et al. (2010) [13]. These models have been popular for modeling air pollution [15] and temperatures [16].

Modeling steps, the selection of orders (p,q,l ₁,l ₂) and f(), and creating counterfactual

We follow the Box–Jenkins three-stage modeling strategy, including identification, estimation, and diagnostic checking [10]. At the identification stage, a visual inspection of the plots for the time series data allows us to check some important features, such as structural changes, outliers, and missing values, and ensure the stationarity of the time series data. Non-stationary time series invalidates the analyses with the ordinary ARIMA model. To check the stationarity, several tests, such as the Augmented Dickey–Fuller (ADF) test, Ljung–Box test, and Kwiatkowski–Philips–Schmidt–Shin test, can be used [9]. If the time series is judged as non-stationary, a common approach is to subtract successive observations—also known as differencing—to stationarize the time series data. By iteratively applying this test-differencing approach, we can select the order d in Eq. (4). Once d is fixed, Y_t in Eq. (4) can be replaced with the differentiated Y_t.

The plots of the autocorrelation function (ACF) and partial autocorrelation function (PACF) should be visually examined to identify the search range of the orders of the AR and MA parameters: 1) to decide an AR(p) model, the ACF should slowly decrease and PACF should cut off after lag p; 2) to decide an MA(q) model, the PACF should slowly decline and the ACF should cut off after lag q; and 3) to decide an ARMA(p,q) model, both the ACF and PACF should tail off.

At the estimation stage, the maximum likelihood method is used to estimate the regression parameters. The Akaike information criterion (AIC) or Bayesian information criterion (BIC) can be used to select the optimal order of (p,q) among the aforementioned search range defined by the ACF or PACF. In addition, regarding the search space defined by the possible combinations of (l₁,l₂) and f(), they can be specified by the researcher based on the literature search or expert opinion; however, the sensitivity analysis for their choice or model fitting statistics such as AIC or BIC can be helpful. In “Application data analysis” section, we check the sensitivity of the results by changing the combination of (l₁,l₂) and f(). The last step is to check the residuals of the selected model by visual inspection of the residual plot and by testing the presence of autocorrelation using methods such as the Ljung–Box test for white noise. If the autocorrelation is still judged to exist in the residuals, different AR and/or MA orders can be chosen.

Finally, once the regression parameters are estimated, the intervention effect can be estimated by calculating the differences between the observed data and counterfactual prediction. The counterfactual (i.e., Y_t in the absence of the intervention) can be created by substituting 0 into all the components in F_t−ks.

Application data analysis

Data description

In response to the global COVID-19 pandemic, governments implemented large-scale public health interventions, such as lockdowns and SoE declarations, to control the spread of the virus. To quantify the effect of such interventions, human mobility has been frequently measured in terms of tracing human contact in several places, such as restaurants and workplaces. Previous studies have illustrated that the SoE declaration can effectively reduce human mobility [17,18,19]. In this study, we use the human mobility index (HMI) in Japan published by Google [20]. The HMI represents the relative percentage changes in the daily number of visitors (or time spent) from the baseline period (i.e., the same day of the 5-week period between January 3 and February 6, 2020) at six locations: retail and recreational places, grocery and pharmacy stores, public transportation (transit) stations, parks and outdoor spaces, workplaces, and residential areas. The HMI is available daily for all 47 prefectures in Japan. For a simple explanation of our method, we mainly use the HMI at workplaces in Tokyo from February 20 to May 19, 2020. In Tokyo, the first SoE was officially declared on April 7, 2020, although the media broadcast the event even before it was declared, which may have impacted human mobility [21, 22]. In this sense, the nominal intervention timing is clear, although the actual intervention (i.e., the first SoE) timing is unclear. In this section, we demonstrate that the proposed method can handle such cases. Lastly, we use the following covariates in X: the daily COVID-19 test-positive rate, daily number of deaths, and daily (average) temperature.

Practical procedure of data analysis

The data are illustrated in Fig. 2A, where the observed HMI is plotted as a black solid line. To reduce the variation in HMI by day of the week, a seven-day rolling average is calculated in advance. The ADF test indicates that the original time series is non-stationary (p=0.571), and thus, the first difference of the time series data (Fig. 2B) is used to induce stationarity (the ADF test shows p=0.049). Therefore, Y_t is the first difference data, and d=0 is fixed. The ACF and PACF of the stationary (i.e., first-differenced) data are plotted in Fig. 2C and D, respectively. In the figures, the black solid bars above or below the blue dashed lines represent statistically significant autocorrelation with p<0.05. In Fig. 2C and D, we can check that autocorrelation does not exist after lag 1. This implies that the search range of the order for p and q might be around lag 1 [11]. We then search over a series of potential models for the best model with the lowest AIC by using the auto.arima() function in the forecast package in R. Each model is optimized using the maximum likelihood method.

Fig. 2

Application results in Tokyo: A Original time series, B First difference of the data, C Autocorrelation function (ACF) plot for the first difference of the data, and D Partial autocorrelation function (PACF) plot for the first difference of the data

To estimate the effect of the SoE on the HMI in Tokyo, the intervention effect is assumed to be distributed uniformly or normally over time, that is, f() in Eq. (5) is set to the pdf of uniform distribution or (truncated) normal distribution (truncated by l₁ and l₂). In addition, for the constraint on the estimation of w_ks, both Eqs. (14) and (15) are used. Consequently, we have four types of ITS models by combining two f() and two constraints, as follows:

• Normal1: f() is a normal distribution N(0,1) that is truncated by l₁=6 and l₂=6 (i.e., the center is the day of the SoE (April 7, 2020) and truncated at April 1 and 13, 2020). In addition, the constraint for w_ks is the polynomial smoothing model Eq. (15), with s=3.

• Normal2: Similar to Normal1, but f() is truncated by different days, that is, f() is a normal distribution N(0,1) that is truncated by l₁=9 and l₂=3 (i.e., the center is the day of the SoE (April 7, 2020) and truncated at March 29 and April 10, 2020). Moreover, the constraint for w_ks is the polynomial smoothing model, Eq. (15), with s=3.

• Normal3: Similar to Normal2, but f() has a different normal distribution, that is, f() is a normal distribution N(−3,3) that is truncated by l₁=9 and l₂=3 (i.e., the center is 3 days before the SoE (April 4, 2020) and truncated at March 29 and April 10, 2020). Further, the constraint for w_ks is the MA model, Eq. (14).

• Uniform: f() is a uniform distribution with l₁=6 and l₂=6 (i.e., the center is the day of the SoE (April 7, 2020) and truncated at April 1 and 13, 2020). In addition, the constraint for w_ks is the MA model, Eq. (14).

The model with the lowest AIC selected by the algorithm is ARIMA(1,0,0) for all the above models. The residuals for all the models have no significant autocorrelation: all p-values of the Ljung–Box test for white noise are more than 0.9 at six lags. As the data comprise the first difference (Fig 2B), we calculate it back into the original time series (Fig. 2A). To compare our method with the other approaches, the conventional ITS model (denoted by cITS in Figures and Tables) and the model proposed by Schaffer et al. (2021) [11], which is an ARIMA-based ITS model (denoted by ARIMA in Figures and Tables), are included.

Lastly, we validate our method in different COVID-19 datasets by using different prefectures (Osaka and Ehime prefectures: Osaka is a populous and large prefecture and Ehime is a local and small prefecture in Japan), different HMI at “grocery and pharmacy store” and “public transportation (transit) stations”, and different timing of SoE (Osaka at April 7 and Ehime at April 16, 2020). The same procedures described above are used. Figure 3 and the Additional files include the detailed results of these validations.

Fig. 3

Application results in Osaka and Ehime: A Upper: Original time series, Lower: First difference of the data in Osaka B Upper: Original time series, Lower: First difference of the data in Ehime

Application results

Figures 2A and B (colored and dashed lines) indicate the observed data and counterfactual predictions in Tokyo by our ARIMAITS-DL models, assuming the absence of the intervention. Dotted lines indicate the comparison groups: cITS (i.e., the conventional ITS model) and ARIMA (i.e., the model by Schaffer et al. (2021)). Table 1 provides more detailed values of the observed and estimated HMIs. It illustrates that the SoE on April 7, 2020, was successful in reducing the mobility of people: Normal1, Normal2, Normal3, and Uniform models show -392.42 (16.0%), -372.33 (15.0%), -388.61 (16.0%), and -392.20 (16.0%) (cumulative) reduction in HMI, respectively, after the SoE became effective in Tokyo. The difference between the observed and counterfactual peaks around the head of May 2020 for all models, Normal1, Normal2, Normal3, and Uniform models, show a reduction of -18.68 (31.1%), -16.66 (27.7%), -17.12 (28.5%), and -18.74 (31.2%), respectively, on May 3, 2020, which corresponds to the middle period of Japan’s long holiday season. In the comparison groups, the similar tendency are observed, while the degree of reduction in HMI are likely to be underestimated: cITS and ARIMA models show -60.15 (3.0%) and -319.21 (15.0%) (cumulative) reduction in HMI, respectively, after the SoE became effective in Tokyo. Importantly, our findings should only be valid for the study period (i.e., February 20 to May 19, 2020, in Tokyo).

Table 1 Observed and estimated human mobility index (HMI) at workplaces in Tokyo

Figure 3A and B indicate the validation results by using different COVID-19 datasets. While only the result of cITS in Osaka doesn’t seem to estimate the counterfactual values well, the tendency is basically the same as in Tokyo (i.e., Fig. 2): compared to the counterfactual predictions made by our method, the conventional methods tend to underestimate the degree of reduction in HMI. The detailed values of the observed and estimated HMIs in the Additional files.

Discussion

ITS analysis is frequently used to quantify the effects of health interventions on health outcomes at the population level. The most popular formulation of the ITS analysis is the (Gauss–Markov type) linear regression model, Eq. (1) [1, 6, 23]. One of the key assumptions is that the residuals are independent and not correlated. However, this assumption is often violated in time series data. By incorporating dependencies between different time points, the ARIMA model is a possible solution to this problem. In addition, the timing of the intervention is assumed to be clear; however, the actual time at which the effect begins is not always clear in practice. For example, as we described in “Application data analysis” section, the SoE itself was issued on a certain date, T₀, but the actual influence was effective long before the issue of the SoE and distributed over time (i.e., [T₀−l₁,T₀+l₂]), as the media, such as TV, announced its issuance in advance. Thus, the actual timing of the SoE is unclear. To address these issues, we proposed the ARIMA ITS model with the distributed lagged functional terms ARIMAITS-DL. The lagged functional term is tailored to represent how the intervention effect is distributed before and after the nominal timing of the intervention (i.e., the date of the SoE). In addition to the distributed effect, another practical feature of our new model is that we can naturally model the cumulative effect with decaying parameter δ, which is explained in detail in “Detailed explanation on distributed lag functional terms” section. One possible difficulty highlighted by this abundance of choice (the orders (p,d,q), the distribution function f() with (l₁,l₂), and the restriction form for w_ks) is how to choose the best combination. Another possible modeling techniques for the unclear timing of intervention would be the combination of fuzzy set theory and time series models, which have been extensively studied in the field of information science such as Chen (1996) and Singh (2021) [24, 25]. In addition, the fuzzy set theory has been frequently used for analyzing the COVID-19 data [26, 27]. These papers point the way to our next research directions.

In “Application data analysis” section, we used AIC to select the orders (p,d,q), and f(), and (l₁,l₂) were varied for the sensitivity check. However, a priori arguments and expert opinions may be helpful for this choice. A previous discussion concluded that the choice of the distributed lag terms should balance between sufficient complexity to capture detail and sufficient simplicity for interpretability from epidemiological or medical perspectives [28]. As we have no consensus on what is an “optimal” ARIMA model, sensitivity analyses and regression diagnoses such as residual analysis are particularly important to assess the robustness of the key conclusions. The R code, “fuzzyARIMAITS”, for the proposed method are provided in a GitHub repository (https://github.com/kingqwert/R/blob/master/ARIMAITS_DS/fuzzyARIMAITS.R) and will be hosted on the CRAN repository (https://www.r-project.org/) shortly, allowing others to apply our method easily.

Real data were examined with a detailed explanation of the analytical procedure to provide practical insights into the effect of the SoE on human mobility in Japan. We used four models consisting of two distributed lag functions (truncated normal and uniform distribution) and two restrictions on w_ks (polynomial and MA) and confirmed that they exhibit almost the same effect for reducing human mobility in several COVID-19 datasets. Our results indicate that the SoE was successful in reducing the movement of people, at least during the study period (i.e., February 20 to May 19, 2020).

A limitation of this study is that our method was examined in only COVID-19 datasets while using four different settings, two comparison methods, and models with potentially different formulations associated with the distributed lag terms to check the sensitivity of the results. We encourage the re-evaluation of our approach using other datasets. Another limitation is that we assume that a stationary time series is available, which is the requirement of the AR model. The stationary assumption may be incompatible with time series data used in ITS analysis, where trends may change at certain points in time. However, in such a case, we can simply use ARIMA(0,0,0), and the idea of simultaneously using the ITS and distributed lag terms can still be valid. In this case, the model reduces to a simple form, such as Eq. (1), in which autocorrelation can be still modeled by including time as a covariate.

Conclusion

The ITS model has been a powerful study design for evaluating health intervention impacts, and its use has been increasing. The most common formulation for ITS, Eq. (1), is not always adequate, especially when the timing of the intervention is unclear. Our model, ARIMAITS-DL, is a useful tool because it can account for such unclear intervention timing and distributed lag effect with autocorrelation and allows for flexible modeling of different types of impacts.