進階已經雙重審閱，尚未人工抽查

群集與重複測量資料

同一個人量四次、同一家醫院收一百個病人——這些列彼此不獨立。把它們當成獨立列會怎樣：組間效果的標準誤太小、組內效果的標準誤太大，方向相反。混合模型、GEE 與 cluster-robust 標準誤各自回答的是不同的問題。

報表看不出來的那一條假設

迴歸模型診斷列了四條假設，前三條都有圖可以看：殘差對預測值、QQ 圖、影響點。第四條——觀察彼此獨立——那一頁只寫了「這一條看不出來」，然後就沒有下文了。Poisson 迴歸兩次提到「要用 GEE 或混合效應模型」，但站上一直沒有那一頁。這一頁把這三處補起來。

獨立性看不出來，是因為它不是資料的性質，是資料怎麼來的性質。同一份 CSV，每一列是一個病人，還是同一個病人的第三次回診，檔案裡長得一模一樣。你必須從研究設計知道這件事，模型不會告訴你。

臨床研究裡違反獨立性的情境比想像中多：

重複測量：同一個病人在術後 30 分鐘、90 分鐘、4 小時、隔天早上各量一次疼痛
群集隨機：以診所或病房為單位隨機分配，同一間診所的病人共享醫師、流程與病人組成
多中心：同一家醫院的病人在照護品質上比跨院的病人相似
配對設計：左右眼、雙胞胎、配對的病例與對照
多階段抽樣：先抽學校再抽學生

共通點是資料有兩層：人（或診所）這一層，以及人裡面的測量這一層。所有的麻煩都來自「同一群裡的觀察比不同群之間的觀察更像」。

這一頁的例子

medicaldata::licorice_gargle 是一個插管前用甘草含漱液預防術後喉嚨痛的隨機對照試驗。 235 位病人隨機分成對照組 117 位與甘草組 118 位，喉嚨痛評分在四個時點各記錄一次。

把寬表攤成長表之後有 940 列，其中 8 個評分是缺的。

個體軌跡圖（spaghetti plot）。橫軸是四個時點（PACU 30 分鐘、PACU 90 分鐘、術後 4 小時、術後第一天早上），縱軸是喉嚨痛評分，觀察到的範圍是 0 到 7。背景是 233 條細線，每條是一位病人的四次評分，依組別上色並加了微小的垂直位移以免重疊；大多數細線貼在 0 附近的低分帶，少數幾條在 4 到 7 之間上下跳動。前景兩條粗線是各組的平均軌跡，附 95% 信賴區間的誤差線：對照組（藍）四個時點的平均從 1.03 緩降到 0.65，整段都在甘草組（紅）之上；甘草組從 0.27 先降到最低的 0.14，在術後 4 小時回升到最高的 0.35，最後是 0.32。兩條平均軌跡在整個時間範圍內沒有交叉。 — 每條細線是一位病人。線與線之間不獨立——同一條線上的四個點來自同一個人，這正是天真模型假裝不存在的東西。產圖腳本 `figures/scripts/B2-07-clustered.R`

這張圖要看兩件事。 第一，同一個人自己跟自己很像——這就是群集。不過這件事不能只靠看圖判斷，而且這份資料剛好示範了為什麼：畫面上佔多數的是貼在 0 分附近那一大群人——233 位病人裡有 122 位四次評分完全一樣，其中 120 位是四次都 0 分，他們的線當然是平的。而分數高的那幾條其實上下跳得很厲害，只是全距達到 4 分以上的總共只有 8 位。所以「線大多是平的」這個印象是那群 0 分的人撐起來的，不是資料真的沒有波動。真正撐起「同一個人很像」這句話的是後面會算出來的 組內相關係數（intraclass correlation coefficient, ICC）0.540，不是視覺印象。個體軌跡圖負責讓你看見變異的形狀，量化這件事要靠變異成分。第二，個體線太擠，看不出組別差異，所以還要疊上平均軌跡——但平均軌跡是用全部資料算的，不是只用畫得出來的那幾條。個體線負責告訴你變異有多大，平均軌跡負責告訴你效果在哪裡，兩者缺一不可。只畫平均軌跡的圖會讓讀者以為每個人都長那樣。

動手跑一次

library(medicaldata)
library(nlme)

lg <- licorice_gargle
lg$id <- seq_len(nrow(lg))

tp <- c("pacu30min_throatPain", "pacu90min_throatPain",
        "postOp4hour_throatPain", "pod1am_throatPain")

# 寬表攤成長表：一列一個「病人 x 時點」
long <- do.call(rbind, lapply(seq_along(tp), function(i)
  data.frame(id = lg$id, treat = lg$treat, t = i, pain = lg[[tp[i]]])))
L <- long[!is.na(long$pain), ]
L$idf <- factor(L$id)

# 先看缺值的形態，再開始分析
table(rowSums(is.na(lg[, tp])))

# 天真模型：把 932 列當成 932 個獨立觀察
summary(lm(pain ~ treat, data = L))

# 混合模型：每個病人一個隨機截距
m <- lme(pain ~ treat, random = ~ 1 | idf, data = L, method = "REML")
summary(m)
VarCorr(m)                      # 變異成分，ICC 從這裡算

# 第三種做法：每人先取平均，再做兩樣本比較
pm <- aggregate(pain ~ id + treat, data = L, FUN = mean)
summary(lm(pain ~ treat, data = pm))

# 組內對比：時間點
summary(lm(pain ~ factor(t) * treat, data = L))                    # 天真
summary(lme(pain ~ factor(t) * treat, random = ~ 1 | idf, data = L))

# cluster-robust（三明治）標準誤，手寫版：
#   (X'X)^-1 [ sum_g X_g' u_g u_g' X_g ] (X'X)^-1
# 一個病人一個 cluster，最後乘上 CR1 小樣本修正。
cluster_robust_se <- function(fit, cluster) {
  X  <- model.matrix(fit)
  u  <- as.numeric(residuals(fit))
  cl <- factor(cluster)
  G <- nlevels(cl); N <- nrow(X); K <- ncol(X)
  bread <- solve(crossprod(X))
  meat  <- matrix(0, K, K)
  for (g in levels(cl)) {
    idx <- which(cl == g)
    sg  <- crossprod(X[idx, , drop = FALSE], u[idx])   # 這個 cluster 的分數向量
    meat <- meat + tcrossprod(sg)
  }
  adj <- (G / (G - 1)) * ((N - 1) / (N - K))           # CR1
  setNames(sqrt(diag(bread %*% meat %*% bread * adj)), colnames(X))
}

cluster_robust_se(lm(pain ~ treat, data = L), L$id)
cluster_robust_se(lm(pain ~ factor(t) * treat, data = L), L$id)
# 實務上用 sandwich::vcovCL(fit, cluster = L$id, type = "HC1")，本站不裝套件所以手寫

驗證環境：R 4.6.0 + nlme 3.1.169 + medicaldata 0.2.0。nlme 隨 R 本體附帶，不必另外裝；lme4 的 lmer() 語法不同但結果相同。

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

lg = sm.datasets.get_rdataset("licorice_gargle", "medicaldata").data
lg["id"] = np.arange(len(lg))

tp = ["pacu30min_throatPain", "pacu90min_throatPain",
      "postOp4hour_throatPain", "pod1am_throatPain"]

L = lg.melt(id_vars=["id", "treat"], value_vars=tp,
            var_name="tp", value_name="pain").dropna(subset=["pain"])
L["t"] = L["tp"].map({v: i + 1 for i, v in enumerate(tp)})

print(smf.ols("pain ~ treat", data=L).fit().summary())        # 天真

md = smf.mixedlm("pain ~ treat", data=L, groups=L["id"])      # 隨機截距
print(md.fit(reml=True).summary())

pm = L.groupby(["id", "treat"], as_index=False)["pain"].mean()
print(smf.ols("pain ~ treat", data=pm).fit().summary())

# cluster-robust：OLS 配以病人為 cluster 的三明治變異數
print(smf.ols("pain ~ treat", data=L)
        .fit(cov_type="cluster", cov_kwds={"groups": L["id"]}).summary())

statsmodels 的 MixedLM 對應 nlme 的 lme()；預設是 REML，與 R 一致。

方向一：組間效果的標準誤太小

先問最簡單的問題：甘草組的喉嚨痛評分比對照組低多少？

天真的做法是把 932 列丟進 lm()。點估計是 -0.582 分（甘草組較低），標準誤 0.0686， t = -8.48，p < 0.001，殘差自由度 930。

混合模型的點估計完全一樣，但標準誤是 0.1113， t = -5.23，p < 0.001，treat 的自由度掉到 231。

天真的標準誤比正確的小了 38.3%。

兩個並排的長條圖，兩邊共用同一條縱軸（標準誤，0 到 0.20）。左圖「組間：治療」有三根長條：天真 lm 最矮（0.0686），cluster-robust（0.1115）與混合模型（0.1113）幾乎一樣高，都比天真的高出約 62%。右圖「組內：時間」有三組長條，對應三個時間點對比；每一組裡天真 lm 都是最高的（0.137），混合模型都是 0.092，明顯比天真的矮——與左圖的方向相反。cluster-robust 在右圖三組之間高低不一，在第一組（PACU 90 分鐘）甚至比混合模型還矮。 — 同一份資料、同一個群集結構。左邊組間效果的標準誤被低估，右邊組內效果的標準誤被高估。方向是相反的。產圖腳本 `figures/scripts/B2-07-clustered.R`

方向二：組內效果的標準誤太大

多數教科書講到上一節就停了，於是留下一個容易記錯的印象：「忽略群集會讓標準誤太小」。 這句話只對一半。

把模型換成 pain ~ factor(時間) * treat，看時間點之間的對比：

項	類別	估計	天真 SE	cluster-robust SE	混合模型 SE	混合 / 天真
截距：對照組在 PACU 30 分鐘	截距	1.026	0.097	0.144	0.097	1.000
PACU 90 分鐘 vs PACU 30 分鐘	組內（時間）	-0.207	0.137	0.064	0.092	0.673
術後 4 小時 vs PACU 30 分鐘	組內（時間）	-0.112	0.137	0.133	0.092	0.673
術後第一天早上 vs PACU 30 分鐘	組內（時間）	-0.379	0.137	0.127	0.092	0.673
甘草含漱 vs 對照，在 PACU 30 分鐘	組間	-0.752	0.137	0.157	0.137	1.000
治療 × PACU 90 分鐘	組內（時間 × 治療）	0.070	0.194	0.076	0.130	0.673
治療 × 術後 4 小時	組內（時間 × 治療）	0.189	0.194	0.159	0.130	0.673
治療 × 術後第一天早上	組內（時間 × 治療）	0.422	0.194	0.155	0.130	0.673

最後一欄告訴你發生了什麼事。treat 與截距那兩列的比值是 1（混合模型與天真模型給出完全相同的標準誤，因為在這個交互作用模型裡 treat 是 PACU 30 分鐘那一個時點的組間對比，每人只貢獻一列，群集沒有作用的空間）。其餘每一列的比值都是 0.673—— 混合模型的標準誤比天真的小了 32.7%。

這才是混合模型真正的價值。

組間比較（甘草 vs 對照）要跨越病人，而病人之間的差異不會因為多量幾次而消失。所以重複測量在這裡幫得上忙，但幫助有限，而且有一個算得出來的天花板。

用同一個模型的變異成分就可以把這條曲線畫出來（平衡設計下，每人平均值的變異數是 $\sigma^2_b + \sigma^2_e/m$ ， $m$ 是每人測量次數）：每人只量一次，組間差異的標準誤是 0.1375；量四次降到 0.1113， 窄了 19.1%；就算量無限多次，也只能降到 0.1010。換句話說，四次已經拿到了可得增益的 71.9%—— 再往上加測量次數擠不出多少，因為天花板是由病人之間的變異決定的，而那一塊不歸重複測量管。

組內比較（PACU 90 分鐘 vs PACU 30 分鐘）在同一個人身上進行，所以「這個人天生比較痛」那部分變異被抵銷掉了。天真模型不知道兩次測量來自同一個人，只好把那份個體差異全部算進誤差項，於是分母被灌水、標準誤被高估。

具體後果：PACU 90 分鐘 vs PACU 30 分鐘這個對比，天真模型的 p 0.132（未達統計顯著），混合模型的 p 0.025。正確處理群集不是「懲罰」，是把資訊還給你。 這也是為什麼 cross-over 試驗與前後測設計比平行組設計省樣本。

混合模型在做什麼

lme(pain ~ treat, random = ~ 1 | 病人) 多了一項：每個病人有自己的截距 $b_i$ ，而這些 $b_i$ 被假設來自一個平均為 0 的常態分布。模型變成

y_{ij} = \beta_0 + \beta_1 \, \text{treat}_i + b_i + \varepsilon_{ij}, \qquad b_i \sim N(0, \sigma^2_b), \quad \varepsilon_{ij} \sim N(0, \sigma^2_e)

$b_i$ 是隨機效果（random effect）：它不是一個要估計的參數，而是一個要估計的分布。這是它與「把病人 id 當成 233 個 level 的固定效果因子」的差別—— 後者要燒掉 232 個自由度，而且組間效果會被整個吃掉（treat 在病人內不變，與病人 id 完全共線）。

那個標準誤，其實就是「每人取平均再比較」

把每位病人的四次評分先平均成一個數字，得到 233 列，再做最普通的兩樣本比較。標準誤是 0.1113—— 與混合模型的 0.1113 相等到小數點後九位（兩者相差 2.4e-10，那是 REML 最佳化的收斂容差，不是估計量的差別）。

這個等式成立的前提是平衡設計，請不要把它記成普遍規律。

這份資料剛好每位病人都有完整的四次測量（缺的兩位是整個人被排除，不是缺其中幾次）。在這種情況下，隨機截距模型對組間效果做的事情，數學上就等於「每人取平均再比較」。

一旦缺值散在各個時點上，等式立刻不成立：混合模型會依每個人有幾次測量給不同的權重（測得多的人權重高），而「每人取平均」把所有人當成一樣可靠。這時候兩者的標準誤會分家，而混合模型的那個才是對的。

這個等式的教學價值不在於它可以拿來取代混合模型，而在於它告訴你混合模型在組間效果上並不神秘：它沒有變出額外的資訊，它只是拒絕把 932 列當成 932 個人。同樣的道理，一篇論文如果重複測量的組間效果報得比「每人取平均」還精確，那是需要解釋的，不是加分。

ICC：多少變異來自「人和人不同」

混合模型把總變異拆成兩塊。病人之間的變異數是 0.5947，同一個病人內部的變異數是 0.5068。 組內相關係數（intraclass correlation coefficient, ICC）是前者佔總和的比例：

\mathrm{ICC} = \frac{\sigma^2_b}{\sigma^2_b + \sigma^2_e}

算出來是 0.540。

這個數字有一個容易漏掉的限定條件。 上面那兩個變異數來自 lme(pain ~ treat, random = ~1 | 病人)，模型裡已經放了 treat，所以它們拆的是扣掉組別差異之後剩下的變異，不是原始評分的總變異。正確的唸法是：在同一組之內，病人之間的差異略大於同一個人在不同時間的波動（0.540 對 0.460）。這個限定條件在這裡是有份量的—— 0.540 只比一半高一點，換一個模型設定就可能落到另一邊。報告 ICC 時一定要說它是從哪個模型的哪個變異拆解來的。

ICC 也可以讀成「同一個病人的任兩次測量之間的相關係數」——這兩種說法在隨機截距模型下是同一件事。

同一條公式還有另一個用途，而論文裡的「ICC = 0.85」多半指的是那一個。 這一頁的 ICC 問的是「一個群裡的兩個人有多像」，量的是集群結構；量測信度的 ICC 問的是「同一個對象被量兩次有多像」，量的是儀器或判讀者。變異拆解是同一條，問題不同，而且信度那一側還分成六種型式（單次 vs 平均量測、consistency vs absolute agreement、rater 視為隨機 vs 固定），不指明是哪一種就無法解讀。細節見量測信度與 ICC。

設計效應：ICC 怎麼換算成「浪費了多少樣本」

樣本數與檢定力提到集群隨機試驗要乘上設計效應 （design effect） $1 + (m-1)\times \mathrm{ICC}$ ，其中 $m$ 是每個集群的人數。這一頁的「集群」就是病人， $m$ 是每人的測量次數 4。

代進去得到設計效應 2.620。也就是說 932 列在組間比較上，只相當於約 356 個獨立觀察。

這個數字要從兩個方向讀，才不會讀成單純的損失。 往下看， 932 列縮水成 356，那是把重複測量當成獨立觀察的代價；往上看， 356 比病人數 233 多，而那個差額正是上一節說的「重複測量確實買到了東西，只是有天花板」。

混合模型、GEE、cluster-robust：三個不同的問題

處理群集有三條主流路線——混合模型、GEE（generalised estimating equations，廣義估計方程）與 cluster-robust 標準誤。它們不是同一件事的三種算法，而是回答三個不同的問題。這一節是本頁最容易被跳過、也最常被誤用的地方。

做法	群集怎麼被使用	估的是什麼效果	本頁的組間 SE
天真 `lm()`	完全忽略	邊際效果，但變異數是錯的	0.0686
cluster-robust（三明治）標準誤	只用在變異數，點估計仍是 OLS	邊際（族群層次）效果	0.1115
GEE（exchangeable 工作相關）	變異數與估計式的加權都用	邊際（族群層次）效果	與 cluster-robust 同量級
混合效應模型（隨機截距）	寫進模型本身	條件（受試者層次）效果	0.1113

條件效果 vs 邊際效果的差別，用一句話問：

混合模型的係數回答「同一個人如果換到另一組，他的評分會差多少」—— 它是在「這個人的隨機截距 $b_i$ 固定住」的條件下定義的。
GEE 與 cluster-robust 的係數回答「整個族群如果全部改用甘草，平均評分會差多少」。

在這一頁的線性模型裡，兩者的點估計相同（-0.582），所以差別看不出來。但那是線性模型的特例。 一旦效果量是不可塌縮的（最典型的就是 logit），條件效果與邊際效果的數值就不一樣，而且條件 OR 一定離 1 更遠——這就是邏輯迴歸那一頁講的 不可塌縮性（non-collapsibility）。但不是每個非線性連結都這樣：log 連結是可塌縮的， Poisson 混合模型與 Poisson GEE 估到的是同一個發生率比，差別只在截距。同一份資料，混合邏輯模型與 GEE 會給出不同大小的 OR，兩個都對，只是回答的問題不同。

怎麼選：

想講「這個病人接受治療會怎樣」（臨床決策、個體預測）→ 混合模型
想講「這個政策推到全族群會怎樣」（公衛、政策評估）→ GEE 或 cluster-robust
群集只是麻煩、不是研究興趣所在（例如多中心試驗的中心效應）→ cluster-robust 最省事，它不要求你把相關結構猜對

cluster-robust 不是「乘上一個固定倍數」

上一節的表格裡，cluster-robust 那一欄的行為值得單獨看一眼。組間效果的 cluster-robust 標準誤 0.1115 與混合模型的 0.1113 幾乎重合。但三個時間點對比的 cluster-robust 標準誤分別是 0.064、 0.133 與 0.127 ——彼此差很多，而且第一個（PACU 90 分鐘 vs PACU 30 分鐘）比混合模型的 0.092 還要小。

原因是三明治估計不強加任何相關結構。隨機截距模型假設同一個人的任兩次測量相關係數都一樣（exchangeable），於是所有組內對比共用同一個修正倍數 0.673。三明治估計直接從資料裡的殘差去估每個對比實際的變異，相鄰時點相關性高的對比就更精確，隔得遠的就沒那麼精確。

這件事的兩面。 三明治估計不需要你猜對相關結構，這是優點；但它是大樣本近似，而且近似的品質取決於群集數而不是總列數。一般的經驗法則是群集數少於 40 左右就要小心，這時候它會低估變異。本頁有 233 個群集，算是充裕；但一個「十二家醫院、每家五百位病人」的多中心研究，群集數是十二不是六千，直接套 cluster-robust 會過度樂觀。

本頁的三明治估計已經帶了 CR1 修正 （乘上 (G/(G-1)) * ((N-1)/(N-K))， $G$ 是群集數、 $N$ 是列數、 $K$ 是參數個數），這是最常見的預設，Stata 與多數套件都用它。但 CR1 只處理自由度、不處理群集數太少——群集少的時候要換成 CR2 或 CR3 這兩種不同的修正，或者乾脆改用混合模型。看到論文只寫「使用 robust standard errors」而沒說哪一種修正、群集數是多少，那個標準誤是查不了帳的。

常見誤用

誤用	為什麼錯
重複測量的資料直接丟進 `lm()` / `glm()`	組間效果的標準誤會太小；自由度是最快的線索
記成「忽略群集一律讓標準誤太小」	組內對比的方向相反，會被高估
只畫平均軌跡不畫個體線	讀者無法判斷變異多大，也看不出有沒有幾條線在拉平均
把病人 id 當成固定效果因子	燒掉大量自由度，而且組間效果與它完全共線、估不出來
直接比較混合模型與 GEE 的 OR 大小	一個是條件效果、一個是邊際效果，在 logit 這種不可塌縮的連結下本來就不同
群集數很少仍套 cluster-robust	三明治估計靠群集數做近似，群集少會低估變異
把設計效應套在組內對比上	那裡重複測量是賺的，不是賠的
缺值形態沒看就開始配模型	平衡與不平衡的分析性質不同，某些簡便做法只在平衡時等價
用「每人取平均」取代混合模型	只在完全平衡時等價；不平衡時會給每個人錯誤的權重
沒有報 ICC	讀者無從判斷群集有多嚴重，後續研究也無法用它估樣本數
隨機效果只放截距就當作處理完了	若效果本身隨人而異（例如各人的時間趨勢不同），要放隨機斜率
把 p 值變小歸功於「用了比較好的模型」	要說清楚是組內對比拿回了被誤算成誤差的個體差異

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-07-clustered.R

讀讀看這張圖

答案取自產生本頁圖表的同一份統計輸出，不是另外打上去的。

同一個治療效果，把 932 列當成互相獨立算出來的標準誤是 0.069。隨機截距混合模型給的標準誤說明了什麼？

看答案與解析

正確答案: 0.111——比 naive 大了將近六成，重複測量帶來的資訊少於列數暗示的

混合模型給 0.111，naive 給 0.069。同一個人的四次測量彼此高度相關，所以它們合起來提供的資訊遠少於四個獨立觀察——naive 的算法把兩者當成一樣，於是標準誤被低估、信賴區間太窄、p 值太小，而點估計三種算法都一樣，從那一欄看不出任何異狀。0.138 是假如每人只測一次會得到的標準誤：混合模型並沒有退化到那個程度，重複測量確實有貢獻，只是不像列數看起來那麼多。

ICC 是 0.540，每人測四次，共 932 列。要說「這些資料實際上相當於多少個獨立觀察」，該用哪一個數字換算？

看答案與解析

正確答案: 用設計效應 2.62 去除列數

設計效應是 2.62，等於一加上（每人次數減一）乘以 ICC；把 932 除以它，有效樣本數只剩三百多。0.54 是 ICC 本身，它是相關的強度而不是折扣倍率；4.00 是每人的測量次數，直接拿它去除等於假設同一個人的四次測量完全相同，那是另一個極端。這也是為什麼「我有將近一千筆資料」在重複測量的研究裡是一句會誤導自己的話。

233 位病人裡，有 122 人四個時間點的分數完全沒有變化。這件事與 ICC 有什麼關係？

看答案與解析

正確答案: 122 人的組內變異是零，人與人之間的差異因此佔了總變異的大部分

ICC 是「人與人之間的變異佔總變異的比例」，而超過一半的病人四次測量完全一樣，等於他們的組內變異是零——分母裡少了這一塊，比例自然高。120 是四次都得零分的人，那是前一群的子集：分數持平在任何一個值都會讓組內變異歸零，不必是零分。233 是分析的病人數，樣本大小本身不會把 ICC 推高。

用到這個方法的章節

隨機對照試驗

素材來源與授權

本頁為原創內容

群集與重複測量資料

報表看不出來的那一條假設

這一頁的例子

動手跑一次

方向一：組間效果的標準誤太小

方向二：組內效果的標準誤太大

混合模型在做什麼

那個標準誤，其實就是「每人取平均再比較」

ICC：多少變異來自「人和人不同」

設計效應：ICC 怎麼換算成「浪費了多少樣本」

混合模型、GEE、cluster-robust：三個不同的問題

cluster-robust 不是「乘上一個固定倍數」

常見誤用

相關頁面

重跑本頁的所有數字

讀讀看這張圖

用到這個方法的章節

素材來源與授權

回報內容問題