Interrupted Time Series (Segmented Regression)

An interrupted time series (ITS) evaluates the effect of an abrupt, well-dated, population-level intervention — a formulary change, a label/boxed-warning revision, a guideline release, a copay redesign, a market withdrawal — by treating the population's own pre-intervention trajectory as the counterfactual.

The outcome is aggregated into equally spaced time points (monthly dispensing rate, weekly hospitalization rate, quarterly cost PMPM) and a segmented regression is fit: Y_t = β0 + β1·time_t + β2·level_t + β3·(time since intervention)_t + ε_t, where `time_t` is the running clock, `level_t` is a 0/1 step that switches on at the intervention, and `(time since intervention)_t` is a ramp that starts counting after it. β1 is the pre-intervention secular trend (slope), β2 is the immediate level change at the interruption, and β3 is the change in slope (the post minus pre trend).

The two estimands together describe whether the intervention produced a sudden jump, a gradual bend, both, or neither. ITS is the strongest single-group quasi-experiment because each unit serves as its own control: time-invariant confounders (stable case-mix, baseline access) are differenced out by the within-series before/after contrast.

Core conceptual distinction

ITS identifies an effect from a discontinuity in time, not from a contrast between treated and untreated people. Three modeling choices are separable and must be pre-specified. (1) Level vs slope: an intervention can shift the intercept (β2, e.g., a one-time stockpiling spike), change the trajectory (β3, e.g., slowly declining uptake), or both — reporting only one when both move misstates the effect. (2) Autocorrelation: consecutive points in a series are correlated, so ordinary least squares understates the standard errors and produces falsely narrow intervals.

Either model the error structure with Newey-West heteroskedasticity-and-autocorrelation-consistent standard errors, fit a Prais-Winsten / Cochrane-Orcutt GLS with an AR(1) term, or estimate an ARIMA error process; always inspect the residual ACF/PACF and the Durbin-Watson statistic. (3) Seasonality: monthly health-services series carry strong annual cycles (respiratory admissions, end-of-year deductible effects), which must be removed with harmonic (Fourier sin/cos) terms, calendar-month indicators, or seasonal differencing, or they will masquerade as level/slope effects.

A fourth, design-level strengthening is the controlled ITS (CITS): add a contemporaneous control series not exposed to the intervention and estimate the difference in segmented parameters, which removes any concurrent shock (a coincident flu season, a payment reform) common to both series. This is the bridge between ITS and difference-in-differences.

Pros, cons, and trade-offs

• vs difference-in-differences (`difference-in-differences-staggered-adoption-rwe`): ITS needs only the treated series and leans on extrapolating the pre-trend; DiD requires a parallel comparison group and assumes parallel trends. Prefer ITS when no untreated comparator exists (a national policy hits everyone) and many pre-intervention points are available; prefer DiD when a credible control group exists and the pre-trend is short or noisy. The controlled-ITS variant is essentially a DiD on segmented-regression parameters and inherits the parallel-trends assumption for the difference.
• vs ecological/aggregate cross-sectional studies (`ecological`): ITS uses the same aggregated data but exploits the longitudinal before/after structure, so it is far less vulnerable to static ecological confounding; its cost is the requirement for enough equally spaced pre/post points (a common rule of thumb is >=8 on each side, more with seasonality) and a sharply dated intervention. Prefer ITS whenever the intervention has a clean date and a usable time series exists.
• vs individual-level cohort analysis: ITS answers a population question cheaply and is robust to stable case-mix, but it cannot estimate individual-level effect modification, cannot adjust for time-varying compositional change (a shifting enrollee mix that coincides with the intervention biases β2/β3), and gives no patient-level effect. Prefer a cohort/PS analysis when the question is who benefits and individual covariates are available.

When NOT to use — and when it is actively misleading or dangerous

• The intervention is gradual or its date is fuzzy. Phased rollouts, anticipation effects, and slow diffusion smear the discontinuity; fitting a sharp step/ramp to a gradual change biases both β2 and β3. Model a transition/phase-in window or abandon ITS.
• A co-timed shock is ignored. If another event (a competing drug's withdrawal, a pandemic, a payment reform) hits the series at nearly the same time, the segmented parameters absorb both and attribute the combined effect to the intervention. A single-group ITS cannot separate them — use a controlled ITS with an unaffected comparator.
• Autocorrelation/seasonality left unmodeled. Reporting OLS standard errors on a serially correlated series produces spuriously significant level/slope changes; un-deseasonalized monthly series routinely manufacture phantom effects. This is the most common ITS error and it is actively misleading.
• Too few points or compositional drift. Fewer than ~8 pre points cannot pin the baseline trend; an enrollee-mix change that coincides with the intervention (e.g., a large employer group entering the plan) is confounding the within-series contrast cannot remove.

Data-source operational depth

• Claims (FFS vs MA): Build the series from a stable, continuously enrolled denominator so the numerator rate is not driven by who is observed. A standing failure mode: Medicare Advantage enrollees generate no fee-for-service claims, so if MA penetration is rising over the study window the FFS-observed population shrinks and changes composition, injecting a spurious trend; restrict to FFS Parts A/B (and D for dispensing outcomes) and hold the denominator construction constant across the whole series. Anchor each point on service/admission/fill dates and account for claims-adjudication lag near the data cut (the last 1-3 points are typically incomplete and should be dropped or the cut moved back).
• EHR: Encounter-driven capture means the at-risk denominator fluctuates with system activity; define an explicit "active in system" denominator (>=1 encounter per period) so apparent rate changes are not artifacts of changing observation, and watch for documentation/coding-policy changes (an ICD-9 to ICD-10 transition, a new EHR module) that create a step in the series unrelated to the intervention.
• Registry: Adjudicated outcomes give clean numerators, but reporting completeness and lag vary by calendar period; verify completeness by period and avoid mistaking a reporting-lag artifact for a level change.

Worked claims example

Question: did a 2021 boxed-warning revision reduce monthly initiation of drug class X among adults in a commercial + Medicare FFS database? (1) Series construction: for each calendar month from 2018-01 to 2023-12, count new initiators (first qualifying fill after a 12-month FFS-observable washout) over the continuously FFS-enrolled adult denominator, giving an initiation rate per 1,000 enrollees per month — 36 pre points and 24 post points.

(2) Model: rate_t = β0 + β1·month_t + β2·post_t + β3·months_since_t + seasonal harmonics (sin/cos at 12- and 6-month periods), fit by Prais-Winsten GLS with an AR(1) error to absorb autocorrelation (residual ACF checked, Durbin-Watson near 2). (3) Estimates: β1 = +0.02/month pre-trend (slowly rising uptake); β2 = -0.45 per 1,000 immediate level drop at the warning (95% CI -0.62 to -0.28); β3 = -0.03/month change in slope (95% CI -0.05 to -0.01), i.e., the rising trend reverses to a decline.

(4) Interpretation: the warning produced both a sudden step-down and a sustained downward bend; the counterfactual is the extrapolated pre-trend, so the 24-month cumulative averted initiations are the area between observed and projected pre-trend lines. (5) Strengthening: add a controlled-ITS comparator — an unaffected therapeutic class with similar baseline trend — and re-estimate the difference in (β2, β3) to rule out a co-timed market shock; run sensitivity analyses dropping the last two (lag-incomplete) months and re-fitting with Newey-West SEs to confirm the inference is not autocorrelation-driven.

Decision diagram

flowchart LR
  Pre[Pre-intervention series<br/>>= 8 equally spaced points] --> Fit[Segmented regression]
  Int((Intervention<br/>sharp date)) --> Fit
  Post[Post-intervention series] --> Fit
  Fit --> B1[beta1: pre-trend slope]
  Fit --> B2[beta2: immediate LEVEL change]
  Fit --> B3[beta3: change in SLOPE]
  Fit --> Corr[Model autocorrelation<br/>Newey-West / Prais-Winsten]
  Fit --> Seas[Remove seasonality<br/>Fourier / month indicators]
  B2 --> CF[Counterfactual = extrapolated pre-trend]
  B3 --> CF

flowchart TD
  Q[Abrupt population-level intervention] --> Date{Sharp, well-dated?}
  Date -->|No, gradual/fuzzy| Stop1[Model phase-in or abandon ITS]
  Date -->|Yes| Pts{>= ~8 equally spaced<br/>pre and post points?}
  Pts -->|No| Stop2[Insufficient data for ITS]
  Pts -->|Yes| Shock{Co-timed shock possible?}
  Shock -->|Yes| CITS[Use controlled ITS<br/>add unaffected comparator series]
  Shock -->|No| Single[Single-group segmented regression]
  CITS --> Diag[Check residual ACF/PACF + Durbin-Watson<br/>deseasonalize]
  Single --> Diag
  Diag --> Report[Report beta2 level + beta3 slope with HAC/AR1 CIs]

Worked example

Scenario

A drug safety team wants to know whether a boxed-warning added to Drug X's label on June 1, 2021 changed prescribing behaviour. Using a commercial claims database, they count new patients starting Drug X each month among continuously enrolled adults, then divide by the number of enrolled adults that month to get an initiation rate per 1,000 enrollees. The table below shows five pre-warning months (January–May 2021) and four post-warning months (June–September 2021). The team wants to estimate (1) the immediate level drop at the warning date and (2) whether the monthly trend also shifted.

Dataset

Monthly Drug X initiation rate (new starters per 1,000 continuously enrolled adults). Month index t runs 1–9; intervention falls between t=5 and t=6.

FIG. 1 — DESIGN TIMELINE

Steps

Result

Label

Level change (β2) = −1.5 per 1,000 (immediate drop at the warning); Slope change (β3) = −0.3 per 1,000 per month (trend reversal from +0.2 to −0.1 per month)

Value

Level Change Beta2

-1.5

Slope Change Beta3

-0.3

Pre Trend Per Month

0.2

Post Trend Per Month

-0.1

Counterfactual June

6.0

Observed June

4.2

Arithmetic Check June

4.8 + 0.2*6 + (-1.5)*1 + (-0.3)*1 = 4.2

Arithmetic Check September

4.8 + 0.2*9 + (-1.5)*1 + (-0.3)*4 = 3.9

Trade-offs

Runnable example

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

def its_segmented(df: pd.DataFrame, period: int = 12, hac_lags: int = 12):
    # df columns: time (0..T-1), rate (outcome per period), intervention_period (int index of the step).
    d = df.sort_values("time").reset_index(drop=True).copy()
    k = d["intervention_period"].iloc[0]
    d["level"] = (d["time"] >= k).astype(int)                 # step: 0 pre, 1 post
    d["trend_post"] = np.where(d["time"] >= k, d["time"] - k + 1, 0)  # ramp after interruption
    # Fourier seasonality (annual + semiannual) to avoid seasonal effects leaking into level/slope.
    for h in (1, 2):
        d[f"sin{h}"] = np.sin(2 * np.pi * h * d["time"] / period)
        d[f"cos{h}"] = np.cos(2 * np.pi * h * d["time"] / period)
    formula = "rate ~ time + level + trend_post + sin1 + cos1 + sin2 + cos2"
    # HAC (Newey-West) covariance corrects SEs for autocorrelation in the residuals.
    model = smf.ols(formula, data=d).fit(cov_type="HAC", cov_kwds={"maxlags": hac_lags})
    return model

# Illustrative series: 36 pre, 24 post; warning causes a level drop and a downward slope change.
rng = np.random.default_rng(7)
t = np.arange(60)
base = 5.0 + 0.02 * t + 0.4 * np.sin(2*np.pi*t/12)            # pre-trend + seasonality
effect = np.where(t >= 36, -0.45 - 0.03 * (t - 36 + 1), 0.0)  # step + ramp after month 36
rate = base + effect + rng.normal(0, 0.1, size=60)
df = pd.DataFrame({"time": t, "rate": rate, "intervention_period": 36})

m = its_segmented(df)
print(m.params[["level", "trend_post"]])   # level change (beta2) and slope change (beta3)
print(m.bse[["level", "trend_post"]])       # HAC standard errors

library(prais)     # Prais-Winsten AR(1) GLS for serially correlated errors
library(sandwich)  # Newey-West HAC covariance
library(lmtest)    # coeftest with a supplied vcov

its_segmented <- function(df, period = 12) {
  df <- df[order(df$time), ]
  k <- df$intervention_period[1]
  df$level      <- as.integer(df$time >= k)                 # step: 0 pre, 1 post
  df$trend_post <- ifelse(df$time >= k, df$time - k + 1, 0)  # ramp after interruption
  df$sin1 <- sin(2*pi*df$time/period); df$cos1 <- cos(2*pi*df$time/period)
  df$sin2 <- sin(4*pi*df$time/period); df$cos2 <- cos(4*pi*df$time/period)
  f <- rate ~ time + level + trend_post + sin1 + cos1 + sin2 + cos2

  # Prais-Winsten GLS with an AR(1) error term (models autocorrelation explicitly).
  pw <- prais_winsten(f, data = df, index = "time")
  # Newey-West HAC cross-check on the OLS fit.
  ols <- lm(f, data = df)
  nw  <- coeftest(ols, vcov = NeweyWest(ols, lag = period, prewhite = FALSE))
  list(prais_winsten = summary(pw)$coefficients[c("level", "trend_post"), ],
       newey_west     = nw[c("level", "trend_post"), ])
}

# Illustrative monthly series: 36 pre + 24 post points.
set.seed(7)
t <- 0:59
base   <- 5.0 + 0.02*t + 0.4*sin(2*pi*t/12)
effect <- ifelse(t >= 36, -0.45 - 0.03*(t - 36 + 1), 0)
df <- data.frame(time = t, rate = base + effect + rnorm(60, 0, 0.1),
                 intervention_period = 36)
print(its_segmented(df))

/* Build the segmented-regression design variables. Intervention starts at period k = 36. */
%let k = 36;
data its;
  set work.series;                 /* time (0..59), rate (outcome per period) */
  level      = (time >= &k);                         /* step: 0 pre, 1 post   */
  if time >= &k then trend_post = time - &k + 1;     /* ramp after interruption */
  else trend_post = 0;
  sin1 = sin(2*constant('pi')*time/12); cos1 = cos(2*constant('pi')*time/12);
  sin2 = sin(4*constant('pi')*time/12); cos2 = cos(4*constant('pi')*time/12);
run;

/* PROC AUTOREG fits OLS plus an autoregressive error structure; NLAG=1 = AR(1),
   BACKSTEP prunes nonsignificant autoregressive lags. DWPROB requests the Durbin-Watson test. */
proc autoreg data=its;
  model rate = time level trend_post sin1 cos1 sin2 cos2
        / nlag=1 backstep method=ml dwprob;
  /* level = immediate level change (beta2); trend_post = change in slope (beta3). */
run;

Citations

1. [1]Bernal JL, Cummins S, Gasparrini A. Interrupted time series regression for the evaluation of public health interventions: a tutorial. International Journal of Epidemiology. 2017;46(1):348-355.
2. [2]Zhang F, Wagner AK, Soumerai SB, Ross-Degnan D. Methods for estimating confidence intervals in interrupted time series analyses of health interventions. Journal of Clinical Epidemiology. 2009;62(2):143-148.
3. [3]Kontopantelis E, Doran T, Springate DA, Buchan I, Reeves D. Regression based quasi-experimental approach when randomisation is not an option: interrupted time series analysis. BMJ. 2015;350:h2750.

1. [4]Bernal JL, Cummins S, Gasparrini A. Corrigendum to: Interrupted time series regression for the evaluation of public health interventions: a tutorial. International Journal of Epidemiology. 2021;50(3):1045.