應用R軟件的robumeta程序包實現基于魯棒方差估計的Meta回歸分析

馮雨嘉，郭梓鑫，曾憲濤，劉小平

本文是基于魯棒方差估計的方法學介紹[1]。魯棒方差估計（robust variance estimation，RVE）是一種處理非獨立效應量的方法，由Hedges，Tipton和Johnson 提出[1]，廣泛的應用于物理學，醫學，生物學等方面[2,3]。在通常使用的多元Meta分析中，要求組內變異是已知的，而使用RVE方法，即使組內變異未知，也可以提供有效的標準誤，點估計和可信區間[4]。

R軟件擁有強大的Meta分析功能，robumeta程序包是R軟件基于魯棒方差估計實現Meta回歸分析的程序包。robumeta程序包由Zachary Fisher，ElizabethTipton，HouZhipeng研發[5]，本文以Oswald等發表的文章中的oswald2013.ex1組數據（附件）為例[6]，演示該程序包的使用方法。

1 R軟件及robumeta程序包的安裝

本文使用R軟件版本為R-3.4.1[7]，可在R語言官方網站（https://www.r-project.org/）下載最新版R軟件，安裝可參閱《R軟件Metafor程序包在Meta分析中的應用》[8]一文。在R軟件安裝完畢后，雙擊桌面的R軟件圖標，即可啟動R軟件。于命令提示符“>”后輸入命令install.packages("robumeta")，在彈出的對話框中選擇某個鏡像安裝(CRAN)，安裝完成后可由library("robumeta") 命令完成加載。至此，軟件及robumeta程序包[5]已安裝完畢。

2 數據加載

使用程序包自帶的數據oswald2013.ex1[6]進行演示，該數據集是oswald2013數據集[6]的一個子集，oswald2013數據集包含從46個個體研究中采集的308個效應量，oswald2013.ex1數據集包含從9個個體研究中采集的32個效應量，該數據集為內隱聯想測驗（Implicit attitude Test，IAT）得出的數據。加載數據的指令為data(oswald2013.ex1)，具體數據參見附表1。

3 原理闡釋

如前所述，本文為基于RVE方法學介紹[1]，RVE是一種處理非獨立效應量的方法，由Hedges,Tipton和Johnson 提出[1]，現被廣泛的運用于醫學，物理學，生物學方面[2,3]。最近，Tipton等對傳統的RVE進行了一些調整[9]，很好的提升了RVE在小樣本中的適用性，當研究的數目小于40時，這樣的調整顯得尤其重要。下面將對RVE的原理進行簡單的介紹。

RVE的基本原理公式[1,3]如下：

其中j=1，……，m,（研究總數為m，），T表示效應量估計，X=（X’1，……，X’m）表示設計的矩陣，ε表示殘差，β表示未知的回歸系數，W表示加權的分塊對角矩陣，b表示用最小二乘法估計的β值。

4 數據分析

在完成數據加載之后，即可開始數據的分析。robumeta程序包的功能主要分為魯棒方差Meta回歸模型擬合，敏感性分析，繪制相應的森林圖。以下將展示具體命令及相關說明。

4.1 魯棒方差Meta回歸模型擬合 魯棒方差Meta回歸模型擬合主要依靠robu()命令實現[1]robu()命令及其結果（表1）如下：data(oswald2013.ex1)

oswald_intercept<-robu(formula=effect.size～1,data=oswald2013.ex1,studynum=Study,var.eff.size=var.eff.size,rho=0.8,small = TRUE)

表1 robu()命令執行后的結果

命令中：data(oswald2013.ex1)表示加載robumeta包自帶的數據oswald2013.ex1。formula可以理解為一個常規格式的線性模型，它的形式大致為y ～ x1 + x2...， y為效應量，x1 + x2...為用戶指定的協變量，在intercept-only模型中，可以指定y～1，這是可選擇的，不是固定的，因而，式中的effect.size ～ 1可以理解為y～1。effect.size為一列從原始數據計算得到的，用戶指定的效應量。formula = effect.size～1表示選用的線性模型為effect.size～1。data = oswald2013.ex1表示使用的數據為oswald2013.ex1，studynum=Study表示選用Study列作為識別研究的特殊值。var.eff.size=var.eff.size表示用戶指定的反映每個效應量組內變異的一列值。rho的缺省值為0.8，用戶指定的值必須在0～1之間。small=TRUE表示使用小樣本調整程序，此命令的缺省值為TURE。

4.2 敏感性分析 敏感性分析主要靠sensitivty()命令實現

sensitivty()命令及其結果（表2）如下：

sensitivity(oswald_intercept)

表2 Sensitivity()命令執行后的結果

若要對Meta回歸模型進行有效的加權，則需要建立組間變異的估計值τ2，而要對τ2值進行點估計，則需要建立ρ，ρ是效應量之間常見的相關關系。靈敏性分析可用來判定ρ對τ2的效應大小。大量的研究表明，τ2值和Meta回歸系數通常對ρ的變化不敏感[10]，但當組間變異明顯小于組內變異時，以上結論可能不成立。

4.3 生成研究內部和研究間的協變量 生成研究內部和研究間的協變量主要靠group.mean()和group.center()命令實現

group.mean()和group.center()命令及其結果（表3～4）如下：

h i e r d a t$f o l l o w u p_m<-g r o u p.mean(hierdat$followup, hierdat$studyid)

h i e r d a t$f o l l o w u p_c<-g r o u p.center(hierdat$followup, hierdat$studyid)

命令中：hierdat$followup表示包含組均值的協變量，hierdat$studyid表示一組用于計算組均值的數據。group.mean()是用于計算組間變異的函數，group.center()是用于計算組內變異的函數。

4.4 Meta回歸分析 進行完以上的數據準備之后，就可以開始進行Meta回歸分析了，具體的命令和結果（表5）如下：

model_hier<-

wup_c+followup_m+binge,data=hierdat,modelw eights="HIER",studynum=studyid, var.eff.size = var,small = TRUE)

命令中：modelweights="HIER"表示選用分層效應的加權模型，另一種模型為“CORR”，表示相關效應模型。命令中其他部分所表示的含義與robu()命令中各項的含義一致。

表3 group.center()命令執行后的結果

表4 group.mean()命令執行后的結果

表5 Meta回歸分析的結果

4.5 輸出Meta分析結果 在以上的前期計算準備完成后，就可以進行森林圖的繪制了。繪制的命令如下：

forest.robu(oswald_intercept,es.lab="Crit.Cat",study.lab="Study","EffectSize"= effect.size,# optional column "Weight" = r.weights) # optional column

命令中：oswald_intercept表示使用的模型為intercept-only模型，es.lab="Crit.Cat"表示用戶指定的個體效應量例如“數學分數”，“物理分數”等，這里指定的個體效應量為"Crit.Cat"。study.lab="Study"表示用戶指定的一列群組效應量例如“作者姓名”，“出版年份”等，這里指定的群組效應量為"Study"。"Effect Size"= effect.size表示一列效應量，這里使用的這列效應量的名稱是effect.size。"Weight" = r.weights表示使用一個名為"Weight"的RVE加權。"Effect Size"= effect.size和"Weight" = r.weights是用來闡明在命令中加入另外的列所需遵循的規則的。此命令的使用需要grid 程序包。所繪制的圖形如圖1。

圖1 robumeta程序包繪制的森林圖

5 結果匯總

在完成上述運算命令后，即可對結果進行匯總和圖形的繪制。執行表6所示命令則可查看表1～表5的結果。若要將結果導出，可采用如下命令(以group.center()命令為例)：

write.csv(hierdat$followup_c,"D:/data/hierdat$followup_c.csv")

上述命令中：hierdat$followup_c表示要輸出的結果的名稱，"D:/data/ hierdat$followup_c.csv"表示文件導出的位置和文件名稱。此命令導出的結果為csv文件，若要保存為xls或xslx格式的文件，則用excel軟件打開后保存為所需格式即可。

表6 執行結果匯總的命令

6 結語

R作為當前最受歡迎的軟件之一，當前國內尚無介紹robumeta程序包的文章，本文以實例對R軟件robumeta程序包實現基于魯棒方差估計的Meta回歸分析的過程及結果進行了全面展示。

魯棒方差估計提供了一種處理非獨立效應量的方法[1]，但傳統的魯棒方差估計只適用于大樣本的分析，而最近對傳統的魯棒方差估計進行改進后，它同時也適用于小樣本的分析[9]。robumeta程序包正是基于魯棒方差估計的原理，采用了改良的魯棒方差估計方法，實現了利用魯棒方差估計進行的Meta回歸分析。

總之，robumeta程序包可以滿足常大樣本和小樣本的魯棒方差估計和Meta回歸分析，且 R 的開放特性也決定了該程序包會在后期不斷完善。

[1] Hedges LV,Tipton E,Johnson MC. Johnson.Robust variance estimation in meta-regression with dependent effect size estimates[J].Res Synth Methods,2010,1(1):39-65.

[2] Mente S,Kuhn M. The use of the R language for medicinal chemistry applications[J]. Curr Top Med Chem,2012,2(18):1957-64.

[3] Tipton E. Robust variance estimation in meta-regression with binary dependent effects[J]. Res Synth Methods,2013,4(2):169-87.

[4] Jackson D,Riley R,White IR. Multivariate meta-analysis:potential and promise[J]. Stat Med, 2011,30(20):2481-98.

[5] Freeman JM,Vining EP,Kossoff EH,et al. A blinded,crossover study of the efficacy of the ketogenic diet[J]. Epilepsia,2009,50(2):322-5.

[6] Oswald FL,Mitchell G,Blanton H,et al. Predicting ethnic and racial discrimination:a meta-analysis of IAT criterion studies[J]. J Pers Soc Psychol,2013,105(2):171-92.

[7] R Core Team (2015). R:A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.URL https://www.R-project.org/.

[8] 董圣杰,曾憲濤,郭毅. R軟件Metafor程序包在Meta分析中的應用[J]. 中國循證醫學雜志,2012(09):1141-7.

[9] Tipton E. Small sample adjustments for robust variance estimation with meta-regression[J]. Psychol Methods,2015,20(3):375-93.

[10] Ishak KJ,Platt RW,Joseph L,et al. Impact of approximating or ignoring within-study covariances in multivariate meta-analyses[J]. Stat Med,2008,27(5):670-86.