在循證實踐和決策中，劑量-反應Meta分析已受到不少學者的關注。它不僅可以反映暴露因素與疾病之間獨特的量效關系，且產生的結果證據等級在觀察性研究中處于較高等級。了解此種類型Meta分析在軟件中如何操作實現，顯得尤為重要。目前，有多種軟件用于劑量-反應Meta分析，不同軟件存在不同的特點。本文將介紹如何使用Stata軟件、R軟件、SAS軟件實現劑量-反應關系Meta分析的計算和繪制劑量趨勢圖，對其進行介紹、歸納及深入探討，并對部分過程進行實例演示。
引用本文： 郭鵬, 徐暢, 周權, 張超, 鄢金柱, 張永剛. 劑量-反應Meta分析之軟件的應用. 中國循證醫學雜志, 2016, 16(3): 363-368. doi: 10.7507/1672-2531.20160056 復制
在流行病學中，劑量-反應關系Meta分析（dose-response meta-analysis）常用于探索暴露（干預）水平的變化與患病風險指標發生風險的潛在關系，通過合并多個相似的劑量-反應研究，得出更全面、更可靠的臨床流行病學證據 [1]。國內已有不少學者介紹了劑量-反應關系Meta分析的實現方法 [2-6]，已發表的方法學文章主要集中于介紹Stata、R、SAS軟件實現。且每種軟件具有其獨特優勢及其局限性。
我們選擇常用Stata、R和SAS軟件，結合實例介紹如何實現劑量-反應Meta分析。盡管實現劑量-反應Meta分析的模型有多種，其中最常用的為限制性立方樣條函數模型，且該模型所有軟件代碼開放，因此本文中所使用的模型均為此模型。本次使用的數據來源于《應用STATA做Meta分析》一書 [7]。此研究共納入原始文獻9篇，其中病例-對照研究6篇、隊列研究3篇。具體數據見表 1。
表1
劑量-反應關系Meta分析實例數據
1 應用Stata軟件實現
1.1 Stata下載與操作
用戶可從http://www.stata.com/stata13/中購買最新版Stata 13.0，本次使用版本為13.0。
1.2 命令加載
需要加載的主要命令如下：
ssc install glst
ssc install xabl
其中代碼來源可參照網址http://nicolaorsini. altervista.org/stata/tutorial/g/ts_glst_bmi_ex.txt或http://nicolaorsini.altervista.org/stata/tutorial/g/alcohol_crc.txt。
1.3 實例分析數據
具體數據見表 1，本文主要使用"glst"、"xabl"命令與限制性立方樣條（restricted cubic spline）函數進行劑量-反應 Meta分析。數據錄入后，首先取對數運算，獲得logrr、loglb、logub、se，命令如下：
gen logadjrr=log(adjrr)
gen loglb=log(lb)
gen logub=log(ub)
gen se=((logub-loglb)/(2*invnorm(.975) ))
1.4 非線性劑量-反應Meta分析
對于量效關系不明確的數據，首先對其非線性趨勢進行檢驗，通過選擇合適的鏈接函數，根據P值判斷暴露劑量與結局指標是否存在趨勢關系，命令如下：
capture drop doses*
pctile dose,percentile(10 50 90)
ret list
mkspline doses = dose,knots(`=r(r1) ' `=r(r2) ' `=r(r3) ') cubic displayknots
glst logadjrr doses*,se(se) cov (n cases) pfirst(id study)
testparm doses2
levelsof dose,local(level)
xblc doses*,c(dose) at(`r(levels)') ref(0) eform
其中第一行命令表示獲取各個參數對應"dose"的劑量值；第二行獲取"dose"劑量的百分位數節點值；第三行列出所獲取節點值；第四行使用mkspline函數獲得3個樣條，并獲得對應新變量"doses"；第五行采用兩步法固定效應劑量-反應模型評估非線性趨勢，進行非線性卡方檢驗（這里test檢驗的P值＞0.05，提示研究間呈線性趨勢；若為之前選取節點數為4，則這里應為testparm doses2 doses3）；第六行列出匯總后各個水平劑量值；第七行以0為參考劑量，獲得不同劑量下對應的效應值及可信區間，呈現非線性點估計值。
1.5 線性劑量-反應Meta分析
命令如下：
glst logadjrr dose,se(se) cov(n cases) pfirst(id study) random eform
lincom dose*3,eform
其中第一行表示采用隨機效應劑量-反應模型評估線性趨勢，此處尚不需"bysort id=gen doser=dose-dose"中心化；第二行則評估增加3個單位牛奶制品對應患卵巢癌的風險。
1.6 劑量-反應曲線的繪制
非線性劑量反應關系曲線的繪制使用的是經典的twoway命令，通過上述結果中的數值，可以用以下命令繪制暴露劑量與對應疾病發生風險的趨勢圖，如下：
twoway(line adjrrwithref dose,sort lp(shortdash) lc(black) ) (scatter adjrr dose [aweight=1/se],msymbol(Oh)),scheme(s1mono) ylabel(0.6 0.8 1 1.2 1.6,angle(horiz) format(%3.2fc)) xlabel (0(1) 4) legend(off)ytitle("Relative Risk (95%CI)",margin(right)) xtitle("Milk consumption,cups/day",margin(top_bottom)) name(figure 1B,replace) yscale(log) plotregion(style(none))
當然，用戶可以使用Graph Editor進行圖形的編輯，調整字體的大小、線條的顏色等，如圖 1及圖 2。
圖1
線性關系趨勢圖
圖2
非線性關系趨勢圖
2 應用R軟件實現
2.1 R下載與操作
作為開源的統計軟件，R軟件的使用是免費的，可以從官網http://www.r-project.org/下載最新R軟件安裝程序和相應軟件包。
2.2 命令包加載
本文主要使用"dosersmeta"軟件包進行劑量-反應Meta分析。在軟件安裝完畢、網絡連接情況下，軟件包加載、調用命令如下：
getwd()
dir()
install.packages("dosresmeta")
install.packages("mvmeta")
install.packages("rms")
install.packages("dosresmeta")
library("Hmisc")
library("survival")
library("SparseM")
library("rms")
library("dosresmeta")
2.3 非線性趨勢劑量-反應Meta分析
具體數據見表 1，可將excel數據格式，轉化為csv或txt格式便于使用。非線性趨勢分析命令如下：
Milk_oc＜- read.table("Milk_oc.txt",header=T)
Milk_oc
knots＜- quantile(Milk_oc$dose,c(0.1,0.5,0.9) )
knots
spl＜- dosresmeta(formula = logadjrr ~ rcs(dose,knots),type = type,id = id,se = se,cases = cases,n = n,data = Milk_oc)
summary(spl)
wald.test(b = coef(spl),Sigma = vcov(spl),Terms = 1:2)
pred＜- predict(spl,data.frame(dose = seq(0,5,0.5) ))
round(pred,3)
第一行用于讀取"Milk_oc.txt"中數據；第二行列出數據"Milk_oc.txt"；第三行則獲取3個百分節點（也可定義為4個節點）；第四行列出節點值；第五行使用dosresmeta程序包進行非線性趨勢分析；第六行將非線性趨勢結果存于spl中；第七行使用wald卡方檢驗判斷是否為非線性關系；第八行以"0"為參考劑量，以劑量"0.5"為間距，以最高劑量為"5"，計算出對應的RR及其95%CI值；第九行獲取上述數據列表。
2.4 線性趨勢劑量-反應Meta分析
探討因變量與自變量自建的線性關系命令如下：
lin＜- dosresmeta(formula = logadjrr ~ dose,type = type,id = id,se = se,cases = cases,n = n,data = Milk_oc)
summary(lin)
predict(lin,delta = 3)
這三行命令的作用分別如下：使用dosresmeta程序包進行線性趨勢分析；將線性趨勢結果存于lin中；計算增加3個單位對應的RR及其95%CI值。
2.5 線性作圖
圖 3所示為線性趨勢圖，R軟件命令如下：
pred＜- predict(lin,data.frame(dose = seq(0,4,1) ),xref = 0)
pred
with(pred,{
plot(dose,pred,type = "l",log = "y",ylab = "Relative risk",las = 1,
xlab = "Milk intake,cups per day",ylim = c (0.8,1.5) ,bty = "l")
lines(dose,ci.lb,lty = "dashed")
lines(dose,ci.ub,lty = "dashed")
})
rug(Milk_oc$dose,quiet = T)
圖3
線性關系趨勢圖
2.6 非線性作圖
圖 4為牛奶與卵巢癌之間的非線性趨勢圖。其繪制與線性圖的繪制類似，主要是第一行命令的改變，命令如下：
newdata = data.frame(dose = seq(0,5,0.5) )
with(predict(spl,newdata),{
plot(get("rcs(dose,knots)dose"),pred,type = "l",log = "y",ylab = "Relative risk",las = 1,
xlab = "Milk intake,grams/day",ylim = c(0.6,2) ,bty = "l")
lines(get("rcs(dose,knots)dose"),ci.lb,lty = "dashed")
lines(get("rcs(dose,knots)dose"),ci.ub,lty = "dashed")
})
rug(Milk_oc$dose,quiet = T)
圖4
非線性關系趨勢圖
3 應用SAS軟件實現
3.1 SAS下載與操作
用戶可以從網址http://www.sas.com/zh_cn/home，購買下載使用SAS軟件。
3.2 metadose宏的獲取
通過網站http://www.hsph.harvard.edu/faculty/spiegelman/metadose，可以下載Metadose宏程序及獲取對應該宏程序講解的pdf文檔（The SAS %METADOSE Macro）。
3.3 具體軟件操作
3.3.1 線性趨勢劑量-反應Meta分析
具體數據見表 1，將上述示例數據錄入SAS軟件，建立名稱為"Milkstudies"數據集，將下載好的%metadose宏導入SAS軟件運行一次，然后運行如下代碼：
%metadose(dat = Milkstudies,ratio = adjrr,UB = ub,Ncase = cases,Ntotal = n,dose = dose,studyname = author,studytype = type,meta = T,wt = 1,unit_wt = servings/day,var_covar = GH,linearCheck = 1,nk = 0,ci = 2,graphtitle = Milk products and ovarian cancer risk)
此處，由于劑量-反應的SAS宏默認是4個節點的非線性立方樣條回歸分析，代碼中常省略"nk=4"，故進行線性分析時，需要重新設定nk值（"nk=0"），然后得到線性分析結果。但是，最終運算得到的結果與非線性結果一致，同時軟件可自動返回線性圖形（圖 5）。
圖5
線性關系圖
3.3.2 非線性趨勢劑量-反應Meta分析
與線性趨勢Meta分析類似，將上述示例數據錄入SAS軟件，建立名稱為"Milkstudies"數據集，將下載好的%metadose宏導入SAS軟件運行一次，然后運行如下代碼：
%metadose(dat = Milkstudies,ratio = adjrr,UB = ub,Ncase = cases,Ntotal = n,dose = dose,studyname = author,studytype = type,meta = T,wt = 1,unit_wt = servings/day,var_covar = GH,linearCheck = 1,ci = 2,graphtitle = Milk products and ovarian cancer risk)
其中，meta＝T/N為是否進行多個研究效應量合并，wt表示單位劑量增量大小，unit_wt表示暴露劑量單位，var_covar表示研究中協方差矩陣估計的方法，包括Greenland method（G）法、Hamling method（H） 法、Greenland and Hamling method（GH）法以及用戶自定義選項（K），其中系統默認為GH；linearCheck用于表示是否進行非線性檢驗；ci可設置為0、1或2，其中0表示不繪制95%CI線條，2表示繪制點線形狀95%CI線條，1表示以灰色區域繪制95%CI線條；graphtitle表示所繪制線條名稱。此處給出的結果是以每增加1個單位暴露劑量合并的RR/OR值。另外，劑量-反應的SAS宏默認是4個節點的非線性立方樣條回歸分析，代碼中常省略"nk=4"這一語句（圖 6）。
圖6
非線性關系圖
4 討論
在流行病學中，常常需要評價暴露因素變化水平與疾病風險的趨勢關系，即是否存在劑量-反應關系，可以建立劑量-反應關系回歸模型，采用廣義最小二乘法（generalized least squares，GLS）對此類數據進行趨勢估計，此方法由Greenland和Longnecker兩位學者開發 [8]，并經Nicola Orsin及其團隊對其進行了良好的完善及推廣 [9-11]。
目前，已知可以用于劑量-反應Meta分析的軟件主要包括Stata、R和SAS軟件，其中Stata軟件利用GLST、XBLC命令實現劑量-反應關系分析 [2, 6-11]，R軟件主要通過dosresmeta、mvmeta程序包實現劑量-反應關系數據分析 [3, 11]，SAS軟件是基于%metadose宏命令實現劑量-反應關系分析 [5]。上述3種軟件進行劑量-反應關系Meta分析時各有特點。在示例展示中，不難發現，運用Stata軟件時，操作界面簡約，代碼的輸入運用靈活，可操作性強。在圖形繪制方面，Graph Editor的運用可使結果的呈現更形象、合理、規范。R代碼為參考GLST命令運算理念設計而來，且R的結果呈現能力更強大，與GLST命令相比，R軟件還可計算出I2值，并增加了擬合優度的內容。但在圖形制作方面，Stata軟件優于R軟件。SAS軟件%metadose宏做劑量-反應關系的優勢在于只需要較少的代碼，就能得到所有的分析結果，特別是繪圖方面，均由宏程序自已生成。這款軟件的局限性在于靈活性不夠，繪圖也不夠美觀。但該三種軟件操作稍顯復雜，需要適當補充相關軟件的基礎知識，合理選用一種或聯合起來使用，充分發揮軟件特點，以達到最佳效果。
在處理劑量-反應數據過程中常遇到的一個問題是，由于GLST命令要求處理的暴露劑量最低值為0，故遇到最低劑量組不可能為0的數據，比如體重指數等時，在處理數據前，應當進行代碼語句中心化或手動中心化 [6]。有興趣的讀者還可參閱相關文獻 [3, 4, 12]對劑量-反應分析中常見問題的處理方式。
在我們熟練運用上述軟件，獲得更多流行病學證據的同時，期望更多學者致力于劑量-反應方法學的研究，為劑量-反應Meta分析的發展打下更堅實的基礎。
在流行病學中，劑量-反應關系Meta分析（dose-response meta-analysis）常用于探索暴露（干預）水平的變化與患病風險指標發生風險的潛在關系，通過合并多個相似的劑量-反應研究，得出更全面、更可靠的臨床流行病學證據 [1]。國內已有不少學者介紹了劑量-反應關系Meta分析的實現方法 [2-6]，已發表的方法學文章主要集中于介紹Stata、R、SAS軟件實現。且每種軟件具有其獨特優勢及其局限性。
我們選擇常用Stata、R和SAS軟件，結合實例介紹如何實現劑量-反應Meta分析。盡管實現劑量-反應Meta分析的模型有多種，其中最常用的為限制性立方樣條函數模型，且該模型所有軟件代碼開放，因此本文中所使用的模型均為此模型。本次使用的數據來源于《應用STATA做Meta分析》一書 [7]。此研究共納入原始文獻9篇，其中病例-對照研究6篇、隊列研究3篇。具體數據見表 1。
表1
劑量-反應關系Meta分析實例數據
1 應用Stata軟件實現
1.1 Stata下載與操作
用戶可從http://www.stata.com/stata13/中購買最新版Stata 13.0，本次使用版本為13.0。
1.2 命令加載
需要加載的主要命令如下：
ssc install glst
ssc install xabl
其中代碼來源可參照網址http://nicolaorsini. altervista.org/stata/tutorial/g/ts_glst_bmi_ex.txt或http://nicolaorsini.altervista.org/stata/tutorial/g/alcohol_crc.txt。
1.3 實例分析數據
具體數據見表 1，本文主要使用"glst"、"xabl"命令與限制性立方樣條（restricted cubic spline）函數進行劑量-反應 Meta分析。數據錄入后，首先取對數運算，獲得logrr、loglb、logub、se，命令如下：
gen logadjrr=log(adjrr)
gen loglb=log(lb)
gen logub=log(ub)
gen se=((logub-loglb)/(2*invnorm(.975) ))
1.4 非線性劑量-反應Meta分析
對于量效關系不明確的數據，首先對其非線性趨勢進行檢驗，通過選擇合適的鏈接函數，根據P值判斷暴露劑量與結局指標是否存在趨勢關系，命令如下：
capture drop doses*
pctile dose,percentile(10 50 90)
ret list
mkspline doses = dose,knots(`=r(r1) ' `=r(r2) ' `=r(r3) ') cubic displayknots
glst logadjrr doses*,se(se) cov (n cases) pfirst(id study)
testparm doses2
levelsof dose,local(level)
xblc doses*,c(dose) at(`r(levels)') ref(0) eform
其中第一行命令表示獲取各個參數對應"dose"的劑量值；第二行獲取"dose"劑量的百分位數節點值；第三行列出所獲取節點值；第四行使用mkspline函數獲得3個樣條，并獲得對應新變量"doses"；第五行采用兩步法固定效應劑量-反應模型評估非線性趨勢，進行非線性卡方檢驗（這里test檢驗的P值＞0.05，提示研究間呈線性趨勢；若為之前選取節點數為4，則這里應為testparm doses2 doses3）；第六行列出匯總后各個水平劑量值；第七行以0為參考劑量，獲得不同劑量下對應的效應值及可信區間，呈現非線性點估計值。
1.5 線性劑量-反應Meta分析
命令如下：
glst logadjrr dose,se(se) cov(n cases) pfirst(id study) random eform
lincom dose*3,eform
其中第一行表示采用隨機效應劑量-反應模型評估線性趨勢，此處尚不需"bysort id=gen doser=dose-dose"中心化；第二行則評估增加3個單位牛奶制品對應患卵巢癌的風險。
1.6 劑量-反應曲線的繪制
非線性劑量反應關系曲線的繪制使用的是經典的twoway命令，通過上述結果中的數值，可以用以下命令繪制暴露劑量與對應疾病發生風險的趨勢圖，如下：
twoway(line adjrrwithref dose,sort lp(shortdash) lc(black) ) (scatter adjrr dose [aweight=1/se],msymbol(Oh)),scheme(s1mono) ylabel(0.6 0.8 1 1.2 1.6,angle(horiz) format(%3.2fc)) xlabel (0(1) 4) legend(off)ytitle("Relative Risk (95%CI)",margin(right)) xtitle("Milk consumption,cups/day",margin(top_bottom)) name(figure 1B,replace) yscale(log) plotregion(style(none))
當然，用戶可以使用Graph Editor進行圖形的編輯，調整字體的大小、線條的顏色等，如圖 1及圖 2。
圖1
線性關系趨勢圖
圖2
非線性關系趨勢圖
2 應用R軟件實現
2.1 R下載與操作
作為開源的統計軟件，R軟件的使用是免費的，可以從官網http://www.r-project.org/下載最新R軟件安裝程序和相應軟件包。
2.2 命令包加載
本文主要使用"dosersmeta"軟件包進行劑量-反應Meta分析。在軟件安裝完畢、網絡連接情況下，軟件包加載、調用命令如下：
getwd()
dir()
install.packages("dosresmeta")
install.packages("mvmeta")
install.packages("rms")
install.packages("dosresmeta")
library("Hmisc")
library("survival")
library("SparseM")
library("rms")
library("dosresmeta")
2.3 非線性趨勢劑量-反應Meta分析
具體數據見表 1，可將excel數據格式，轉化為csv或txt格式便于使用。非線性趨勢分析命令如下：
Milk_oc＜- read.table("Milk_oc.txt",header=T)
Milk_oc
knots＜- quantile(Milk_oc$dose,c(0.1,0.5,0.9) )
knots
spl＜- dosresmeta(formula = logadjrr ~ rcs(dose,knots),type = type,id = id,se = se,cases = cases,n = n,data = Milk_oc)
summary(spl)
wald.test(b = coef(spl),Sigma = vcov(spl),Terms = 1:2)
pred＜- predict(spl,data.frame(dose = seq(0,5,0.5) ))
round(pred,3)
第一行用于讀取"Milk_oc.txt"中數據；第二行列出數據"Milk_oc.txt"；第三行則獲取3個百分節點（也可定義為4個節點）；第四行列出節點值；第五行使用dosresmeta程序包進行非線性趨勢分析；第六行將非線性趨勢結果存于spl中；第七行使用wald卡方檢驗判斷是否為非線性關系；第八行以"0"為參考劑量，以劑量"0.5"為間距，以最高劑量為"5"，計算出對應的RR及其95%CI值；第九行獲取上述數據列表。
2.4 線性趨勢劑量-反應Meta分析
探討因變量與自變量自建的線性關系命令如下：
lin＜- dosresmeta(formula = logadjrr ~ dose,type = type,id = id,se = se,cases = cases,n = n,data = Milk_oc)
summary(lin)
predict(lin,delta = 3)
這三行命令的作用分別如下：使用dosresmeta程序包進行線性趨勢分析；將線性趨勢結果存于lin中；計算增加3個單位對應的RR及其95%CI值。
2.5 線性作圖
圖 3所示為線性趨勢圖，R軟件命令如下：
pred＜- predict(lin,data.frame(dose = seq(0,4,1) ),xref = 0)
pred
with(pred,{
plot(dose,pred,type = "l",log = "y",ylab = "Relative risk",las = 1,
xlab = "Milk intake,cups per day",ylim = c (0.8,1.5) ,bty = "l")
lines(dose,ci.lb,lty = "dashed")
lines(dose,ci.ub,lty = "dashed")
})
rug(Milk_oc$dose,quiet = T)
圖3
線性關系趨勢圖
2.6 非線性作圖
圖 4為牛奶與卵巢癌之間的非線性趨勢圖。其繪制與線性圖的繪制類似，主要是第一行命令的改變，命令如下：
newdata = data.frame(dose = seq(0,5,0.5) )
with(predict(spl,newdata),{
plot(get("rcs(dose,knots)dose"),pred,type = "l",log = "y",ylab = "Relative risk",las = 1,
xlab = "Milk intake,grams/day",ylim = c(0.6,2) ,bty = "l")
lines(get("rcs(dose,knots)dose"),ci.lb,lty = "dashed")
lines(get("rcs(dose,knots)dose"),ci.ub,lty = "dashed")
})
rug(Milk_oc$dose,quiet = T)
圖4
非線性關系趨勢圖
3 應用SAS軟件實現
3.1 SAS下載與操作
用戶可以從網址http://www.sas.com/zh_cn/home，購買下載使用SAS軟件。
3.2 metadose宏的獲取
通過網站http://www.hsph.harvard.edu/faculty/spiegelman/metadose，可以下載Metadose宏程序及獲取對應該宏程序講解的pdf文檔（The SAS %METADOSE Macro）。
3.3 具體軟件操作
3.3.1 線性趨勢劑量-反應Meta分析
具體數據見表 1，將上述示例數據錄入SAS軟件，建立名稱為"Milkstudies"數據集，將下載好的%metadose宏導入SAS軟件運行一次，然后運行如下代碼：
%metadose(dat = Milkstudies,ratio = adjrr,UB = ub,Ncase = cases,Ntotal = n,dose = dose,studyname = author,studytype = type,meta = T,wt = 1,unit_wt = servings/day,var_covar = GH,linearCheck = 1,nk = 0,ci = 2,graphtitle = Milk products and ovarian cancer risk)
此處，由于劑量-反應的SAS宏默認是4個節點的非線性立方樣條回歸分析，代碼中常省略"nk=4"，故進行線性分析時，需要重新設定nk值（"nk=0"），然后得到線性分析結果。但是，最終運算得到的結果與非線性結果一致，同時軟件可自動返回線性圖形（圖 5）。
圖5
線性關系圖
3.3.2 非線性趨勢劑量-反應Meta分析
與線性趨勢Meta分析類似，將上述示例數據錄入SAS軟件，建立名稱為"Milkstudies"數據集，將下載好的%metadose宏導入SAS軟件運行一次，然后運行如下代碼：
%metadose(dat = Milkstudies,ratio = adjrr,UB = ub,Ncase = cases,Ntotal = n,dose = dose,studyname = author,studytype = type,meta = T,wt = 1,unit_wt = servings/day,var_covar = GH,linearCheck = 1,ci = 2,graphtitle = Milk products and ovarian cancer risk)
其中，meta＝T/N為是否進行多個研究效應量合并，wt表示單位劑量增量大小，unit_wt表示暴露劑量單位，var_covar表示研究中協方差矩陣估計的方法，包括Greenland method（G）法、Hamling method（H） 法、Greenland and Hamling method（GH）法以及用戶自定義選項（K），其中系統默認為GH；linearCheck用于表示是否進行非線性檢驗；ci可設置為0、1或2，其中0表示不繪制95%CI線條，2表示繪制點線形狀95%CI線條，1表示以灰色區域繪制95%CI線條；graphtitle表示所繪制線條名稱。此處給出的結果是以每增加1個單位暴露劑量合并的RR/OR值。另外，劑量-反應的SAS宏默認是4個節點的非線性立方樣條回歸分析，代碼中常省略"nk=4"這一語句（圖 6）。
圖6
非線性關系圖
4 討論
在流行病學中，常常需要評價暴露因素變化水平與疾病風險的趨勢關系，即是否存在劑量-反應關系，可以建立劑量-反應關系回歸模型，采用廣義最小二乘法（generalized least squares，GLS）對此類數據進行趨勢估計，此方法由Greenland和Longnecker兩位學者開發 [8]，并經Nicola Orsin及其團隊對其進行了良好的完善及推廣 [9-11]。
目前，已知可以用于劑量-反應Meta分析的軟件主要包括Stata、R和SAS軟件，其中Stata軟件利用GLST、XBLC命令實現劑量-反應關系分析 [2, 6-11]，R軟件主要通過dosresmeta、mvmeta程序包實現劑量-反應關系數據分析 [3, 11]，SAS軟件是基于%metadose宏命令實現劑量-反應關系分析 [5]。上述3種軟件進行劑量-反應關系Meta分析時各有特點。在示例展示中，不難發現，運用Stata軟件時，操作界面簡約，代碼的輸入運用靈活，可操作性強。在圖形繪制方面，Graph Editor的運用可使結果的呈現更形象、合理、規范。R代碼為參考GLST命令運算理念設計而來，且R的結果呈現能力更強大，與GLST命令相比，R軟件還可計算出I2值，并增加了擬合優度的內容。但在圖形制作方面，Stata軟件優于R軟件。SAS軟件%metadose宏做劑量-反應關系的優勢在于只需要較少的代碼，就能得到所有的分析結果，特別是繪圖方面，均由宏程序自已生成。這款軟件的局限性在于靈活性不夠，繪圖也不夠美觀。但該三種軟件操作稍顯復雜，需要適當補充相關軟件的基礎知識，合理選用一種或聯合起來使用，充分發揮軟件特點，以達到最佳效果。
在處理劑量-反應數據過程中常遇到的一個問題是，由于GLST命令要求處理的暴露劑量最低值為0，故遇到最低劑量組不可能為0的數據，比如體重指數等時，在處理數據前，應當進行代碼語句中心化或手動中心化 [6]。有興趣的讀者還可參閱相關文獻 [3, 4, 12]對劑量-反應分析中常見問題的處理方式。
在我們熟練運用上述軟件，獲得更多流行病學證據的同時，期望更多學者致力于劑量-反應方法學的研究，為劑量-反應Meta分析的發展打下更堅實的基礎。
