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.