盡管 Meta 分析技術得到了快速發展，但目前還缺乏針對縱向數據的合并技術。基于廣義線性混合效應模型的 Meta 分析模型可以充分考慮各個時間點之間的相關性，并準確地估計出最終的合并效應值，是縱向數據進行 Meta 分析的理想模型。本文通過實例數據，介紹如何利用 SAS 軟件實現縱向數據 Meta 分析，并提供 SAS 編程代碼。
引用本文： 陳櫻, 黃碧芬, 鄭建清, 肖麗華, 吳敏. 基于 SAS 軟件的混合效應模型實現重復測量數據的 Meta 分析. 中國循證醫學雜志, 2019, 19(8): 998-1006. doi: 10.7507/1672-2531.201902004 復制
縱向數據是一種特殊形式的數據，主要來源于醫學和社會學等領域。在醫學研究領域，縱向研究通常用于臨床和流行病學研究，縱向數據通常也被稱為重復測量數據。重復測量數據主要來自于同一受試對象在不同時間點上對同一個觀測指標進行多次測量獲得的數據[1]。被納入 Meta 分析的縱向研究通常具有以下特點：① 不同時間點上的測量對象是同一個研究內的研究對象集合，而非單一個體，因此研究對象存在失訪、脫落等現象；② 不同研究針對的觀測指標，即結局指標是相同的，但受到失訪、脫落等顯著影響；③ 不同研究的時間間隔不同，而且并非報告了全部時間點上的結局指標。對縱向研究數據進行 Meta 分析，系統評價員不僅關心主要干預措施對研究對象的整體效應，也會關心在某些特殊時間點的測量指標效應量、主要測量效應的時間趨勢以及不同時間點之間的變化值等[2]。
盡管 Meta 分析技術得到了快速發展，但目前針對重復測量數據的最佳合并技術還在探索中。縱向數據 Meta 分析的難點是解決不同時間點測量結局指標的內在相關性，這意味著此類數據的合并分析不僅僅需要考慮研究內變異和研究間變異，還需要考慮結局變量的出現時間點以及不同時間點的變化等[3]。Ishak 等[4]提出一種基于廣義線性混合效應模型來合并縱向數據的 SAS 程序，其不僅考慮了固定效應，而且也考慮了隨機效應，同時也可以獲得不同時間點的結局變量的效應值及其變化趨勢。此外，亦可以在模型中引入協變量，分析協變量對結局測量的影響等。該模型是進行縱向數據 Meta 分析的有力工具。本文基于實例，詳細介紹了廣義線性混合效應模型 Meta 分析技術在 SAS 軟件中的實現方法，并提供相關的 SAS 代碼。
1 基于廣義線性混合效應模型 Meta 分析的原理及方法簡介
1.1 縱向數據的一般線性混合模型
在介紹基于廣義線性混合效應模型 Meta 分析之前，先簡單介紹一般線性混合效應模型的基本原理，該模型是縱向數據的重要分析方法。
假設對納入分析的 N 個單位（units，此處單位是指觀測結果來源于單個個體；而 Meta 分析所指個體是 study，即單個研究）進行 K 個時間點上的結局測量；可以用 yi 表示來自第 i 個單位的觀測值，也就是 K×1 個結局向量，并用 yij 表示該單位的第 j 個觀測值。根據上述變量，可以擬合出一個可以解釋觀測值之間相關性的一般線性混合模型：
 |
在上述模型中，Xiβ 是模型的固定效應部分，其中 Xi 是包含潛在時間相關的解釋變量的 K×p 的矩陣；而 β 是一個包含固定效應的 p×1 向量，是待估計的固定效應參數向量，β 反映了解釋變量 Xi 對響應變量 yi 的影響程度。Ziγi 是模型的隨機效應部分。Zi 是某些協變量的 K×q 的設計矩陣，其構成方式與 Xi 相同，但更通常的情況是 Xi 中的一部分；而 γi 是一個包含隨機效應的 q×1 向量，是待估計的隨機效應參數向量，服從均值向量為 0、方差協方差矩陣為 G 的正態分布，表示為 γi～N（0，G）。εi 是 K×1 的殘差的設計矩陣，是隨機誤差向量。在混合效應模型中，放寬了對 εi 的限制條件，其元素不必為獨立分布，即對 εi 沒有 Var（εi）=σe2 及 cov（εi，εj）=0，cov（Ei，Ej）=0 的假定。用符號表示隨機誤差向量 εi～N（0，R）。不要求 εi 的方差、協方差矩陣 R 的主對角元素為 σe2、非主對角元素為 0。該模型假設不同單位（研究）的觀察結果是獨立的，因此 cov（εij，εml）=0，其中 i≠m，在任意一個觀測結局中，j，l=1···K。模型的另外一個假設是殘差和隨機效應是獨立的，即此 cov（γi，εi）=Cov（G，R）=0，即 G 與 R 間無相關關系。
模型中，γi 的期望值 E（γi）=0，εi 的期望值 E（εi）=0，或者 E
=0。
 |
 |
因此，響應變量 yi 的邊際分布是一個多元正態分布（multivariate normal distribution，MVN）；可以被指定為 MVNq（Xiβ，Vi），其中 Vi=var（yi）=Zi'GZi+Ri。
G 和 Ri 的協方差結構決定了如何處理個體間和個體內（研究間和研究內）的相關性。這種可以通過以下幾種方法解決。例如，可以選擇不包括任何隨機效應，并將 Ri 設置為對于所有研究都恒定的一般非結構化正定矩陣：即 Ri=R。另一種方法是引入隨機截距（也就是將 Zi 設置為 1S 的列）具有（標量）方差 G，并且設置 Ri=σ2IK，其中 IK 是 K×K 的單位矩陣。這樣就形成了一個復合對稱相關結構，其中當 j≠l 時，corr（yij，yil）=G/G+σ2。當涉及多個隨機效應時，必須給 G 指定一個結構來描述每個研究內隨機效應之間的關系。一些最常用的結構是：復合對稱（compound-symmetry，CS）或恒定相關（constant correlations，CC）、一階自回歸結構（first-order auto-regressive，AR（1））。
1.2 縱向數據 Meta 分析的混合模型
上面描述的一般模型稍微修改即可適用于縱向數據的 Meta 分析。實際上，在實踐中常用的單變量隨機效應 Meta 分析模型就是一個指定 K=1 和 Zi=1（即隨機截距）的特例。研究內方差（是單變量情況下的標量）被設置為每個研究中報告的觀察到的估計值的方差，并且假設已知且無誤差。通過該方差，可以通過倒方差來加權每個觀察的效應。
在縱向數據的 Meta 分析中，觀察的單位是在 K 個時點上收集的治療效應估計值的研究。這些時間的選擇通常是基于可用的數據和分析目標。研究之間可能會因預設的測量的時間而異，并且可能無法報告某些測量時間點的結果。如果報告的結果與效應的大小明顯相關，則可能是一種發表偏倚，類似于研究中可能發生的選擇偏倚。
Meta 分析數據的縱向模型也不同于協方差矩陣規范中的一般情況；實際上，上述一些方法在這種情況下是不合適的。例如，隨機效應模型越來越被認為更適合 Meta 分析。因此，省略隨機效應并僅考慮 Ri 的相關性（即采用固定效應方法）不是理想的。此外，研究之間可能會因樣本量、測量方法和其他方面而有所不同，這些方面將導致不同研究之間的估計值存在或多或少精確性差異。例如，由于失訪，研究中的差異也可能隨著時間的推移而變化。因此，將研究內協方差矩陣 Ri 設置為 σ2IK 是不夠的，因為它假設所有研究中的所有效應估計的方差是不變的。至少應該將 Ri 設置為是具有不同值的對角矩陣。并且，為了確保對數據進行適當的加權，我們將 Ri 的對角元素設置為研究中報告的值，正如單變量 Meta 分析。
模型的擬合始于設置固定效應，這將指導隨機效應的選擇。由于 Meta 分析的目的是總結治療隨時間的變化，因此 Xi 包括將響應變量與時間相關聯的協變量。這可能包括一些時間參數化（例如：線性、二次、對數等）。在這里，我們考慮在每個 K 測量時刻總結效應的情況。這可以通過省略截距并包括 K 個時間指標，或通過包括一個截距和 K?1 個時間指標來完成。例如，當 K=4 時，我們可以設置：
 |
在第一種情況下，匯總估計值是在每個時間點獲得的。而在第二種情況下，每個時間點的匯總效應估計值的協變量測量差異與第一種情況相關，其可以通過截距估計。研究水平的協變量（例如：男性的比例、平均年齡）也可以包括在 Xi 中。這些協變量通常具有固定效應，因此不會影響協方差的相關結構。
在下文中，除非另有說明，否則我們假設 Ri 是對角矩陣，其值設置為研究中報告的估計的方差（我們用研究 i 的 Rij2 表示第 j 個對角線元素）。為了便于記憶，我們假設 Xi 只包括時間指標（即沒有截距或其他協變量）；將其他固定效應納入到模型中不會改變相關結構。
1.2.1 隨機研究效應模型
考慮觀測之間相關性的最簡單方法是隨機效應模型，它是給定研究的所有觀測的共性結果。這可以被認為是隨機截距模型。例如 K=4，我們會設置 
，這樣 γi 是一個標量；然后可以將模型寫成：
 |
其中 Var（γi）=G。對所有的 i 和 j 來說，殘差向量的協方差是 Cov（εi）=Ri 和 Cov（γi，εij）=0。yij 的邊際分布的方差是 G+Rij2。在研究 i 中的 j 和 l 時刻收集的兩個觀測值之間的協方差由 Cov（yij，yil）=G 給出；不同研究的觀察結果是獨立的。因此，觀測之間的相關性是 corr（yij，yil）=G/
，在這種情況下，觀測之間的不同相關性是被允許的，因為研究內方差在每個時點都可以不同。
1.2.2 隨機時間效應模型
隨機研究效應模型有一定的限制性，因為它假設研究之間的異質性以同樣的方式影響給定研究中每個測量時刻的觀察結果。然而，這可能不是真的。例如，在每項研究中，病例丟失可能會有所不同，導致后期估計的影響具有更大的異質性。 我們可以通過在每個測量時刻引入隨機效應來擴展模型。也就是說，我們設置 Zi=Xi；然后通過以下方式給出模型：
 |
或者
 |
因此，γi 現在是具有協方差矩陣 G 的隨機效應 K? 向量。假設這些是獨立的（即，D 是對角線）等同于假設觀測本身是獨立的，因為：
 |
因為 G 和 Ri 都是對角線，殘差與隨機效應無關。這相當于將每個觀察視為來自不同的研究，或者分別對每個時點的數據進行 Meta 分析。
因此，必須對 G 指定結構以考慮觀察之間的相關性。理想情況下，G 保持非結構化，以允許隨機效應對之間的不同相關性。但是，這會增加 K（K?1）/2 參數到模型中。考慮到數據通常在 Meta 分析中受到限制，這可能使得無法可靠地估計模型的參數。更復雜的結構，如復合對稱或 AR（1）更實用，因為它們只涉及一個附加參數。然而，這是以隨機效應之間更受約束的關系為代價的。
1.2.3 多變量模型：研究內和研究間相關性
到目前為止，上面所描述的模型都是假設在研究內變異中在不同的測量時刻獨立地發生（對角線 Ri），并且通過隨機效應考慮觀察之間的相關性。然而，研究內部和研究之間都存在相關性。例如，抽樣變異和測量誤差可能以類似的方式影響研究中不同時間點的估計。一般而言，導致研究之間變異的因素（例如：盲法、隨機化等）也可能導致研究內的相關性，因為所有觀測均受到這些因素的影響。另一方面，特定研究某些背景因素會在某個時點高估效果，也可能會對其他時點產生類似的影響。
到目前為止所考慮的方法僅僅獲得了觀測的邊際分布的總相關性。為了單獨考慮這些，我們通過放寬研究內協方差矩陣是對角線的假設來擴展最后一節中的模型。兩種方法是可能的：如果研究報告了估計之間的協方差，則可以使用這些值來固定 Ri 的非對角線元素，如同方差一樣。然而，這種情況很少發生。或者，可以采用混合方法，其中 Ri 的對角元素保持固定，并且從數據估計非對角線元素。然而，這需要規定 Ri 的結構，這為模型引入了更多的參數。
由于 Meta 分析中的數據通常存在限制性，因此一些簡化是為了限制添加到模型中的參數。首先，我們假設研究內的相關性在研究中是恒定的（即，研究 i 和 m 的 corr（yij，yil）=corr（ymj，yml），但是仍然允許方差不同。此外，我們采用單參數相關結構，如復合對稱或 AR（1）。
用于 Meta 分析的混合模型通常通過最大似然估計或廣義最小二乘法來估計。雖然廣義估計方程有時用于重復測量患者數據，但它們很少應用于 Meta 分析，因為它們不允許估計研究級隨機效應。這些方法很有用，例如，可用于推導每個研究中效應的經驗貝葉斯置信區間。
我們示例中的模型使用 SAS/STAT PROC MIXED 軟件通過最大似然估計。Van Houwelingen 等描述了該程序如何用于具有固定的研究內協方差的一般多變量（多結局）Meta 分析。然而，如果只有研究中的方差是固定的，并且要從數據中估計協方差分量，則需要額外的操作。這些在附錄中有描述。Nan 等用 BUGS 在貝葉斯方法中實現了一個類似的模型。
本文介紹的模型來源于 Ishak 等發表的文章[4]，張天嵩等在《實用循證醫學方法學》一書[5]中做了簡單的中文介紹，具體的模型介紹，可以參考原始文獻。
2 資料與方法
2.1 示例數據介紹
為便于說明 Ishak 等提出的基于廣義線性混合效應模型來合并縱向數據具體操作方法，我們采用《高級 Meta 分析方法—基于 Stata 實現》第十一章實例 1 的數據為實例進行演示[5]。該數據集的詳細資料見表 1 所示，其中 study 表示研究；study_id 表示研究代碼；Meas_time 代表 3 個測量時間點；numt、meant、sdt 分別代表試驗組的樣本量、均數、標準差；numc、meanc、sdc 分別代表對照組的樣本量、均數、標準差；“.”表示數據缺失。表 2 是經過轉換的橫向數據，eff_est1、eff_est2、eff_est3 分別表示時間點 1、2、3 的效應值，因此時間點 1 只有四個研究報道了效應值，時間點 2 共有 5 個研究報道了效應值，而時間點 3 只有 3 個研究報道了效應值。
表1
待分析數據集
表2
待分析數據集
2.2 SAS 代碼介紹
我們在下面描述通用的 SAS 代碼，以實現縱向效應估計的多變量 Meta 分析模型，該模型可以估計研究內和研究間的相關性。模型擬合的難度在于指定 SAS 代碼來估計研究間協方差矩陣（G）的所有參數和研究內相關性 ρR，同時保持研究內方差固定于其觀察值。
將表 1 的數據錄入 SAS 軟件，并存儲在名為 Example_data 的數據集中（因篇幅限制，數據步的 SAS 代碼簡要輸入）。
數據分析的第一步是利用試驗組和對照組的均數計算每個時間點（K=3）的效應估計值（eff_est，yij），并計算每個估計值的方差（var_est 或 Rij2）；同時利用估計值的方差權重（weight），w=1/var_est。正常情況下，不同時間點之間的協方差不會被報道，本示例同樣未提供協方差。
本文使用 PROC MIXED 擬合多變量模型，該模型必須指定研究內和研究之間協方差矩陣的結構（在我們的示例中為 AR（1））。研究間協方差矩陣對所有研究都是通用的，其所有參數都是從數據中估算出來的。另一方面，研究內協方差矩陣的規范需要特別考慮，因為必須允許方差分量在不同研究中變化，并且它們的值固定在觀察到的估計方差上，而相關參數對所有研究都是共同的（以確保可識別性），并這些參數從數據中估算。我們可以表示研究內協方差矩陣的結構如下：
 |
其中 C 是具有未知相關參數的矩陣，對于所有研究是共同的，而 Wi 對于每個研究是特定的并且必須固定于它們的觀察值。這可以使用 PROC MIXED 中的 parms 和 weight 命令完成。parms 語句首先接受協方差矩陣參數的起始值（存儲在數據集中），并允許將特定參數固定為起始值（使用 hold 選項）。我們設置了 G 和 C 的初始值，并在整個估計過程中將 C 中的方差參數保持為 1。研究內方差的倒數被指定為 PROC MIXED 中的權重，其定義為常數，并且用作上述公式中描述的研究內相關矩陣的乘數。適合多變量模型的 SAS 代碼如附錄所述。
3 結果
運行文后附錄代碼后，一共可以獲得 13 個計算表格。但其中只有 4 個表格里的內容是系統評價者關心的結果（表 3～6）。
表3
協方差參數估計結果
表4
模型估計的固定效應結果
表5
模型估計隨機效應的結果
表6
固定效應的 3 型檢驗結果
3.1 協方差參數估計結果
協方差參數估計結果見表 3。由于模型在整個估計過程中將 C 中的方差參數保持為 1。模型迭代計算后獲得了 G 的結果，反映的是隨機效應的參數向量。
3.2 固定效應的估計結果
表 4 是模型估計的固定效應結果。可見兩種干預方法在 6、12、18 個月等不同時間點的合并效應量分別為：［β=1.649，95%CI（?0.922，4.221），P=0.208］、［β=1.007，95%CI（?1.414，3.427），P=0.415］和［β=3.329，95%CI（?1.002，7.66），P=0.132］。根據相應P值，我們可以認為兩種干預措施在不同時間點的處理效應并無統計學差異。
3.3 隨機效應的估計結果
表 5 是模型估計的隨機效應結果。
3.4 固定效應的 3 型檢驗結果
結果見表 6。根據表 6 結果，Meas_time 對主要干預措施并無顯著效應，即兩組干預效應不受時間影響。
4 討論
對于縱向數據的 Meta 分析，大多數研究可能只報告了一些但不是所有時間點的效應值。在原始研究中，“未報告”的效應值被視為缺失。為了保證數據分析的有效性，縱向數據中的缺失值被認為具有時間特異性，其缺失機制被假設為完全隨機缺失（missing completely at random，MCAR）[4]。然而現實中，由于研究人員預設的隨訪時間點存在差異，在納入這些數據時，應該考慮數據之間的相關性。此時缺失數據發生的概率與所觀察到的變量是有關的，缺失值的缺失機制具有隨機丟失（missing at random，MAR）的特點。但即使如此，多變量模型也可以通過觀察到的數據信息，來考慮數據之間的相關性，總體并不影響模型的實現。在多數情況下，MCAR 與 MAR 均被稱為可忽略的缺失形式。
這些效應值的 Meta 分析，最常用的方法是將單個時間點上的數據進行單獨匯總分析。但這種方法忽略了各個時間點間的內在關系，從而忽略測量結局之間的相關結構[6]。統計推斷和模擬研究顯示：研究內相關會影響終點結局，忽略研究內相關可能會使 Meta 分析的結果存在偏差，如合并效應量均方和標準誤的估計過高[7]。多變量 Meta 分析作為單變量 Meta 分析的一個拓展，該模型可以考慮內部數據的相關性，多變量模型更好地擬合并產生更精確的匯總效應值，并通過協變量分析，對研究內及研究間的異質性進行估計。多變量模型特別適用于數據有限的資料，并可以通過多個主流統計軟件實施。Musekiwa 等[8]利用 R 語言，基于一般線性混合模型實現了多個時間點的效應大小的 Meta 分析。Trikalinos 等[9]利用 Stata 實現多個時間點的效應大小的 Meta 分析。Jansen[10]則利用網絡 Meta 分析的思想實現多個時間點的效應大小的 Meta 分析。
在本文中，筆者通過實例演示了基于廣義線性混合效應模型重復測量數據 Meta 分析的原理及方法。在筆者引用的示例研究中，共納入 6 個研究，涉及 3 個時間點和 12 組效應值，其中第三個時間點只有 3 個研究報道了效應值，因此屬于小樣本數據集。和多數研究者一樣，系統評價者不僅關系每個時間點上的匯集效應值，也會關心時間對總效應值的影響。筆者通過附錄的 SAS 軟件代碼實現了最終的 Meta 分析結果。張天嵩等主編的《高級 Meta 分析—基于 Stata 實現》一書中介紹了多變量模型在 Stata 軟件中的實現方法，可以通過 Stata 中的 mvmeta 命令來實現。該書假設所有研究、所有干預措施、所有時間點的相關系數相同，并通過假設相關系數等于 0、0.4、0.8 進行敏感性分析。最終的分析結果顯示，模型假設的相關系數對最終的合并效應值影響甚小。組內相關系數是未知的，通常也很難獲得。因此需要提前假設，然后進行敏感性分析。筆者通過 SAS 軟件，同樣實現了多變量模型的 Meta 分析，而筆者采用的代碼更為簡潔。我們假設在整個估計過程中將 C 中的方差參數保持為 1（研究內協方差參數），并將研究間協方差參數設置為 0.1，計算結果與張天嵩等的方法相似，證明筆者的代碼是行之有效的。總體而言，研究中的協方差可以使用預設的近似相關性來固定，而最終匯總效應估計并沒有太大誤差[11]。
在本文提供的實例中，各研究的隨訪時間點相同，因此將時間變量作為分類變量進行處理。但臨床實踐中，如果各研究的隨訪時間點不同時，可以將時間變量作為連續性變量處理，并通過合理的假設時間與效應關系即可以處理時間-效應。時間與效應量之間最簡單的關系可以假設為線性關系，在 SAS 語法中，可以引入“時間變量*效應值”項來考慮時間與效應變量之間的交互效應，Arends[12]提出了這種線性關系的可行性。
雖然本文的例子是基于試驗組內獲得的效應平均值及標準差，但本文描述的模型可應用于不同的效果測量，例如，計算在不同時間點的對數優勢比（odds ratio，OR），當然相對危險度（relative risk，RR）也是可用的。關于效應尺度的選擇，只要它們符合正態分布即可。
值得注意的是，相關的隨機時間效應和多變量模型需要指定合適的隨機效應（G）和/或殘差（Ri）的協方差結構[13]。在小樣本數據集中，如本文采用的示例數據以及多數 Meta 分析，多參數協方差結構（例如，非結構化協方差結構或 Toeplitz 協方差結構）可能使最終的估計復雜化或甚至無法估計（無法收斂）[14]。事實上，在本文例子中，使用這些結構會使程序不收斂或產生只有點估計的非正定的協方差矩陣估計（讀者可以在附錄代碼中指定“UN”來嘗試）。因此，我們將協方差結構設置為 AR（1）和復合對稱結構。理論上，AR（1）非常適合于縱向數據，因為引入異質性的背景因素可能隨時間變化，并且對于時間接近情況下的測量效應具有更相似的影響[15]。在多變量模型中，我們根據相關參數對 Ri 的協方差分量進行參數化，我們假設這些參數對所有研究都是共同的。也就是說，對于任何研究 i 和觀察 j 和 l，
，其中 ρjlR 是所有研究中在時間 j 和 l 的估計之間的相關性。另一方面，復合對稱結構假設相關性是恒定的，因此模型也可以獲得收斂。本文的示例數據采用這兩種協方差結構獲得的研究結果相似。但因為篇幅有限，筆者只報告基于 AR（1）結構模型的研究結果。
綜上所述，基于廣義線性混合效應模型的 Meta 分析模型可以充分考慮各個時間點之間的相關性，并準確地估計出最終的合并效應值，是縱向數據 Meta 分析的理想模型。
附錄代碼：
data example_data;
input study_id Meas_timenumt meantsdtnumc meancsdc @@;
cards;
11......
1 2 61 ?1.6 5.2 59 ?2.3 5.2
……
;
run;
data example_data;/*計算效應值及方差、權重*/
set example_data;
eff_est=meant-meanc;
var_est=sqrt(((numt?1)*sdt**2+(numc?1)*sdc**2)/(numt+numc?2));
w=1/var_est;
keep study_idMeas_time eff_est var_est w;
run;
Proc print data=example_data noobs;run;
/*設置研究內協方差*/
data within_study_cov;
param='Var?1'; est=1; output;
param='Var?2'; est=1; output;
param='Var?3'; est=1; output;
param='Corr'; est=0.5; output;
run;
/*設置研究間協方差*/
data btw_study_cov;
param='Var?1'; est=0.1; output;
param='Var?2'; est=0.1; output;
param='Var?3'; est=0.1; output;
param='Corr'; est=0.5; output;
run;
/*合并研究內協方差和研究間協方差*/
data initial_values;
set btw_study_cov within_study_cov;
keep param est;
run;
proc mixed method=REML cl data=example_data;
class study_id meas_time;
model eff_est=meas_time /noint s cl ddf=1000,1000,1000;
random meas_time /subject=study_id type=arh(1) s;
repeated meas_time /subject= study_id type=arh(1);
parms/parmsdata=initial_values hold=5 to 8;
weight w;
run;
縱向數據是一種特殊形式的數據，主要來源于醫學和社會學等領域。在醫學研究領域，縱向研究通常用于臨床和流行病學研究，縱向數據通常也被稱為重復測量數據。重復測量數據主要來自于同一受試對象在不同時間點上對同一個觀測指標進行多次測量獲得的數據[1]。被納入 Meta 分析的縱向研究通常具有以下特點：① 不同時間點上的測量對象是同一個研究內的研究對象集合，而非單一個體，因此研究對象存在失訪、脫落等現象；② 不同研究針對的觀測指標，即結局指標是相同的，但受到失訪、脫落等顯著影響；③ 不同研究的時間間隔不同，而且并非報告了全部時間點上的結局指標。對縱向研究數據進行 Meta 分析，系統評價員不僅關心主要干預措施對研究對象的整體效應，也會關心在某些特殊時間點的測量指標效應量、主要測量效應的時間趨勢以及不同時間點之間的變化值等[2]。
盡管 Meta 分析技術得到了快速發展，但目前針對重復測量數據的最佳合并技術還在探索中。縱向數據 Meta 分析的難點是解決不同時間點測量結局指標的內在相關性，這意味著此類數據的合并分析不僅僅需要考慮研究內變異和研究間變異，還需要考慮結局變量的出現時間點以及不同時間點的變化等[3]。Ishak 等[4]提出一種基于廣義線性混合效應模型來合并縱向數據的 SAS 程序，其不僅考慮了固定效應，而且也考慮了隨機效應，同時也可以獲得不同時間點的結局變量的效應值及其變化趨勢。此外，亦可以在模型中引入協變量，分析協變量對結局測量的影響等。該模型是進行縱向數據 Meta 分析的有力工具。本文基于實例，詳細介紹了廣義線性混合效應模型 Meta 分析技術在 SAS 軟件中的實現方法，并提供相關的 SAS 代碼。
1 基于廣義線性混合效應模型 Meta 分析的原理及方法簡介
1.1 縱向數據的一般線性混合模型
在介紹基于廣義線性混合效應模型 Meta 分析之前，先簡單介紹一般線性混合效應模型的基本原理，該模型是縱向數據的重要分析方法。
假設對納入分析的 N 個單位（units，此處單位是指觀測結果來源于單個個體；而 Meta 分析所指個體是 study，即單個研究）進行 K 個時間點上的結局測量；可以用 yi 表示來自第 i 個單位的觀測值，也就是 K×1 個結局向量，并用 yij 表示該單位的第 j 個觀測值。根據上述變量，可以擬合出一個可以解釋觀測值之間相關性的一般線性混合模型：
|
在上述模型中，Xiβ 是模型的固定效應部分，其中 Xi 是包含潛在時間相關的解釋變量的 K×p 的矩陣；而 β 是一個包含固定效應的 p×1 向量，是待估計的固定效應參數向量，β 反映了解釋變量 Xi 對響應變量 yi 的影響程度。Ziγi 是模型的隨機效應部分。Zi 是某些協變量的 K×q 的設計矩陣，其構成方式與 Xi 相同，但更通常的情況是 Xi 中的一部分；而 γi 是一個包含隨機效應的 q×1 向量，是待估計的隨機效應參數向量，服從均值向量為 0、方差協方差矩陣為 G 的正態分布，表示為 γi～N（0，G）。εi 是 K×1 的殘差的設計矩陣，是隨機誤差向量。在混合效應模型中，放寬了對 εi 的限制條件，其元素不必為獨立分布，即對 εi 沒有 Var（εi）=σe2 及 cov（εi，εj）=0，cov（Ei，Ej）=0 的假定。用符號表示隨機誤差向量 εi～N（0，R）。不要求 εi 的方差、協方差矩陣 R 的主對角元素為 σe2、非主對角元素為 0。該模型假設不同單位（研究）的觀察結果是獨立的，因此 cov（εij，εml）=0，其中 i≠m，在任意一個觀測結局中，j，l=1···K。模型的另外一個假設是殘差和隨機效應是獨立的，即此 cov（γi，εi）=Cov（G，R）=0，即 G 與 R 間無相關關系。
模型中，γi 的期望值 E（γi）=0，εi 的期望值 E（εi）=0，或者 E=0。
|
|
因此，響應變量 yi 的邊際分布是一個多元正態分布（multivariate normal distribution，MVN）；可以被指定為 MVNq（Xiβ，Vi），其中 Vi=var（yi）=Zi'GZi+Ri。
G 和 Ri 的協方差結構決定了如何處理個體間和個體內（研究間和研究內）的相關性。這種可以通過以下幾種方法解決。例如，可以選擇不包括任何隨機效應，并將 Ri 設置為對于所有研究都恒定的一般非結構化正定矩陣：即 Ri=R。另一種方法是引入隨機截距（也就是將 Zi 設置為 1S 的列）具有（標量）方差 G，并且設置 Ri=σ2IK，其中 IK 是 K×K 的單位矩陣。這樣就形成了一個復合對稱相關結構，其中當 j≠l 時，corr（yij，yil）=G/G+σ2。當涉及多個隨機效應時，必須給 G 指定一個結構來描述每個研究內隨機效應之間的關系。一些最常用的結構是：復合對稱（compound-symmetry，CS）或恒定相關（constant correlations，CC）、一階自回歸結構（first-order auto-regressive，AR（1））。
1.2 縱向數據 Meta 分析的混合模型
上面描述的一般模型稍微修改即可適用于縱向數據的 Meta 分析。實際上，在實踐中常用的單變量隨機效應 Meta 分析模型就是一個指定 K=1 和 Zi=1（即隨機截距）的特例。研究內方差（是單變量情況下的標量）被設置為每個研究中報告的觀察到的估計值的方差，并且假設已知且無誤差。通過該方差，可以通過倒方差來加權每個觀察的效應。
在縱向數據的 Meta 分析中，觀察的單位是在 K 個時點上收集的治療效應估計值的研究。這些時間的選擇通常是基于可用的數據和分析目標。研究之間可能會因預設的測量的時間而異，并且可能無法報告某些測量時間點的結果。如果報告的結果與效應的大小明顯相關，則可能是一種發表偏倚，類似于研究中可能發生的選擇偏倚。
Meta 分析數據的縱向模型也不同于協方差矩陣規范中的一般情況；實際上，上述一些方法在這種情況下是不合適的。例如，隨機效應模型越來越被認為更適合 Meta 分析。因此，省略隨機效應并僅考慮 Ri 的相關性（即采用固定效應方法）不是理想的。此外，研究之間可能會因樣本量、測量方法和其他方面而有所不同，這些方面將導致不同研究之間的估計值存在或多或少精確性差異。例如，由于失訪，研究中的差異也可能隨著時間的推移而變化。因此，將研究內協方差矩陣 Ri 設置為 σ2IK 是不夠的，因為它假設所有研究中的所有效應估計的方差是不變的。至少應該將 Ri 設置為是具有不同值的對角矩陣。并且，為了確保對數據進行適當的加權，我們將 Ri 的對角元素設置為研究中報告的值，正如單變量 Meta 分析。
模型的擬合始于設置固定效應，這將指導隨機效應的選擇。由于 Meta 分析的目的是總結治療隨時間的變化，因此 Xi 包括將響應變量與時間相關聯的協變量。這可能包括一些時間參數化（例如：線性、二次、對數等）。在這里，我們考慮在每個 K 測量時刻總結效應的情況。這可以通過省略截距并包括 K 個時間指標，或通過包括一個截距和 K?1 個時間指標來完成。例如，當 K=4 時，我們可以設置：
|
在第一種情況下，匯總估計值是在每個時間點獲得的。而在第二種情況下，每個時間點的匯總效應估計值的協變量測量差異與第一種情況相關，其可以通過截距估計。研究水平的協變量（例如：男性的比例、平均年齡）也可以包括在 Xi 中。這些協變量通常具有固定效應，因此不會影響協方差的相關結構。
在下文中，除非另有說明，否則我們假設 Ri 是對角矩陣，其值設置為研究中報告的估計的方差（我們用研究 i 的 Rij2 表示第 j 個對角線元素）。為了便于記憶，我們假設 Xi 只包括時間指標（即沒有截距或其他協變量）；將其他固定效應納入到模型中不會改變相關結構。
1.2.1 隨機研究效應模型
考慮觀測之間相關性的最簡單方法是隨機效應模型，它是給定研究的所有觀測的共性結果。這可以被認為是隨機截距模型。例如 K=4，我們會設置 ，這樣 γi 是一個標量；然后可以將模型寫成：
|
其中 Var（γi）=G。對所有的 i 和 j 來說，殘差向量的協方差是 Cov（εi）=Ri 和 Cov（γi，εij）=0。yij 的邊際分布的方差是 G+Rij2。在研究 i 中的 j 和 l 時刻收集的兩個觀測值之間的協方差由 Cov（yij，yil）=G 給出；不同研究的觀察結果是獨立的。因此，觀測之間的相關性是 corr（yij，yil）=G/，在這種情況下，觀測之間的不同相關性是被允許的，因為研究內方差在每個時點都可以不同。
1.2.2 隨機時間效應模型
隨機研究效應模型有一定的限制性，因為它假設研究之間的異質性以同樣的方式影響給定研究中每個測量時刻的觀察結果。然而，這可能不是真的。例如，在每項研究中，病例丟失可能會有所不同，導致后期估計的影響具有更大的異質性。 我們可以通過在每個測量時刻引入隨機效應來擴展模型。也就是說，我們設置 Zi=Xi；然后通過以下方式給出模型：
|
或者
|
因此，γi 現在是具有協方差矩陣 G 的隨機效應 K? 向量。假設這些是獨立的（即，D 是對角線）等同于假設觀測本身是獨立的，因為：
|
因為 G 和 Ri 都是對角線，殘差與隨機效應無關。這相當于將每個觀察視為來自不同的研究，或者分別對每個時點的數據進行 Meta 分析。
因此，必須對 G 指定結構以考慮觀察之間的相關性。理想情況下，G 保持非結構化，以允許隨機效應對之間的不同相關性。但是，這會增加 K（K?1）/2 參數到模型中。考慮到數據通常在 Meta 分析中受到限制，這可能使得無法可靠地估計模型的參數。更復雜的結構，如復合對稱或 AR（1）更實用，因為它們只涉及一個附加參數。然而，這是以隨機效應之間更受約束的關系為代價的。
1.2.3 多變量模型：研究內和研究間相關性
到目前為止，上面所描述的模型都是假設在研究內變異中在不同的測量時刻獨立地發生（對角線 Ri），并且通過隨機效應考慮觀察之間的相關性。然而，研究內部和研究之間都存在相關性。例如，抽樣變異和測量誤差可能以類似的方式影響研究中不同時間點的估計。一般而言，導致研究之間變異的因素（例如：盲法、隨機化等）也可能導致研究內的相關性，因為所有觀測均受到這些因素的影響。另一方面，特定研究某些背景因素會在某個時點高估效果，也可能會對其他時點產生類似的影響。
到目前為止所考慮的方法僅僅獲得了觀測的邊際分布的總相關性。為了單獨考慮這些，我們通過放寬研究內協方差矩陣是對角線的假設來擴展最后一節中的模型。兩種方法是可能的：如果研究報告了估計之間的協方差，則可以使用這些值來固定 Ri 的非對角線元素，如同方差一樣。然而，這種情況很少發生。或者，可以采用混合方法，其中 Ri 的對角元素保持固定，并且從數據估計非對角線元素。然而，這需要規定 Ri 的結構，這為模型引入了更多的參數。
由于 Meta 分析中的數據通常存在限制性，因此一些簡化是為了限制添加到模型中的參數。首先，我們假設研究內的相關性在研究中是恒定的（即，研究 i 和 m 的 corr（yij，yil）=corr（ymj，yml），但是仍然允許方差不同。此外，我們采用單參數相關結構，如復合對稱或 AR（1）。
用于 Meta 分析的混合模型通常通過最大似然估計或廣義最小二乘法來估計。雖然廣義估計方程有時用于重復測量患者數據，但它們很少應用于 Meta 分析，因為它們不允許估計研究級隨機效應。這些方法很有用，例如，可用于推導每個研究中效應的經驗貝葉斯置信區間。
我們示例中的模型使用 SAS/STAT PROC MIXED 軟件通過最大似然估計。Van Houwelingen 等描述了該程序如何用于具有固定的研究內協方差的一般多變量（多結局）Meta 分析。然而，如果只有研究中的方差是固定的，并且要從數據中估計協方差分量，則需要額外的操作。這些在附錄中有描述。Nan 等用 BUGS 在貝葉斯方法中實現了一個類似的模型。
本文介紹的模型來源于 Ishak 等發表的文章[4]，張天嵩等在《實用循證醫學方法學》一書[5]中做了簡單的中文介紹，具體的模型介紹，可以參考原始文獻。
2 資料與方法
2.1 示例數據介紹
為便于說明 Ishak 等提出的基于廣義線性混合效應模型來合并縱向數據具體操作方法，我們采用《高級 Meta 分析方法—基于 Stata 實現》第十一章實例 1 的數據為實例進行演示[5]。該數據集的詳細資料見表 1 所示，其中 study 表示研究；study_id 表示研究代碼；Meas_time 代表 3 個測量時間點；numt、meant、sdt 分別代表試驗組的樣本量、均數、標準差；numc、meanc、sdc 分別代表對照組的樣本量、均數、標準差；“.”表示數據缺失。表 2 是經過轉換的橫向數據，eff_est1、eff_est2、eff_est3 分別表示時間點 1、2、3 的效應值，因此時間點 1 只有四個研究報道了效應值，時間點 2 共有 5 個研究報道了效應值，而時間點 3 只有 3 個研究報道了效應值。
表1
待分析數據集
表2
待分析數據集
2.2 SAS 代碼介紹
我們在下面描述通用的 SAS 代碼，以實現縱向效應估計的多變量 Meta 分析模型，該模型可以估計研究內和研究間的相關性。模型擬合的難度在于指定 SAS 代碼來估計研究間協方差矩陣（G）的所有參數和研究內相關性 ρR，同時保持研究內方差固定于其觀察值。
將表 1 的數據錄入 SAS 軟件，并存儲在名為 Example_data 的數據集中（因篇幅限制，數據步的 SAS 代碼簡要輸入）。
數據分析的第一步是利用試驗組和對照組的均數計算每個時間點（K=3）的效應估計值（eff_est，yij），并計算每個估計值的方差（var_est 或 Rij2）；同時利用估計值的方差權重（weight），w=1/var_est。正常情況下，不同時間點之間的協方差不會被報道，本示例同樣未提供協方差。
本文使用 PROC MIXED 擬合多變量模型，該模型必須指定研究內和研究之間協方差矩陣的結構（在我們的示例中為 AR（1））。研究間協方差矩陣對所有研究都是通用的，其所有參數都是從數據中估算出來的。另一方面，研究內協方差矩陣的規范需要特別考慮，因為必須允許方差分量在不同研究中變化，并且它們的值固定在觀察到的估計方差上，而相關參數對所有研究都是共同的（以確保可識別性），并這些參數從數據中估算。我們可以表示研究內協方差矩陣的結構如下：
|
其中 C 是具有未知相關參數的矩陣，對于所有研究是共同的，而 Wi 對于每個研究是特定的并且必須固定于它們的觀察值。這可以使用 PROC MIXED 中的 parms 和 weight 命令完成。parms 語句首先接受協方差矩陣參數的起始值（存儲在數據集中），并允許將特定參數固定為起始值（使用 hold 選項）。我們設置了 G 和 C 的初始值，并在整個估計過程中將 C 中的方差參數保持為 1。研究內方差的倒數被指定為 PROC MIXED 中的權重，其定義為常數，并且用作上述公式中描述的研究內相關矩陣的乘數。適合多變量模型的 SAS 代碼如附錄所述。
3 結果
運行文后附錄代碼后，一共可以獲得 13 個計算表格。但其中只有 4 個表格里的內容是系統評價者關心的結果（表 3～6）。
表3
協方差參數估計結果
表4
模型估計的固定效應結果
表5
模型估計隨機效應的結果
表6
固定效應的 3 型檢驗結果
3.1 協方差參數估計結果
協方差參數估計結果見表 3。由于模型在整個估計過程中將 C 中的方差參數保持為 1。模型迭代計算后獲得了 G 的結果，反映的是隨機效應的參數向量。
3.2 固定效應的估計結果
表 4 是模型估計的固定效應結果。可見兩種干預方法在 6、12、18 個月等不同時間點的合并效應量分別為：［β=1.649，95%CI（?0.922，4.221），P=0.208］、［β=1.007，95%CI（?1.414，3.427），P=0.415］和［β=3.329，95%CI（?1.002，7.66），P=0.132］。根據相應P值，我們可以認為兩種干預措施在不同時間點的處理效應并無統計學差異。
3.3 隨機效應的估計結果
表 5 是模型估計的隨機效應結果。
3.4 固定效應的 3 型檢驗結果
結果見表 6。根據表 6 結果，Meas_time 對主要干預措施并無顯著效應，即兩組干預效應不受時間影響。
4 討論
對于縱向數據的 Meta 分析，大多數研究可能只報告了一些但不是所有時間點的效應值。在原始研究中，“未報告”的效應值被視為缺失。為了保證數據分析的有效性，縱向數據中的缺失值被認為具有時間特異性，其缺失機制被假設為完全隨機缺失（missing completely at random，MCAR）[4]。然而現實中，由于研究人員預設的隨訪時間點存在差異，在納入這些數據時，應該考慮數據之間的相關性。此時缺失數據發生的概率與所觀察到的變量是有關的，缺失值的缺失機制具有隨機丟失（missing at random，MAR）的特點。但即使如此，多變量模型也可以通過觀察到的數據信息，來考慮數據之間的相關性，總體并不影響模型的實現。在多數情況下，MCAR 與 MAR 均被稱為可忽略的缺失形式。
這些效應值的 Meta 分析，最常用的方法是將單個時間點上的數據進行單獨匯總分析。但這種方法忽略了各個時間點間的內在關系，從而忽略測量結局之間的相關結構[6]。統計推斷和模擬研究顯示：研究內相關會影響終點結局，忽略研究內相關可能會使 Meta 分析的結果存在偏差，如合并效應量均方和標準誤的估計過高[7]。多變量 Meta 分析作為單變量 Meta 分析的一個拓展，該模型可以考慮內部數據的相關性，多變量模型更好地擬合并產生更精確的匯總效應值，并通過協變量分析，對研究內及研究間的異質性進行估計。多變量模型特別適用于數據有限的資料，并可以通過多個主流統計軟件實施。Musekiwa 等[8]利用 R 語言，基于一般線性混合模型實現了多個時間點的效應大小的 Meta 分析。Trikalinos 等[9]利用 Stata 實現多個時間點的效應大小的 Meta 分析。Jansen[10]則利用網絡 Meta 分析的思想實現多個時間點的效應大小的 Meta 分析。
在本文中，筆者通過實例演示了基于廣義線性混合效應模型重復測量數據 Meta 分析的原理及方法。在筆者引用的示例研究中，共納入 6 個研究，涉及 3 個時間點和 12 組效應值，其中第三個時間點只有 3 個研究報道了效應值，因此屬于小樣本數據集。和多數研究者一樣，系統評價者不僅關系每個時間點上的匯集效應值，也會關心時間對總效應值的影響。筆者通過附錄的 SAS 軟件代碼實現了最終的 Meta 分析結果。張天嵩等主編的《高級 Meta 分析—基于 Stata 實現》一書中介紹了多變量模型在 Stata 軟件中的實現方法，可以通過 Stata 中的 mvmeta 命令來實現。該書假設所有研究、所有干預措施、所有時間點的相關系數相同，并通過假設相關系數等于 0、0.4、0.8 進行敏感性分析。最終的分析結果顯示，模型假設的相關系數對最終的合并效應值影響甚小。組內相關系數是未知的，通常也很難獲得。因此需要提前假設，然后進行敏感性分析。筆者通過 SAS 軟件，同樣實現了多變量模型的 Meta 分析，而筆者采用的代碼更為簡潔。我們假設在整個估計過程中將 C 中的方差參數保持為 1（研究內協方差參數），并將研究間協方差參數設置為 0.1，計算結果與張天嵩等的方法相似，證明筆者的代碼是行之有效的。總體而言，研究中的協方差可以使用預設的近似相關性來固定，而最終匯總效應估計并沒有太大誤差[11]。
在本文提供的實例中，各研究的隨訪時間點相同，因此將時間變量作為分類變量進行處理。但臨床實踐中，如果各研究的隨訪時間點不同時，可以將時間變量作為連續性變量處理，并通過合理的假設時間與效應關系即可以處理時間-效應。時間與效應量之間最簡單的關系可以假設為線性關系，在 SAS 語法中，可以引入“時間變量*效應值”項來考慮時間與效應變量之間的交互效應，Arends[12]提出了這種線性關系的可行性。
雖然本文的例子是基于試驗組內獲得的效應平均值及標準差，但本文描述的模型可應用于不同的效果測量，例如，計算在不同時間點的對數優勢比（odds ratio，OR），當然相對危險度（relative risk，RR）也是可用的。關于效應尺度的選擇，只要它們符合正態分布即可。
值得注意的是，相關的隨機時間效應和多變量模型需要指定合適的隨機效應（G）和/或殘差（Ri）的協方差結構[13]。在小樣本數據集中，如本文采用的示例數據以及多數 Meta 分析，多參數協方差結構（例如，非結構化協方差結構或 Toeplitz 協方差結構）可能使最終的估計復雜化或甚至無法估計（無法收斂）[14]。事實上，在本文例子中，使用這些結構會使程序不收斂或產生只有點估計的非正定的協方差矩陣估計（讀者可以在附錄代碼中指定“UN”來嘗試）。因此，我們將協方差結構設置為 AR（1）和復合對稱結構。理論上，AR（1）非常適合于縱向數據，因為引入異質性的背景因素可能隨時間變化，并且對于時間接近情況下的測量效應具有更相似的影響[15]。在多變量模型中，我們根據相關參數對 Ri 的協方差分量進行參數化，我們假設這些參數對所有研究都是共同的。也就是說，對于任何研究 i 和觀察 j 和 l，，其中 ρjlR 是所有研究中在時間 j 和 l 的估計之間的相關性。另一方面，復合對稱結構假設相關性是恒定的，因此模型也可以獲得收斂。本文的示例數據采用這兩種協方差結構獲得的研究結果相似。但因為篇幅有限，筆者只報告基于 AR（1）結構模型的研究結果。
綜上所述，基于廣義線性混合效應模型的 Meta 分析模型可以充分考慮各個時間點之間的相關性，并準確地估計出最終的合并效應值，是縱向數據 Meta 分析的理想模型。
附錄代碼：
data example_data;
input study_id Meas_timenumt meantsdtnumc meancsdc @@;
cards;
11......
1 2 61 ?1.6 5.2 59 ?2.3 5.2
……
;
run;
data example_data;/*計算效應值及方差、權重*/
set example_data;
eff_est=meant-meanc;
var_est=sqrt(((numt?1)*sdt**2+(numc?1)*sdc**2)/(numt+numc?2));
w=1/var_est;
keep study_idMeas_time eff_est var_est w;
run;
Proc print data=example_data noobs;run;
/*設置研究內協方差*/
data within_study_cov;
param='Var?1'; est=1; output;
param='Var?2'; est=1; output;
param='Var?3'; est=1; output;
param='Corr'; est=0.5; output;
run;
/*設置研究間協方差*/
data btw_study_cov;
param='Var?1'; est=0.1; output;
param='Var?2'; est=0.1; output;
param='Var?3'; est=0.1; output;
param='Corr'; est=0.5; output;
run;
/*合并研究內協方差和研究間協方差*/
data initial_values;
set btw_study_cov within_study_cov;
keep param est;
run;
proc mixed method=REML cl data=example_data;
class study_id meas_time;
model eff_est=meas_time /noint s cl ddf=1000,1000,1000;
random meas_time /subject=study_id type=arh(1) s;
repeated meas_time /subject= study_id type=arh(1);
parms/parmsdata=initial_values hold=5 to 8;
weight w;
run;
