實驗介入成效之各種統計方法比較（下篇，GEE與LMM）

文章作者：林星帆
本文為重新編輯 2013/11起 發表在晨晰統計部落格的文章，這裡、這裡與這裡



圖片來源：在此

處理重複測量資料的新一代統計方法

如本系列的上篇文章所述，關於實驗（介入）效果的檢驗，目前我們已經學習到了 t 檢定（獨立 t 與配對 t）、差異中的差異分析（DID）、共變數分析（ANCOVA）以及多因子變異數分析（Multi-factorial ANOVA）。

基本上這些方法都是至少在 1950 年代以前就發展完畢的工具，如果用 Google scholar 搜尋「Analysis of covariance」或「Analysis of variance」，發現從 1930 年代就有文獻在報導，而在 1930 到 1980 年代其間，其實上述這些方法在應用於實驗介入效果上並沒有太多的改變。

而一直到 1986 年，由現時陽明大學校長梁賡義教授（當時任教於約翰霍普金斯大學生物統計系）發展出廣義估計方程式（Generalized estimating equation, GEE）應用於縱貫性研究（Longitudinal data analysis），並發表於生物統計重量級期刊 Biometrika，此後縱貫資料分析就常常見到 GEE 分析的身影。

而之後我們即將要介紹的線性混合模式（Linear mixed model, LMM）則是稍微早一點應用於縱貫資料，為哈佛大學生物統計系 Nan Laird 與 James Ware 教授首先於 1982 年提出以隨機效果模式（Random effect model）處理縱貫資料，並提出一個更廣泛的廣義線性模式（Generalized linear model）。

因此，從 1990 年代之後，GEE 與隨機效果模式就成為縱貫資料分析的最主流方式，而實驗型（介入）研究由於至少會有一次後測資料，加上前測則至少有2次時間點的資料，可算是縱貫型資料的最簡化形式，既然屬於縱貫資料那麼當然可以用 GEE 與隨機效果模式這兩個方法來驗證介入效果。

雖然這兩種方法在許多面向都比過去發展的方法更為彈性且需要較少的假設，但是這兩種方法在初始概念上卻差異極大，主要是在看待「重複測量」這件事情上，它們採取完全不一樣的角度（與假設），接下來我們就開始介紹這兩種方法的原理與應用。


廣義估計方程式（Generalized estimating equation, GEE）

筆者曾分別在 2009 及 2012 年寫過兩篇 GEE 的介紹文章（文章連結在這裡與這裡），但皆是以敘述其概念為主，這一篇文章我即將以簡單的（但不能說是精確）迴歸方程式來表示 GEE，讓各位讀者一窺其奧妙。

為了簡化，我們仍然以 2×2 的設計來討論，也就是組別為實驗組與對照組，測驗的時間點只有前測與後測，此時 GEE 若以迴歸方程式來表示，如以下所示：

其中

y_i1 = 某人的前測分數 pretest score;

y_i2 = 某人的後測分數 posttest score;

x₁ = 1 代表實驗組、x₁ = 0 代表對照組;

x₂ = 1 代表後測、x₂ = 0 代表前測;

CORR 為工作相關矩陣;

下標 i 代表某一位研究對象

讀者可立刻發現這個公式跟多因子變異數分析的公式非常雷同，但是差在多了「CORR」這個項目（自變項），這部分我們後面會再詳細介紹，讓我們先把焦點放在方程式。

與多因子變異數分析相同的是，交互作用（x1x2）的存在就是假設實驗組與對照組的斜率不同，若 β₃ 達顯著差異就表示兩組的斜率確實不同，這表示實驗組的前後測差值與對照組的前後測差值之間有顯著不同，因此表示介入確實有效（如果是實驗組比對照組多 10 分，表示實驗組的進步幅度比對照組還多 10 分，而這 10 分就是實驗介入的淨效果），在這個部分的解釋方式其實與多因子變異數分析是相同的。

GEE 的工作相關矩陣

接著開始討論方程式多出的「CORR」此項目，稱之為工作相關矩陣（Working correlation matrix），這是 GEE 兩大特色之一，而這個工作相關矩陣就是 GEE 如何「看待」重複測量的精神所在。

簡單地來說，工作相關矩陣允許同一受測者的不同時間點之間的依變項是具有相關的（正相關），例如前測分數越高者其後測分數通常也會高（反之亦然）。

不僅如此，GEE 的工作相關矩陣還可以選擇不同的相關類型，我們就先簡介幾種最常見的，包括獨立（independent）、未結構化（unstructured）、可交換（Exchangeable）以及 AR1（Auto-regressive first order），下圖。

上圖列出幾種最常見的工作相關矩陣類型，我們假設每個樣本都有五個時間點的資料（假設前測 1 次，然後連續收集 4 次後測）。接著，我們依序討論各種矩陣的意涵：

獨立矩陣（Independent）

我們可以看到非對角線（off-diagonal）的相關係數全部都被規定是 0，這表示同一個受試者在不同時間點的依變項得分是完全沒有關係的，很明顯的如果是重複測量（縱貫型）的研究，獨立矩陣是絕對不合理的。

其實，除非樣本數非常的小，否則獨立矩陣皆不適於在任何情況之下採用（然而，如果是使用 SPSS 所發表的論文，常常看到用獨立矩陣，這是因為 SPSS 的預設是獨立矩陣，請各位讀者一定要注意）

AR(1) 矩陣（Auto-regressive first order）

我們可以看到 t1 跟 t2 的相關是「ρ」，這個相關係數會自動由實際的觀察資料中計算出來（透過 GEE 的平均數與共變異數矩陣所計算的），不過可發現 t1 跟 t3 的相關是「ρ²」，也就是說，如果 ρ 等於 0.70，那麼 ρ² 就等於0.49。

如果我們的 t1 指的是前測，t2 是介入後 1 個月，t3 是介入後 2 個月，此時使用 AR(1) 是非常合理的選擇，因為如果前測跟剛介入完 1 個月的相關是 0.7，那麼可以預期前測與介入後 2 個月的相關應該會比較低，因為 t1 跟 t2 只距離了一個月，而 t1 跟 t3 卻距離了兩個月，AR(1) 矩陣假設距離越久的時間點之間的相關越低，而且此相關係數會等於 ρk（k 為距離幾個時間點）。

通常 AR(1) 適用於重複測量的間隔（interval）是相同長度的研究，例如不同時間點之間的時間間隔是相同的，例如 t1 跟 t2 離一個月，t2 跟 t3 也是離一個月，以此類推。

反之，如果 t3 是介入後 6 個月，那麼此時選擇 AR(1) 或許就不是這麼適合，t1 跟 t2 的相關是 ρ（例如 0.70），雖然 t1 跟 t3 的相關是 ρ²（0.49）非常合理，但是此時 t2 跟 t3 由於只距離一個時間點因此相關係數也是 ρ（0.70）就變得不合理了，因為 t2 跟 t3 之間距離了五個月之久。

可交換矩陣（Exchangeable）

可交換矩陣比較容易理解，因為大家可以觀察到五個時間點（一共 10 個相關係數）全部都被設定一樣（全部都是 ρ），也就是說我們假設不同時間點之間的相關係數是相同的。

通常在縱貫型分析，這也是極為常用的矩陣之一，特別是重複測量的間隔不是相同長度時。

然而在非縱貫型研究的資料中，一般我們稱之為集群資料（Clustered data），例如同一個學校之內的學生（的依變項）比較容易有相關、同一個主治醫師有偏好的治療方針也導致底下病人的預後比較容易有相關（另外一種實際狀況，病情比較嚴重的病人會被送到「名醫」那邊，所以該名醫底下病人的預後狀況有相關），此時可交換是最適合的工作相關矩陣。

值得注意的是，在多因子變異數分析中（例如重複測量變異數分析、混合設計變異數分析），所使用的就是可交換矩陣，不過在 ANOVA 中稱之為複合對稱（Compound symmetry）。

未結構化矩陣（Unstructured）

未結構化工作相關矩陣不假設各時間點之間的相關係數為多少，而是以實際觀察資料作估計，因此十個相關係數可能都不一樣，聽起來好像最「準確」，不是嗎？

確實如此，但是必須瞭解到一個事實，在 GEE 的方程式中，除了迴歸係數需要估計之外（β₀, β₁, β₂, β₃），工作相關矩陣也需要被「估計」，如果我們採用獨立矩陣，那麼不需要再額外估計（因為是 0），若是可交換或 AR(1) 都只需要額外估計一個相關係數（就是 ρ），但是在未結構化矩陣中，竟有「10 個」相關矩陣需要被估計，也就是說，雖然從方程式中看起來只有 4 個參數，但是其實是 14 個參數（加上 10 個需要估計的相關係數），這在統計術語來說是很沒有效率（Efficient）的模型。

因此除非我們樣本數非常的大（例如上千、上萬），否則一般我們是不考慮用未結構化的矩陣的。

還有一種狀況，即使樣本數非常大，但時間點如果很多個，那麼也不會考慮使用未結構化，例如 10 個時間點就有 45 個相關係數要被估計（C10 取 2），等於方程式中竟然多達 45 個參數需要被估計，這會使得整個模型其他變項很可能都不顯著。

目前為止，我們已經可以理解到 GEE 是如何「看待及處理」同一個個案的重複測量資料，就是透過工作相關矩陣正視「不同時間點之間的依變項有正相關」這個事實，並將此工作相關矩陣加到方程式之中加以估計。接著，我們要介紹角度完全不同的另一個主流方法，稱之為混合效果模式（LMM）或隨機效果模式（Random effect model）。

如果對於 GEE 計算過程有興趣，但卻又不嫻熟於代數或矩陣的讀者，可以參考 Hanley et al.（2003）的文章，裡頭非常詳盡地介紹了 GEE 的計算過程與其代表的意義，而且可以透過手算完成（放心，只需要用到加減乘除跟一點點的矩陣），非常值得各位花時間閱讀。

混合線性模式（Linear mixed model, LMM）

目前為止，對於檢驗實驗介入成效的方法，我們已經知道非常多種方式，而在這其中目前又以 GEE 及現在即將要介紹的 LMM 蔚為主流。

不過很有趣的是，這兩種方法皆克服了傳統統計方法（例如 t-test, ANCOVA, ANOVA）的某些限制，因此才會廣受研究者的歡迎，但這兩種當代主流方法卻在「看待及處理」同一個個案的重複測量資料時，採用完全不同的角度，以下就讓我們開始介紹 LMM。

LMM 全名為混合線性模式（Linear mixed model）或稱為 Mixed effect model，不過它在各個領域有不同的名稱，在生物統計領域習慣稱作LMM，在應用統計領域則常被稱為多層次模型或多層次迴歸（Multilevel model / Multilevel regression），但在教育或心理領域則常以階層線性模型（Hierarchical linear modeling, HLM）稱呼它，而在經濟或財金領域最可能稱之隨機效果模式（Random effect model），不過無論如何稱呼，其背後的原理大致是差不多的。

筆者曾在 2009 年寫過一篇簡介文章，有興趣的讀者可以參考。

LMM 的基本想法

為了簡化，我們仍然以 2×2 的設計來討論，也就是組別為實驗組與對照組，測驗的時間點只有前測與後測，此時 LMM 若以迴歸方程式來表示，如以下方程式所示：

Level 1（第一層方程式）

Y_ij = β_0i + β_1ix₂ + r_ij

Level 2（第二層方程式）

β_0i = γ₀₀ + γ₀₁ x₁ + μ_0i

β_1i = γ₁₀ + γ₁₁ x₁ + μ_1i

其中

y_ij = 第 i 個個案的第j次的依變項得分，j = 1 or 2

x₂ = 1 代表後測、x₂ = 0 代表前測;

x₁ = 1 代表實驗組、x₁ = 0 代表對照組;
下標 i 代表某一位研究對象

各位讀者可以從很多的 HLM/LMM 教科書發現，在說明 LMM 的方程式時，常用這種「多層次」的方式來教學，其實並不是一定要用這樣來呈現，但是若不拆成這種格式呈現，大多數非統計背景的讀者是不容易理解的。

讓我們先看一下下圖，LMM 在看待同一位個案的多次資料點，是以「鑲套或鑲嵌」（nested）的方式來看待，此時重複測量資料點是巢套在個人之下，如下圖。

讓我們再回到第一層方程式，「β₁」代表的是時間的效果（因為 x₂ 是時間），可以說時間每增加一個單位則 Y 增加多少個單位，不過由於我們現在的例子是只有 2 個時間點，因此就是後測（1）減前測（0）的差值。

但是眼尖的讀者是否有看到迴歸係數為「β_1i」而非「β₁」，很重要的是，這個「i」事關重大，可以說是 LMM 的核心概念，「i」代表的是第 k 位個案，若 k 等於 100 則 i 就介於 1 到 100，β_1i 代表的是「每一位個案都有自己的一條線」，統計術語是每一位個案的時間趨勢（線性效果）被允許可以不同。

但讀者可能會納悶，每一個個案都有自己的成長趨勢，這不是很合理嗎？是的，確實是非常合理，但問題是傳統統計方法只能假設/強迫每一位個案的成長趨勢是相同的，而 LMM 就是可以突破這個假設，此為 LMM 受歡迎的主要原因之一。

然而除了時間趨勢（線性效果）可以隨個案而改變之外，大家可以注意到連截距項（intercept）也有下標 i，因此「β_0i」代表的是當 x₂ 為 0 的時候（注意我們的編碼方式，x₂ = 0 代表是前測），此時 Y 等於多少，也就是說「β_0i」就是第 i 位個案的前測分數，而下標 i 表示 LMM 可以允許每一位個案的前測分數是不一樣的，相同的傳統統計方法只能假設/強迫每一位個案的截距項是相同的。

LMM 受到歡迎的原因

所以我們可以得知，LMM 廣受歡迎的最主要原因就是可以設定「隨機效果」（random effect），例如允許每一位個案的初始值（在我們這個例子中，就是前測分數）可以不同。

換句話說，每個人的初始值具有變異性，我們可以從前測分數（β_0i）的第二層迴歸係數中再拆解為3個項目，其中第3項的「μ_0i」就是很關鍵的隨機效果，有這一項就代表我們允許前測分數（β_0i）具有變異數，到時候統計報表會提供這個變異數的顯著性，此變異數的統計符號以「τ₀₀」表示（τ 念做 tau）。

然而此時的「γ₀₁」代表的是實驗組與對照組在前測的分數是否有顯著差異，也就是兩組前測分數的同質性檢驗，值得注意的是，雖然 LMM 是以兩層方程式來表示，但是在之前介紹的多因子變異數分析或 GEE 的方程式中也都有這個項目，只是 LMM 看起來會有一點不太一樣。

另外一個不是很具有實務意義的就是「γ₀₀」，它代表的是對照組的前測分數是否不同於 0，一般不去解釋這個項目。

接著讓我們再關注第二層迴歸方程式，「β_1i」代表是每一個人的成長曲線，相同的道理，「μ_1i」的存在就是允許每一個人都有自己的線，也因此斜率也會有一個變異數，代號為「τ₁₁」。

此時的「γ₁₁」代表的是實驗組與對照組在成長幅度是否有顯著差異，這個項目非常重要，它就是我們用來驗證介入效果是否存在的效應項，也就是實驗組的前後測差值減掉對照組的前後測差值。

承上，其實我們在 GEE 模式中也是用相同的效應項來驗證介入成效，只不過 GEE 並沒有允許每一個個案可以有自己成長的趨勢，換句話說就是 GEE 的時間效果是設為固定效果（Fixed effect），因此 GEE 算是邊際模型（Marginal model）的一類。

附帶一提，此時的「γ₁₀」具有實務意義，因為它代表對照組的線性效果是否為零，白話一點來說，就是對照組的前後測差值是否顯著不同於 0，也就是對照組本身的配對 t 檢定。實務上有時候我們會將實驗組的編碼設為 0，對照組設為 1，此時的「γ₁₀」就變成實驗組本身的前後測分數是否具有差異，具有報導的價值。

Level 1（第一層方程式）

Y_ij = β_0i + β_1ix₂ + r_ij

Level 2（第二層方程式）

β_0i = γ₀₀ + γ₀₁ x₁ + μ_0i

β_1i = γ₁₀ + γ₁₁ x₁ + μ_1i

其中

y_ij = 第 i 個個案的第 j 次的依變項得分，j = 1 or 2

x₂ = 1 代表後測、x₂ = 0 代表前測;

x₁ = 1 代表實驗組、x₁ = 0 代表對照組;

下標 i 代表某一位研究對象

總結

結合如本系列的上篇文章，我們可以這麼總結，如果我們的研究是隨機分組，然後又只有一個後測的時間點，那麼應該可以放心地使用共變數分析（ANCOVA），因為它在隨機分派時具有最強的檢定力（power）。

當不是隨機分組的時候，那我們最好看兩組的進步幅度是否具有差異，再來說實驗介入是否有效，然而之前提到 DID 分析會有「向平均數迴歸」的問題，所以可以建議採用 GEE 或 LMM（而且可以控制其他混淆因子的效果）。

根據許多模擬研究指出，如果依變項是連續的，那麼 GEE 跟 LMM 其實結論幾乎都是相同的，不過如果樣本數很小（例如小於 30），LMM 很有可能模式無法收斂，因為隨機效果的設定是很耗費估計參數的作法，所以如果非隨機分派又是小樣本的話，可以優先考慮 GEE，大樣本的話兩者皆可考慮。

無庸置疑的是，在當代的統計方法中，目前已經是以 GEE 及 LMM 為主要的方法，特別是時間點 3 次以上的實驗設計，也已有許多方法論及統計學的專家在研究這兩種方法的表現，之後有機會我也會再分享這些學者的研究精華。

參考文獻

1. Hanley, J. A., Negassa, A., & Forrester, J. E. (2003). Statistical analysis of correlated data using generalized estimating equations: An orientation. American journal of epidemiology, 157(4), 364-375. 
2. Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963-974. 
3. Liang, K.-Y., & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13-22. 

更多文章  或  獲得協助