引用本文： 張云權, 馬露, 馮仁杰, 朱耀輝, 李存祿. 模型回歸系數的合并分析在R軟件metafor包中的實現. 中國循證醫學雜志, 2015, 15(3): 367-372. doi: 10.7507/1672-2531.20150062 復制
Meta分析是對以往的研究結果進行系統定量綜合的統計分析方法 [1]，在循證醫學當中被公認是最高級別的循證決策證據 [2]。伴隨在醫學領域中的廣泛應用，Meta分析的類型也不斷豐富，目前主要包括常規Meta分析（二分類變量和連續型變量的Meta分析等）、單組率的Meta分析、累積Meta分析、診斷性Meta分析和網狀Meta分析等 [3]。而針對不同研究中通過構建相近模型（如多元線性回歸、logistic回歸、Cox回歸模型等）進行分析獲得的關于某個特定解釋變量與結局變量關聯的回歸系數的效應合并可以看作Meta分析中的一種特殊類型，較少受到國內研究者關注。
Meta分析的軟件十分多樣，除了傳統的Review Manager和Stata外，免費、開源的R軟件在Meta分析和統計分析中的應用也已逐漸受到越來越多的科研工作者的青睞 [4]。R軟件中可用于Meta分析的程序包主要有metafor、meta、Rmeta等，其中metafor包功能更加強大，應用更靈活 [5, 6]。本文主要通過一個已發表的Meta分析實例，介紹模型回歸系數合并的Meta分析在metafor包中的實現。
1 研究實例
1.1 數據資料
以Stieb等 [7]發表在《Environmental Research》上的《Ambient air pollution，birth weight and preterm birth：A systematic review and meta-analysis》一文中關于孕婦整個孕期的大氣顆粒物PM10暴露對新生兒低出生體重（出生體重＜2 500 g）影響的研究數據（森林圖）為例，數據整理如表 1所示。
表1
顆粒物PM10暴露與新生兒低出生體重發生風險關系的14項研究數據
表2
轉換和整理后的研究數據
表 1中location表示研究地點，year為文獻發表時間，OR為孕婦孕期PM10暴露每升高20 μg/m3，新生兒發生低出生體重的優勢比，ci.lb、ci.ub分別為相應OR值的可信區間（CI）下限和上限。
1.2 安裝并加載metafor程序包
＜install.packages(“metafor”)
＜library(metafor)
1.3 數據集的準備
由于上述文獻中不同研究均采用logistic回歸模型分析PM10暴露與低出生體重結局的關聯，研究中獲取的數據均為各研究模型分析得到的經回歸系數轉換而得的優勢比OR及其CI，其95%CI=exp(ln(OR)±1.96se)（se為回歸系數的標準誤）。由于R軟件metafor包中如下基本命令可進行Meta分析的計算：rma(yi,sei,data)，其中yi為各研究的效應值，sei為各效應值相應的標準誤。而優勢比OR由于為logistic回歸模型的回歸系數經指數轉換而來，因而需對研究數據按上述公式進行轉換和整理。
1.4 數據讀取和模型構建
通過read.table將上述csv數據文件導入R數據框metadat中（excel另存為csv格式）（read.table后引號內為csv文件存放的磁盤路徑，#后為R程序注釋語句，不在程序中運行）。
＜metadat＜- read.table ("F:/SkyDrive/Learning/Rlearning/meta.csv",header=T,sep=",")
＜metamod＜-rma(yi=lnOR,sei=se,data=metadat,method="DL") #利用rma函數構建模型
＜summary(metamod) #查看模型結果
由上述summary( )命令的結果（圖 1）總結可知，納入研究數為13，采用隨機效應模型，DL法方差分量tau2=0.002 2，se=0.047 0，I2=69.08%，H2=3.23，異質性檢驗Q=42.04，P＜0.01，說明納入研究間存在較大異質性。效應值合并結果（lnOR）為0.093 3，標準誤為0.022 9，P＜0.01。OR=exp(0.093 3) =1.10，95%CI下限為exp(0.048 5) =1.05，上限為exp(0.138 2) =1.15，這與文獻報道結果一致。上述對數轉換結果也可通過predict(metamod,transf=exp)獲得。
圖1
利用rma函數構建模型結果
1.5 森林圖
R中metafor包可利用forest函數繪制森林圖，本研究可通過代碼forest(metamod,transf=exp,refline=1) 粗略地實現（圖 2）。為使森林圖反映的信息更全面、可讀性更強，在forest中添加了一些額外設置，并通過text語句分別在森林圖區域中相應坐標添加了研究的地點信息和數據列標簽（圖 3）。
圖2
粗略森林圖
圖3
信息更全面、可讀性更強的森林圖
＜forest(metamod,transf=exp,refline=1,#對效應值進行指數轉換并指定參考線
mlab="Random-effect Model for All Studies",#設定合并效應的標簽
slab=paste(metadat$author,metadat$year,sep=","),#設置研究標簽study label
xlim=c(-3,5) ,alim=c(0,3) ,steps=7,xlab="OR",showweights=T) #設置整個森林圖x軸范圍、實際x軸（包括范圍、軸上觸須數和標簽）、顯示各研究權重
＜text(-1,14:1,pos=4,metadat$location) #在坐標(-1,14:1) 往右添加各研究的地點
＜text(c(-3,-1,4,5) ,16,pos=c(4,4,2,2) ,c("Author(s) and Year","Location","Weight","OR [95% CI]"),cex = 1,font = 2) #在相應坐標、方向添加研究相關信息（作者、發表年份、研究地點）、研究權重和效應值，并設置文本字號、字體
1.6 亞組分析
亞組分析是探索不同研究異質性來源的方式之一，R metafor包中可通過森林圖中的坐標設置（具體坐標須通過不斷調整達到最佳效果）將研究分為不同亞組（代碼1），并利用addpoly函數在森林圖上添加相應亞組模塊的效應值合并結果（代碼2），從而實現亞組分析（圖 4）。本研究將發表年份在2008年以前和以后的研究分為兩個亞組分別進行效應合并，僅做示例（并無具體實際意義），實際應用中應依據異質性來源分析情況進行合理的亞組分組后進行亞組分析。
代碼1：
＜forest(metamod,transf=exp,refline=1,
mlab="RE Model for All Studies",slab=paste(metadat$author,metadat$year,sep=","),xlim=c(-1,4) ,ylim=c(-1,22) ,alim=c(0,3) ,steps=7,cex=.75,xlab="OR",
order=order(-metadat$year),rows=c(2:11,15:18) ) #以研究發表年份由大到小排序，坐標圖中2～11行共10項研究為一個亞組，15～18行共4項研究為另一亞組
＜text(c(-1,4) ,21,pos=c(4,2) ,c("Author(s) and Year","OR [95% CI]"),font = 2) #在相應坐標、方向添加研究相關信息（作者、發表年份）和效應值（包含CI）
＜text(-1,c(19,12) ,c("Publication Years＜=2008","Publication Years＞2008"),pos=4,cex=.75,font=2) #分別設置兩個亞組的分類標簽及坐標位置
圖4
亞組分析的森林圖
代碼2：
＜subgroup1 ＜- rma(yi=lnOR,sei=se,data=metadat,method="DL",subset=(year＜=2008) )
＜subgroup2＜- rma(yi=lnOR,sei=se,data=metadat,method="DL",subset=(year＞2008) ) #將發表年份在2008年以前和以后的研究分為兩個亞組分別進行效應合并
＜addpoly(subgroup1,cex=.75,font=4,#設置亞組1字號、字體
mlab="RE Model for Subgroup1",row=14,#設置合并效應標簽及所在行位置
transf=exp,atransf=exp) #將研究效應值和坐標軸均做指數轉換為OR
＜addpoly(subgroup2,cex=.75,font=4,mlab="RE Model for Subgroup2",row=1,
transf=exp,atransf=exp) #同亞組1
1.7 漏斗圖及對稱性的檢驗
可簡單利用funnel(metamod)實現漏斗圖的繪制，亦可具體設置坐標軸的范圍和標簽。
＜funnel(metamod,xlim=c(-1.5,2) ,ylim=c(0,0.8) ,
xlab="Log Odds Ratio",ylab="Standard Error of logOR",steps=4,at=c(-1,0,1,2) ,digits=1)
漏斗圖的對稱性可在一定程度上反映是否存在發表偏倚，圖 5直觀上判斷對稱性欠佳。metafor包中可采用regtest和ranktest函數分別進行Egger’s檢驗和Begg’s 檢驗，檢驗漏斗圖的對稱性。不同的資料類型，采用Egger’s檢驗和Begg’s檢驗的檢驗效能并不相同，甚至出現檢驗結果矛盾的情況 [8]。目前一般認為，Egger’s檢驗較Begg’s檢驗更高 [9]。本研究采用regtest檢驗結果見圖 6，可以看出對稱性檢驗Z=5.33，P＜0.01，可以認為漏斗圖不對稱，提示可能存在一定的發表偏倚。
圖5
漏斗圖
圖6
regtest檢驗結果
1.8 敏感性分析
leave1out函數可以進行敏感性分析，計算分別剔除單個入選研究后合并的效應值、標準誤以及95%CI等。本例顯示，剔除單個研究對合并效應值影響較小，說明敏感性較低，結果較為穩健可信（圖 7）。
圖7
采用leavleout進行敏感性分析結果
2 討論
本文以效應值OR為例，比較簡略地介紹了模型回歸系數合并的Meta分析在R軟件metafor包中的實現方法。通過利用公開發表的Meta分析文獻數據進行實例演示，我們發現，metafor包進行模型回歸系數合并的Meta分析的結果比較可靠，具有較高的應用價值。
在利用繪制森林圖時，metafor包可十分靈活地通過坐標和方向設置添加研究相關信息和自定義文本信息，從而增加可讀性，同時可根據實際需要設置森林圖中的元素顯示方式（文本字體、字號和坐標軸參數）和位置，并輕松實現亞組分析。與同樣免費的Review Manager軟件相比，Metafor包提供的regtest和ranktest函數可分別對漏斗圖對稱性進行Egger’s檢驗和Begg’s檢檢驗 [5]，leave1out函數則可輕松實現敏感性分析，應用起來更為全面、高效。
此外，metafor包也可進行多變量Meta回歸分析，探討Meta分析中不同研究間異質性的來源，還能進行混合效應模型（包括單個、多個分類或連續性變量）擬合運算，檢驗模型系數并獲得CI，以及對參數進行精確檢驗如置換檢驗（permutation tests） [6]。
綜上所述，R軟件metafor包不僅可準確靈活地實現模型回歸系數合并的Meta分析，同時在發表偏倚的識別、檢驗以及異質性來源分析等方面的功能也十分強大，在Meta分析中應用前景十分廣泛。
Meta分析是對以往的研究結果進行系統定量綜合的統計分析方法 [1]，在循證醫學當中被公認是最高級別的循證決策證據 [2]。伴隨在醫學領域中的廣泛應用，Meta分析的類型也不斷豐富，目前主要包括常規Meta分析（二分類變量和連續型變量的Meta分析等）、單組率的Meta分析、累積Meta分析、診斷性Meta分析和網狀Meta分析等 [3]。而針對不同研究中通過構建相近模型（如多元線性回歸、logistic回歸、Cox回歸模型等）進行分析獲得的關于某個特定解釋變量與結局變量關聯的回歸系數的效應合并可以看作Meta分析中的一種特殊類型，較少受到國內研究者關注。
Meta分析的軟件十分多樣，除了傳統的Review Manager和Stata外，免費、開源的R軟件在Meta分析和統計分析中的應用也已逐漸受到越來越多的科研工作者的青睞 [4]。R軟件中可用于Meta分析的程序包主要有metafor、meta、Rmeta等，其中metafor包功能更加強大，應用更靈活 [5, 6]。本文主要通過一個已發表的Meta分析實例，介紹模型回歸系數合并的Meta分析在metafor包中的實現。
1 研究實例
1.1 數據資料
以Stieb等 [7]發表在《Environmental Research》上的《Ambient air pollution，birth weight and preterm birth：A systematic review and meta-analysis》一文中關于孕婦整個孕期的大氣顆粒物PM10暴露對新生兒低出生體重（出生體重＜2 500 g）影響的研究數據（森林圖）為例，數據整理如表 1所示。
表1
顆粒物PM10暴露與新生兒低出生體重發生風險關系的14項研究數據
表2
轉換和整理后的研究數據
表 1中location表示研究地點，year為文獻發表時間，OR為孕婦孕期PM10暴露每升高20 μg/m3，新生兒發生低出生體重的優勢比，ci.lb、ci.ub分別為相應OR值的可信區間（CI）下限和上限。
1.2 安裝并加載metafor程序包
＜install.packages(“metafor”)
＜library(metafor)
1.3 數據集的準備
由于上述文獻中不同研究均采用logistic回歸模型分析PM10暴露與低出生體重結局的關聯，研究中獲取的數據均為各研究模型分析得到的經回歸系數轉換而得的優勢比OR及其CI，其95%CI=exp(ln(OR)±1.96se)（se為回歸系數的標準誤）。由于R軟件metafor包中如下基本命令可進行Meta分析的計算：rma(yi,sei,data)，其中yi為各研究的效應值，sei為各效應值相應的標準誤。而優勢比OR由于為logistic回歸模型的回歸系數經指數轉換而來，因而需對研究數據按上述公式進行轉換和整理。
1.4 數據讀取和模型構建
通過read.table將上述csv數據文件導入R數據框metadat中（excel另存為csv格式）（read.table后引號內為csv文件存放的磁盤路徑，#后為R程序注釋語句，不在程序中運行）。
＜metadat＜- read.table ("F:/SkyDrive/Learning/Rlearning/meta.csv",header=T,sep=",")
＜metamod＜-rma(yi=lnOR,sei=se,data=metadat,method="DL") #利用rma函數構建模型
＜summary(metamod) #查看模型結果
由上述summary( )命令的結果（圖 1）總結可知，納入研究數為13，采用隨機效應模型，DL法方差分量tau2=0.002 2，se=0.047 0，I2=69.08%，H2=3.23，異質性檢驗Q=42.04，P＜0.01，說明納入研究間存在較大異質性。效應值合并結果（lnOR）為0.093 3，標準誤為0.022 9，P＜0.01。OR=exp(0.093 3) =1.10，95%CI下限為exp(0.048 5) =1.05，上限為exp(0.138 2) =1.15，這與文獻報道結果一致。上述對數轉換結果也可通過predict(metamod,transf=exp)獲得。
圖1
利用rma函數構建模型結果
1.5 森林圖
R中metafor包可利用forest函數繪制森林圖，本研究可通過代碼forest(metamod,transf=exp,refline=1) 粗略地實現（圖 2）。為使森林圖反映的信息更全面、可讀性更強，在forest中添加了一些額外設置，并通過text語句分別在森林圖區域中相應坐標添加了研究的地點信息和數據列標簽（圖 3）。
圖2
粗略森林圖
圖3
信息更全面、可讀性更強的森林圖
＜forest(metamod,transf=exp,refline=1,#對效應值進行指數轉換并指定參考線
mlab="Random-effect Model for All Studies",#設定合并效應的標簽
slab=paste(metadat$author,metadat$year,sep=","),#設置研究標簽study label
xlim=c(-3,5) ,alim=c(0,3) ,steps=7,xlab="OR",showweights=T) #設置整個森林圖x軸范圍、實際x軸（包括范圍、軸上觸須數和標簽）、顯示各研究權重
＜text(-1,14:1,pos=4,metadat$location) #在坐標(-1,14:1) 往右添加各研究的地點
＜text(c(-3,-1,4,5) ,16,pos=c(4,4,2,2) ,c("Author(s) and Year","Location","Weight","OR [95% CI]"),cex = 1,font = 2) #在相應坐標、方向添加研究相關信息（作者、發表年份、研究地點）、研究權重和效應值，并設置文本字號、字體
1.6 亞組分析
亞組分析是探索不同研究異質性來源的方式之一，R metafor包中可通過森林圖中的坐標設置（具體坐標須通過不斷調整達到最佳效果）將研究分為不同亞組（代碼1），并利用addpoly函數在森林圖上添加相應亞組模塊的效應值合并結果（代碼2），從而實現亞組分析（圖 4）。本研究將發表年份在2008年以前和以后的研究分為兩個亞組分別進行效應合并，僅做示例（并無具體實際意義），實際應用中應依據異質性來源分析情況進行合理的亞組分組后進行亞組分析。
代碼1：
＜forest(metamod,transf=exp,refline=1,
mlab="RE Model for All Studies",slab=paste(metadat$author,metadat$year,sep=","),xlim=c(-1,4) ,ylim=c(-1,22) ,alim=c(0,3) ,steps=7,cex=.75,xlab="OR",
order=order(-metadat$year),rows=c(2:11,15:18) ) #以研究發表年份由大到小排序，坐標圖中2～11行共10項研究為一個亞組，15～18行共4項研究為另一亞組
＜text(c(-1,4) ,21,pos=c(4,2) ,c("Author(s) and Year","OR [95% CI]"),font = 2) #在相應坐標、方向添加研究相關信息（作者、發表年份）和效應值（包含CI）
＜text(-1,c(19,12) ,c("Publication Years＜=2008","Publication Years＞2008"),pos=4,cex=.75,font=2) #分別設置兩個亞組的分類標簽及坐標位置
圖4
亞組分析的森林圖
代碼2：
＜subgroup1 ＜- rma(yi=lnOR,sei=se,data=metadat,method="DL",subset=(year＜=2008) )
＜subgroup2＜- rma(yi=lnOR,sei=se,data=metadat,method="DL",subset=(year＞2008) ) #將發表年份在2008年以前和以后的研究分為兩個亞組分別進行效應合并
＜addpoly(subgroup1,cex=.75,font=4,#設置亞組1字號、字體
mlab="RE Model for Subgroup1",row=14,#設置合并效應標簽及所在行位置
transf=exp,atransf=exp) #將研究效應值和坐標軸均做指數轉換為OR
＜addpoly(subgroup2,cex=.75,font=4,mlab="RE Model for Subgroup2",row=1,
transf=exp,atransf=exp) #同亞組1
1.7 漏斗圖及對稱性的檢驗
可簡單利用funnel(metamod)實現漏斗圖的繪制，亦可具體設置坐標軸的范圍和標簽。
＜funnel(metamod,xlim=c(-1.5,2) ,ylim=c(0,0.8) ,
xlab="Log Odds Ratio",ylab="Standard Error of logOR",steps=4,at=c(-1,0,1,2) ,digits=1)
漏斗圖的對稱性可在一定程度上反映是否存在發表偏倚，圖 5直觀上判斷對稱性欠佳。metafor包中可采用regtest和ranktest函數分別進行Egger’s檢驗和Begg’s 檢驗，檢驗漏斗圖的對稱性。不同的資料類型，采用Egger’s檢驗和Begg’s檢驗的檢驗效能并不相同，甚至出現檢驗結果矛盾的情況 [8]。目前一般認為，Egger’s檢驗較Begg’s檢驗更高 [9]。本研究采用regtest檢驗結果見圖 6，可以看出對稱性檢驗Z=5.33，P＜0.01，可以認為漏斗圖不對稱，提示可能存在一定的發表偏倚。
圖5
漏斗圖
圖6
regtest檢驗結果
1.8 敏感性分析
leave1out函數可以進行敏感性分析，計算分別剔除單個入選研究后合并的效應值、標準誤以及95%CI等。本例顯示，剔除單個研究對合并效應值影響較小，說明敏感性較低，結果較為穩健可信（圖 7）。
圖7
采用leavleout進行敏感性分析結果
2 討論
本文以效應值OR為例，比較簡略地介紹了模型回歸系數合并的Meta分析在R軟件metafor包中的實現方法。通過利用公開發表的Meta分析文獻數據進行實例演示，我們發現，metafor包進行模型回歸系數合并的Meta分析的結果比較可靠，具有較高的應用價值。
在利用繪制森林圖時，metafor包可十分靈活地通過坐標和方向設置添加研究相關信息和自定義文本信息，從而增加可讀性，同時可根據實際需要設置森林圖中的元素顯示方式（文本字體、字號和坐標軸參數）和位置，并輕松實現亞組分析。與同樣免費的Review Manager軟件相比，Metafor包提供的regtest和ranktest函數可分別對漏斗圖對稱性進行Egger’s檢驗和Begg’s檢檢驗 [5]，leave1out函數則可輕松實現敏感性分析，應用起來更為全面、高效。
此外，metafor包也可進行多變量Meta回歸分析，探討Meta分析中不同研究間異質性的來源，還能進行混合效應模型（包括單個、多個分類或連續性變量）擬合運算，檢驗模型系數并獲得CI，以及對參數進行精確檢驗如置換檢驗（permutation tests） [6]。
綜上所述，R軟件metafor包不僅可準確靈活地實現模型回歸系數合并的Meta分析，同時在發表偏倚的識別、檢驗以及異質性來源分析等方面的功能也十分強大，在Meta分析中應用前景十分廣泛。
