Interrupted Time Series (Segmented Regression)

In plain language

An interrupted time series study watches a population's rate of some outcome — such as how many patients per month start a new drug — over many months, then asks whether a specific event (a safety warning, a policy change, a formulary switch) caused that rate to change. The key idea is that the population's own past trend, if the intervention had never happened, acts as the stand-in for 'what would have happened anyway.' The method separately estimates two things: did the rate jump or drop instantly the moment the event occurred, and did the long-run direction of the trend bend after that? One honest limitation: if another unrelated event happened at almost the same time, the method cannot separate its effect from the intervention's without adding a second comparison series.

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 to use

A single, abrupt, sharply dated population-level intervention; an outcome that can be aggregated into a regularly spaced series with enough pre-intervention points to characterize the baseline trend and seasonality; no clean individual-level comparator (so a cohort/PS design is infeasible); decision-makers who need a transparent before/after policy evaluation (FDA safety-label impact, PQA/CMS quality-measure or formulary changes, HTA coverage-policy assessments). ITS is the default design for drug-utilization and safety-communication impact studies in claims.

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.

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.

Steps

• Fit the pre-period data (t = 1 to 5) to a straight line. The line starts at 5.0 in January and rises by 0.2 per month, so the formula is: expected rate = 4.8 + 0.2 × t.

• Extend that line forward to predict what June (t = 6) would have been without the warning: 4.8 + 0.2 × 6 = 6.0 per 1,000.

• The observed rate in June is 4.2. The gap between the counterfactual (6.0) and what was actually observed at t = 6 breaks into two parts: the level change and the first month of slope change.

• The level change (β2) is −1.5 per 1,000 — a sudden, one-time drop attributable to the warning announcement itself.

• The slope change (β3) is −0.3 per 1,000 per month. This means the trend shifts from +0.2/month (rising) to −0.1/month (slowly falling) after the warning.

• Verify June: counterfactual 6.0 + level change (−1.5) + slope change × months_since_warning 1 (−0.3 × 1 = −0.3) = 6.0 − 1.5 − 0.3 = 4.2. Matches the observed rate.

• Verify September (t = 9, months_since_warning = 4): counterfactual 4.8 + 0.2×9 = 6.6; + (−1.5) + (−0.3×4 = −1.2) = 6.6 − 1.5 − 1.2 = 3.9. Matches the observed rate.

• The upward pre-period trend (+0.2/month) has reversed to a downward post-period trend (−0.1/month), meaning the warning changed not just the level but the direction of prescribing.

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.26 + (-1.5)1 + (-0.3)*1 = 4.2

Arithmetic Check September

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

Timeline Spec

Title

Drug X monthly initiation rate: boxed-warning interruption, Jan–Sep 2021

Caption

Monthly initiation rate (per 1,000 enrollees) before and after the June 2021 boxed warning. The dashed line shows the counterfactual trend if the warning had never been issued. The solid orange line shows observed post-warning rates. The gap between them widens each month because the slope also changed.

Alt Text

Segmented timeline showing Drug X initiation rate rising from 5.0 to 5.8 per 1,000 during January–May 2021, a vertical intervention marker at June 2021 labeled Boxed Warning Added, then the observed rate dropping to 4.2 in June and declining further to 3.9 by September, while a dashed counterfactual line continues rising from 5.8 toward 6.6, illustrating both a level drop and a slope reversal.

Window

Start

2021-01

End

2021-09

Label

9-month observation window (5 pre, 4 post)

Events

• Label

  Boxed Warning Added

  Date

  2021-06

  Kind

  intervention

  Description

  FDA boxed warning added to Drug X label — the 'interruption' point in the series

Spans

• Kind

  pre_intervention

  Start

  2021-01

  End

  2021-05

  Label

  Pre-intervention: rising trend (+0.2/month), rate 5.0 → 5.8

• Kind

  post_intervention

  Start

  2021-06

  End

  2021-09

  Label

  Post-intervention: reversed trend (−0.1/month), rate 4.2 → 3.9

Result

Label

Level drop = −1.5 per 1,000 at warning; slope reversed from +0.2 to −0.1 per month

Value

Level Change

-1.5

Slope Change

-0.3

Pre Slope

0.2

Post Slope

-0.1

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))