劑量-反應 Meta 分析在循證證據產生及臨床決策方面的應用日益廣泛。這種研究方法是通過將針對相同臨床或公共衛生問題的不同研究的某種劑量-反應關系加權合并以獲得“平均”劑量效應，以全面反映某種暴露與結局之間的劑量-效應關系。當前使用最為廣泛的劑量-反應 Meta 分析模型是基于經典的二階段法，該方法的優勢在于可采用固定或隨機效應模型按照異質性大小進行效應量合并。而實現二階段隨機效應模型下劑量-反應 Meta 分析有兩種不同的方式，即普通模型和考慮參數間相關性的隨機效應模型。本文對這兩種隨機模型進行介紹，并通過 Stata 軟件演示相關操作過程，以期為循證實踐提供理論參考。
引用本文： 張超, 高崢巖, 黃靜宇, 李勝, 牛玉明, 徐暢, 韓芳芳. 劑量-反應 Meta 分析之兩種不同隨機效應模型的應用. 中國循證醫學雜志, 2017, 17(5): 616-620. doi: 10.7507/1672-2531.201703042 復制
流行病學研究的重要任務是探討或建立暴露/干預與疾病結局的因果關系（cause-and-effect relationship）[1]。探討不同水平下的暴露/干預與結局指標的動態關系稱為劑量-反應關系（dose-response relationship）；而將針對相同臨床或公共衛生問題的不同研究的這種劑量-反應關系加權合并以獲得“平均”劑量效應的方法稱為劑量-反應 Meta 分析（dose-response meta-analysis，DRMA）[2]。
當前使用最為廣泛的 DRMA 模型是基于經典的二階段法（two-stage approach）下實現的，即首先通過線性或非線性模型擬合出單篇研究的回歸系數（第一階段），再采用固定或隨機效應模型將第一階段的回歸系數進行合并（第二階段）[3,4]。固定效應模型的 DRMA 因無需考慮研究間異質性，在模型開發及實現上較為簡單；而對于隨機效應模型的 DRMA，特別是存在多個參數時（即非線性模型），軟件程序開發及實現較為困難。
本文將介紹兩種基于 Stata 軟件、實現多參數隨機效應模型的 DRMA 方法，并代入睡眠與全因死亡率的數據分別進行 DRMA，以期為循證實踐提供相關參考。
1 資料來源
1.1 軟件代碼來源
本文中涉及的部分軟件代碼由 Suhail A.R Doi 教授和徐暢博士提供，并已獲得他們同意。
1.2 數據來源和基本特征
以 Liu 等[5]2017 年發表的關于睡眠時間與全因死亡率的劑量-反應 Meta 分析為例，采用兩種不同方法重新進行 DRMA。該研究從 51 篇關于睡眠時間與全因死亡率的前瞻性隊列研究中，納入 40 篇符合劑量-反應 Meta 分析的原始研究數據進行了分析。我們重新進行數據提取，其中 4 個原始研究的基本特征如表 1。
表1
部分納入研究的基本特征
1.3 數據整理
按照劑量-反應 Meta 分析模型的格式要求，數據整理見表 2。由于原始研究給出的睡眠的時間多為區間（如 8～9 h），而不是確切的點估計數值，為了方便進行運算，將數值區間的中位數（如：8.5）當作“平均”睡眠時間，對于開放性區間（如>9 h），可用相鄰區間的寬度作為當前開放區間的寬度，再取中位數作為“平均”劑量。也即本例中“>9 h”這個開放性區間的另一端可估算為 10 h，因為其相鄰區間（8～9 h）的寬度為 1。另一種估算辦法是將 9 h 乘以 1.2 或 1.5 當作開放區間的另一端，然后取估算的區間的中位數當作“平均”劑量。值得注意的是，本例最短的睡眠時間是 4.5 h，并且每個納入研究參照劑量水平并非為最短的睡眠時間，因此我們需對數據進行中心化（通常做法為將所有劑量減去均值），并將參照劑量水平轉化為最低的睡眠時間。
表2
數據整理后的納入研究情況
2 兩種模型的應用
2.1 基于普通 GLST 模型的應用
GLST 模型由 Orsini N 教授提出，是目前使用最廣泛的劑量-反應 Meta 分析方法[6]。該方法在 Stata 軟件的使用需要安裝兩個程序：glst、xblc。安裝命令如下：
ssc install glst
findit xblc
方法一：將以上整理后的數據輸入或復制到 Stata 軟件數據窗口，并使用以下命令（由 Suhail A.R Doi 教授和徐暢博士提供，可直接復制以下代碼并使用）進行劑量-反應 Meta 分析。
對每篇研究劑量數據進行中心化并計算權重：
bysort id：gen dosec=dose–dose[1]
gen wt=1/（se^2）
使用限制性立方樣條函數插值：
mkspline dosecs=dosec，cubic nk（3）disp
隨機效應模型下非線性劑量-反應 Meta 分析：
glst logrr dosecs*，cov（n case）se（se）pfirst（id
studyt）r eform
predict phatgr，xb
產生預測置信區間：
predict segr，stdp
gen ucigr=phatgr+（1.96*segr）
gen lcigr=phatgr–（1.96*segr）
隨機效應模型下分段線性劑量-反應 Meta 分析：
lincom（（（6–2）*dosecs1）+（（5.25–0.32）*
dosecs2））/（6–2），eform
lincom（（（2–0）*dosecs1）+（（0.32–0）*
dosecs2））/（–2），eform
非線性檢驗：
test dosecs1=dosecs2
圖形繪制：
quietly levelsof dosec，local（levels）
xblc dosecs*，covname（ dosec ）at（‘r（levels）’）
ref （2）eform line
本命令中使用兩個回歸系數的平均值計算線性趨勢（參見命令中的*提示），并且該代碼可直接計算分段線性趨勢，該部分第一行命令為睡眠時間超過 7 h 的分段，第二行命令為小于 7 h 的分段。分段線性結果為：睡眠時間在 7 h 及以上的人群，每增加 1 h 的睡眠時間，全因死亡相對風險增加 10%，95%CI（8%，11%）。本例非線性檢驗的P 值為 0.00，拒絕睡眠與全因死亡率為線性相關的零假設，提示本例中更適合使用非線性關系。根據上述代碼可獲取非線性劑量-反應關系圖，見圖 1。
圖1
睡眠時間與全因死亡率的非線性劑量-反應關系
將 glst 所在行的命令，后面的 r 改為 f，即可計算固定效應模型下的結果，r 和 f 是 random 及 fixed 的簡寫，分別表示隨機和固定效應模型，此處不再贅述。
實際上，命令的最后一行，xblc 命令，會在 Stata 結果輸出窗口（圖 2）列出不同劑量對應的相對危險度。其意義表示，相對于每天睡眠 7 h，5 h 的全因死亡率相對危險度為 1.12，95%CI（1.08，1.16），以此類推。值得注意的是，這里提到的 7 h，是需要將（圖 2）顯示的時間加上 5 h，因為在我們分析數據之前，對數據進行了中心化，將每個納入研究的睡眠時間平均減去了 5 h。
圖2
睡眠時間的 Stata 軟件顯示
2.2 基于 Mvmeta 模型的應用
使用 Mvmeta 代碼，程序運算思路是首先使用 glst 命令計算單個研究的回歸參數，然后采用 Mvmeta 將這些回歸參數按照固定或隨機效應模型進行合并。與上一方法的不同之處在于，該代碼考慮了研究內及研究間回歸參數的相關性，并對相關性進行了校正。非線性模型代碼如下：
中心化并計算權重：
bysort id：gen dosec=dose–dose[1] if id !=9 & id
!=12 & id !=13 & id !=14 & id !=27
gen wt=1/（se^2）
使用限制性立方樣條函數插值：
mkspline dosecs=dosec，cubic nk（3）disp
mkspline doses=dose，nk（3）cubic
計算單篇研究的回歸系數并將數據轉換為 mvmeta 命令所需的結構：
mvmeta_make glst logrr dosecs1 dosecs2，cov
（n case）se（se）pfirst（id studyt）
saving（ssest_mks）replace by（id）names（b V）
隨機效應模型劑量-反應 Meta 分析：
preserve
use ssest_mks，clear
mvmeta b V，reml i2 nc
非線性檢驗：
test bdosecs1=bdosecs2
預測非線性結果：
capture estimates save streg，replace
restore
estimates use streg
predictnl rr_pl=exp（_b[bdosecs1]*（doses1-5）+
_b[bdosecs2]*（doses2-0）），ci（lci uci）
繪制劑量-反應 Meta 分析曲線：
twoway （line rr_pl dose，sort lc（black）lp（–））
（scatter rr dose [aweight=1/se]，
msymbol（Oh）），scheme（s1mono）
ytitle（“Relative Risk”，margin（right））xtitle
（“Sleep duration，（h）”，margin（top_bottom））
通過上述命令，可獲得如下結果，見圖 3。對比圖 1 我們可以發現，該結果跟第一種方法有一定區別，原因在于這種方法是使用了多變量隨機效應 Meta 回歸命令 mvmeta，對回歸系數間的相關性進行了調整。但第一種命令能更方便地獲取線性及分段線性的結果。
圖3
睡眠時間與全因死亡率的劑量-反應曲線（Mvmeta 軟件結果）
3 討論
本文對當前兩種不同實現隨機效應模型下 DRMA 的方法進行了介紹，并采用睡眠于全因死亡率的數據展示了兩種方法在 Stata 軟件下的實現。我們發現，兩種不同方法計算的劑量-效應結果實際上存在較大差異，盡管結果均提示睡眠時間過長或過短均會增加全因死亡率的風險。
這兩種實現 DRMA 的隨機效應模型方法的相同點在于均基于二階段模型理論，并且均考慮了研究內劑量-效應之間的相關性。不同點在于，第一種方法未考慮參數間的相關性，即單個研究內多參數（回歸系數）之間的相關性，而第二種方法采用 Mvmeta 框架將這種相關性考慮在內。實際上，在隨機模型下，一個基本的前提是假設單個研究的效應量服從正態分布。盡管被廣泛使用，但這種假設一直未被證實，所得的結果應批判性對待。
近年來，國內劑量反應 Meta 分析的發展較為迅速，中國作者發表了大量相關方法學、軟件應用及循證實踐的學術論文。據統計，自 2011 年至 2015 年，我國學者公開發表了約 256 篇劑量-反應 Meta 分析，約占此期間國際范圍內劑量-反應 Meta 分析總數的 56%[2]。盡管如此，這些國內作者的劑量-反應 Meta 分析報告質量參差不齊，需要引起足夠重視。
總之，實現劑量-反應 Meta 分析的方法較多，但不同方法存在各自的優劣性。且不同方法的適用范圍并非完全一致。劑量-反應 Meta 分析是一種“年輕”的統計學方法，仍存在許多缺陷，相關方法仍需后續進一步研究并完善。
流行病學研究的重要任務是探討或建立暴露/干預與疾病結局的因果關系（cause-and-effect relationship）[1]。探討不同水平下的暴露/干預與結局指標的動態關系稱為劑量-反應關系（dose-response relationship）；而將針對相同臨床或公共衛生問題的不同研究的這種劑量-反應關系加權合并以獲得“平均”劑量效應的方法稱為劑量-反應 Meta 分析（dose-response meta-analysis，DRMA）[2]。
當前使用最為廣泛的 DRMA 模型是基于經典的二階段法（two-stage approach）下實現的，即首先通過線性或非線性模型擬合出單篇研究的回歸系數（第一階段），再采用固定或隨機效應模型將第一階段的回歸系數進行合并（第二階段）[3,4]。固定效應模型的 DRMA 因無需考慮研究間異質性，在模型開發及實現上較為簡單；而對于隨機效應模型的 DRMA，特別是存在多個參數時（即非線性模型），軟件程序開發及實現較為困難。
本文將介紹兩種基于 Stata 軟件、實現多參數隨機效應模型的 DRMA 方法，并代入睡眠與全因死亡率的數據分別進行 DRMA，以期為循證實踐提供相關參考。
1 資料來源
1.1 軟件代碼來源
本文中涉及的部分軟件代碼由 Suhail A.R Doi 教授和徐暢博士提供，并已獲得他們同意。
1.2 數據來源和基本特征
以 Liu 等[5]2017 年發表的關于睡眠時間與全因死亡率的劑量-反應 Meta 分析為例，采用兩種不同方法重新進行 DRMA。該研究從 51 篇關于睡眠時間與全因死亡率的前瞻性隊列研究中，納入 40 篇符合劑量-反應 Meta 分析的原始研究數據進行了分析。我們重新進行數據提取，其中 4 個原始研究的基本特征如表 1。
表1
部分納入研究的基本特征
1.3 數據整理
按照劑量-反應 Meta 分析模型的格式要求，數據整理見表 2。由于原始研究給出的睡眠的時間多為區間（如 8～9 h），而不是確切的點估計數值，為了方便進行運算，將數值區間的中位數（如：8.5）當作“平均”睡眠時間，對于開放性區間（如>9 h），可用相鄰區間的寬度作為當前開放區間的寬度，再取中位數作為“平均”劑量。也即本例中“>9 h”這個開放性區間的另一端可估算為 10 h，因為其相鄰區間（8～9 h）的寬度為 1。另一種估算辦法是將 9 h 乘以 1.2 或 1.5 當作開放區間的另一端，然后取估算的區間的中位數當作“平均”劑量。值得注意的是，本例最短的睡眠時間是 4.5 h，并且每個納入研究參照劑量水平并非為最短的睡眠時間，因此我們需對數據進行中心化（通常做法為將所有劑量減去均值），并將參照劑量水平轉化為最低的睡眠時間。
表2
數據整理后的納入研究情況
2 兩種模型的應用
2.1 基于普通 GLST 模型的應用
GLST 模型由 Orsini N 教授提出，是目前使用最廣泛的劑量-反應 Meta 分析方法[6]。該方法在 Stata 軟件的使用需要安裝兩個程序：glst、xblc。安裝命令如下：
ssc install glst
findit xblc
方法一：將以上整理后的數據輸入或復制到 Stata 軟件數據窗口，并使用以下命令（由 Suhail A.R Doi 教授和徐暢博士提供，可直接復制以下代碼并使用）進行劑量-反應 Meta 分析。
對每篇研究劑量數據進行中心化并計算權重：
bysort id：gen dosec=dose–dose[1]
gen wt=1/（se^2）
使用限制性立方樣條函數插值：
mkspline dosecs=dosec，cubic nk（3）disp
隨機效應模型下非線性劑量-反應 Meta 分析：
glst logrr dosecs*，cov（n case）se（se）pfirst（id
studyt）r eform
predict phatgr，xb
產生預測置信區間：
predict segr，stdp
gen ucigr=phatgr+（1.96*segr）
gen lcigr=phatgr–（1.96*segr）
隨機效應模型下分段線性劑量-反應 Meta 分析：
lincom（（（6–2）*dosecs1）+（（5.25–0.32）*
dosecs2））/（6–2），eform
lincom（（（2–0）*dosecs1）+（（0.32–0）*
dosecs2））/（–2），eform
非線性檢驗：
test dosecs1=dosecs2
圖形繪制：
quietly levelsof dosec，local（levels）
xblc dosecs*，covname（ dosec ）at（‘r（levels）’）
ref （2）eform line
本命令中使用兩個回歸系數的平均值計算線性趨勢（參見命令中的*提示），并且該代碼可直接計算分段線性趨勢，該部分第一行命令為睡眠時間超過 7 h 的分段，第二行命令為小于 7 h 的分段。分段線性結果為：睡眠時間在 7 h 及以上的人群，每增加 1 h 的睡眠時間，全因死亡相對風險增加 10%，95%CI（8%，11%）。本例非線性檢驗的P 值為 0.00，拒絕睡眠與全因死亡率為線性相關的零假設，提示本例中更適合使用非線性關系。根據上述代碼可獲取非線性劑量-反應關系圖，見圖 1。
圖1
睡眠時間與全因死亡率的非線性劑量-反應關系
將 glst 所在行的命令，后面的 r 改為 f，即可計算固定效應模型下的結果，r 和 f 是 random 及 fixed 的簡寫，分別表示隨機和固定效應模型，此處不再贅述。
實際上，命令的最后一行，xblc 命令，會在 Stata 結果輸出窗口（圖 2）列出不同劑量對應的相對危險度。其意義表示，相對于每天睡眠 7 h，5 h 的全因死亡率相對危險度為 1.12，95%CI（1.08，1.16），以此類推。值得注意的是，這里提到的 7 h，是需要將（圖 2）顯示的時間加上 5 h，因為在我們分析數據之前，對數據進行了中心化，將每個納入研究的睡眠時間平均減去了 5 h。
圖2
睡眠時間的 Stata 軟件顯示
2.2 基于 Mvmeta 模型的應用
使用 Mvmeta 代碼，程序運算思路是首先使用 glst 命令計算單個研究的回歸參數，然后采用 Mvmeta 將這些回歸參數按照固定或隨機效應模型進行合并。與上一方法的不同之處在于，該代碼考慮了研究內及研究間回歸參數的相關性，并對相關性進行了校正。非線性模型代碼如下：
中心化并計算權重：
bysort id：gen dosec=dose–dose[1] if id !=9 & id
!=12 & id !=13 & id !=14 & id !=27
gen wt=1/（se^2）
使用限制性立方樣條函數插值：
mkspline dosecs=dosec，cubic nk（3）disp
mkspline doses=dose，nk（3）cubic
計算單篇研究的回歸系數并將數據轉換為 mvmeta 命令所需的結構：
mvmeta_make glst logrr dosecs1 dosecs2，cov
（n case）se（se）pfirst（id studyt）
saving（ssest_mks）replace by（id）names（b V）
隨機效應模型劑量-反應 Meta 分析：
preserve
use ssest_mks，clear
mvmeta b V，reml i2 nc
非線性檢驗：
test bdosecs1=bdosecs2
預測非線性結果：
capture estimates save streg，replace
restore
estimates use streg
predictnl rr_pl=exp（_b[bdosecs1]*（doses1-5）+
_b[bdosecs2]*（doses2-0）），ci（lci uci）
繪制劑量-反應 Meta 分析曲線：
twoway （line rr_pl dose，sort lc（black）lp（–））
（scatter rr dose [aweight=1/se]，
msymbol（Oh）），scheme（s1mono）
ytitle（“Relative Risk”，margin（right））xtitle
（“Sleep duration，（h）”，margin（top_bottom））
通過上述命令，可獲得如下結果，見圖 3。對比圖 1 我們可以發現，該結果跟第一種方法有一定區別，原因在于這種方法是使用了多變量隨機效應 Meta 回歸命令 mvmeta，對回歸系數間的相關性進行了調整。但第一種命令能更方便地獲取線性及分段線性的結果。
圖3
睡眠時間與全因死亡率的劑量-反應曲線（Mvmeta 軟件結果）
3 討論
本文對當前兩種不同實現隨機效應模型下 DRMA 的方法進行了介紹，并采用睡眠于全因死亡率的數據展示了兩種方法在 Stata 軟件下的實現。我們發現，兩種不同方法計算的劑量-效應結果實際上存在較大差異，盡管結果均提示睡眠時間過長或過短均會增加全因死亡率的風險。
這兩種實現 DRMA 的隨機效應模型方法的相同點在于均基于二階段模型理論，并且均考慮了研究內劑量-效應之間的相關性。不同點在于，第一種方法未考慮參數間的相關性，即單個研究內多參數（回歸系數）之間的相關性，而第二種方法采用 Mvmeta 框架將這種相關性考慮在內。實際上，在隨機模型下，一個基本的前提是假設單個研究的效應量服從正態分布。盡管被廣泛使用，但這種假設一直未被證實，所得的結果應批判性對待。
近年來，國內劑量反應 Meta 分析的發展較為迅速，中國作者發表了大量相關方法學、軟件應用及循證實踐的學術論文。據統計，自 2011 年至 2015 年，我國學者公開發表了約 256 篇劑量-反應 Meta 分析，約占此期間國際范圍內劑量-反應 Meta 分析總數的 56%[2]。盡管如此，這些國內作者的劑量-反應 Meta 分析報告質量參差不齊，需要引起足夠重視。
總之，實現劑量-反應 Meta 分析的方法較多，但不同方法存在各自的優劣性。且不同方法的適用范圍并非完全一致。劑量-反應 Meta 分析是一種“年輕”的統計學方法，仍存在許多缺陷，相關方法仍需后續進一步研究并完善。
