邊際模型與廣義估計方程

12.1 一條迴歸式裡藏著兩個問題

想像一個戒菸方案，以及兩個表面上一模一樣的問題。第一個問的是：如果整個母體都戒菸一個月，復發的風險會降低多少；第二個問的是：如果某一位特定的吸菸者戒菸，他的風險會降低多少。第一個是母體平均 (population-average) 或邊際 (marginal) 的問題，比較的是母體在不同狀態下的平均結果；第二個是個體特定 (subject-specific) 或條件 (conditional) 的問題，比較的是一個人與他自己的反事實 (counterfactual)。對於用恆等連結 (identity link) 分析的連續結果變項，兩個問題的答案相同，因為取平均是線性運算，個體改變量的平均等於平均數的改變量。對於用羅吉斯連結 (logistic link) 分析的二元結果變項，兩者則不相同，而理由是幾何上的。圖 12.1 畫出一群個體羅吉斯曲線，每一條都很陡，彼此只差在隨機截距 (random intercept)，並把它們在母體上的平均一併畫出。平均曲線明顯比任何一條個體曲線都平緩，因為把一群水平位移的 S 形曲線平均起來，得到的是一條較淺的 S 形，而較平緩的曲線意味著較小的斜率。母體平均的羅吉斯係數因此相對於個體特定的係數被衰減 (attenuated) 了。

註：每一條細線是一個人的結果機率隨某個預測變項變化的羅吉斯曲線，各條之間只差在隨機截距 (random intercept)。它們在母體上的平均，也就是紅色的邊際曲線 (marginal curve)，比任何一條個體曲線都淺。母體平均斜率因而相對於個體特定斜率被衰減；對非線性的連結函數而言，兩種估計標的並不相同。

這個衰減是可以量化的。對一個隨機截距變異數為 \(\sigma_u^2\) 的羅吉斯模型而言，母體平均斜率大約等於個體特定斜率除以 \(\sqrt{1 + c^2\sigma_u^2}\)，其中 \(c = 16\sqrt{3}/(15\pi) \approx 0.588\)，來自以機率單位 (probit) 近似羅吉斯函數。圖 12.2 以模擬確認了這個關係：把兩種模型都配適到隨機截距變異數已知的資料上，母體平均斜率與個體特定斜率的比值緊貼著公式，並隨個體間變異數增大而穩定下降，由「所有人共用同一個截距」時的相等，一路降到個體間變異數很大時的大約一半。實務上的後果是，一個羅吉斯迴歸係數的意義取決於它出自哪一種模型，而兩者不可以當成同一個量來比較或合併。表 12.1 陳述了這項區分，而這一節真正要交付的，是它與研究問題之間的對應：方案評估、流行病學與政策方面的問題關心的是母體會怎麼樣，需要的是母體平均的估計標的；關於個體歷程與機制的問題關心的是一個人身上會怎麼樣，需要的則是個體特定的那一個。

註：母體平均與個體特定羅吉斯斜率的比值，由模擬得到（點），對隨機截距變異數作圖，並附上近似式 \(1/\sqrt{1 + c^2\sigma_u^2}\)（曲線）。個體間的變異數越大，邊際斜率被壓得比個體特定斜率越低。兩者只有在完全沒有個體間變異時才會相等。

12.2 廣義估計方程這套做法

一個廣義估計方程由三個彼此獨立的部件指定 (Liang & Zeger, 1986)。平均數模型 (mean model) 把平均結果的某個函數連結到預測變項上，與一般的廣義線性模型 (generalized linear model) 完全相同，所以對治療試驗的二元緩解 (remission) 結果而言，一個羅吉斯平均數模型把緩解的對數勝算 (log-odds) 表示成組別與週次的函數。變異數函數 (variance function) 說明結果變項的變異數如何隨它的平均數而變，同樣與廣義線性模型一致。而工作相關 \(R(\alpha)\) 描述重複測量之間的個體內依賴，這正是一般廣義線性模型所缺的那個部件。圖 12.3 顯示了幾個標準選項：獨立 (independence)，完全不理個體內相關；可交換 (exchangeable)，假定任意兩個時點之間都是同一個相關；自我迴歸 (autoregressive)，讓相關隨時間間隔衰退；以及無結構 (unstructured)，把每一對相關都自由估計。這套方法最值得注意也最具定義性的性質，在基礎概念方塊中推導，是平均數模型的係數估計值即使工作相關設錯了仍然是一致的，因為平均數的估計方程不論假定什麼依賴結構都是不偏的。工作相關只影響效率 (efficiency)，不影響推論的有效性，這也就是它被稱為「工作」相關的原因：它是一個用來提高精確度的裝置，而不是結論所倚賴的假設。這一點在治療試驗中得到實徵上的確認：不論假定獨立、可交換或自我迴歸相關，估計出來的組別與週次交互作用對緩解的效果在對數勝算量尺上都在 \(0.14\) 附近，只有標準誤略有差異。表 12.2 依設計提供選擇指引。

註：四種工作相關結構的熱圖 (heatmap)：獨立（沒有個體內相關）、可交換（單一個共同相關）、自我迴歸（隨時間間隔衰退）與無結構（每一對自由估計）。在其中任何一種之下係數估計值都是一致的；選擇只影響效率，而結構應該與設計相稱。

三明治估計式的穩健不是免費的，它的代價是一項小樣本脆弱性 (small-sample fragility)：三明治的夾餡是殘差乘積在各群集上的平均，而當群集數很少時，那個平均是真實變異數的一個很差、而且向下偏誤的估計，於是標準誤太小、信賴區間太窄。圖 12.4 以一項涵蓋率模擬顯示了後果：有一百個群集時，百分之九十五的區間涵蓋率幾乎達到名目水準，但只有十個或十五個群集時，涵蓋率只有百分之九十二，這樣的真實錯誤率膨脹會讓許多已發表的小樣本廣義估計方程分析變得過度寬鬆。解藥是偏誤修正過的三明治估計式，它們調整殘差以移除向下的偏誤；實務上的建議是，只要群集數低於大約四十到五十就套用一種修正 (Kauermann & Carroll, 2001; Mancl & DeRouen, 2001)。表 12.3 列出這些修正以及各自何時是必要的。

註：一項模擬中，某個時間斜率的名目百分之九十五信賴區間的實徵涵蓋率，對群集數作圖。群集數少時穩健標準誤向下偏誤，區間的涵蓋率因而低於名目水準；到大約五十個群集時，落差已可忽略。低於大約四十個群集就需要小樣本修正。

12.3 模型評估與選擇

因為廣義估計方程不是靠最大化概似 (likelihood) 配適出來的，那些熟悉的以概似為基礎的工具都不能用，最接近的替代品是 準概似訊息準則 (quasi-likelihood information criterion, QIC) (Pan, 2001)，它把 Akaike 訊息準則 (Akaike information criterion, AIC) 移植到準概似的情境，可以用來比較工作相關結構，其變形也可用來比較平均數模型。這個準則必須審慎使用，絕不可奉為圭臬，因為它不是真正的概似，行為也不如它所模仿的 Akaike 準則那樣為人所理解；它是挑選工作相關的一個粗略指引，不是實質模型的裁判。沒有概似同時也意味著沒有自然的決定係數 (coefficient of determination)，而誠實傳達一個廣義估計方程配適結果的方式，是透過模型所隱含的、在關心的量尺上的預測值，也就是母體平均的機率或平均數，而不是透過一個這套方法本來就給不出來的摘要統計量。

12.4 遺漏資料：這套方法的罩門

廣義估計方程帶著一項很容易被忽略的弱點，而且它與一個常見的誤解正好相反。因為這套方法不作完整的分配假設，人們常以為它處理遺漏資料很優雅，但事實正好相反：未加權的估計方程只有在資料為完全隨機遺漏 (missing completely at random) 時才有效，而這個要求遠比第 6 章那些概似方法所需的隨機遺漏 (missing at random) 條件嚴格得多。在縱貫研究典型的隨機遺漏退出之下，也就是狀況較差的人比較可能離開，未加權的方法是有偏誤的，因為它實際上是在留下來的那個被選擇過的樣本上取平均。圖 12.5 是本章的招牌展示，一項真值已知的模擬：症狀隨時間下降，而退出與否取決於前一次觀察到的值，於是病況較重的人先離開、留下來的人系統性地較健康。未加權的廣義估計方程跟著觀察到的存活者走，報告出來的平均軌跡遠低於真值，因而高估了改善；而以完全訊息最大概似法 (full-information maximum likelihood) 估計的混合模型幾乎精確地還原了真實軌跡，因為概似在隨機遺漏之下仍然有效。在邊際架構之內的解藥是逆機率加權 (inverse-probability-weighted) 的廣義估計方程 (Robins et al., 1995)，它把每一筆觀察到的紀錄依「仍然留在研究中的估計機率」的倒數加權，把容易退出的那類人加權放大，使觀察到的樣本得以代表完整的母體。圖 12.5 顯示這樣的加權修補了一部分偏誤，但也只有一部分，而理由很有啟發性：加權只用到邊際的觀察資料，丟掉了「一個已經離開的人，他的軌跡本來要往哪裡去」這種個體內訊息，而概似方法把這份訊息用滿了。這張展示因此一次交付兩個教訓：邊際方法在現實的遺漏機制之下是脆弱的，而即使是它的解藥，也仍然被本書其餘部分所偏好的概似取向所支配。

註：一項真值已知的模擬：症狀下降，而病況較重的人先退出（隨機遺漏，取決於前一次觀察到的值）。未加權的廣義估計方程（橘）落在真實世代平均數（虛線）之下，高估了改善；以最大概似估計的混合模型（藍）還原了真值；逆機率加權的估計方程（綠）修補了一部分偏誤但不是全部，因為它丟掉了概似方法所用的個體內訊息。

12.5 廣義估計方程的實作

主要的實作範例，是把治療試驗中母體平均的緩解機率配適成組別與週次的函數。模型報告的是各組的平均緩解率如何演變，而圖 12.6 把這些模型所隱含的母體平均機率畫了出來：兩組都改善，實驗組的緩解機率上升得比較快，而組別與週次的交互作用就是母體平均的治療效果。這正是一項方案評估想要的量，也就是母體比率的改變；它以機率而不是以對數勝算係數報告，恰恰是因為母體平均機率是可解讀的，而被衰減的係數則會誘發上面所警告的那種個體特定式的誤讀。第二個例子，也就是把經驗取樣資料中的日層次二元結果對個體內預測變項建模，可以用來說明這套取向誠實的界線：當關心的真的是個體內的歷程，例如一個人在比自己平常更有壓力的日子裡是不是更可能有某種行為，母體平均的估計標的回答的是一個微妙地不同的問題，而第 15 章與第 23 章的個體特定混合模型才是比較好的工具。表 12.4 是報告檢核表。

註：由一個羅吉斯廣義估計方程得到的、治療試驗中兩組跨週次的母體平均緩解機率。兩組都改善，而實驗組改善得比較快；組別與週次的交互作用就是母體平均的治療效果，也就是一項方案評估要問的那個量。

12.6 在邊際模型與混合模型之間選擇

在廣義估計方程與混合模型之間的選擇，不是統計品味的問題，而是估計標的與假設的問題，圖 12.7 把這個決策攤開來。第一個、也是主導性的問題是研究想要哪一種估計標的：要母體平均的答案就指向邊際模型，要個體特定的答案就指向混合模型，而對線性結果變項而言兩者一致，用哪一個都行。第二個問題是遺漏機制：常見的隨機遺漏，以概似為基礎的混合模型可以有效處理，邊際模型則需要加權。第三個問題是分析必須交付什麼：只有混合模型能給出隨機效果變異數 (random-effect variance) 的估計，也就是改變上的個別差異，而那往往正是實質上關心的量；邊際模型沒有它，既是一種簡化，也是一項真實的損失。其餘的考量，也就是群集數與計算上的穩健性，在群集很少、混合模型的隨機效果結構難以估計時偏向邊際模型。表 12.5 摘要了這項比較。兩種方法是協調的，不是對立的：第 15 章的廣義線性混合模型就是這裡所發展的邊際模型的個體特定對應項，而第 12.1 節的衰減關係正是它們的係數之間的那座橋。

註：估計標的主導這項選擇：母體平均的問題指向廣義估計方程，個體特定的問題指向廣義線性混合模型。超出完全隨機遺漏的遺漏機制需要加權或改用概似方法；需要隨機效果變異數則必須用混合模型。對線性結果變項而言兩種估計標的一致。

12.7 在 R 中執行廣義估計方程

geepack 套件可以配適廣義估計方程，群集識別碼由 id 引數傳入，工作相關由 corstr 指定，並預設報告穩健標準誤。母體平均的緩解模型，以及跨工作相關的比較，只要幾行。

library(geepack); library(dplyr)
rct <- readRDS("Examples/data/therapy_rct.rds") |>
  mutate(arm = factor(arm, levels = c("Control","Treatment")),
         remit = as.integer(hdrs <= 7)) |>
  filter(!is.na(hdrs)) |> arrange(patient_id, week)

# --- 母體平均的羅吉斯模型；穩健標準誤是預設 ---
fit <- geeglm(remit ~ week * arm, id = patient_id, data = rct,
              family = binomial, corstr = "exchangeable")
summary(fit)                                   # week:armTreatment 就是母體平均效果

# --- 係數在各種工作相關之下都一致 ---
sapply(c("independence","exchangeable","ar1"), function(cs)
  coef(geeglm(remit ~ week * arm, id = patient_id, data = rct,
              family = binomial, corstr = cs))["week:armTreatment"])

邊際圖所用的母體平均機率，來自對已配適模型呼叫 predict，而逆機率加權的分析則是把權重傳給同一個配適函式。完整的程式碼，包括衰減與涵蓋率兩項模擬，以及帶累積退出權重的隨機遺漏脆弱性展示，都在隨附的腳本 ch12_analysis_V01.R 之中。圖形則由中文版的 ch12_figures_zh_V01.R 繪出。

# --- 依組別與週次的母體平均預測機率 ---
newd <- expand.grid(week = 0:11, arm = factor(c("Control","Treatment")))
newd$p <- predict(fit, newdata = newd, type = "response")

# --- 逆機率加權：把 MCAR 放寬到 MAR ---
# 權重 w = 1 /（仍留在研究中的累積機率），由一個退出模型算出；
# 傳給 geeglm，並使用獨立的工作相關
gee_ipw <- geeglm(y ~ factor(t), id = id, data = observed, weights = w,
                  corstr = "independence")

12.8 解讀與報告結果

一項廣義估計方程的結果，報告時要點名估計標的、描述三個模型部件，並在關心的量尺上解讀效果。這項治療試驗的示範段落可以這樣寫：「母體平均的緩解機率以一個羅吉斯廣義估計方程模型化，採可交換的工作相關與穩健標準誤，並以病人作為群集（\(N = 240\)）。組別與週次的交互作用，也就是母體平均的治療效果，在對數勝算量尺上為正（\(\hat\beta = 0.13\)，穩健 \(SE = 0.09\)），而模型所隱含的母體平均緩解機率在實驗組上升較快，到第 12 週時達到較高的比率（圖 12.6）。在獨立、可交換與自我迴歸三種工作相關之下估計值幾乎沒有變動，這與係數的一致性相符。由於退出比較可能是隨機遺漏而不是完全隨機遺漏，主要的推論倚賴第 13 章的混合模型，廣義估計方程則作為母體平均的互補結果報告。」估計標的被點名，穩健標準誤與工作相關都寫出來了，而遺漏的但書也很誠實。

12.9 常見的迷思

關於廣義估計方程，有幾種說法會誤導讀者。第一種是廣義估計方程不過是用另一套軟體配適的混合模型；兩者回答的是不同的問題，邊際的與條件的，而對非線性連結而言那是不同的量，兩者也在不同的遺漏假設之下才有效。第二種是穩健標準誤可以救所有的問題；三明治估計式保護的是工作相關被設錯 (misspecified)，而不是平均數模型被設錯，而且它自己在群集少時也很脆弱。第三種是可交換相關對縱貫資料是安全的預設；縱貫的相關會隨時間衰退 (decay)，所以自我迴歸結構通常更說得通，儘管這個選擇只影響效率。第四種、也是後果最嚴重的一種，是因為這套方法不作分配假設，所以它處理遺漏資料很好；事實正好相反，未加權的方法需要完全隨機遺漏這個很強的條件，而在概似方法可以妥善處理的隨機遺漏退出之下它是有偏誤的。還有一個反覆被問到的問題，是審稿人為什麼會要求用某一種方法，答案要回到估計標的：把研究關心的那個量點名出來，方法自然就跟著定了。

接下來讀哪裡

邊際模型是走出古典方法的兩條路之一，另一條是混合模型，由第四部整個展開。第 13 章引進線性混合模型，它的隨機截距與隨機斜率 (random slope) 把邊際模型略去的個別差異還原回來，而它的概似估計在隨機遺漏之下有效；這裡學到的穩健標準誤會在那裡以群集推論的選項再次出現。第 15 章完成本章開啟的那一對，把廣義線性混合模型發展成邊際模型的個體特定對應項，屆時第 12.1 節的衰減關係會成為兩組係數之間明確的翻譯。逆機率加權的邏輯會在第 29 章的存活模型中回來，而「先點名估計標的」這項紀律，則會在任何需要把因果或描述目標講精確的地方反覆出現，其中最尖銳的是第 21 章的交叉延宕追蹤模型爭論。要帶走的教訓是：問題選定估計標的，估計標的選定模型，而再多的穩健性都替代不了把問題問對。

習題

1. 12.1 驗證衰減。模擬隨機截距變異數已知的個體特定羅吉斯資料，同時配適一個廣義線性混合模型與一個廣義估計方程，確認兩者斜率的比值符合衰減公式。
2. 12.2 工作相關。在獨立、可交換與自我迴歸三種工作相關之下重配緩解的例子，報告係數與它的穩健標準誤如何變化，並解釋係數為什麼那麼穩定。
3. 12.3 把脆弱性造出來。在一個模擬資料集上建構一個隨機遺漏的退出機制，對照混合模型顯示未加權廣義估計方程的偏誤，再用逆機率加權把它修補回來。
4. 12.4 小群集的涵蓋率。模擬一項有十五個群集的研究，估計樸素與偏誤修正三明治標準誤的實徵涵蓋率。
5. 12.5 估計標的短論。針對三個研究情境，判斷要的是母體平均還是個體特定的估計標的，並以一小段文字為每一個選擇辯護。