縱向數據在不同時間點存在內在相關性問題，而傳統的 Meta 分析技術無法處理這種問題。基于多水平模型的回歸系數可以充分考慮縱向數據各個時間點的相關性。本文利用 SAS 軟件實現多水平回歸系數模型的 Meta 分析，并提供簡潔易操作的 SAS 代碼。
引用本文： 肖麗華, 鄭建清, 黃碧芬, 蘇菁菁, 吳敏. 利用 SAS 軟件實現基于多水平模型縱向數據的 Meta 分析. 中國循證醫學雜志, 2019, 19(5): 614-621. doi: 10.7507/1672-2531.201810044 復制
在涉及隨訪資料的臨床研究或流行病學研究中，縱向研究是極為常見的研究類型，這種研究設計的目的在于評估治療（干預措施）或暴露隨時間變化的影響[1]。這種類型的研究通常會在預定的時間點（time-point）進行重復的效應值測量，例如惡性腫瘤患者的生存率等[2]。縱向數據的原始研究，經常使用綜合測量（summary measure）來表示治療效應量（例如，單個試驗組內所有參與者效應值的平均值及其標準差等）。在這種情況下，治療效果是以綜合測量的絕對差異或相對差異表示（例如，治療組和對照組的平均值之差）[3]。現有針對此類數據的 Meta 分析方法，通常是提取所關心的時間點的估計值，繼而通過固定效應模型或隨機效應模型來獲取特定時間點的合并效應值，從而描述所關心的治療效果或暴露效果，但此類方法只能綜合縱向研究的單個時點上的效應量，忽略不同時間點效應量的內在關系，從而導致統計學上的分析偏倚。近年來方法學家開始重視 Meta 分析過程中的分析偏倚（analytical bias），即統計學偏倚[4]。在縱向研究中，各個時間點的綜合測量是通過對相同研究人群進行多次測量獲得的，這些研究對象可能存在失訪、脫落，從而表現為各個時間點之間特殊的內在關系，這種內在相關可以是生物學的（如疾病的自然病程）、結構性的（如相對于共同基線在不同時間計算的效應）或統計學的（如受到相似“錯誤”的測量或者來自相同方程的計算結果，如截距和斜率）[3]。由于時間相關關系的存在，通過單個時間點獲得組間的 Meta 分析結果可能會使后續 Meta 分析的結果發生混淆，甚至得出的結果互相矛盾、難以解釋。
近年來，越來越多的學者關注了縱向數據的 Meta 分析。Lopes 等[5]描述了一種基于混合多元正態分布的患者水平數據的貝葉斯模型來合并分析縱向研究效應量。Maas 等[6]描述了一種混合模型來合并分析縱向效應量，這種方法通過隨機截距和線性時間效應模型來處理觀測之間的相關性。Ishak 等[3]提出廣義線性緩和效應模型，該模型既考慮固定效應，也考慮隨機效應，可獲得更多信息，并通過 SAS 軟件的 MIXED 過程步來實現，但缺點是計算過程復雜，且無法使用森林圖提供客觀結果。Idris[7]基于三水平模型，提出一種利用回歸系數從聚合數據中獲得各個試驗的協變量回歸系數，并對回歸系數進行 Meta 分析，從而估計各個協變量的效應，這種方法可以充分利用所有時間點的數據信息，不僅可以充分考慮治療效應（干預效應），也可以分析時間效應，還可以分析治療-時間的交互效應。更重要的是，如果可以獲得各個研究的基線變量，則可以在聚合數據水平上分析協變量對治療效應的噪音影響，從而分析臨床異質性。張天嵩等[8]在《高級 Meta 分析方法－基于 Stata 實現》一書中介紹了該模型在 STATA 中的實現方法，但目前尚無統計學家利用 SAS 軟件來實現該模型的計算。在本文中，筆者利用 SAS 軟件的 REG 過程逐步獲得協變量的回歸系數，再利用 MIXED 過程逐步實現回歸系數的加權 Meta 分析，最后通過 SGPLOT 繪制森林圖，使結果清晰可觀。
1 三水平模型簡介
縱向研究數據可以被拓展為聚合水平的三水平模型。最簡單的模型可被認為是帶有一個截距和三個協變量［其中一個協變量代表干預措施（X）與時間（Time）的交互作用］的廣義線性模型，即：
 |
其中，yijk 代表第 k 個研究中干預措施 j 在時間點 i 上的觀測效應值。在評價兩個干預措施的情況下，當 j=0（通常代表對照組）時，則 Xijk 取值為 0；當 j=1（通常代表研究組）時，則 Xijk 取值為 1；在多臂試驗時，j 可擴展至 n 個試驗組。時間點 i（水平 1）嵌套于干預組別 j（水平 2），并進一步嵌套于研究水平 k（水平 3）。在該模型中，重復測量數據的層次結構變成時間點包含在干預組內，并被認為是一個潛在的協變量。該模型在水平 1 的隨機誤差項為
。
在此情況下，時間點i觀測效應量的相應方差為：
 |
影響該方差大小的因素包括研究間方差、研究內方差、研究 k 中干預措施 j 在 i 時間點上的樣本量。
假設觀測對象在縱向隨訪過程中的所有時間點上沒有失訪、脫落，那么兩個不同時間點上觀測效應的協方差為
；但實際情況研究對象隨著隨訪時間延長存在脫落可能，樣本量會逐漸減小。因此，水平 1 誤差項的方差會增大，則協方差 
，其中 n1jk 為兩個時間點的最大觀測對象人數。
值得注意的是，公式 1 是最簡單的模型，實際上該模型可以擴展至更多的協變量，讀者如感興趣可參閱相關文獻。
2 資料與方法
2.1 示例數據介紹
在介紹 SAS 程序代碼之前，筆者采用 Tan 等[9]發表的系統評價數據進行分析。該研究共納入 4 個研究比較經尿道鈥激光前列腺剜除術（HoLEP）與經尿道前列腺切除術（TURP）治療癥狀性前列腺梗阻的效果。原作者分析了 4 個觀測指標，但本文僅選擇其中一個觀測指標作為介紹。該研究報道了不同時點（分別為基線、6 個月、12 個月）的最大尿流速度（Qmax）。Qmax 為一個典型的縱向數據，其綜合測量基于聚合水平的平均值及其標準差。各個研究的原始資料見表 1。
表1
待分析數據集
2.2 SAS 代碼介紹
筆者開發的 SAS 代碼包括兩部分：第一部分利用 SAS 中的 PROC REG 對擬分析的協變量進行回歸分析，獲得回歸系數的參數估計值，并將其輸出到數據集 Pest；第二部分利用 SAS 中的 PROC MIXED 對第一部分獲得的回歸系數進行 Meta 分析，并將 Meta 分析的結果輸出到數據集 solF。最后利用 PROC SGPLOT 繪制森林圖。第一、二部分計算見代碼 1，第三部分繪制森林圖見代碼 2。
代碼 1
/*第一部分：回歸系數的估計*/
%macro reg(isize);
%do i=1 %to &isize.;
data origin&i.;
set origin;
if trial = &i.; w=1/(sd/sqrt(n)); timetrt=time*treat;
run;
proc reg data= origin&i.;
model mean=treat time timetrt;
weight w;
ods output ParameterEstimates=PEst&i.;
run;
%end;
%mend;
%reg(&trialnum);
/*輸出回歸系數數據集*/
data PEst;
set PEst1-PEst4;
if Variable = "Intercept" then do; var=1; end;
else if Variable ="treat" then do; var=2; end;
else if Variable ="time" then do; var=3; end;
else do;var=4; end;
keep Variable var Estimate StdErr;
run;
/*第二部分：對每一個回歸系數進行 meta 分析*/
data start;/*設置初始值*/
input est;
cards;
0.01;
run;
%macro meta(isize);
%do i=1 %to &isize.;
data PEst&i.;
set PEst;
if var = &i.;
run;
data PEst&i.;
set PEst&i.;
study=_N_;
run;
data covvars&i.;
set PEst&i.;
est=StdErr**2; keep est;
run;
data covvars&i.;
set start covvars&i.;
run;
proc mixed data=PEst&i. method = ml NOINFO NOITPRINT NOPROFILE;
class study;
model Estimate=/s;
random int /subject = study;
repeated/group=study;
parms/parmsdata= covvars&i. eqcons=2 to 5;
ods output SolutionF=solF&i.;
run;
%end;
%mend;
%meta(®cofnum);
data solF;
set solF1-solF4;
format Estimate 6.2;
format StdErr 6.2;
format Probt 6.3;
format Lower 6.2;
format Upper 6.2;
Lower=round(Estimate-1.96*StdErr,.01);
Upper=round(Estimate+1.96*StdErr,.01);
run;
proc print data=solF noobs; run;
代碼 2
/*第三部分：繪制森林圖*/
%macro forestdata(isize);
%do i=1 %to &isize.;
data PEst&i.;
set PEst&i.;
format Estimate 6.2;
keep Estimate StdErr study;
run;
data solF&i.;
set solF&i.;
keep Estimate StdErr Probt;
run;
data forest&i.(rename=(study=studyNUM));
set PEst&i. solF&i.;
run;
%end;
%mend;
%forestdata(®cofnum);
%macro Anno_data(isize);
%do i=1 %to &isize.;
data forest&i.;
set forest&i.;
if studyNUM=. then do;wi=0;end;
else do;wi=1/(StdErr)**2;end;
run;
proc sql;
create table forest&i. as select StudyNUM,Estimate,StdErr,wi,Probt, sum(wi) as sumwi from forest&i.;
quit;
data forest&i.(rename=(Estimate=Beta));
set forest&i.;
format weight percent5.;
Estimate=round(Estimate,.01);
LowerCL=round(Estimate-1.96*StdErr,.01);
UpperCL=round(Estimate+1.96*StdErr,.01);
Weight=round(wi/sumwi,.0001);
if Weight=0 then do;Weight=.; end;
run;
data forest&i.;
set forest&i.;
length Study $ 16;
format Q1 Q3 6.2;
beta=round(beta,.01);
if studyNUM=1 then do;Study="Gupta";end;
else if studyNUM=2 then do;Study="Tan";end;
else if studyNUM=3 then do;Study="Kuntz";end;
else if studyNUM=4 then do;Study="Montorsi"; end;
else do;Study="";end;
lcl2=lowercl;
ucl2=uppercl;
Q1=min(round(Beta-weight,.01),round(Beta+weight,.01));
Q3=max(round(Beta-weight,.01),round(Beta+weight,.01));
ES='ES'; LCL='LCL'; UCL='UCL'; WT='Weight';Pvalue='Pvalue';
drop wi studyNUM StdErr sumwi;
run;
%end;
%mend;
%Anno_data(®cofnum);
data forest;
set forest1-forest4;
arm=ceil(_n_/5);
if Study ="" then do; Overall="Overall"; end;
else do; Overall=""; end;
run;
proc print data=forest;run;
%macro Meta_forest(isize);
%do i=1 %to &isize.;
data forest&i.;
set forest;
beta_signal="-β";
if arm=&i. and Study ^="" then do; Study&i.=cats(Study,beta_signal,&i.-1); end;
else do; study&i.=""; end;
if arm=&i. and Overall ^="" then do;Overall&i.=cats(Overall,beta_signal,&i.-1);end;
else do;Overall&i.="";end;
drop beta_signal;
run;
%end;
%mend;
%Meta_forest(®cofnum);
data forest;
merge forest1-forest4;
if Overall1^="" then do;studyref=Overall1;end;
if Overall2^="" then do;studyref=Overall2;end;
if Overall3^="" then do;studyref=Overall3;end;
if Overall4^="" then do;studyref=Overall4;end;
run;
proc print data=forest;run;
data _null_;
pct=0.75/nobs;
height=6;
dpi=100;
call symputx("pct", pct);
call symputx("pct2", 2*pct);
call symputx("dpi", dpi);
call symputx("height", height);
call symputx("heightin", height || "in");
call symputx("thickness", floor(height*dpi*pct));
set forest nobs=nobs;
run;
/*創建自定義森林圖樣式*/
proc template;
define style styles.ForestColor93;
parent = Styles.htmlBlue;
style GraphFonts from GraphFonts /
'GraphDataFont' = ("<sans-serif>, <MTsans-serif>",7pt)
'GraphValueFont' = ("<sans-serif>, <MTsans-serif>",7pt);
style GraphData1 from GraphData1 /
contrastcolor = GraphColors('gcdata2')
color = GraphColors('gdata2');
style GraphData2 from GraphData2 /
contrastcolor = GraphColors('gcdata1')
color = GraphColors('gdata1');
end;
run;
/*設置森林圖樣式，DPI，圖像參數和標題等*/
ods html;/*closeods html close 可以只生成圖片到文件夾;*/
ods graphics / reset width=6in height=&heightin imagename="ForestPlotColor";
ods listing sge=off style=Styles.ForestColor93 image_dpi=&dpi;
title h=12pt "Impact of Covariates on Survival by Study";
title2 h=8pt 'The estimated regression coefficients and 95% CL';
%macro sgplot(isize);
proc sgplot data=forest noautolegend nocycleattrs;
refline studyref / lineattrs=(thickness=&thickness) transparency=0.85;
%do i=1 %to &isize.;
highlow y=study&i. low=lcl2 high=ucl2 / type=line;
highlow y=study&i. low=q1 high=q3 / type=bar;
scatter y=study&i. x=ES / markerchar=beta x2axis;
scatter y=study&i. x=lcl / markerchar=lowercl x2axis;
scatter y=study&i. x=ucl / markerchar=uppercl x2axis;
scatter y=study&i. x=wt / markerchar=weight x2axis;
scatter y=Overall&i. x=beta / markerattrs=(symbol=diamondfilled size=10);
scatter y=Overall&i. x=ES / markerchar=beta x2axis;
scatter y=Overall&i. x=lcl / markerchar=lowercl x2axis;
scatter y=Overall&i. x=ucl / markerchar=uppercl x2axis;
scatter y=Overall&i. x=wt / markerchar=weight x2axis;
scatter y=Overall&i. x=Pvalue / markerchar=Probt x2axis;
%end;
refline -10 10 / axis=x lineattrs=(pattern=shortdash) transparency=0.5;
inset 'Unfavorable factor' / position=bottomleft;
inset 'Favorable factor' / position=bottom;
xaxis type=LINEAR offsetmin=0 offsetmax=0.35 min=-30 max=30 minor display=(nolabel) ;
x2axis offsetmin=0.7 display=(noticks nolabel);
yaxis display=(noticks nolabel) offsetmin=&pct offsetmax=&pct2 reverse;
run;
%mend;
%sgplot(4);
title;
3 結果
3.1 經典 Meta 分析結果
傳統重復測量的 Meta 分析通常采用重要時間點（relevant time-point meta-analysis，RTM）、最初/最終時間點（first/final time-point meta-analysis，FTM）或所有時間點 Meta 分析（all time-point meta-analysis，ATM）。作為對比，筆者也提供了傳統重復測量的 Meta 分析結果（表 2）。該結果可采用 RevMan 軟件、STATA 軟件或 R 軟件實現。在表 2 中，0 時點代表基線，1 為最初時間點 Meta 分析結果，2 代表最終時間點 Meta 分析結果，All 代表所有時間點 Meta 分析結果。最初時間點 Meta 分析結果顯示，與 TURP 相比，HoLEP 并不能顯著提高 Qmax（P>0.05）；而最終時間點 Meta 分析和所有時間點 Meta 分析結果卻顯示，與 TURP 相比，HoLEP 可以顯著提高 Qmax（P<0.05），不同時間點的 Meta 分析結果存在矛盾。
表2
經典 Meta 分析結果
3.2 多水平回歸系數 Meta 分析模型結果
運行代碼 1 后，即可獲得表 3 結果。表 3 中，截距 β0 代表各個試驗中不同干預組基線 Qmax 的 Meta 分析結果，β1 代表干預措施（我們采用 0 代表 HoLEP，1 代表 TURP），β2 代表不同時間點（采用 0 代表基線時點，1 代表 6 個月時點、2 代表 12 個月時點），β3 代表干預與時間的交互作用（干預*時間）。三水平模型的 Meta 分析結果顯示，Qmax 與各干預組基線顯著相關，且基線 Qmax 越高，后期 Qmax 越高［β0=8.93，95%CI（1.83，16.03），P=0.028］。與 TURP 相比，HoLEP 并不能顯著提高 Qmax［β1=–0.44，95%CI（–10.55，9.67），P=0.899］。不同時間點的 Qmax 存在顯著差異（P=0.02），且隨著時間的增加，Qmax 逐漸增加［β2=9.42，95%CI（2.82，16.02）］。干預措施與時間點間無顯著的交互作用［β3=–0.39，95%CI（–10.29，9.5），P=0.907］。利用代碼 2 可以獲得多水平回歸系數 Meta 分析模型的森林圖（圖 1）。
表3
多水平回歸系數 Meta 分析模型結果
圖1
多水平回歸系數 Meta 分析的森林圖
4 討論
我們對 Tan 等[9]的 Meta 分析數據進行的再次分析顯示縱向數據的 Meta 分析不能忽略各個時間點的內在相關性。經典 Meta 分析結果顯示：在最初時間點，與 TURP 相比，HoLEP 可提高 Qmax，但差異并無統計學差異（P=0.528）；6 個月時，HoLEP 也不能顯著提高 Qmax；然而最終時間點和所有時間點的 Meta 分析結果卻顯示，與 TURP 相比，HoLEP 可顯著提高 Qmax（P<0.05）。不同時間點的 Meta 分析結果存在矛盾，因此難以說明與 TURP 相比，HoLEP 是否可以提高 Qmax。經典 Meta 分析結果更無法得知，第 12 個月的 Qmax 是因為治療因素所致，還是時間因素所致。而我們采用 Idris[7]提出的三水平模型，則使計算變得簡單且易于解釋。將治療、時間和治療-時間作為自變量，引入三水平模型得到的結果表明：兩種干預方法對前列腺梗阻患者 Qmax 的影響差異沒有統計學意義，而時間則是影響 Qmax 的重要因素。另一方面，三水平模型 Meta 分析結果表明，Qmax 與各干預組基線顯著相關，基線 Qmax 越高，后期 Qmax 越高。
上述結果表明，對于縱向數據，基于回歸系數的三水平模型的 Meta 分析結果比傳統的 Meta 分析結果更容易獲得可靠、易懂、易解釋的統計學結果。利用三水平模型，可通過擬合多協變量效應模型來計算總體估計，并獲得各個協變量的效應估計，甚至分析基線變量對總體效應的影響。與基于單水平數據產生的 Meta 估計相比，對回歸系數執行多響應模型提供了非常好的效應估計。雖然基于單水平數據的 Meta 分析方法因其非常簡單更容易被系統評價研究者選擇，但其只能提供固定參數的估計值。而多水平模型可以獲得 3 級水平的隨機參數的估計值，可能是更優選的方法。然而，多水平模型的應用需要對 2 級協方差矩陣進行初始假設，并且涉及更復雜的統計基礎。盡管如此，筆者認為三水平模型是縱向數據 Meta 分析的最優被擇方案。如果單個研究樣本量越大，則回歸系數估計值預計會更精確，因為回歸系數的估計值將具有較小的標準誤。反過來，該研究在 Meta 分析中的所占權重會更大，從而提供更好的總體估計。該方法的另一個優點是即使只有少數的原始研究被納入 Meta 分析，也可使用本文的方法進行 Meta 分析，而這恰恰是 Meta 分析的常見問題之一。
近年來，多變量模型亦被提出用于縱向數據的 Meta 分析[1]，但該模型需要對原始數據的協方差進行估計，因為多數原始研究不會報道協方差。而且當研究涉及很多的時間點，那么多變量模型將會失去很大功效，因為得到的多變量模型具有很高的維度。因此，所關心時間點的數量以及不同時間點之間的時間跨度都可能受到可用數據的影響[10]。因此，時間點太多或太少，均無法通過多變量模型獲得精確的效應估計。此外，如果原始研究在完全不同的時間點給出綜合測量數據的情況下，多變量模型將完全不可行[1]。而本文提出的基于回歸系數的三水平模型的 Meta 分析方法則不需要考慮時間點數量和間隔的影響，可將時間點直接作為連續性變量引入模型。
值得注意的是，基于回歸系數的三水平模型的 Meta 分析方法尚存在以下局限性：① 該模型需要完成兩個階段的參數估計，即回歸系數的估計和總體參數的估計，因此可能帶來 Ⅱ 類錯誤風險；② 在獲得回歸系數的估計時，需要假設各個系數之間具有獨立性，而且假設不同時間段內的平均效應具有線性；③ 在估計回歸系數時，設置的 H0 假設為 β=0，即回歸系數=0。最后，該方法無法估計研究內相關性，這是在 Meta 分析中對聚合水平重復測量數據進行建模的主要問題之一，也是目前縱向數據 Meta 分析懸而未決的重要問題之一。
總之，多水平回歸系數模型可以充分考慮并解釋縱向研究中不同時點效應值之間的相關性，可以提供更好的數據擬合和更精確的合并效應估計，另外該方法還計算了時間的影響以及治療-時間的相互作用，故對于縱向數據的 Meta 分析，多水平分析模型值得推廣應用。
在涉及隨訪資料的臨床研究或流行病學研究中，縱向研究是極為常見的研究類型，這種研究設計的目的在于評估治療（干預措施）或暴露隨時間變化的影響[1]。這種類型的研究通常會在預定的時間點（time-point）進行重復的效應值測量，例如惡性腫瘤患者的生存率等[2]。縱向數據的原始研究，經常使用綜合測量（summary measure）來表示治療效應量（例如，單個試驗組內所有參與者效應值的平均值及其標準差等）。在這種情況下，治療效果是以綜合測量的絕對差異或相對差異表示（例如，治療組和對照組的平均值之差）[3]。現有針對此類數據的 Meta 分析方法，通常是提取所關心的時間點的估計值，繼而通過固定效應模型或隨機效應模型來獲取特定時間點的合并效應值，從而描述所關心的治療效果或暴露效果，但此類方法只能綜合縱向研究的單個時點上的效應量，忽略不同時間點效應量的內在關系，從而導致統計學上的分析偏倚。近年來方法學家開始重視 Meta 分析過程中的分析偏倚（analytical bias），即統計學偏倚[4]。在縱向研究中，各個時間點的綜合測量是通過對相同研究人群進行多次測量獲得的，這些研究對象可能存在失訪、脫落，從而表現為各個時間點之間特殊的內在關系，這種內在相關可以是生物學的（如疾病的自然病程）、結構性的（如相對于共同基線在不同時間計算的效應）或統計學的（如受到相似“錯誤”的測量或者來自相同方程的計算結果，如截距和斜率）[3]。由于時間相關關系的存在，通過單個時間點獲得組間的 Meta 分析結果可能會使后續 Meta 分析的結果發生混淆，甚至得出的結果互相矛盾、難以解釋。
近年來，越來越多的學者關注了縱向數據的 Meta 分析。Lopes 等[5]描述了一種基于混合多元正態分布的患者水平數據的貝葉斯模型來合并分析縱向研究效應量。Maas 等[6]描述了一種混合模型來合并分析縱向效應量，這種方法通過隨機截距和線性時間效應模型來處理觀測之間的相關性。Ishak 等[3]提出廣義線性緩和效應模型，該模型既考慮固定效應，也考慮隨機效應，可獲得更多信息，并通過 SAS 軟件的 MIXED 過程步來實現，但缺點是計算過程復雜，且無法使用森林圖提供客觀結果。Idris[7]基于三水平模型，提出一種利用回歸系數從聚合數據中獲得各個試驗的協變量回歸系數，并對回歸系數進行 Meta 分析，從而估計各個協變量的效應，這種方法可以充分利用所有時間點的數據信息，不僅可以充分考慮治療效應（干預效應），也可以分析時間效應，還可以分析治療-時間的交互效應。更重要的是，如果可以獲得各個研究的基線變量，則可以在聚合數據水平上分析協變量對治療效應的噪音影響，從而分析臨床異質性。張天嵩等[8]在《高級 Meta 分析方法－基于 Stata 實現》一書中介紹了該模型在 STATA 中的實現方法，但目前尚無統計學家利用 SAS 軟件來實現該模型的計算。在本文中，筆者利用 SAS 軟件的 REG 過程逐步獲得協變量的回歸系數，再利用 MIXED 過程逐步實現回歸系數的加權 Meta 分析，最后通過 SGPLOT 繪制森林圖，使結果清晰可觀。
1 三水平模型簡介
縱向研究數據可以被拓展為聚合水平的三水平模型。最簡單的模型可被認為是帶有一個截距和三個協變量［其中一個協變量代表干預措施（X）與時間（Time）的交互作用］的廣義線性模型，即：
|
其中，yijk 代表第 k 個研究中干預措施 j 在時間點 i 上的觀測效應值。在評價兩個干預措施的情況下，當 j=0（通常代表對照組）時，則 Xijk 取值為 0；當 j=1（通常代表研究組）時，則 Xijk 取值為 1；在多臂試驗時，j 可擴展至 n 個試驗組。時間點 i（水平 1）嵌套于干預組別 j（水平 2），并進一步嵌套于研究水平 k（水平 3）。在該模型中，重復測量數據的層次結構變成時間點包含在干預組內，并被認為是一個潛在的協變量。該模型在水平 1 的隨機誤差項為。
在此情況下，時間點i觀測效應量的相應方差為：
|
影響該方差大小的因素包括研究間方差、研究內方差、研究 k 中干預措施 j 在 i 時間點上的樣本量。
假設觀測對象在縱向隨訪過程中的所有時間點上沒有失訪、脫落，那么兩個不同時間點上觀測效應的協方差為；但實際情況研究對象隨著隨訪時間延長存在脫落可能，樣本量會逐漸減小。因此，水平 1 誤差項的方差會增大，則協方差 ，其中 n1jk 為兩個時間點的最大觀測對象人數。
值得注意的是，公式 1 是最簡單的模型，實際上該模型可以擴展至更多的協變量，讀者如感興趣可參閱相關文獻。
2 資料與方法
2.1 示例數據介紹
在介紹 SAS 程序代碼之前，筆者采用 Tan 等[9]發表的系統評價數據進行分析。該研究共納入 4 個研究比較經尿道鈥激光前列腺剜除術（HoLEP）與經尿道前列腺切除術（TURP）治療癥狀性前列腺梗阻的效果。原作者分析了 4 個觀測指標，但本文僅選擇其中一個觀測指標作為介紹。該研究報道了不同時點（分別為基線、6 個月、12 個月）的最大尿流速度（Qmax）。Qmax 為一個典型的縱向數據，其綜合測量基于聚合水平的平均值及其標準差。各個研究的原始資料見表 1。
表1
待分析數據集
2.2 SAS 代碼介紹
筆者開發的 SAS 代碼包括兩部分：第一部分利用 SAS 中的 PROC REG 對擬分析的協變量進行回歸分析，獲得回歸系數的參數估計值，并將其輸出到數據集 Pest；第二部分利用 SAS 中的 PROC MIXED 對第一部分獲得的回歸系數進行 Meta 分析，并將 Meta 分析的結果輸出到數據集 solF。最后利用 PROC SGPLOT 繪制森林圖。第一、二部分計算見代碼 1，第三部分繪制森林圖見代碼 2。
代碼 1
/*第一部分：回歸系數的估計*/
%macro reg(isize);
%do i=1 %to &isize.;
data origin&i.;
set origin;
if trial = &i.; w=1/(sd/sqrt(n)); timetrt=time*treat;
run;
proc reg data= origin&i.;
model mean=treat time timetrt;
weight w;
ods output ParameterEstimates=PEst&i.;
run;
%end;
%mend;
%reg(&trialnum);
/*輸出回歸系數數據集*/
data PEst;
set PEst1-PEst4;
if Variable = "Intercept" then do; var=1; end;
else if Variable ="treat" then do; var=2; end;
else if Variable ="time" then do; var=3; end;
else do;var=4; end;
keep Variable var Estimate StdErr;
run;
/*第二部分：對每一個回歸系數進行 meta 分析*/
data start;/*設置初始值*/
input est;
cards;
0.01;
run;
%macro meta(isize);
%do i=1 %to &isize.;
data PEst&i.;
set PEst;
if var = &i.;
run;
data PEst&i.;
set PEst&i.;
study=_N_;
run;
data covvars&i.;
set PEst&i.;
est=StdErr**2; keep est;
run;
data covvars&i.;
set start covvars&i.;
run;
proc mixed data=PEst&i. method = ml NOINFO NOITPRINT NOPROFILE;
class study;
model Estimate=/s;
random int /subject = study;
repeated/group=study;
parms/parmsdata= covvars&i. eqcons=2 to 5;
ods output SolutionF=solF&i.;
run;
%end;
%mend;
%meta(®cofnum);
data solF;
set solF1-solF4;
format Estimate 6.2;
format StdErr 6.2;
format Probt 6.3;
format Lower 6.2;
format Upper 6.2;
Lower=round(Estimate-1.96*StdErr,.01);
Upper=round(Estimate+1.96*StdErr,.01);
run;
proc print data=solF noobs; run;
代碼 2
/*第三部分：繪制森林圖*/
%macro forestdata(isize);
%do i=1 %to &isize.;
data PEst&i.;
set PEst&i.;
format Estimate 6.2;
keep Estimate StdErr study;
run;
data solF&i.;
set solF&i.;
keep Estimate StdErr Probt;
run;
data forest&i.(rename=(study=studyNUM));
set PEst&i. solF&i.;
run;
%end;
%mend;
%forestdata(®cofnum);
%macro Anno_data(isize);
%do i=1 %to &isize.;
data forest&i.;
set forest&i.;
if studyNUM=. then do;wi=0;end;
else do;wi=1/(StdErr)**2;end;
run;
proc sql;
create table forest&i. as select StudyNUM,Estimate,StdErr,wi,Probt, sum(wi) as sumwi from forest&i.;
quit;
data forest&i.(rename=(Estimate=Beta));
set forest&i.;
format weight percent5.;
Estimate=round(Estimate,.01);
LowerCL=round(Estimate-1.96*StdErr,.01);
UpperCL=round(Estimate+1.96*StdErr,.01);
Weight=round(wi/sumwi,.0001);
if Weight=0 then do;Weight=.; end;
run;
data forest&i.;
set forest&i.;
length Study $ 16;
format Q1 Q3 6.2;
beta=round(beta,.01);
if studyNUM=1 then do;Study="Gupta";end;
else if studyNUM=2 then do;Study="Tan";end;
else if studyNUM=3 then do;Study="Kuntz";end;
else if studyNUM=4 then do;Study="Montorsi"; end;
else do;Study="";end;
lcl2=lowercl;
ucl2=uppercl;
Q1=min(round(Beta-weight,.01),round(Beta+weight,.01));
Q3=max(round(Beta-weight,.01),round(Beta+weight,.01));
ES='ES'; LCL='LCL'; UCL='UCL'; WT='Weight';Pvalue='Pvalue';
drop wi studyNUM StdErr sumwi;
run;
%end;
%mend;
%Anno_data(®cofnum);
data forest;
set forest1-forest4;
arm=ceil(_n_/5);
if Study ="" then do; Overall="Overall"; end;
else do; Overall=""; end;
run;
proc print data=forest;run;
%macro Meta_forest(isize);
%do i=1 %to &isize.;
data forest&i.;
set forest;
beta_signal="-β";
if arm=&i. and Study ^="" then do; Study&i.=cats(Study,beta_signal,&i.-1); end;
else do; study&i.=""; end;
if arm=&i. and Overall ^="" then do;Overall&i.=cats(Overall,beta_signal,&i.-1);end;
else do;Overall&i.="";end;
drop beta_signal;
run;
%end;
%mend;
%Meta_forest(®cofnum);
data forest;
merge forest1-forest4;
if Overall1^="" then do;studyref=Overall1;end;
if Overall2^="" then do;studyref=Overall2;end;
if Overall3^="" then do;studyref=Overall3;end;
if Overall4^="" then do;studyref=Overall4;end;
run;
proc print data=forest;run;
data _null_;
pct=0.75/nobs;
height=6;
dpi=100;
call symputx("pct", pct);
call symputx("pct2", 2*pct);
call symputx("dpi", dpi);
call symputx("height", height);
call symputx("heightin", height || "in");
call symputx("thickness", floor(height*dpi*pct));
set forest nobs=nobs;
run;
/*創建自定義森林圖樣式*/
proc template;
define style styles.ForestColor93;
parent = Styles.htmlBlue;
style GraphFonts from GraphFonts /
'GraphDataFont' = ("<sans-serif>, <MTsans-serif>",7pt)
'GraphValueFont' = ("<sans-serif>, <MTsans-serif>",7pt);
style GraphData1 from GraphData1 /
contrastcolor = GraphColors('gcdata2')
color = GraphColors('gdata2');
style GraphData2 from GraphData2 /
contrastcolor = GraphColors('gcdata1')
color = GraphColors('gdata1');
end;
run;
/*設置森林圖樣式，DPI，圖像參數和標題等*/
ods html;/*closeods html close 可以只生成圖片到文件夾;*/
ods graphics / reset width=6in height=&heightin imagename="ForestPlotColor";
ods listing sge=off style=Styles.ForestColor93 image_dpi=&dpi;
title h=12pt "Impact of Covariates on Survival by Study";
title2 h=8pt 'The estimated regression coefficients and 95% CL';
%macro sgplot(isize);
proc sgplot data=forest noautolegend nocycleattrs;
refline studyref / lineattrs=(thickness=&thickness) transparency=0.85;
%do i=1 %to &isize.;
highlow y=study&i. low=lcl2 high=ucl2 / type=line;
highlow y=study&i. low=q1 high=q3 / type=bar;
scatter y=study&i. x=ES / markerchar=beta x2axis;
scatter y=study&i. x=lcl / markerchar=lowercl x2axis;
scatter y=study&i. x=ucl / markerchar=uppercl x2axis;
scatter y=study&i. x=wt / markerchar=weight x2axis;
scatter y=Overall&i. x=beta / markerattrs=(symbol=diamondfilled size=10);
scatter y=Overall&i. x=ES / markerchar=beta x2axis;
scatter y=Overall&i. x=lcl / markerchar=lowercl x2axis;
scatter y=Overall&i. x=ucl / markerchar=uppercl x2axis;
scatter y=Overall&i. x=wt / markerchar=weight x2axis;
scatter y=Overall&i. x=Pvalue / markerchar=Probt x2axis;
%end;
refline -10 10 / axis=x lineattrs=(pattern=shortdash) transparency=0.5;
inset 'Unfavorable factor' / position=bottomleft;
inset 'Favorable factor' / position=bottom;
xaxis type=LINEAR offsetmin=0 offsetmax=0.35 min=-30 max=30 minor display=(nolabel) ;
x2axis offsetmin=0.7 display=(noticks nolabel);
yaxis display=(noticks nolabel) offsetmin=&pct offsetmax=&pct2 reverse;
run;
%mend;
%sgplot(4);
title;
3 結果
3.1 經典 Meta 分析結果
傳統重復測量的 Meta 分析通常采用重要時間點（relevant time-point meta-analysis，RTM）、最初/最終時間點（first/final time-point meta-analysis，FTM）或所有時間點 Meta 分析（all time-point meta-analysis，ATM）。作為對比，筆者也提供了傳統重復測量的 Meta 分析結果（表 2）。該結果可采用 RevMan 軟件、STATA 軟件或 R 軟件實現。在表 2 中，0 時點代表基線，1 為最初時間點 Meta 分析結果，2 代表最終時間點 Meta 分析結果，All 代表所有時間點 Meta 分析結果。最初時間點 Meta 分析結果顯示，與 TURP 相比，HoLEP 并不能顯著提高 Qmax（P>0.05）；而最終時間點 Meta 分析和所有時間點 Meta 分析結果卻顯示，與 TURP 相比，HoLEP 可以顯著提高 Qmax（P<0.05），不同時間點的 Meta 分析結果存在矛盾。
表2
經典 Meta 分析結果
3.2 多水平回歸系數 Meta 分析模型結果
運行代碼 1 后，即可獲得表 3 結果。表 3 中，截距 β0 代表各個試驗中不同干預組基線 Qmax 的 Meta 分析結果，β1 代表干預措施（我們采用 0 代表 HoLEP，1 代表 TURP），β2 代表不同時間點（采用 0 代表基線時點，1 代表 6 個月時點、2 代表 12 個月時點），β3 代表干預與時間的交互作用（干預*時間）。三水平模型的 Meta 分析結果顯示，Qmax 與各干預組基線顯著相關，且基線 Qmax 越高，后期 Qmax 越高［β0=8.93，95%CI（1.83，16.03），P=0.028］。與 TURP 相比，HoLEP 并不能顯著提高 Qmax［β1=–0.44，95%CI（–10.55，9.67），P=0.899］。不同時間點的 Qmax 存在顯著差異（P=0.02），且隨著時間的增加，Qmax 逐漸增加［β2=9.42，95%CI（2.82，16.02）］。干預措施與時間點間無顯著的交互作用［β3=–0.39，95%CI（–10.29，9.5），P=0.907］。利用代碼 2 可以獲得多水平回歸系數 Meta 分析模型的森林圖（圖 1）。
表3
多水平回歸系數 Meta 分析模型結果
圖1
多水平回歸系數 Meta 分析的森林圖
4 討論
我們對 Tan 等[9]的 Meta 分析數據進行的再次分析顯示縱向數據的 Meta 分析不能忽略各個時間點的內在相關性。經典 Meta 分析結果顯示：在最初時間點，與 TURP 相比，HoLEP 可提高 Qmax，但差異并無統計學差異（P=0.528）；6 個月時，HoLEP 也不能顯著提高 Qmax；然而最終時間點和所有時間點的 Meta 分析結果卻顯示，與 TURP 相比，HoLEP 可顯著提高 Qmax（P<0.05）。不同時間點的 Meta 分析結果存在矛盾，因此難以說明與 TURP 相比，HoLEP 是否可以提高 Qmax。經典 Meta 分析結果更無法得知，第 12 個月的 Qmax 是因為治療因素所致，還是時間因素所致。而我們采用 Idris[7]提出的三水平模型，則使計算變得簡單且易于解釋。將治療、時間和治療-時間作為自變量，引入三水平模型得到的結果表明：兩種干預方法對前列腺梗阻患者 Qmax 的影響差異沒有統計學意義，而時間則是影響 Qmax 的重要因素。另一方面，三水平模型 Meta 分析結果表明，Qmax 與各干預組基線顯著相關，基線 Qmax 越高，后期 Qmax 越高。
上述結果表明，對于縱向數據，基于回歸系數的三水平模型的 Meta 分析結果比傳統的 Meta 分析結果更容易獲得可靠、易懂、易解釋的統計學結果。利用三水平模型，可通過擬合多協變量效應模型來計算總體估計，并獲得各個協變量的效應估計，甚至分析基線變量對總體效應的影響。與基于單水平數據產生的 Meta 估計相比，對回歸系數執行多響應模型提供了非常好的效應估計。雖然基于單水平數據的 Meta 分析方法因其非常簡單更容易被系統評價研究者選擇，但其只能提供固定參數的估計值。而多水平模型可以獲得 3 級水平的隨機參數的估計值，可能是更優選的方法。然而，多水平模型的應用需要對 2 級協方差矩陣進行初始假設，并且涉及更復雜的統計基礎。盡管如此，筆者認為三水平模型是縱向數據 Meta 分析的最優被擇方案。如果單個研究樣本量越大，則回歸系數估計值預計會更精確，因為回歸系數的估計值將具有較小的標準誤。反過來，該研究在 Meta 分析中的所占權重會更大，從而提供更好的總體估計。該方法的另一個優點是即使只有少數的原始研究被納入 Meta 分析，也可使用本文的方法進行 Meta 分析，而這恰恰是 Meta 分析的常見問題之一。
近年來，多變量模型亦被提出用于縱向數據的 Meta 分析[1]，但該模型需要對原始數據的協方差進行估計，因為多數原始研究不會報道協方差。而且當研究涉及很多的時間點，那么多變量模型將會失去很大功效，因為得到的多變量模型具有很高的維度。因此，所關心時間點的數量以及不同時間點之間的時間跨度都可能受到可用數據的影響[10]。因此，時間點太多或太少，均無法通過多變量模型獲得精確的效應估計。此外，如果原始研究在完全不同的時間點給出綜合測量數據的情況下，多變量模型將完全不可行[1]。而本文提出的基于回歸系數的三水平模型的 Meta 分析方法則不需要考慮時間點數量和間隔的影響，可將時間點直接作為連續性變量引入模型。
值得注意的是，基于回歸系數的三水平模型的 Meta 分析方法尚存在以下局限性：① 該模型需要完成兩個階段的參數估計，即回歸系數的估計和總體參數的估計，因此可能帶來 Ⅱ 類錯誤風險；② 在獲得回歸系數的估計時，需要假設各個系數之間具有獨立性，而且假設不同時間段內的平均效應具有線性；③ 在估計回歸系數時，設置的 H0 假設為 β=0，即回歸系數=0。最后，該方法無法估計研究內相關性，這是在 Meta 分析中對聚合水平重復測量數據進行建模的主要問題之一，也是目前縱向數據 Meta 分析懸而未決的重要問題之一。
總之，多水平回歸系數模型可以充分考慮并解釋縱向研究中不同時點效應值之間的相關性，可以提供更好的數據擬合和更精確的合并效應估計，另外該方法還計算了時間的影響以及治療-時間的相互作用，故對于縱向數據的 Meta 分析，多水平分析模型值得推廣應用。
