R 軟件是一款具備強大的數據處理和圖形繪制功能的免費軟件,其在 Meta 分析中的應用日趨廣泛。R 軟件中有 meta、metafor、metaplus 等程序包可實現 Meta 分析。本文以實例對 R 軟件 metaplus 程序包實現 Meta 分析的過程和結果進行展示。
引用本文: 郭梓鑫, 馮雨嘉, 王清華, 曾憲濤, 劉小平. 應用 R 軟件 metaplus 程序包實現 Meta 分析. 中國循證醫學雜志, 2018, 18(7): 763-768. doi: 10.7507/1672-2531.201708086 復制
隨著循證醫學的發展,Meta 分析理論與方法學也隨著需求而不斷發展。目前 Meta 分析被公認為客觀評價和合成針對某一特定問題研究證據的最佳手段,是循證決策的重要依據[1-3]。R 軟件是一款功能強大、免費開放的統計分析及繪圖軟件[4],metaplus 程序包是 R 專用于實現 robust meta 分析和 Meta 回歸的程序包[5]。metaplus 程序包由 Beath 等[6]研發,版本為 Version 0.7~9,該程序包是通過標準 Meta 分析法和基于輪廓似然置信區間的穩健法來實現的。本文以 Lawlor 等[7]《運動干預在抑郁癥治療中的有效性:基于隨機對照試驗的系統評價與 Meta 回歸分析》中的數據為例,介紹 metaplus 程序包的使用方法。
1 R 軟件及 metaplus 程序包的安裝及加載
R 軟件的簡介及安裝詳見前文[8]。本文以 R 3.4.1 版為例,待安裝完畢后,雙擊桌面的 R 軟件圖標,即可啟動 R 的交互式窗口(R-GUI)。于命令提示符“>”后輸入命令 install.packages("metaplus"),在彈出的對話框中選擇某個鏡像安裝(CRAN),安裝完成后可由 library("metaplus")命令完成加載。至此,R 軟件及 metaplus 程序包就已安裝及加載完畢。
2 數據加載
在進行數據加載之前須對數據進行處理,整理后的數據排列格式見表 1。該研究數據為運動與不運動比較、運動與認知治療比較的貝克抑郁量表評分的效應量和加權平均差的標準化均值差,并采用 Meta 分析和 Meta 回歸方法來驗證運動干預在抑郁癥治療中的有效性。

在 Excel 軟件中整理數據并以 csv 格式保存。導入數據語句為 mydata<-read.csv(file.choose()),輸入到 R 軟件后按回車鍵即可選擇文件夾位置,導入要分析的 csv 數據。數據導入后可使用命令 edit(mydata)對數據進行編輯和修改,還可使用命令 summary(madata)輸出基本的數據描述性信息。
3 數據分析
metaplus 程序包總體上可以分為 Meta 分析和 Meta 回歸兩大功能。其主要執行函數命令有 metaplus 函數、testOutliers 函數和 outlierProbs 函數。其中,metaplus 函數用于模型擬合、選擇合適模型及繪制森林圖等;testOutliers 及 outlierProbs 函數用于檢測數據中的離群值。下面以實例說明數據分析的過程。
3.1 模型的擬合與選擇
命令:mydata<-mydata[order(mydata$duration),]
3.1.1 標準正態隨機效應模型
命令:mydata.meta<-metaplus(smd, sesmd, mods=duration, slab=study, random="normal", data=mydata)
summary(mydata.meta)
命令中 smd 為觀察的效應量,sesmd 為效應量的標準誤,mods 為各研究相對應的協變量數據框架,slab 為與每個研究相對應的字符串向量(Vector of character strings corresponding to each study),random 為隨機效應的類型(normal 表示正態分布,t-dist 表示 t 分布,mixture 表示混合分布),data 為數據。上述數據的標準正態隨機效應模型 Meta 分析結果見表 2 和表 3。


3.1.2 t 分布隨機效應模型
metaplus 軟件包進行 t 分布隨機效應模型 Meta 分析與標準正態隨機效應模型類似,只需要將參數 random 賦值為“t-dist”即可,具體命令如下:
mydata.meta<-metaplus(smd, sesmd, mods=duration, slab=study, random="t-dist", data=mydata)
summary(mydata.meta)
上述數據的 t 分布隨機效應模型 Meta 分析結果見表 4 和表 5。


3.1.3 混合隨機效應模型
metaplus 進行混合隨機效應模型合并只需要參數 random 賦值為“mixture”即可,具體命令如下:
mydata.meta<-metaplus(smd, sesmd, mods=duration, slab=study, random="mixture", data=mydata)
summary(mydata.meta)
上述數據的混合隨機效應模型 Meta 分析結果見表 6 和表 7。


3.2 離群值檢驗
對于 t 分布隨機效應模型,離群值的檢驗需要自由度趨近為無窮大;對于混合隨機效應模型,離群值的檢驗需要離群研究概率密度函數的權重趨近為零。由于這兩個檢驗都涉及參數空間邊界上的參數,漸近理論不適用,所以在零假設下運用參數自舉的方法來檢驗統計量的經驗分布。然后將觀察到的統計量似然比與此分布進行比較,以確定 P 值。由于參數自舉中使用大量隨機生成的數據,P 值是不同的[9]。metaplus 提供了 testOutliers 函數進行離群值檢驗,命令如下:
mydata.testOutliers<- testOutliers(mydata.mix)
summary(mydata.testOutliers)
離群值檢驗的結果見表 8。

此外,metaplus 還提供 outliersProbs 用于計算離群概率,對于混合隨機效應模型和 t 分布隨機效應模型可以計算出每個研究的離群概率,并使用 plot 命令繪制離群概率圖像。命令如下:
mydata.outlierProbs<- outlierProbs(mydata.mix)
plot(mydata.outlierProbs)
離群概率圖結果見圖 1。

3.3 模型選擇
對于相關模型的擬合與選擇,總體上依據 AIC(赤池信息量準則)和 BIC(貝葉斯信息規則)的數值。AIC 和 BIC 都是對模型的擬合效果進行評價的指標,若 AIC 和 BIC 值越小,則模型對數據的擬合越好。當 t 分布隨機效應模型和混合隨機效應模型比標準正態隨機效應模型的 AIC 和 BIC 的值低,那么傾向于選擇 t 分布隨機效應模型和混合隨機效應模型。結合 testOutliers 和 outlierProbs 函數,可以對研究數據中的離群值進行檢測,若數據中存在明顯的離群值,則傾向于選擇混合隨機效應模型[9]。因本例中數據存在明顯的離群值,故選擇應用混合隨機效應模型。
3.4 輸出 Meta 分析結果
經過前期的模型擬合與選擇,并完成離群值檢驗后,metaplus 程序包可以調用 metafor 程序包的繪制森林圖功能,完成森林圖的繪制。命令如下:
plot(exercise.mix,extrameta=NULL)
森林圖結果見圖 2。

3.5 Meta 回歸
目前,metaplus 程序包尚未有預測方法,計算 4、8 和 12 周效應的替代方法是在這些時間將數據置于中心,并對每個數據進行 Meta 回歸[10]。命令如下:
mydata$duration4<-mydata$duration-4
mydata$duration8<-mydata$duration-8
mydata$duration12<-mydata$duration-12
mydata.nodurn<-metaplus(smd, sesmd, label = "Random Mixture (No Duration)", slab=study, random="mixture", data = mydata)
mydata.wk4<-metaplus (smd, sesmd, mods=duration4, label = "Random Mixture (Week 4)", slab= study, random = "mixture", data=mydata)
mydata.wk8<-metaplus (smd, sesmd, mods=duration8, label="Random Mixture (Week 8)", slab = study, random = "mixture", data = mydata)
mydata.wk12<-metaplus (smd, sesmd, mods=duration8, label="Random Mixture (Week 12)", slab=study, random="mixture", data=mydata)
plot (mydata.nodurn, extrameta=list (mydata.wk4, mydata.wk8, mydata.wk12), xlab="Standardised mean difference")
Meta 回歸森林圖見圖 3,Meta 回歸的截距可以估計對應的平均效應。圖 3 結果表明,研究效應在持續時間延長時會下降,可能與安慰劑效應減低或消失有關。

4 結果匯總
在完成上述運算命令后,即可對結果進行匯總和圖形的繪制。執行表 9 中的命令后,即可展示以上圖表結果。

5 模型的比較
在 Meta 分析中,最常見的兩種模型是固定效應模型和隨機效應模型。固定效應模型假設所有納入的研究擁有共同的真實效應量,而隨機效應模型中的真實效應量隨研究的不同而改變[11]。基于不同模型的運算,所得到的合并后的效應量的均數值也不相同[12]。盡管不一定適用于實踐,但為了數據處理的方便,傳統上假定隨機效應模型的真實效應量呈正態分布[13]。而 t 分布可以更靈活地來替代這個假設,通過馬爾科夫鏈蒙特卡羅方法實現正態分布和 t 分布的偏態擴張[14, 15],比常規正態分布更加靈活。混合隨機效應模型可以通過比較具有和不具有離群值的模型來確定是否存在異常值,并且使用后驗概率來確定,其優點在于 Meta 分析時考慮了含離群值的研究,且具有適當的權重[16]。
6 結語
隨著循證醫學的發展,為了適應不同類型 Meta 分析需求,多種 Meta 分析軟件被研發[17]。目前,國內有多篇介紹 Meta 分析相關軟件[18]的文章,包括 STATA 軟件[19]、RevMan 軟件[20]、SAS 軟件[21]、R 軟件[8]等。R 軟件作為一款免費的統計分析及繪圖軟件,受到廣大的歡迎。R 軟件中有 meta、metafor 等程序包可以實現 Meta 分析。本文以實例對 R 軟件 metaplus 程序包實現 Meta 分析的過程和結果進行了展示。
metaplus 程序包是 R 軟件中實現 Meta 分析和 Meta 回歸的程序包。metaplus 運用的是標準 Meta 分析法和基于輪廓似然置信區間的 robust 方法。metaplus 程序包的目的不是為了取代一般的 Meta 分析程序包,如 metafor 程序包[22],而是為了提供額外的專業分析,如離群值的檢測和離群值的概率計算等。metaplus 程序包可以繪制離群概率圖,并調用 Metafor 程序包繪制森林圖等結果圖形。然而,值得注意的是,metaplus 程序包的不足在于進行大量數據分析計算時會延長計算時間,特別是進行 testOutliers() 命令時。另一個不足之處是 metaplus 程序包只適用于標準正態隨機效應、t 分布隨機效應和混合正態隨機效應 3 種模型,沒有擴展到效應量的其他分布類型[9]。
總之,metaplus 程序包可以在實現基本 Meta 分析功能的基礎上進行離群值的檢測和概率計算。雖然 metaplus 程序包有些許不足,但是 R 軟件的開放特性也會使該程序包不斷地完善。
隨著循證醫學的發展,Meta 分析理論與方法學也隨著需求而不斷發展。目前 Meta 分析被公認為客觀評價和合成針對某一特定問題研究證據的最佳手段,是循證決策的重要依據[1-3]。R 軟件是一款功能強大、免費開放的統計分析及繪圖軟件[4],metaplus 程序包是 R 專用于實現 robust meta 分析和 Meta 回歸的程序包[5]。metaplus 程序包由 Beath 等[6]研發,版本為 Version 0.7~9,該程序包是通過標準 Meta 分析法和基于輪廓似然置信區間的穩健法來實現的。本文以 Lawlor 等[7]《運動干預在抑郁癥治療中的有效性:基于隨機對照試驗的系統評價與 Meta 回歸分析》中的數據為例,介紹 metaplus 程序包的使用方法。
1 R 軟件及 metaplus 程序包的安裝及加載
R 軟件的簡介及安裝詳見前文[8]。本文以 R 3.4.1 版為例,待安裝完畢后,雙擊桌面的 R 軟件圖標,即可啟動 R 的交互式窗口(R-GUI)。于命令提示符“>”后輸入命令 install.packages("metaplus"),在彈出的對話框中選擇某個鏡像安裝(CRAN),安裝完成后可由 library("metaplus")命令完成加載。至此,R 軟件及 metaplus 程序包就已安裝及加載完畢。
2 數據加載
在進行數據加載之前須對數據進行處理,整理后的數據排列格式見表 1。該研究數據為運動與不運動比較、運動與認知治療比較的貝克抑郁量表評分的效應量和加權平均差的標準化均值差,并采用 Meta 分析和 Meta 回歸方法來驗證運動干預在抑郁癥治療中的有效性。

在 Excel 軟件中整理數據并以 csv 格式保存。導入數據語句為 mydata<-read.csv(file.choose()),輸入到 R 軟件后按回車鍵即可選擇文件夾位置,導入要分析的 csv 數據。數據導入后可使用命令 edit(mydata)對數據進行編輯和修改,還可使用命令 summary(madata)輸出基本的數據描述性信息。
3 數據分析
metaplus 程序包總體上可以分為 Meta 分析和 Meta 回歸兩大功能。其主要執行函數命令有 metaplus 函數、testOutliers 函數和 outlierProbs 函數。其中,metaplus 函數用于模型擬合、選擇合適模型及繪制森林圖等;testOutliers 及 outlierProbs 函數用于檢測數據中的離群值。下面以實例說明數據分析的過程。
3.1 模型的擬合與選擇
命令:mydata<-mydata[order(mydata$duration),]
3.1.1 標準正態隨機效應模型
命令:mydata.meta<-metaplus(smd, sesmd, mods=duration, slab=study, random="normal", data=mydata)
summary(mydata.meta)
命令中 smd 為觀察的效應量,sesmd 為效應量的標準誤,mods 為各研究相對應的協變量數據框架,slab 為與每個研究相對應的字符串向量(Vector of character strings corresponding to each study),random 為隨機效應的類型(normal 表示正態分布,t-dist 表示 t 分布,mixture 表示混合分布),data 為數據。上述數據的標準正態隨機效應模型 Meta 分析結果見表 2 和表 3。


3.1.2 t 分布隨機效應模型
metaplus 軟件包進行 t 分布隨機效應模型 Meta 分析與標準正態隨機效應模型類似,只需要將參數 random 賦值為“t-dist”即可,具體命令如下:
mydata.meta<-metaplus(smd, sesmd, mods=duration, slab=study, random="t-dist", data=mydata)
summary(mydata.meta)
上述數據的 t 分布隨機效應模型 Meta 分析結果見表 4 和表 5。


3.1.3 混合隨機效應模型
metaplus 進行混合隨機效應模型合并只需要參數 random 賦值為“mixture”即可,具體命令如下:
mydata.meta<-metaplus(smd, sesmd, mods=duration, slab=study, random="mixture", data=mydata)
summary(mydata.meta)
上述數據的混合隨機效應模型 Meta 分析結果見表 6 和表 7。


3.2 離群值檢驗
對于 t 分布隨機效應模型,離群值的檢驗需要自由度趨近為無窮大;對于混合隨機效應模型,離群值的檢驗需要離群研究概率密度函數的權重趨近為零。由于這兩個檢驗都涉及參數空間邊界上的參數,漸近理論不適用,所以在零假設下運用參數自舉的方法來檢驗統計量的經驗分布。然后將觀察到的統計量似然比與此分布進行比較,以確定 P 值。由于參數自舉中使用大量隨機生成的數據,P 值是不同的[9]。metaplus 提供了 testOutliers 函數進行離群值檢驗,命令如下:
mydata.testOutliers<- testOutliers(mydata.mix)
summary(mydata.testOutliers)
離群值檢驗的結果見表 8。

此外,metaplus 還提供 outliersProbs 用于計算離群概率,對于混合隨機效應模型和 t 分布隨機效應模型可以計算出每個研究的離群概率,并使用 plot 命令繪制離群概率圖像。命令如下:
mydata.outlierProbs<- outlierProbs(mydata.mix)
plot(mydata.outlierProbs)
離群概率圖結果見圖 1。

3.3 模型選擇
對于相關模型的擬合與選擇,總體上依據 AIC(赤池信息量準則)和 BIC(貝葉斯信息規則)的數值。AIC 和 BIC 都是對模型的擬合效果進行評價的指標,若 AIC 和 BIC 值越小,則模型對數據的擬合越好。當 t 分布隨機效應模型和混合隨機效應模型比標準正態隨機效應模型的 AIC 和 BIC 的值低,那么傾向于選擇 t 分布隨機效應模型和混合隨機效應模型。結合 testOutliers 和 outlierProbs 函數,可以對研究數據中的離群值進行檢測,若數據中存在明顯的離群值,則傾向于選擇混合隨機效應模型[9]。因本例中數據存在明顯的離群值,故選擇應用混合隨機效應模型。
3.4 輸出 Meta 分析結果
經過前期的模型擬合與選擇,并完成離群值檢驗后,metaplus 程序包可以調用 metafor 程序包的繪制森林圖功能,完成森林圖的繪制。命令如下:
plot(exercise.mix,extrameta=NULL)
森林圖結果見圖 2。

3.5 Meta 回歸
目前,metaplus 程序包尚未有預測方法,計算 4、8 和 12 周效應的替代方法是在這些時間將數據置于中心,并對每個數據進行 Meta 回歸[10]。命令如下:
mydata$duration4<-mydata$duration-4
mydata$duration8<-mydata$duration-8
mydata$duration12<-mydata$duration-12
mydata.nodurn<-metaplus(smd, sesmd, label = "Random Mixture (No Duration)", slab=study, random="mixture", data = mydata)
mydata.wk4<-metaplus (smd, sesmd, mods=duration4, label = "Random Mixture (Week 4)", slab= study, random = "mixture", data=mydata)
mydata.wk8<-metaplus (smd, sesmd, mods=duration8, label="Random Mixture (Week 8)", slab = study, random = "mixture", data = mydata)
mydata.wk12<-metaplus (smd, sesmd, mods=duration8, label="Random Mixture (Week 12)", slab=study, random="mixture", data=mydata)
plot (mydata.nodurn, extrameta=list (mydata.wk4, mydata.wk8, mydata.wk12), xlab="Standardised mean difference")
Meta 回歸森林圖見圖 3,Meta 回歸的截距可以估計對應的平均效應。圖 3 結果表明,研究效應在持續時間延長時會下降,可能與安慰劑效應減低或消失有關。

4 結果匯總
在完成上述運算命令后,即可對結果進行匯總和圖形的繪制。執行表 9 中的命令后,即可展示以上圖表結果。

5 模型的比較
在 Meta 分析中,最常見的兩種模型是固定效應模型和隨機效應模型。固定效應模型假設所有納入的研究擁有共同的真實效應量,而隨機效應模型中的真實效應量隨研究的不同而改變[11]。基于不同模型的運算,所得到的合并后的效應量的均數值也不相同[12]。盡管不一定適用于實踐,但為了數據處理的方便,傳統上假定隨機效應模型的真實效應量呈正態分布[13]。而 t 分布可以更靈活地來替代這個假設,通過馬爾科夫鏈蒙特卡羅方法實現正態分布和 t 分布的偏態擴張[14, 15],比常規正態分布更加靈活。混合隨機效應模型可以通過比較具有和不具有離群值的模型來確定是否存在異常值,并且使用后驗概率來確定,其優點在于 Meta 分析時考慮了含離群值的研究,且具有適當的權重[16]。
6 結語
隨著循證醫學的發展,為了適應不同類型 Meta 分析需求,多種 Meta 分析軟件被研發[17]。目前,國內有多篇介紹 Meta 分析相關軟件[18]的文章,包括 STATA 軟件[19]、RevMan 軟件[20]、SAS 軟件[21]、R 軟件[8]等。R 軟件作為一款免費的統計分析及繪圖軟件,受到廣大的歡迎。R 軟件中有 meta、metafor 等程序包可以實現 Meta 分析。本文以實例對 R 軟件 metaplus 程序包實現 Meta 分析的過程和結果進行了展示。
metaplus 程序包是 R 軟件中實現 Meta 分析和 Meta 回歸的程序包。metaplus 運用的是標準 Meta 分析法和基于輪廓似然置信區間的 robust 方法。metaplus 程序包的目的不是為了取代一般的 Meta 分析程序包,如 metafor 程序包[22],而是為了提供額外的專業分析,如離群值的檢測和離群值的概率計算等。metaplus 程序包可以繪制離群概率圖,并調用 Metafor 程序包繪制森林圖等結果圖形。然而,值得注意的是,metaplus 程序包的不足在于進行大量數據分析計算時會延長計算時間,特別是進行 testOutliers() 命令時。另一個不足之處是 metaplus 程序包只適用于標準正態隨機效應、t 分布隨機效應和混合正態隨機效應 3 種模型,沒有擴展到效應量的其他分布類型[9]。
總之,metaplus 程序包可以在實現基本 Meta 分析功能的基礎上進行離群值的檢測和概率計算。雖然 metaplus 程序包有些許不足,但是 R 軟件的開放特性也會使該程序包不斷地完善。