中斷時間序列與差異中的差異（ITS 與 DiD）

一個政策上路了，然後呢

台灣讀者最容易撞到這個問題形狀的地方有三個：健保把某個藥或某項檢查從自費改成給付、醫院導入一套感染管制 bundle、以及 COVID 期間就醫行為整個位移。三件事都長成同一個樣子——一個時點、一個介入、一條逐月的指標，然後有人問「所以降了多少」。

最省事的答法是前後比：介入前平均這麼多，介入後平均那麼多，相除。這份資料上做出來是每月 873 對 571，比值 0.654（95% CI 0.642–0.665），看起來乾淨俐落。

它的問題不是精度不夠，是它把本來就在走的趨勢整段記在政策頭上。指標如果已經在下降，前後比會把那段下降算成政策的功勞；如果指標有季節性、而介入剛好落在旺季與淡季的交界，前後比連方向都可能講反。

這一頁講的兩個方法，差別只在一句話：反事實從哪裡來。

• 中斷時間序列（interrupted time series, ITS）把介入前的趨勢延伸出去，當作「如果政策沒上路，指標會走到哪裡」。反事實是外插來的。
• 差異中的差異（difference-in-differences, DiD）找一條政策打不到的對照序列，用它在同一段時間裡的變化當反事實。反事實是借來的。

兩者都不是憑空變出對照組。ITS 賭的是「趨勢會延續」，DiD 賭的是「兩條線本來會平行地走下去」。這一頁把兩個賭注放在同一份資料上，讓你看到多一條對照序列到底買到了什麼、又要付什麼。

這一頁用的資料

datasets::Seatbelts 是 base R 內建的，data(Seatbelts) 一行就有：英國 192 個月的道路傷亡月資料，從 1969-01 到 1984-12。它剛好同時提供介入序列與天然對照序列：

• front——前座傷亡。1983 年的強制繫安全帶法規範的就是前座。
• rear——後座傷亡。同一批車、同一條路、同一個經濟環境，但法律管不到後座。

還有一件事必須先處理：介入前後的行駛里程從每月 14463 變成 18890，多了 30.6%。如果不把它放進模型，「開得比較多所以撞得比較多」會被折進「法律的效果」裡。所以本頁的每個模型都帶 offset(log(kms))，把結果從「每月幾個人」換成「每單位里程幾個人」——這就是 Poisson 迴歸那一頁講的 offset。

本頁最後一節的 empirical calibration 用的是模擬資料，不是 Seatbelts；那一節開頭會再說明一次。

中斷時間序列：把介入前的趨勢延伸出去當反事實

ITS 就是一個帶了兩個特別項的迴歸。用 Poisson 寫成：

• $t$ 是月份序號，$\beta_1$ 是介入前的長期趨勢。這一項就是前後比缺的東西。
• $\mathrm{law}_t$ 是介入後為 1 的指示變項，$\beta_2$ 是階躍（level change）：介入當下那一跳。
• $\mathrm{tsl}_t$ 是介入後累積的月數（介入前為 0），$\beta_3$ 是斜率變化（slope change）：介入之後每個月比介入前多走多少。本頁程式碼裡的 tsl 在介入當月已經是 1，所以 $\exp(\beta_2)$ 是把不連續點外推到 $\mathrm{tsl} = 0$ 的那個值，介入後第一個月實際觀察到的比值是 $\exp(\beta_2 + \beta_3)$。把 $\exp(\beta_2)$ 報成階躍是慣例，但編碼方式會在階躍與斜率之間搬動一小塊，Methods 要寫清楚用的是哪一種。
• $X_t$ 是隨時間變動的共變項，本頁放的是油價。

階躍與斜率回答的是兩個不同的臨床問題。階躍問的是「這件事一發生就改變了什麼」——法規、給付、開關型的介入通常長這樣。斜率問的是「這件事改變了後續的走向」——衛教、訓練、逐步擴散的介入通常長這樣。一個介入可以只有其中一種，也可以兩種都有。

配適出來的階躍是 RR（rate ratio，發生率比；本頁每一個取過指數的係數都是這個量）0.680（95% CI 0.654–0.707，p < 0.0001），也就是法律上路當下，前座傷亡相對於反事實下降了約 32.0%。斜率變化是每月 RR 1.0123（95% CI 1.010–1.015）——大於 1，代表法律之後的月趨勢比法律之前往上翹。

前後比不是中斷時間序列

把上面的模型跟本頁開頭那個前後比並排：

兩個數字沒有差很多，這其實是這份資料的運氣——介入前的趨勢本來就相對平緩。真正的差別在最右邊那一欄：前後比連「這是一跳還是一個轉向」都答不出來，而這正是臨床上最想知道的那件事。政策如果只是把指標往下推一格、之後照舊往上爬，跟政策把整條走勢扳彎，後續的決策完全不同。

還有一個看不見的代價：前後比沒有辦法被自己的圖抓錯。ITS 的圖上有一條反事實虛線，讀者可以自己判斷那條線延伸得合不合理；前後比只給兩個平均數，沒有任何東西可以被質疑。

沒先扣掉的東西，會被記在政策頭上

這一節是這一頁最實際的一段。同一個階躍，把校正項一項一項拿掉：

第一件值得注意的事跟直覺不合：把油價與里程 offset 拿掉，估計只在 0.660 到 0.682 之間動了 3.4%，信賴區間大幅重疊。 也就是說，在這份資料上，「你有沒有校正油價」根本改不了結論。這件事本來被寫成本頁的主要教學點，是實際跑過之後被自己的數字推翻的。

真正會改變結論的是季節性。加進日曆月之後，階躍從 0.680 變成 0.739（95% CI 0.711–0.768）——與完整校正模型的區間 0.654–0.707 不重疊。AIC 同時從 5204 掉到 2941，跌了 2263，代表季節性是這份資料裡被忽略掉的最大一塊結構。

原因不神秘：道路傷亡有明顯的年內週期（冬季日照短、天候差，傷亡高），而法律在 1983-02 上路，介入後那 23 個月覆蓋的月份組合跟介入前 169 個月的不一樣。沒有把日曆月扣掉，那塊季節差就被記在法律頭上。

安慰劑序列：法律管不到的那條線

ITS 最大的弱點是它沒有對照組——它的對照是同一群人的過去。所以只要有另一個原因剛好在同一個月改變了指標（另一項政策、一次油價衝擊、統計口徑換了），ITS 分不出來。

Seatbelts 給了一個很乾淨的檢查：把完全一樣的模型跑在 rear 上。後座傷亡受同樣的季節、同樣的油價、同樣的里程影響，但法律不管後座。所以它是一個負對照結果（negative control outcome）：如果模型是對的，這裡不應該看到階躍。

完整校正的安慰劑序列給出階躍 RR 0.963，95% CI 0.919–1.008 跨過 1，未達統計顯著。這是這一頁最有價值的一列：法律真的只打到前座，而這件事現在是讀者自己在表上看得到的證據，不是作者的斷言。

措辭要準確：這一列說的是「在這份資料上未偵測到後座的階躍」，不是「後座完全沒有變化」。CI 的下界到 0.919，也就是一個約 8.1% 的真實下降仍然與資料相容。要主張「後座沒有受影響」，需要的是事先設定的等效界線，不是一個跨過 1 的區間。

第三列是附帶的教訓：同一個安慰劑序列，不校正就「顯著」了（RR 0.953，95% CI 0.910–0.998，p = 0.0400）。安慰劑檢定本身也吃模型設定——它不是一個免費的保險，一個設定錯誤的安慰劑檢定會給你假的安心感，也會給你假的警報。

負對照的完整分類（負對照暴露與負對照結果各能反駁什麼、一個負對照不為零時怎麼解讀）在 工具變數那一頁；本頁最後一節談的是同一個想法擴大到幾十個結果之後能做什麼。

差異中的差異：反事實用借的

有了 rear 這條線，可以做的就不只是安慰劑檢定。把兩條序列疊成 long format，加一個組別指示變項，反事實就從「延伸出去」換成「借過來」：

$\alpha_3$——那個交互作用項——就是 DiD 估計。它的意思是：治療組在介入前後的變化，減掉對照組在同一段時間裡的變化。任何同時打到兩組的東西（油價、天氣、經濟、統計口徑）在相減時被消掉，這就是多一條對照序列買到的東西。

三個係數各自的意思值得逐一拆開，因為論文表格常常只印交互作用那一列，讀者不知道另外兩列被拿去做了什麼：

組別主效果 2.182 說的是介入前前座傷亡本來就是後座的兩倍多。DiD 不要求兩組的水準相同，它只要求兩組的變化趨勢平行——這是 DiD 與配對、加權那一類方法最不一樣的地方，也是它常被誤解的地方。

平行趨勢：DiD 唯一的識別假設

DiD 用對照序列換掉了外插，代價是換來一個新的假設：平行趨勢（parallel trends）——如果政策沒有上路，兩條序列的差距會維持不變。

這個假設沒有辦法被驗證，因為它講的是一個沒有發生的世界。能做的是事前檢查：把介入前那段切成若干期，各期估一個「治療組相對於對照組」的值，看它們是不是都貼在 1 附近。這就是 event-study 圖。

結果很明確，而且對這一頁是好消息：13 個非基準年裡有 10 個的區間不含 1，全部為零的聯合檢定是卡方 129.6（13 自由度，p < 0.0001）。1969 到 1972 那四年，前座對後座的比值是基準年的 1.143 到 1.174 倍；要到 1975 之後才大致趨平。

所以本頁的 DiD 估計 0.642 不能被寫成比 ITS 的 0.680 更可信的數字。 誠實的寫法是：事前檢查否定了全期的平行趨勢，只有介入前的後段勉強站得住；如果要用 DiD，應該把樣本限縮到後段、或改用允許趨勢差異的設定，並且把這個決定寫進 Methods。

怎麼算

四個模型全部用 base R 的 glm()，沒有任何額外套件。

data(Seatbelts)
d <- as.data.frame(Seatbelts)
d$time <- seq_len(nrow(d))
d$mon  <- (d$time - 1) %% 12 + 1
d$month <- factor(d$mon)

# 介入 index 從資料推導，不要寫死月份。公文上的生效日與資料庫裡
# 真正開始變化的月份常常差一到兩個月。
law_idx <- which(d$law == 1)[1]

# 介入後累積的月數，介入前為 0。這個編碼讓 law 的係數是「介入當下
# 那一跳」、tsl 的係數是「之後每個月的斜率變化」。原點挪一個月會
# 在兩個係數之間搬動數值，所以要固定寫法。
d$tsl <- pmax(0, d$time - law_idx + 1)

# --- 1. ITS：前座 ---------------------------------------------------
its <- glm(front ~ time + law + tsl + PetrolPrice + offset(log(kms)),
           data = d, family = poisson())
exp(cbind(RR = coef(its), confint.default(its)))["law", ]

# --- 2. 季節性：加日曆月。本頁真正會改變結論的那一項 -----------------
its_season <- glm(front ~ time + law + tsl + PetrolPrice + month +
                    offset(log(kms)), data = d, family = poisson())
exp(coef(its_season)["law"])

# --- 3. 安慰劑序列：一模一樣的模型跑在 rear 上 -----------------------
# 法律不管後座。這裡出現階躍就代表模型抓到的不只是法律。
placebo <- glm(rear ~ time + law + tsl + PetrolPrice + offset(log(kms)),
               data = d, family = poisson())
exp(cbind(RR = coef(placebo), confint.default(placebo)))["law", ]

# --- 4. DiD：兩條序列疊成 long format -------------------------------
long <- rbind(
  data.frame(y = d$front, grp = 1L, d[, c("time", "law", "kms", "PetrolPrice")]),
  data.frame(y = d$rear,  grp = 0L, d[, c("time", "law", "kms", "PetrolPrice")])
)
did <- glm(y ~ grp * law + time + PetrolPrice + offset(log(kms)),
           data = long, family = poisson())
exp(cbind(RR = coef(did), confint.default(did)))["grp:law", ]

# --- 5. 平行趨勢的事前檢查（event study）----------------------------
# 只用介入前的資料，逐年估一個 grp x year 交互作用。
pre <- subset(long, law == 0 & time < law_idx)
pre$year <- 1969 + (pre$time - 1) %/% 12
pre$yf <- relevel(factor(pre$year), ref = "1982")   # 最後一個完整的介入前年
es <- glm(y ~ grp + yf + grp:yf + PetrolPrice + offset(log(kms)),
          data = pre, family = poisson())

# 聯合檢定：所有介入前交互作用同時為零。拒絕就是否定平行趨勢；
# 不拒絕「不」等於平行趨勢成立。
es_null <- glm(y ~ grp + yf + PetrolPrice + offset(log(kms)),
               data = pre, family = poisson())
anova(es_null, es, test = "LRT")

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

d = sm.datasets.get_rdataset("Seatbelts", "datasets").data.copy()
d["time"] = np.arange(1, len(d) + 1)
d["mon"] = (d["time"] - 1) % 12 + 1

law_idx = int(d.loc[d["law"] == 1, "time"].iloc[0])
d["tsl"] = np.maximum(0, d["time"] - law_idx + 1)

# offset 在 statsmodels 是 fit 的參數，不是公式的一部分。
its = smf.glm("front ~ time + law + tsl + PetrolPrice",
              data=d, family=sm.families.Poisson(),
              offset=np.log(d["kms"])).fit()
print(np.exp(its.params["law"]), np.exp(its.conf_int().loc["law"]))

# 季節性：C(mon) 是 patsy 的因子寫法，等同 R 的 factor()
its_season = smf.glm("front ~ time + law + tsl + PetrolPrice + C(mon)",
                     data=d, family=sm.families.Poisson(),
                     offset=np.log(d["kms"])).fit()

# DiD：long format，grp * law 在 patsy 裡與 R 的意思相同
long = pd.concat([
    d.assign(y=d["front"], grp=1),
    d.assign(y=d["rear"], grp=0),
], ignore_index=True)
did = smf.glm("y ~ grp * law + time + PetrolPrice",
              data=long, family=sm.families.Poisson(),
              offset=np.log(long["kms"])).fit()
print(np.exp(did.params["grp:law"]))

ITS 與 DiD 並排：多一條對照序列買到什麼

兩個估計差了 0.038，DiD 那個離 1 更遠，也就是降幅更大（0.642 對 0.680）。差距的來源看得到：DiD 的時期主效果 1.190 大於 1，也就是借來的反事實說「介入後這段期間本來應該上升」，於是前座相對於它就要多降一點。

這不代表 DiD 比較好。 它代表兩個方法對「本來會怎樣」給了不同的答案，而在這份資料上，DiD 那個答案背後的假設沒有通過事前檢查。多一條對照序列買到的是「共同衝擊被消掉」，付出的是「兩條線的相對走勢必須穩定」——而後者在這裡站不住。

實務上常見的做法是兩個都報：ITS 當主分析、DiD 當敏感度分析（或反過來），並且說明兩者的假設不同，一致才是支持，不一致就要說明差在哪。只報其中一個而不提另一個，讀者沒有辦法判斷結論對方法的選擇有多敏感。

負對照與 empirical calibration

單一個負對照結果（像上面的 rear）能回答「有沒有偏誤」。一批負對照可以回答更進一步的問題：偏誤有多大、往哪個方向。

左圖那團點是這一節的全部論點。39 個結果的真值都是 1，如果分析是無偏的，它們應該散布在 1 的兩側；實際上整團偏到右邊，配適出來的虛無分布中心是 RR 1.316，偏離 1 達 31.6%，離散程度（log 尺度的 tau）是 0.117。

更刺眼的是這一個數字：39 個負對照裡有 34 個的 95% CI 不含 1——換句話說，如果照常規讀 p 值，這 34 個「顯著關聯」你會全部相信，而它們每一個的真值都是 1。名目上 95% 的區間，實際涵蓋率只有 12.8%。

校正的做法本身很簡單：把每個估計在 log 尺度上減掉那個偏移量 mu，把標準誤改成原本的 se 與虛無分布離散度 tau 的平方和開根號。結果是三件事同時發生：

1. 負對照的涵蓋率從 12.8% 回到 94.9%，也就是名目值該有的樣子。
2. 目標結果從 2.99（95% CI 2.72–3.29，不含真值 2.0）變成 2.27（95% CI 1.77–2.91，蓋住真值）。
3. 代價是區間寬了 2.6 倍。校正不會讓你更有把握，它讓你的把握變誠實。

讀一篇政策評估論文要問的問題

看到「某政策上路後 X 下降了多少」這種標題時，四個問題可以問掉大部分的疑慮：

1. 對照組是誰？ ITS 的對照是同一群人的過去，DiD 的對照是另一群人的同期。沒有講清楚是哪一種，那個數字就沒有辦法評價。用了 DiD 的話，對照組憑什麼被政策打不到，要寫得出理由。
2. 平行趨勢查了沒？ 有沒有 event-study 圖，還是只有一句「我們檢查過平行趨勢且未達顯著」。後者是把檢定力不足讀成證據。
3. 季節性與長期趨勢扣掉了沒？ 本頁的例子裡，只差一個日曆月項就讓階躍從 0.680 走到 0.739，兩個區間不重疊。
4. 有沒有負對照或安慰劑序列？ 一個「法律管不到」的結果或族群，是最便宜也最有說服力的檢查。

另外兩個常被跳過的細節：介入時點是怎麼定的（公文日期與資料庫裡的轉折常常差一到兩個月），以及有沒有處理自相關（沒有的話標準誤偏窄，區間比實際窄）。

常見誤用

這一頁與其他頁的關係

• 臨床上會用到它的地方——見資料庫與健保申報資料研究。跨年度的趨勢分析必須查給付政策沿革，而政策變更本身就是一個可以用本頁方法估計的介入。
• 底層的模型——見 Poisson 迴歸與發生率。本頁的每一個模型都是帶 offset 的 Poisson GLM，係數的指數就是 RR。
• 另一種「自己當自己對照」的設計——見自我對照設計。它把對照放在同一個人的其他時間窗內，與 ITS 用整條序列的過去當對照是同一個家族的兩種切法。
• 反事實與識別假設的語言——見因果圖與後門準則與目標試驗模擬。平行趨勢是一個關於反事實的假設，跟可交換性一樣不可驗證。
• 負對照的完整分類——見工具變數。
• 時間相關的陷阱——見時間相依共變項。政策評估的資料常常同時有時間相依的暴露與時間相依的干擾。

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B6-10-its-did.R

回報這一頁的問題