28.3 Combining both RDiT and ITS
Combining Regression Discontinuity in Time and Interrupted Time Series in a single regression framework is not standard, but hybrid or augmented approaches can capture both:
- A sharp local discontinuity at the cutoff (\(T^*\)).
- Longer-term trend changes before and after the intervention.
Below are some conceptual strategies for merging RDiT and ITS ideas.
28.3.1 Augment an ITS Model with a Local Discontinuity Term
A more comprehensive ITS can be written using the following form:
\[ Y_t = \beta_0 + \beta_1 T_t + \beta_2 D_t + \beta_3 \bigl(T_t \times D_t\bigr) + \beta_4 P_t + \epsilon_t, \]
where:
- \(T_t\) is the continuous time index
- \(\beta_1\) is the baseline slope (pre-intervention).
- \(D_t = 1\) if \(t \geq T^*\) (0 otherwise)
- \(\beta_2\) is the immediate jump (level change) at \(T^*\).
- \(\bigl(T_t \times D_t\bigr)\) is the interaction term that allows the slope to differ after \(T^*\)
- \(\beta_3\) is the difference in slope post-\(T^*\).
- \(P_t\) is the time elapsed since \(T^*\) (0 before \(T^*\))
- \(\beta_4\) captures the sustained effect of the intervention.
- \(\epsilon_t\) is the error term.
To embed an RDiT-style local discontinuity, define a bandwidth \(h\) around \(T^*\) and add local polynomial terms (e.g., linear) inside that window:
\[ \begin{aligned} Y_t =\;& \beta_0 + \beta_1 T_t + \beta_2 D_t + \beta_3 \bigl(T_t \times D_t\bigr) + \beta_4 P_t \\ &\quad + \alpha_1 \bigl(T_t - T^*\bigr)\mathbf{1}\bigl(\lvert T_t - T^* \rvert < h\bigr) + \alpha_2 \bigl(T_t - T^*\bigr) D_t \,\mathbf{1}\bigl(\lvert T_t - T^* \rvert < h\bigr) + \epsilon_t. \end{aligned} \]
- \(\mathbf{1}(\cdot)\) is an indicator function that is 1 if the time index is within \(h\) of \(T^*\) (and 0 otherwise).
- The \(\beta\) terms capture the global pre/post trends, level changes, and sustained effects.
- The \(\alpha\) terms capture local curvature or a sharper jump within the \(\pm h\) region of \(T^*\).
This approach leverages the entire time series for trend estimation (ITS), while also including local (RDiT-like) estimation of an immediate discontinuity around \(T^*\).
Cautions:
- You need sufficient data globally (for the ITS) and locally near \(T^*\) (for the RDiT terms).
- Extra parameters can lead to complexity and potential overfitting.
28.3.2 Two-Stage (or Multi-Stage) Modeling
An ad hoc but sometimes practical workflow is:
Stage 1: Local RDiT
- Focus on a chosen bandwidth around \(T^*\).
- Fit a local polynomial to capture the immediate jump at \(T^*\).
- Estimate \(\hat{\tau}\) as the jump.
Stage 2: ITS on the Full Series
- Use the entire time series in a segmented (interrupted) regression framework.
- Incorporate \(\hat{\tau}\) explicitly as a known offset or treat it as a prior estimate for the jump.
In practice, you could adjust the outcome by the estimated jump:
\[ Y_t^* = Y_t - \hat{\tau}\,D_t, \]
and then fit the augmented model:
\[ Y_t^* = \beta_0 + \beta_1 T_t + \beta_2 D_t + \beta_3 \bigl(T_t \times D_t\bigr) + \beta_4 P_t + \nu_t. \]
Here, \(Y_t^*\) is the outcome after removing the locally estimated jump from Stage 1. This is not a single unified regression, but a two-step approach combining local (RDiT) and global (ITS) analyses.
28.3.3 Hierarchical or Multi-Level Modeling
A more unified, Bayesian hierarchical approach can also combine RDiT and ITS:
\[ \begin{aligned} Y_{t} &= \underbrace{\beta_0 + \beta_1 T_t + \beta_2 D_t + \beta_3 \bigl(T_t \times D_t\bigr) + \beta_4 P_t}_{\text{Global ITS component}} \\ &\quad\;+\; \underbrace{\alpha_1 \bigl(T_t - T^*\bigr)\mathbf{1}\bigl(\lvert T_t - T^* \rvert < h\bigr) + \alpha_2 \bigl(T_t - T^*\bigr) D_t\,\mathbf{1}\bigl(\lvert T_t - T^* \rvert < h\bigr)}_{\text{Local RDiT component}} \;+\; \epsilon_t. \end{aligned} \]
- Global ITS component captures overall level shifts, slope changes, and sustained effects.
- Local RDiT component captures a sharp jump or local polynomial shape around \(T^*\).
- Additional hierarchical layers (e.g., group-level effects) can handle multiple groups or multiple cutoffs.
Using one of these strategies ensures you capture both the global pre/post-intervention trends (ITS) and any local discontinuities near \(T^*\) (RDiT).
28.3.4 Empirical Example
- Generating synthetic data with a global ITS pattern (immediate jump, slope change, and sustained effect) plus a localized “extra” jump around \(T^* \pm h\) (RDiT-style).
- Fitting a single unified model that includes:
- Global ITS terms (beta coefficients)
- Local RDiT terms (alpha coefficients) within a bandwidth h
- Showing a two-stage approach where we first estimate a local jump and then adjust the outcome.
- Visualizing the results.
# -------------------------------------------------------------------
# 0. Libraries
# -------------------------------------------------------------------
if(!require("sandwich")) install.packages("sandwich", quiet=TRUE)
if(!require("lmtest")) install.packages("lmtest", quiet=TRUE)
library(sandwich)
library(lmtest)
# -------------------------------------------------------------------
# 1. Generate Synthetic Data
# -------------------------------------------------------------------
set.seed(111)
n <- 150 # total number of time points
T_star <- 80 # cutoff/intervention time
t_vals <- seq_len(n) # time index: 1, 2, ..., n
We’ll create a dataset that has:
Baseline slope (pre-intervention)
Immediate jump at \(T^*\)
Post-intervention slope change (ITS style)
Sustained effect
PLUS a local polynomial “extra jump” around \(T^* \pm h=5\) (RDiT style)
h <- 5 # bandwidth for local discontinuity
within_h <- abs(t_vals - T_star) < h # indicator for local region
Create outcome \(Y\) with:
baseline slope: \(0.3 * t\)
immediate jump at \(T^*\): +5
slope change after \(T^*\): +0.2 per unit time (beyond the baseline)
sustained effect: \(0.1*( \text{time since } T^* )\)
local polynomial “extra jump” in the region \([T^* \pm h]\): \((t - T^*) * 2\) only active within \(\pm h\) of \(T^*\)
random noise
Y <- 0.3 * t_vals + # baseline slope
ifelse(t_vals >= T_star, 5, 0) + # immediate jump
ifelse(t_vals >= T_star, 0.2*(t_vals - T_star), 0) + # slope change
ifelse(t_vals >= T_star, 0.1*(t_vals - T_star), 0) + # sustained effect
ifelse(within_h, 2*(t_vals - T_star), 0) + # local polynomial jump
rnorm(n, sd=2) # noise
# Put it in a data frame
df <- data.frame(t = t_vals, Y = Y)
# -------------------------------------------------------------------
# 2. Define Variables for a Single Unified Model
# -------------------------------------------------------------------
df$D <- ifelse(df$t >= T_star, 1, 0) # intervention dummy
df$T_c <- df$t # rename time to T_c for clarity
df$P <- ifelse(df$t >= T_star, df$t - T_star, 0) # time since intervention
# Indicator for local region (± h around T_star)
df$local_indicator <- ifelse(abs(df$t - T_star) < h, 1, 0)
# Center time around T_star for local polynomial
df$t_centered <- df$t - T_star
# -------------------------------------------------------------------
# 3. Fit the Unified "Augmented ITS + RDiT" Model
# -------------------------------------------------------------------
# Y_t = beta_0
# + beta_1 * T_c
# + beta_2 * D
# + beta_3 * (T_c * D)
# + beta_4 * P
# + alpha_1 * (t_centered)* local_indicator
# + alpha_2 * (t_centered)* D * local_indicator
# + epsilon_t
#
mod_hybrid <- lm(
Y ~ T_c + D + I(T_c*D) + P +
I(t_centered * local_indicator) +
I(t_centered * D * local_indicator),
data = df
)
# Robust standard errors
res_hybrid <- coeftest(mod_hybrid, vcov = vcovHC(mod_hybrid, type = "HC1"))
# -------------------------------------------------------------------
# 4. A Two-Stage Approach for Illustrative Purposes
# -------------------------------------------------------------------
# STAGE 1: Local RDiT to estimate extra local jump around T_star ± h
df_local <- subset(df, abs(t - T_star) < h)
mod_local <- lm(Y ~ t_centered*D, data = df_local)
# estimate of the jump at T_star from local model
tau_hat <- coef(mod_local)["D"]
# Adjust outcome by subtracting local jump * D
df$Y_star <- df$Y - tau_hat * df$D
# STAGE 2: Fit a standard ITS on the adjusted outcome Y_star
mod_its_adjusted <- lm(Y_star ~ T_c + D + I(T_c*D) + P, data = df)
res_its_adjusted <-
coeftest(mod_its_adjusted, vcov = vcovHC(mod_its_adjusted, type = "HC1"))
# -------------------------------------------------------------------
# 5. Plot the Data and Fitted Lines
# -------------------------------------------------------------------
par(mfrow=c(1,2))
# Plot 1: Observed data vs. fitted "Hybrid Model"
plot(
df$t,
df$Y,
pch = 16,
xlab = "Time (t)",
ylab = "Outcome (Y)",
main = "Hybrid RDiT+ITS Model"
)
abline(v = T_star, lwd = 2, lty = 2) # cutoff
lines(df$t,
predict(mod_hybrid),
col = "blue",
lwd = 2)
# Plot 2: Adjusted Data (Two-Stage) vs. Fitted ITS
plot(
df$t,
df$Y_star,
pch = 16,
xlab = "Time (t)",
ylab = "Outcome (Y*)",
main = "Two-Stage Approach: Adjusted Y*"
)
abline(v = T_star, lwd = 2, lty = 2)
lines(df$t,
predict(mod_its_adjusted),
col = "red",
lwd = 2)
print(res_hybrid)
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -1.066743 0.453019 -2.3547 0.01989
#> T_c 0.328522 0.009090 36.1409 < 2.2e-16
#> D -19.072974 1.557761 -12.2438 < 2.2e-16
#> I(T_c * D) 0.282215 0.015636 18.0488 < 2.2e-16
#> I(t_centered * local_indicator) 2.197553 0.323646 6.7900 2.726e-10
#> I(t_centered * D * local_indicator) -0.317783 0.399128 -0.7962 0.42723
#>
#> (Intercept) *
#> T_c ***
#> D ***
#> I(T_c * D) ***
#> I(t_centered * local_indicator) ***
#> I(t_centered * D * local_indicator)
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n Two-Stage Approach:\n")
#>
#> Two-Stage Approach:
cat(" (1) Local RDiT => estimated local jump = ", round(tau_hat,2), "\n")
#> (1) Local RDiT => estimated local jump = 3.27
cat(" (2) Standard ITS on adjusted outcome:\n")
#> (2) Standard ITS on adjusted outcome:
print(res_its_adjusted)
#>
#> t test of coefficients:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.553196 0.507846 -1.0893 0.2778
#> T_c 0.308729 0.012701 24.3075 <2e-16 ***
#> D -20.268345 1.889371 -10.7276 <2e-16 ***
#> I(T_c * D) 0.281836 0.019730 14.2844 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
- Beta terms (\(T_c, D, T_c:D, P\)) capture the global ITS components:
- Baseline slope (\(\beta_1\))
- Immediate jump (\(\beta_2\))
- Post-intervention slope change (\(\beta_3\))
- Sustained effect (\(\beta_4\))
- Alpha terms (t_centered * local_indicator, etc.) capture a local RDiT effect.
- In the two-stage approach, \(\hat{\tau}\) is the local jump. Subtracting it (Stage 1) yields a ‘cleaned’ \(Y^*\), and we then fit a simpler ITS (Stage 2).
28.3.5 Practical Guidance
- Data Requirements: You must have enough data in the local window around \(T^*\) for RDiT and enough pre- and post-intervention observations for ITS.
- Avoid Over-Parameterization: Adding many interaction and polynomial terms can quickly increase complexity.
- Interpretation: A single unified model that merges local jumps with global trends may be difficult to interpret. Distinguish clearly which parameters capture the sharp discontinuity versus the long-run trend.
In practice, a hybrid approach can work if you truly believe there is both a local immediate jump and a longer-run trend change that standard ITS alone may not fully capture. However, weigh the added complexity against the quality and richness of your data.