劑量-反應關系Meta分析在探討潛在風險因素暴露水平與對應疾病發生風險方面具有重要作用,且應用日益廣泛。隨著各種函數模型的不斷優化,劑量-反應關系函數模型的建立逐漸成熟,已由單純的線性關系發展成集線性及非線性關系于一體的函數模型,由單純的一次回歸函數發展為多次回歸函數。基于R軟件實現這一分析的軟件程序也逐漸研發出來,本文介紹了如何使用R軟件mvmeta及dosresmeta進行劑量-反應關系Meta分析。
引用本文： 徐暢, 曾憲濤, 張超, 李勝, 張永剛, 劉同族. 應用R軟件dosresmeta程序包和mvmeta程序包實現劑量-反應關系Meta分析. 中國循證醫學雜志, 2015, 15(4): 479-483. doi: 10.7507/1672-2531.20150079 復制
劑量-反應關系（dose-response relations）是指對應結局指標隨著暴露因素/水平的增加或降低而對應產生變化的趨勢，基于劑量-反應關系數據的Meta分析正日益得到廣泛應用 [1, 2]。隨著數學方法學理論及編程手段的逐步深入，已有越來越多的模型及軟件可用于實現劑量-反應關系的Meta分析 [3]。R軟件進行劑量-反應Meta分析目前較常用的方法有基于二次函數模型下的程序及限制性立方樣條函數下的程序這兩種。2013年發布的dosresmeta程序包為首款用于劑量-反應分析的程序包 [4]，該程序包需要使用其他方法對數據進行處理；由White及其同事研發的mvmeta程序包是針對多元Meta分析的一種數據分析程序 [5, 6]。mvmeta程序包可合并多個研究的多種參數，為劑量-反應關系Meta分析模型提供了很好的平臺：不僅可進行線性劑量-反應關系中一個回歸參數的合并，也可用于非線性劑量-反應關系中多個回歸參數的合并 [5, 6]。本文介紹如何使用R軟件dosresmeta程序包和mvmeta程序包進行劑量-反應關系Meta分析。
1 方法與材料
1.1 程序包的準備
R軟件實現劑量-反應Meta分析需先行加載mvmeta、Hmisc、rms和dosresmeta這4個程序包。此外，還可加載專用于制圖的ggplot2程序包。R軟件的安裝及相關軟件包的安裝、加載方法請參閱《R軟件Metafor程序包在Meta分析中的應用》 [7]一文。
1.2 二次函數回歸模型
二次函數回歸模型是基于最高次數為二次的多項式函數發展而來，該模型既包含傳統線性-劑量反應關系模型（即將二次項參數設為零），也包含非線性-劑量反應模型，于2009年被Liu等首次應用于制作劑量-反應關系Meta分析 [8]。
該模型函數表達式如下：
| $In\lambda \left( x,z \right)=a+{{b}_{1}}x+{{b}_{2}}{{x}^{2}}+\theta 'z$ |
其中x表示暴露劑量（例如，脂肪攝入量），Inλ為對應的Risk ratio（RR）或Odds ratio（OR）或Hazard risk（HR）的自然對數，b1、b2分別為x及x2的參數，θ'z為協方差。公式（1）為單篇研究的劑量-反應關系表達式，對于Meta分析則需對公式進行變換。方法為：假設i為研究數，x0為參照劑量，xj為第j分層的暴露劑量，參照劑量對應的RR=1，對數轉換后為0，代入（1）則公式可變換為公式（2）：
| $InR{{R}_{ij}}={{b}_{1}}i({{x}_{ij}}-{{x}_{i0}})\text{ }+{{b}_{2i}}({{x}_{ij}}^{2}~-{{x}_{i0}}^{2}~)\text{ }+{{\xi }_{ij}}~$ |
其中ξij表示研究內及研究間協方差矩陣。該函數擬合的趨勢噪音（波動）較小，但對局部的變動不靈敏 [9]。
1.3 限制性立方樣條函數模型
限制性立方樣條函數由樣條函數演變而來，Orsini等于2012年首先發布了應用于劑量-反應Meta分析的程序包 [10]。該函數擬合效果介于參數法及非參數法之間，函數擬合較為接近真實趨勢，但局部波動較大，可能存在過度擬合 [9]。
該模型函數可表達如下：
| $log\left( RR|x \right)={{b}_{1}}x+{{b}_{2}}{{x}^{2}}+{{b}_{3}}{{x}^{3}}+\sum {{~}_{i=1}}^{c-1}{{b}_{i3}}(x-{{k}_{i}}){{~}_{+}}^{3}$ |
前面幾項函數意義同公式（1），此外c表示將暴露劑量x分為c個分層及c-1個節點，記為k1…ki，(x-ki)+為函數的陽性部分，其中：
| $\begin{array}{*{35}{l}} {{\left( x-{{k}_{i}} \right)}_{+}}=0Ifx\text{ }{{k}_{i}} \\ {{\left( x-{{k}_{i}} \right)}_{+}}=x-{{k}_{i}}Ifx{{k}_{i}} \\ \end{array}$ |
此模型亦已嵌入線性劑量-反應關系模型，即設二次和三次項為零。限制性立方樣條函數中的“限制”原意是將原立方樣條函數中第一個節點以前、第三個節點以后的趨勢強制性設置為線性，以避免出現龍格現象（Runge's phenomenon）[11]。
1.4 參數估計
參數估計一般分為最小二乘法（least square method）、廣義最小二乘法（generalized least squares method）、最大似然估計法（maximum likelihood estimate）和貝葉斯法（Bayes method）[12]，目前使用最廣泛的仍然是廣義最小二乘法及最大似然估計法，該法對目標參數進行求導或微分積分，然后令變換后的表達式為零，其解即為目標結果。
2 軟件操作及數據分析
2.1 數據資料
本文數據來自于《應用STATA做Meta分析》書中所引用的實例（表 1）[2]，數據導入方法請參閱前期文章 [7, 13]，導入數據格式如表 1所示。本文將表 1中數據轉為“.txt”格式并用read.table命令導入，命令如下:
data＜-read.table("D:\Rdata\data.txt",header= TRUE)
data
導入后，計算adjrr對數及相應se；也可先行計算后直接放入數據表格中再導入。計算命令如下：
logrr＞-log(adjrr)
se＞-(lb-ub)/ (2*invnormal(0.975) )
需要注意的是：計算后參照劑量的se改為“NA”，而不是“0”。
表1
導入數據格式
2.2 二次函數回歸模型的劑量-反應Meta分析
二次函數回歸模型強制設定合并模型為隨機效應模型（random effect model），且使用的參數估計方法為最大似然估計法。導入數據后，輸入如下命令以執行數據分析：
quadr＜- dosresmeta(formula = logrr ~ dose + I(dose^2) ,id = id,
type = type,se = se,cases = case,n = n,
data = data,method = "mm")
summary(quadr)
上述命令中quadr是二次函數（quadratic）的縮寫；“mm”為參數估計方法，即最大似然估計法。用如下命令執行整體趨勢預測、曲線繪制、及結果計算：
newdata＜- data.frame(dose = seq(0,7,1) )
prediction＜- predict(quadr,newdata,expo = T)
prediction$ci.ub[prediction$ci.ub＞=8]＜- 8
prediction$ci.lb[prediction$ci.lb＜=.3]＜- .3
ggplot(prediction,aes(dose,pred)) + geom_line() + theme_classic() +geom_ribbon(aes(ymin = ci.lb,ymax = i.ub),alpha = .05) +scale_y_continuous("RR of ovarian cancer",trans="log",breaks = c(1,2,3,4,),limits = (c(1,4) )) +scale_x_continuous(",grams/day",breaks = seq(0,7,1) )
newdata＜- data.frame(dose = seq(0,7,1) )
predict(quadr,newdata,expo = T)
上述命令中seq設定劑量范圍及增加量，如seq（0，7，1）表示劑量范圍為0至7，每次增加1個單位；breaks=c是用于設定Y軸的顯示范圍，limits則限制Y軸的值。
本例運行結果見表 2。
表2
運行結果
本例單篇研究的劑量趨勢圖與最終合并的劑量趨勢圖如下分別見圖 1和圖 2。
圖1
單篇研究的劑量圖匯總
圖2
合并后的劑量趨勢圖
上述結果表明，雖然卵巢癌發病風險隨著劑量的增加而增加，但95%置信區間均跨過1，提示無統計學意義。更多細節方面的研究，如計算單篇研究下一次及二次項參數、將研究類型等可能影響結果的因素當作協變量的計算方法，請參閱Liu等 [8]的研究。
2.3 限制性立方樣條模型的劑量-反應Meta分析
限制性立方樣條模型提供了固定及隨機效應兩種合并模型，使用的參數估計方法為廣義最小二乘法。本文為與二次函數法保持一致并比較結果，僅采用隨機效應模型進行展示。數據導入方法同上述，本處直接進入數據分析及圖形繪制步驟。具體如下：
① 線性關系的命令如下：
lin＜- dosresmeta(formula = logrr ~ dose,type = type,id = id,se = se,cases = case,n = n,data = data)
summary(lin)
predict(lin,delta = 12)
② 非線性關系的命令如下：
knots＜- quantile(data$dose,c(.1,.5,.9) )
spl＜- dosresmeta(formula = logrr ~ rcs(dose,knots),type = type,id = id,se = se,cases = case,n = n,data = data)
summary(spl)
其中knots設置節點數；c(0.1,0.5,0.9) 表示在10%、50%、90%的dose分布上各取一節點，若原始數據分層最大為3層（當然最小必須為3層），則最多只能取3個節點。
③ 非線性檢驗的作用是判斷趨勢是否為非線性，本例中假設第二條樣條系數為0，命令如下：
wald.test(b = coef(spl),Sigma = vcov(spl),Terms = 1:2)
命令執行后得出結果為：
Chi-squared test:
X2 = 1.3,df = 2,P(＞X2) = 0.53
結果表明接受零假設，即本例數據乳糖攝入與卵巢癌為線性關系，因此需要使用上述線性關系結果。本例線性劑量-反應關系結果為每天每增加1 mL乳糖攝入量，對應卵巢癌風險增量為1.02[95%CI（0.94，1.10）]，即兩者并無明顯關聯性。
④ 非線性劑量-反應結果
返回各暴露劑量對應卵巢癌發病風險值，命令如下：
pred＜- predict(spl,data.frame(dose = seq(0,7,1) ),xref = 0)
print(pred,digits = 2)
命令執行后本例結果如下（表 3）：
表3
執行命令所得結果
因本例數據呈線性關系，因此建議先行將每篇文章數據單獨進行線性劑量反應Meta分析，然后將結果用常規Meta分析程序 [2, 3, 7]繪制森林圖。
此處仍給出非線性關系劑量趨勢圖（圖 3）。
圖3
非線性關系劑量趨勢圖
3 討論
從結果部分來看，上述兩種方法結果在乳糖攝入量小于5 mL/d時較為接近，當超過該值時，結果相差很大，但最終得出的結論基本一致。究其原因，或許與兩種模型的特性有關。二次函數模型的CI較限制性立方樣條函數寬，并且并未考慮趨勢是否為線性，加之函數本身擬合效果缺陷，在實際應用中并不多見。非線性檢驗表明該趨勢為線性相關，提示線性函數模型可能為最佳擬合模型。
除了上述兩種方法外，還有多種曲線逼近的方法，如二次B樣條、分段多項式等方法。二次B樣條屬于非參數回歸，逼近效果較一般參數回歸好，但因過度擬合等問題只被用于評估新函數模型的擬合效果；分段多項式法也是一種較好的方法，與限制性立方樣條回歸為目前兩大最常用劑量-反應Meta分析的方法。這兩者擬合效果相近，目前公開用于分段多項式法分析程序僅限于SAS軟件 [2, 14-16]，這將會在本團隊后期文章中進行詳細介紹。
當然，每種曲線逼近（擬合）模型均有其優缺點，比如線性或二次函數模型能較好的反應整體變化趨勢，但對于局部細微的變化的察覺靈敏度較低；高次（三次及以上）函數模型對局部變化反應較真實，但可能由于局部變化并不真實等而存在過度擬合問題，導致結果產生偏倚。因此使用者需要對每種方法進行對比并結合臨床背景，選出最適合的模型。
劑量-反應關系（dose-response relations）是指對應結局指標隨著暴露因素/水平的增加或降低而對應產生變化的趨勢，基于劑量-反應關系數據的Meta分析正日益得到廣泛應用 [1, 2]。隨著數學方法學理論及編程手段的逐步深入，已有越來越多的模型及軟件可用于實現劑量-反應關系的Meta分析 [3]。R軟件進行劑量-反應Meta分析目前較常用的方法有基于二次函數模型下的程序及限制性立方樣條函數下的程序這兩種。2013年發布的dosresmeta程序包為首款用于劑量-反應分析的程序包 [4]，該程序包需要使用其他方法對數據進行處理；由White及其同事研發的mvmeta程序包是針對多元Meta分析的一種數據分析程序 [5, 6]。mvmeta程序包可合并多個研究的多種參數，為劑量-反應關系Meta分析模型提供了很好的平臺：不僅可進行線性劑量-反應關系中一個回歸參數的合并，也可用于非線性劑量-反應關系中多個回歸參數的合并 [5, 6]。本文介紹如何使用R軟件dosresmeta程序包和mvmeta程序包進行劑量-反應關系Meta分析。
1 方法與材料
1.1 程序包的準備
R軟件實現劑量-反應Meta分析需先行加載mvmeta、Hmisc、rms和dosresmeta這4個程序包。此外，還可加載專用于制圖的ggplot2程序包。R軟件的安裝及相關軟件包的安裝、加載方法請參閱《R軟件Metafor程序包在Meta分析中的應用》 [7]一文。
1.2 二次函數回歸模型
二次函數回歸模型是基于最高次數為二次的多項式函數發展而來，該模型既包含傳統線性-劑量反應關系模型（即將二次項參數設為零），也包含非線性-劑量反應模型，于2009年被Liu等首次應用于制作劑量-反應關系Meta分析 [8]。
該模型函數表達式如下：
| $In\lambda \left( x,z \right)=a+{{b}_{1}}x+{{b}_{2}}{{x}^{2}}+\theta 'z$ |
其中x表示暴露劑量（例如，脂肪攝入量），Inλ為對應的Risk ratio（RR）或Odds ratio（OR）或Hazard risk（HR）的自然對數，b1、b2分別為x及x2的參數，θ'z為協方差。公式（1）為單篇研究的劑量-反應關系表達式，對于Meta分析則需對公式進行變換。方法為：假設i為研究數，x0為參照劑量，xj為第j分層的暴露劑量，參照劑量對應的RR=1，對數轉換后為0，代入（1）則公式可變換為公式（2）：
| $InR{{R}_{ij}}={{b}_{1}}i({{x}_{ij}}-{{x}_{i0}})\text{ }+{{b}_{2i}}({{x}_{ij}}^{2}~-{{x}_{i0}}^{2}~)\text{ }+{{\xi }_{ij}}~$ |
其中ξij表示研究內及研究間協方差矩陣。該函數擬合的趨勢噪音（波動）較小，但對局部的變動不靈敏 [9]。
1.3 限制性立方樣條函數模型
限制性立方樣條函數由樣條函數演變而來，Orsini等于2012年首先發布了應用于劑量-反應Meta分析的程序包 [10]。該函數擬合效果介于參數法及非參數法之間，函數擬合較為接近真實趨勢，但局部波動較大，可能存在過度擬合 [9]。
該模型函數可表達如下：
| $log\left( RR|x \right)={{b}_{1}}x+{{b}_{2}}{{x}^{2}}+{{b}_{3}}{{x}^{3}}+\sum {{~}_{i=1}}^{c-1}{{b}_{i3}}(x-{{k}_{i}}){{~}_{+}}^{3}$ |
前面幾項函數意義同公式（1），此外c表示將暴露劑量x分為c個分層及c-1個節點，記為k1…ki，(x-ki)+為函數的陽性部分，其中：
| $\begin{array}{*{35}{l}} {{\left( x-{{k}_{i}} \right)}_{+}}=0Ifx\text{ }{{k}_{i}} \\ {{\left( x-{{k}_{i}} \right)}_{+}}=x-{{k}_{i}}Ifx{{k}_{i}} \\ \end{array}$ |
此模型亦已嵌入線性劑量-反應關系模型，即設二次和三次項為零。限制性立方樣條函數中的“限制”原意是將原立方樣條函數中第一個節點以前、第三個節點以后的趨勢強制性設置為線性，以避免出現龍格現象（Runge's phenomenon）[11]。
1.4 參數估計
參數估計一般分為最小二乘法（least square method）、廣義最小二乘法（generalized least squares method）、最大似然估計法（maximum likelihood estimate）和貝葉斯法（Bayes method）[12]，目前使用最廣泛的仍然是廣義最小二乘法及最大似然估計法，該法對目標參數進行求導或微分積分，然后令變換后的表達式為零，其解即為目標結果。
2 軟件操作及數據分析
2.1 數據資料
本文數據來自于《應用STATA做Meta分析》書中所引用的實例（表 1）[2]，數據導入方法請參閱前期文章 [7, 13]，導入數據格式如表 1所示。本文將表 1中數據轉為“.txt”格式并用read.table命令導入，命令如下:
data＜-read.table("D:\Rdata\data.txt",header= TRUE)
data
導入后，計算adjrr對數及相應se；也可先行計算后直接放入數據表格中再導入。計算命令如下：
logrr＞-log(adjrr)
se＞-(lb-ub)/ (2*invnormal(0.975) )
需要注意的是：計算后參照劑量的se改為“NA”，而不是“0”。
表1
導入數據格式
2.2 二次函數回歸模型的劑量-反應Meta分析
二次函數回歸模型強制設定合并模型為隨機效應模型（random effect model），且使用的參數估計方法為最大似然估計法。導入數據后，輸入如下命令以執行數據分析：
quadr＜- dosresmeta(formula = logrr ~ dose + I(dose^2) ,id = id,
type = type,se = se,cases = case,n = n,
data = data,method = "mm")
summary(quadr)
上述命令中quadr是二次函數（quadratic）的縮寫；“mm”為參數估計方法，即最大似然估計法。用如下命令執行整體趨勢預測、曲線繪制、及結果計算：
newdata＜- data.frame(dose = seq(0,7,1) )
prediction＜- predict(quadr,newdata,expo = T)
prediction$ci.ub[prediction$ci.ub＞=8]＜- 8
prediction$ci.lb[prediction$ci.lb＜=.3]＜- .3
ggplot(prediction,aes(dose,pred)) + geom_line() + theme_classic() +geom_ribbon(aes(ymin = ci.lb,ymax = i.ub),alpha = .05) +scale_y_continuous("RR of ovarian cancer",trans="log",breaks = c(1,2,3,4,),limits = (c(1,4) )) +scale_x_continuous(",grams/day",breaks = seq(0,7,1) )
newdata＜- data.frame(dose = seq(0,7,1) )
predict(quadr,newdata,expo = T)
上述命令中seq設定劑量范圍及增加量，如seq（0，7，1）表示劑量范圍為0至7，每次增加1個單位；breaks=c是用于設定Y軸的顯示范圍，limits則限制Y軸的值。
本例運行結果見表 2。
表2
運行結果
本例單篇研究的劑量趨勢圖與最終合并的劑量趨勢圖如下分別見圖 1和圖 2。
圖1
單篇研究的劑量圖匯總
圖2
合并后的劑量趨勢圖
上述結果表明，雖然卵巢癌發病風險隨著劑量的增加而增加，但95%置信區間均跨過1，提示無統計學意義。更多細節方面的研究，如計算單篇研究下一次及二次項參數、將研究類型等可能影響結果的因素當作協變量的計算方法，請參閱Liu等 [8]的研究。
2.3 限制性立方樣條模型的劑量-反應Meta分析
限制性立方樣條模型提供了固定及隨機效應兩種合并模型，使用的參數估計方法為廣義最小二乘法。本文為與二次函數法保持一致并比較結果，僅采用隨機效應模型進行展示。數據導入方法同上述，本處直接進入數據分析及圖形繪制步驟。具體如下：
① 線性關系的命令如下：
lin＜- dosresmeta(formula = logrr ~ dose,type = type,id = id,se = se,cases = case,n = n,data = data)
summary(lin)
predict(lin,delta = 12)
② 非線性關系的命令如下：
knots＜- quantile(data$dose,c(.1,.5,.9) )
spl＜- dosresmeta(formula = logrr ~ rcs(dose,knots),type = type,id = id,se = se,cases = case,n = n,data = data)
summary(spl)
其中knots設置節點數；c(0.1,0.5,0.9) 表示在10%、50%、90%的dose分布上各取一節點，若原始數據分層最大為3層（當然最小必須為3層），則最多只能取3個節點。
③ 非線性檢驗的作用是判斷趨勢是否為非線性，本例中假設第二條樣條系數為0，命令如下：
wald.test(b = coef(spl),Sigma = vcov(spl),Terms = 1:2)
命令執行后得出結果為：
Chi-squared test:
X2 = 1.3,df = 2,P(＞X2) = 0.53
結果表明接受零假設，即本例數據乳糖攝入與卵巢癌為線性關系，因此需要使用上述線性關系結果。本例線性劑量-反應關系結果為每天每增加1 mL乳糖攝入量，對應卵巢癌風險增量為1.02[95%CI（0.94，1.10）]，即兩者并無明顯關聯性。
④ 非線性劑量-反應結果
返回各暴露劑量對應卵巢癌發病風險值，命令如下：
pred＜- predict(spl,data.frame(dose = seq(0,7,1) ),xref = 0)
print(pred,digits = 2)
命令執行后本例結果如下（表 3）：
表3
執行命令所得結果
因本例數據呈線性關系，因此建議先行將每篇文章數據單獨進行線性劑量反應Meta分析，然后將結果用常規Meta分析程序 [2, 3, 7]繪制森林圖。
此處仍給出非線性關系劑量趨勢圖（圖 3）。
圖3
非線性關系劑量趨勢圖
3 討論
從結果部分來看，上述兩種方法結果在乳糖攝入量小于5 mL/d時較為接近，當超過該值時，結果相差很大，但最終得出的結論基本一致。究其原因，或許與兩種模型的特性有關。二次函數模型的CI較限制性立方樣條函數寬，并且并未考慮趨勢是否為線性，加之函數本身擬合效果缺陷，在實際應用中并不多見。非線性檢驗表明該趨勢為線性相關，提示線性函數模型可能為最佳擬合模型。
除了上述兩種方法外，還有多種曲線逼近的方法，如二次B樣條、分段多項式等方法。二次B樣條屬于非參數回歸，逼近效果較一般參數回歸好，但因過度擬合等問題只被用于評估新函數模型的擬合效果；分段多項式法也是一種較好的方法，與限制性立方樣條回歸為目前兩大最常用劑量-反應Meta分析的方法。這兩者擬合效果相近，目前公開用于分段多項式法分析程序僅限于SAS軟件 [2, 14-16]，這將會在本團隊后期文章中進行詳細介紹。
當然，每種曲線逼近（擬合）模型均有其優缺點，比如線性或二次函數模型能較好的反應整體變化趨勢，但對于局部細微的變化的察覺靈敏度較低；高次（三次及以上）函數模型對局部變化反應較真實，但可能由于局部變化并不真實等而存在過度擬合問題，導致結果產生偏倚。因此使用者需要對每種方法進行對比并結合臨床背景，選出最適合的模型。
