基于不同研究劑量反應數據之間是否存在異質性，以及分析過程中是否需要考慮非線性成分，劑量反應關系的Meta分析過程中有多種模型可供選擇。鑒于目前國內外并沒有相對明確的指南來指導如何規范化劑量反應關系的Meta分析流程，本文在系統總結既往研究的基礎上，詳細闡述如何結合不同研究間劑量反應數據之間的異質性、擬合優度檢驗結果及模型參數部分總體無效假設檢驗的結果，選擇正確的分析模型，并在此基礎上初步歸納出劑量反應關系Meta分析的分析流程。
引用本文： 黃育北, 李衛芹, 席波, 榮瑩, 宋豐舉, 陳可欣. 劑量反應關系Meta分析的模型選擇及分析流程. 中國循證醫學雜志, 2016, 16(2): 223-228. doi: 10.7507/1672-2531.20160035 復制
Meta分析方法，由于能夠將同類研究結果進行合并，以達到增加樣本量，提高檢驗效能的目的，從而成為循證醫學領域重要的定量分析方法 [1-3]。但是對于暴露為分類變量的資料，傳統的Meta分析方法只能簡單合并二分類變量對應的效應值，無法利用多分類變量對應的劑量反應關系數據，從而在一定程度上浪費了大量有效數據。在此基礎上，尋求可以分析劑量反應關系的Meta分析方法隨之而生。
劑量反應關系的Meta分析方法最早于1992年由Greenland和Longnecker提出 [4]。其原理是運用廣義最小二乘法來估計暴露因素每增加單位劑量時，結局效應的變化值。這種估計的基本假設為線性假設，即構建的劑量反應關系為線性的劑量反應關系 [4, 5]。隨著線性劑量反應關系研究的深入，更多的研究者發現了很多特殊類型的非線性劑量反應類型 [6-9]，如表現為J型或U型曲線（或倒J型或倒U型曲線）等。這種特殊類型的劑量反應關系，提示暴露引起的人群疾病風險分布中存在一類風險最低或最高的人群，尋找這類人群相應的暴露水平，具有非常重要的公共衛生學意義。基于此類需求，非線性劑量反應關系的Meta分析方法隨之迅速發展。目前Stata中已有相應的模塊——GLST模塊（Generalized least squares for trend estimation of summarized dose-response data），用于Meta分析中構建線性及非線性劑量反應關系模型 [10, 11]。但是由于GLST本身構建的多種模型之間存在較大的相似性，同時模塊本身并沒有提供一個相對明確的指南來指導如何規范化劑量反應關系的Meta分析流程，因此，國內外研究者在運用GLST模塊進行劑量反應關系的Meta分析時，在模型的選擇以及結果的解釋上并不統一。因此，本文在系統總結既往研究的基礎上，探討如何結合不同研究間劑量反應數據之間的異質性、擬合優度檢驗結果，及模型參數部分總體無效假設檢驗的結果，選擇正確的分析模型，并在此基礎上初步歸納出劑量反應關系Meta分析的分析流程。
1 數據簡介
本文采用的數據來源于GLST方法開發者所采用的示例數據之一，即飲酒與肺癌發病風險劑量反應關系的Meta分析數據 [10]，如表 1所示。該項Meta分析共納入4個原始研究（由于原始文獻采用作者名的縮寫來標注不同研究，這里予以保留，規范的分析中建議采用第一作者+發表年份表示）。每一個研究中，飲酒量（dose）均存在5個不同的水平。每一個水平上均報告了對應的觀察人年數（peryears）和相應的肺癌發生數（cases），以及調整了其他混雜變量之后，相對于參考組水平，其他飲酒組肺癌發病的相對風險（logrr）及其相應的標準誤（se）。由于納入的原始研究均為基于人時的發病率數據，因此，數據類型（type）均為GLST模塊中的第二類。
表1
4個原始研究不同飲酒攝入量對應的肺癌發病風險及相應的其他基本信息
2 運用GLST模塊進行劑量反應關系Meta分析的規范化流程
關于GLST模塊各項參數的意義，本文不加以贅述，詳細參考GLST幫助（具體可以在Stata中輸入help glst來了解詳細內容） [11]。本文將在系統總結既往研究的基礎上 [12-15]，詳細闡述進行劑量反應關系Meta分析的模型選擇及分析流程。
2.1 第一步：根據線性模型擬合優度檢驗結果判斷模型是否需要考慮不同研究間線性趨勢的異質性
模型選擇：兩階段線性模型，先估計單個研究的線性趨勢，然后再合并各個研究的趨勢估計。固定或隨機效應模型均可以，因為線性模型中，不論是固定效應模型還是隨機效應模型，并不會改變模型擬合優度檢驗（意義等價于異質性檢驗）的P值（原理等價于方差分析中將固定效應模型的組間變異分解成隨機效應模型中的組間變異及隨機變異之和，但因變量的總體變異不變）。
詳解：
（1）構建兩階段固定效應線性模型
glst logrr dose,se(se) cov(peryears cases) pfirst(study type) ts(f)
模型總體分析結果如下：
Two-stage fixed-effects dose-response model Number of studies=4
Generalized least-squares regression Number of obs=4
Goodness-of-fit chi2(3) =6.13 Model chi2(1) =9.50
Prob>chi2=0.1055 Prob>chi2=0.0021
結果解釋：模型擬合優度檢驗的P值為0.105?5>0.05。因此，可以認為該模擬的擬合結果沒有問題，也即線性模型中可以不用考慮不同研究間線性劑量反應關系之間的異質性。
（2）構建兩階段隨機效應線性模型
glst logrr dose,se(se) cov(peryears cases) pfirst(study type) ts(r)
模型總體分析結果如下：
Two-stage random-effects dose-response model Number of studies=4
Generalized least-squares regression Number of obs=4
Goodness-of-fit chi2(3) =6.13 Model chi2(1) =4.40
Prob>chi2=0.1055 Prob>chi2=0.0358
結果解釋：模型擬合優度檢驗的P值與兩階段固定效應線性模型的擬合優度檢驗的P值是一致的（均為0.105?5），意義同上。因此，在第一步的分析，只需選取任意一種兩階段線性模型，來判斷異質性即可。
2.2 第二步：根據異質性判斷結果，選擇正確線性模型，分析線性趨勢
模型選擇：根據異質性的判斷結果，選擇兩階段隨機（有異質性）或固定（無異質性）效應線性模型。
詳解：根據第一步分析結果，可以認為線性模型中可以不用考慮不同研究間線性劑量反應之間的異質性，所以選擇兩階段固定效應線性模型。相應模型參數估計結果如下：
logrr | Coef. Std. Err. z P>|z| [95% Conf. Interval]
dose | .005548.0018002 3.08 0.002 .0020197 .0090762
結果解釋：模型參數dose無效假設檢驗的P值為0.002，與模型參數部分總體無效假設檢驗P值一致。因為當模型只有一個參數時，模型參數部分總體無效假設檢驗結果的P值等價于模型唯一參數無效假設檢驗的P值。進一步可以通過lincom函數功能是將線性模型參數估計值β轉換為RR值或OR值。
lincom dose,eform
結果如下：
(1) dose=0
logrr | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
(1) | 1.005563.0018102 3.08 0.002 1.002022 1.009118
從上述結果可以看到：飲酒量每增加一個單位（g/d），肺癌的發病風險增加0.56%。此時，線性劑量反應關系分析目的已經達到。但從表 1可以看到，納入的4個原始研究的飲酒量與肺癌劑量反應關系似乎并不完全服從線性，初步提示應該探討非線性的劑量反應關系。
表2
不同研究線性趨勢點估計值及95%置信區間
2.3 第三步：綜合模型擬合優度檢驗及模型參數部分總體無效假設檢驗結果，選擇正確非線性模型，分析非線性趨勢
模型選擇：綜合模型擬合優度檢驗及模型參數部分總體無效假設檢驗結果，選擇固定或隨機效應非線性模型。
詳解：非線性模型中，由于不符合線性方差分析的前提條件，因此，固定效應模型與隨機效應模型總體擬合優度檢驗的P值是不一致的[10, 11]，因此，不能簡單根據模型擬合優度檢驗的P值來選擇固定效應還是隨機效應的非線性模型，此時需同時參考兩模型的參數部分總體無效假設檢驗的P值。由于非線性模型通常采用限制性立方樣條（restricted cubic splines，RCS）方法來構建[10, 16]，因此，本文也采用RCS來構建模型的非線性成分，并采用常用的4節點RCS方法（具體可以在Stata輸入help mkspline來了解詳細內容）來分別構建固定或隨機效應非線性模型。
（1）構建固定效應非線性模型
glst logrr doses*,se(se) cov(peryears cases)pfirst(study type)
模型總體分析結果如下：
Fixed-effects dose-response model Number of studies=4
Generalized least-squares regression Number of obs=16
Goodness-of-fit chi2(13) =22.10 Model chi2(3) =12.50
Prob>chi2=0.0539 Prob>chi2=0.0059
（2）構建隨機效應非線性模型
glst logrr doses*,se(se) cov(peryears cases)pfirst(study type) random
模型總體分析結果如下：
Random-effects dose-response model Number of studies=4
Iterative Generalized least-squares regression Number of obs=16
Goodness-of-fit chi2(13) =14.78 Model chi2(3) =4.25
Prob>chi2=0.3215 Prob>chi2=0.2357
可以看到，兩模型擬合優度檢驗的P值均>0.05，進一步說明非線性模型中也可以不用考慮不同研究間非線性劑量反應關系之間的異質性。同時，固定效應非線性模型參數部分的總體無效假設檢驗的P值為0.005?9<0.05，而隨機效應非線性模型參數部分的總體無效假設檢驗的P值為0.235?7>0.05，說明固定效應非線性模型中納入參數總體上拒絕無效假設，也即具有統計學意義。相反，隨機效應非線性模型中納入參數總體上是沒有統計學意義的。綜合說明本例中采用固定效應的非線性模型是合適的。如果兩模型參數部分的總體無效假設檢驗的P值均<0.05，則說明兩種模型均適合，結合異質性檢驗（也即模型擬合優度檢驗）結果選擇其中一個作為最終模型。反之，如果兩模型參數部分的總體無效假設檢驗的P值均>0.05，則提示非線性模型是不合適的。
在此基礎上，可進一步檢驗模型非線性成分是否有統計學意義。根據RCS模型構建非線性成分的原理 [10, 16]，RCS模型中第一個變量默認為線性成分。此外，所有其他非線性變量的綜合即為RCS模型的非線性成分。其中，非線性變量的數目等于節點數（n）減2，即n-2。由于本例中采用常用的4節點固定效應非線性模型中，因此，dose1即為模型的線性成分，dose2和dose3共同構成模型的非線性成分。
logrr | Coef. Std. Err. z P>|z| [95% Conf. Interval]
doses1 |.0166507 .0174593 0.95 0.340 -.0175689 .0508703
doses2 | -.4920428 .4698667 -1.05 0.295 -1.412965 .428879
doses3 | .5779587 .5406277 1.07 0.285 -.4816523 1.63757
因此，檢驗本例中非線性成分是否有統計學意義，即是計算dose2和dose3綜合的無效假設檢驗的P值，方法如下：
testparm doses2 doses3
結果如下：
(1) doses2=0
(2) doses3=0
chi2(2) =3.00
Prob>chi2=0.2229
可以看到，模型非線性成分總體無效假設檢驗的P值為0.222?9>0.05，說明模型非線性成分不具有統計學意義。同時模型線性成分的P值為0.340>0.05，說明模型線性成分也不具有統計學意義。但從前面模型參數部分（包括線性成分及非線性成分）的總體無效假設檢驗P值（0.005?9）可以看到，說明模型參數部分總體上是有意義的。因此，似乎出現一個矛盾，模型參數部分總體上有意義，但參數部分的線性成分和非線性成分均不具有統計學意義。此時，模型的結果只能解釋為：模型因變量的變異是由線性成分與非線性成分來共同決定，而不取決于其中任一成分。這樣的解釋似乎難以理解，而且研究者通常也會忽略模型參數部分總體無效假設檢驗的P值，而將該結果解釋成不存在非線性趨勢。然而，這種解釋是不合適的。因為這種現象在非線性模型中很常見。此時，繪制非線性劑量反應圖形顯得尤為重要。
2.4 第四步：繪制非線性圖形
從圖 1（具體編碼省略）可以看出飲酒量與肺癌發病風險之間確實存在非線性關系。
圖1
飲酒量與肺癌發病風險關系
3 運用GLST模塊進行劑量反應關系Meta分析的常見誤用
誤區1：未采用兩階段固定（或隨機）效應線性模型，而直接采用簡單的固定（或隨機）效應線性模型來判斷線性模型的異質性。如下兩個模型：
方法1：glst logrr dose,se(se) cov(peryears cases) pfirst(study type) ts(f)
方法2：glst logrr dose,se(se) cov(peryears cases) pfirst(study type)
兩種方法的輸出結果在參數估計上并沒有差別，唯一差別是模型擬合優度檢驗的P值。第一種方法擬合優度檢驗的P值為0.105?5（表明模型擬合較好，可以不用考慮不同研究線性趨勢之間的異質性），第二種方法（構建的是固定效應的一般線性模型）模型擬合優度檢驗的P值為0.048?6（表明模型擬合不好，提示應該考慮不同研究間的異質性）。出現這種差異的原因是：第二種方法人為假定不同研究來自于“同一總體”，而忽略了不同研究之間可能存在的異質性。
誤區2：簡單地先采用GLST單獨計算每篇研究的線性趨勢RR值，再用metan模塊合并不同研究的線性趨勢RR值，得到綜合的線性趨勢RR值。這種做法等同于分解了GLST的兩階段模型，如果確實存在線性關系，其計算結果與GLST兩階段線性模型結果是一致的。但該方法無法用于處理非線性關系。即使是采用上述方法用于處理線性關系，由于對方法不熟悉，同樣也容易出錯。具體可以采用mvmeta_make（在stata中輸入help mvmeta_make可獲得詳細幫助）程序來實現GLST兩階段線性模型構建原理的詳細分解計算過程：
mvmeta_make glst logrr dose,se(se) cov(peryears cases) pfirst(study type) saving(test) replace by(study) names(b V)
use test,clear
metan bdose bdoselci bdoseuci,eform effect(RR) lcols(study)
由于mvmeta_make模塊只能保存線性點估計值（b）和斜方差陣（V），不能保存置信區間。因此metan模塊采用的點估計值的95%上限（bdoselci）和下限（bdoseuci）分別來自mvmeta_make語句的中間輸出結果，整理如下：
最終metan分析得到的4個研究的合并趨勢［RR=1.006，（1.002，1.009）］，同時異質性檢驗的P值和總體模型的P值分別為0.106和0.002（如下），這與前述兩階段固定效應線性模型結果一致，進一步說明兩者本質是相同的。但由于該方法同時能輸出異質性定量評價指標I2的大小，這自然是對GLST模塊一個重要補充（GLST模塊本身并不能提供I2）。
Heterogeneity chi-squared=6.13(d.f.=3) p=0.106
I-squared (variation in ES attributable to heterogeneity)=51.0%
Test of ES=1 : z=3.08 p=0.002
誤區3：GLST分析前，將原始研究暴露因素的參考組水平進行中心化處理，也即參考組的暴露水平在中心化處理后變為0。這種處理帶來的后果包括：① 人為改變了原始研究的暴露劑量與結局之間的真實關聯。如原始研究的膳食纖維攝入量分別為11.5、14.3、16.4、18.8和22.9 g/d，中心化處理后，膳食纖維攝入量分別為0、2.8、4.9、7.3和11.4，但對應的RR其實并沒有修改 [11]。由此，可以看到真實的劑量反應關系已經被歪曲。② 人為增加了零組附近的觀察數，同時暴露劑量的有效變化范圍明顯縮小。如上例中的膳食纖維攝入量在中心化處理后，暴露劑量的有效變化范圍由11.5～22.9變成0～11.4。在納入研究數很多，而且每篇研究的參考組暴露水平分布很廣的話，這種處理會大大縮小暴露劑量的有效變化范圍。因此，在GLST分析中，對原始研究的參考暴露水平進行中心化處理會導致研究結果嚴重偏離真實的觀察結果。
綜上所述，隨著劑量反應關系原始研究的增加，未來基于劑量反應關系的Meta分析也會增加。目前雖然有很好的方法（如GLST）來進行該項分析，但為了得到更加準確的分析結果，以及更加合理地解釋結果，非常需要規范化劑量反應關系的Meta分析流程。本文在系統總結既往研究的基礎上初步歸納出該方法的規范化流程，未來需要更多的證據來完善該流程。
Meta分析方法，由于能夠將同類研究結果進行合并，以達到增加樣本量，提高檢驗效能的目的，從而成為循證醫學領域重要的定量分析方法 [1-3]。但是對于暴露為分類變量的資料，傳統的Meta分析方法只能簡單合并二分類變量對應的效應值，無法利用多分類變量對應的劑量反應關系數據，從而在一定程度上浪費了大量有效數據。在此基礎上，尋求可以分析劑量反應關系的Meta分析方法隨之而生。
劑量反應關系的Meta分析方法最早于1992年由Greenland和Longnecker提出 [4]。其原理是運用廣義最小二乘法來估計暴露因素每增加單位劑量時，結局效應的變化值。這種估計的基本假設為線性假設，即構建的劑量反應關系為線性的劑量反應關系 [4, 5]。隨著線性劑量反應關系研究的深入，更多的研究者發現了很多特殊類型的非線性劑量反應類型 [6-9]，如表現為J型或U型曲線（或倒J型或倒U型曲線）等。這種特殊類型的劑量反應關系，提示暴露引起的人群疾病風險分布中存在一類風險最低或最高的人群，尋找這類人群相應的暴露水平，具有非常重要的公共衛生學意義。基于此類需求，非線性劑量反應關系的Meta分析方法隨之迅速發展。目前Stata中已有相應的模塊——GLST模塊（Generalized least squares for trend estimation of summarized dose-response data），用于Meta分析中構建線性及非線性劑量反應關系模型 [10, 11]。但是由于GLST本身構建的多種模型之間存在較大的相似性，同時模塊本身并沒有提供一個相對明確的指南來指導如何規范化劑量反應關系的Meta分析流程，因此，國內外研究者在運用GLST模塊進行劑量反應關系的Meta分析時，在模型的選擇以及結果的解釋上并不統一。因此，本文在系統總結既往研究的基礎上，探討如何結合不同研究間劑量反應數據之間的異質性、擬合優度檢驗結果，及模型參數部分總體無效假設檢驗的結果，選擇正確的分析模型，并在此基礎上初步歸納出劑量反應關系Meta分析的分析流程。
1 數據簡介
本文采用的數據來源于GLST方法開發者所采用的示例數據之一，即飲酒與肺癌發病風險劑量反應關系的Meta分析數據 [10]，如表 1所示。該項Meta分析共納入4個原始研究（由于原始文獻采用作者名的縮寫來標注不同研究，這里予以保留，規范的分析中建議采用第一作者+發表年份表示）。每一個研究中，飲酒量（dose）均存在5個不同的水平。每一個水平上均報告了對應的觀察人年數（peryears）和相應的肺癌發生數（cases），以及調整了其他混雜變量之后，相對于參考組水平，其他飲酒組肺癌發病的相對風險（logrr）及其相應的標準誤（se）。由于納入的原始研究均為基于人時的發病率數據，因此，數據類型（type）均為GLST模塊中的第二類。
表1
4個原始研究不同飲酒攝入量對應的肺癌發病風險及相應的其他基本信息
2 運用GLST模塊進行劑量反應關系Meta分析的規范化流程
關于GLST模塊各項參數的意義，本文不加以贅述，詳細參考GLST幫助（具體可以在Stata中輸入help glst來了解詳細內容） [11]。本文將在系統總結既往研究的基礎上 [12-15]，詳細闡述進行劑量反應關系Meta分析的模型選擇及分析流程。
2.1 第一步：根據線性模型擬合優度檢驗結果判斷模型是否需要考慮不同研究間線性趨勢的異質性
模型選擇：兩階段線性模型，先估計單個研究的線性趨勢，然后再合并各個研究的趨勢估計。固定或隨機效應模型均可以，因為線性模型中，不論是固定效應模型還是隨機效應模型，并不會改變模型擬合優度檢驗（意義等價于異質性檢驗）的P值（原理等價于方差分析中將固定效應模型的組間變異分解成隨機效應模型中的組間變異及隨機變異之和，但因變量的總體變異不變）。
詳解：
（1）構建兩階段固定效應線性模型
glst logrr dose,se(se) cov(peryears cases) pfirst(study type) ts(f)
模型總體分析結果如下：
Two-stage fixed-effects dose-response model Number of studies=4
Generalized least-squares regression Number of obs=4
Goodness-of-fit chi2(3) =6.13 Model chi2(1) =9.50
Prob>chi2=0.1055 Prob>chi2=0.0021
結果解釋：模型擬合優度檢驗的P值為0.105?5>0.05。因此，可以認為該模擬的擬合結果沒有問題，也即線性模型中可以不用考慮不同研究間線性劑量反應關系之間的異質性。
（2）構建兩階段隨機效應線性模型
glst logrr dose,se(se) cov(peryears cases) pfirst(study type) ts(r)
模型總體分析結果如下：
Two-stage random-effects dose-response model Number of studies=4
Generalized least-squares regression Number of obs=4
Goodness-of-fit chi2(3) =6.13 Model chi2(1) =4.40
Prob>chi2=0.1055 Prob>chi2=0.0358
結果解釋：模型擬合優度檢驗的P值與兩階段固定效應線性模型的擬合優度檢驗的P值是一致的（均為0.105?5），意義同上。因此，在第一步的分析，只需選取任意一種兩階段線性模型，來判斷異質性即可。
2.2 第二步：根據異質性判斷結果，選擇正確線性模型，分析線性趨勢
模型選擇：根據異質性的判斷結果，選擇兩階段隨機（有異質性）或固定（無異質性）效應線性模型。
詳解：根據第一步分析結果，可以認為線性模型中可以不用考慮不同研究間線性劑量反應之間的異質性，所以選擇兩階段固定效應線性模型。相應模型參數估計結果如下：
logrr | Coef. Std. Err. z P>|z| [95% Conf. Interval]
dose | .005548.0018002 3.08 0.002 .0020197 .0090762
結果解釋：模型參數dose無效假設檢驗的P值為0.002，與模型參數部分總體無效假設檢驗P值一致。因為當模型只有一個參數時，模型參數部分總體無效假設檢驗結果的P值等價于模型唯一參數無效假設檢驗的P值。進一步可以通過lincom函數功能是將線性模型參數估計值β轉換為RR值或OR值。
lincom dose,eform
結果如下：
(1) dose=0
logrr | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
(1) | 1.005563.0018102 3.08 0.002 1.002022 1.009118
從上述結果可以看到：飲酒量每增加一個單位（g/d），肺癌的發病風險增加0.56%。此時，線性劑量反應關系分析目的已經達到。但從表 1可以看到，納入的4個原始研究的飲酒量與肺癌劑量反應關系似乎并不完全服從線性，初步提示應該探討非線性的劑量反應關系。
表2
不同研究線性趨勢點估計值及95%置信區間
2.3 第三步：綜合模型擬合優度檢驗及模型參數部分總體無效假設檢驗結果，選擇正確非線性模型，分析非線性趨勢
模型選擇：綜合模型擬合優度檢驗及模型參數部分總體無效假設檢驗結果，選擇固定或隨機效應非線性模型。
詳解：非線性模型中，由于不符合線性方差分析的前提條件，因此，固定效應模型與隨機效應模型總體擬合優度檢驗的P值是不一致的[10, 11]，因此，不能簡單根據模型擬合優度檢驗的P值來選擇固定效應還是隨機效應的非線性模型，此時需同時參考兩模型的參數部分總體無效假設檢驗的P值。由于非線性模型通常采用限制性立方樣條（restricted cubic splines，RCS）方法來構建[10, 16]，因此，本文也采用RCS來構建模型的非線性成分，并采用常用的4節點RCS方法（具體可以在Stata輸入help mkspline來了解詳細內容）來分別構建固定或隨機效應非線性模型。
（1）構建固定效應非線性模型
glst logrr doses*,se(se) cov(peryears cases)pfirst(study type)
模型總體分析結果如下：
Fixed-effects dose-response model Number of studies=4
Generalized least-squares regression Number of obs=16
Goodness-of-fit chi2(13) =22.10 Model chi2(3) =12.50
Prob>chi2=0.0539 Prob>chi2=0.0059
（2）構建隨機效應非線性模型
glst logrr doses*,se(se) cov(peryears cases)pfirst(study type) random
模型總體分析結果如下：
Random-effects dose-response model Number of studies=4
Iterative Generalized least-squares regression Number of obs=16
Goodness-of-fit chi2(13) =14.78 Model chi2(3) =4.25
Prob>chi2=0.3215 Prob>chi2=0.2357
可以看到，兩模型擬合優度檢驗的P值均>0.05，進一步說明非線性模型中也可以不用考慮不同研究間非線性劑量反應關系之間的異質性。同時，固定效應非線性模型參數部分的總體無效假設檢驗的P值為0.005?9<0.05，而隨機效應非線性模型參數部分的總體無效假設檢驗的P值為0.235?7>0.05，說明固定效應非線性模型中納入參數總體上拒絕無效假設，也即具有統計學意義。相反，隨機效應非線性模型中納入參數總體上是沒有統計學意義的。綜合說明本例中采用固定效應的非線性模型是合適的。如果兩模型參數部分的總體無效假設檢驗的P值均<0.05，則說明兩種模型均適合，結合異質性檢驗（也即模型擬合優度檢驗）結果選擇其中一個作為最終模型。反之，如果兩模型參數部分的總體無效假設檢驗的P值均>0.05，則提示非線性模型是不合適的。
在此基礎上，可進一步檢驗模型非線性成分是否有統計學意義。根據RCS模型構建非線性成分的原理 [10, 16]，RCS模型中第一個變量默認為線性成分。此外，所有其他非線性變量的綜合即為RCS模型的非線性成分。其中，非線性變量的數目等于節點數（n）減2，即n-2。由于本例中采用常用的4節點固定效應非線性模型中，因此，dose1即為模型的線性成分，dose2和dose3共同構成模型的非線性成分。
logrr | Coef. Std. Err. z P>|z| [95% Conf. Interval]
doses1 |.0166507 .0174593 0.95 0.340 -.0175689 .0508703
doses2 | -.4920428 .4698667 -1.05 0.295 -1.412965 .428879
doses3 | .5779587 .5406277 1.07 0.285 -.4816523 1.63757
因此，檢驗本例中非線性成分是否有統計學意義，即是計算dose2和dose3綜合的無效假設檢驗的P值，方法如下：
testparm doses2 doses3
結果如下：
(1) doses2=0
(2) doses3=0
chi2(2) =3.00
Prob>chi2=0.2229
可以看到，模型非線性成分總體無效假設檢驗的P值為0.222?9>0.05，說明模型非線性成分不具有統計學意義。同時模型線性成分的P值為0.340>0.05，說明模型線性成分也不具有統計學意義。但從前面模型參數部分（包括線性成分及非線性成分）的總體無效假設檢驗P值（0.005?9）可以看到，說明模型參數部分總體上是有意義的。因此，似乎出現一個矛盾，模型參數部分總體上有意義，但參數部分的線性成分和非線性成分均不具有統計學意義。此時，模型的結果只能解釋為：模型因變量的變異是由線性成分與非線性成分來共同決定，而不取決于其中任一成分。這樣的解釋似乎難以理解，而且研究者通常也會忽略模型參數部分總體無效假設檢驗的P值，而將該結果解釋成不存在非線性趨勢。然而，這種解釋是不合適的。因為這種現象在非線性模型中很常見。此時，繪制非線性劑量反應圖形顯得尤為重要。
2.4 第四步：繪制非線性圖形
從圖 1（具體編碼省略）可以看出飲酒量與肺癌發病風險之間確實存在非線性關系。
圖1
飲酒量與肺癌發病風險關系
3 運用GLST模塊進行劑量反應關系Meta分析的常見誤用
誤區1：未采用兩階段固定（或隨機）效應線性模型，而直接采用簡單的固定（或隨機）效應線性模型來判斷線性模型的異質性。如下兩個模型：
方法1：glst logrr dose,se(se) cov(peryears cases) pfirst(study type) ts(f)
方法2：glst logrr dose,se(se) cov(peryears cases) pfirst(study type)
兩種方法的輸出結果在參數估計上并沒有差別，唯一差別是模型擬合優度檢驗的P值。第一種方法擬合優度檢驗的P值為0.105?5（表明模型擬合較好，可以不用考慮不同研究線性趨勢之間的異質性），第二種方法（構建的是固定效應的一般線性模型）模型擬合優度檢驗的P值為0.048?6（表明模型擬合不好，提示應該考慮不同研究間的異質性）。出現這種差異的原因是：第二種方法人為假定不同研究來自于“同一總體”，而忽略了不同研究之間可能存在的異質性。
誤區2：簡單地先采用GLST單獨計算每篇研究的線性趨勢RR值，再用metan模塊合并不同研究的線性趨勢RR值，得到綜合的線性趨勢RR值。這種做法等同于分解了GLST的兩階段模型，如果確實存在線性關系，其計算結果與GLST兩階段線性模型結果是一致的。但該方法無法用于處理非線性關系。即使是采用上述方法用于處理線性關系，由于對方法不熟悉，同樣也容易出錯。具體可以采用mvmeta_make（在stata中輸入help mvmeta_make可獲得詳細幫助）程序來實現GLST兩階段線性模型構建原理的詳細分解計算過程：
mvmeta_make glst logrr dose,se(se) cov(peryears cases) pfirst(study type) saving(test) replace by(study) names(b V)
use test,clear
metan bdose bdoselci bdoseuci,eform effect(RR) lcols(study)
由于mvmeta_make模塊只能保存線性點估計值（b）和斜方差陣（V），不能保存置信區間。因此metan模塊采用的點估計值的95%上限（bdoselci）和下限（bdoseuci）分別來自mvmeta_make語句的中間輸出結果，整理如下：
最終metan分析得到的4個研究的合并趨勢［RR=1.006，（1.002，1.009）］，同時異質性檢驗的P值和總體模型的P值分別為0.106和0.002（如下），這與前述兩階段固定效應線性模型結果一致，進一步說明兩者本質是相同的。但由于該方法同時能輸出異質性定量評價指標I2的大小，這自然是對GLST模塊一個重要補充（GLST模塊本身并不能提供I2）。
Heterogeneity chi-squared=6.13(d.f.=3) p=0.106
I-squared (variation in ES attributable to heterogeneity)=51.0%
Test of ES=1 : z=3.08 p=0.002
誤區3：GLST分析前，將原始研究暴露因素的參考組水平進行中心化處理，也即參考組的暴露水平在中心化處理后變為0。這種處理帶來的后果包括：① 人為改變了原始研究的暴露劑量與結局之間的真實關聯。如原始研究的膳食纖維攝入量分別為11.5、14.3、16.4、18.8和22.9 g/d，中心化處理后，膳食纖維攝入量分別為0、2.8、4.9、7.3和11.4，但對應的RR其實并沒有修改 [11]。由此，可以看到真實的劑量反應關系已經被歪曲。② 人為增加了零組附近的觀察數，同時暴露劑量的有效變化范圍明顯縮小。如上例中的膳食纖維攝入量在中心化處理后，暴露劑量的有效變化范圍由11.5～22.9變成0～11.4。在納入研究數很多，而且每篇研究的參考組暴露水平分布很廣的話，這種處理會大大縮小暴露劑量的有效變化范圍。因此，在GLST分析中，對原始研究的參考暴露水平進行中心化處理會導致研究結果嚴重偏離真實的觀察結果。
綜上所述，隨著劑量反應關系原始研究的增加，未來基于劑量反應關系的Meta分析也會增加。目前雖然有很好的方法（如GLST）來進行該項分析，但為了得到更加準確的分析結果，以及更加合理地解釋結果，非常需要規范化劑量反應關系的Meta分析流程。本文在系統總結既往研究的基礎上初步歸納出該方法的規范化流程，未來需要更多的證據來完善該流程。
