多变量Meta分析在R语言mvmeta包的实现

引用本文

姚宁宁, 陈炳为, 黄灏, 钱刘兰. 多变量Meta分析在R语言mvmeta包的实现.循证医学, 2015,15(6): 376-379
YAO Ning-ning, CHEN Bing-wei, HUANG Hao, QIAN Liu-lan. The Implementation of Multivariate Meta-Analysis by mvmeta Package in R Language. Journal of Evidence-based Medicine, 2015,15(6): 376-379 复制到剪切板

Permissions

多变量Meta分析在R语言mvmeta包的实现

姚宁宁, 陈炳为, 黄灏, 钱刘兰

东南大学公共卫生学院流行病与卫生统计学系, 南京 210009

通信作者：陈炳为,Tel:025-83272562; E-mail:drchenbw @126.com

作者简介：姚宁宁(1988-),女,山东临沂人,在读硕士研究生,研究方向为卫生统计学。

摘要

目的介绍多变量Meta分析的原理及其在R语言mvmeta包的实现。方法以β-受体阻滞剂、硬化疗法及对照组三种非手术治疗方式预防肝硬化出血为例,分别计算两种疗法与对照组的初次出血危险性的危险度(odds ratio, OR)值,将ln(OR)值视为多变量,然后用R语言mvmeta包实现多变量Meta分析。结果多变量Meta分析的结果为,β-受体阻滞剂。硬化疗法与对照组的OR值及95%可信区间分别为0.508 9(0.270 5,0.957 1)、0.552 2(0.318 4,0.957 7),说明两种疗法与对照组比较的初次出血危险性均小于1,差异具有统计学意义。研究间相关系数为0.436 5。结论当Meta分析的结局为多变量或多组时,多变量Meta分析可同时考虑多个结局变量以及研究间结局变量之间的相关,使得分析更加合理,结果更加准确和可信。

关键词: 多变量Meta分析; R语言; mvmeta包

中图分类号:R195.1 文献标识码:A 收稿日期:2015-01-06

The Implementation of Multivariate Meta-Analysis by mvmeta Package in R Language

YAO Ning-ning, CHEN Bing-wei, HUANG Hao, QIAN Liu-lan

Department of Epidemiology and Statistics, School of Public Health, Southeast University, Nanjing 210009, China

Abstract

Objective To introduce the principle of the multivariate meta-analysis and its’ implementation by the mvmeta package in R language.Methods The randomized trials of non-surgical treatment of first bleeding in cirrhosis using beta-blockers, sclerotherapy and control group was used as an example. The log-scale of first bleeding risk of OR values that two therapies compared control group were calculated separately and regarded as multivariate, and then multivariate meta-analysis was used by mvmeta package in R language.Results Multivariate meta-analysis showed: OR and their 95% confidence interval of beta-blockers, sclerotherapy compared the control group were 0.508 9 (0.270 5, 0.957 1) and 0.552 2(0.318 4, 0.957 7) respectively. The first bleeding risk of the two therapies are less than 1. The difference is statistically significant. The correlation of between-study is 0.436 5.Conclusion When meta-analysis has the multiple endpoints or groups, multivariate meta-analysis takes multiple endpoints and the between-study correlation into account simultaneously, making the progress more reasonable and the result of multivariate meta-analysis more accurate and credible.

Key words: multivariate meta-analysis; R language; mvmeta package

Show Figures

Meta分析是对相同研究目的的多个独立研究的结果进行综合统计分析。定量综合的一种研究方法^[1]。通常Meta分析假定效应量来自于独立的研究, 因此统计结果也是独立的。然而, 许多研究不能满足独立性的假设, 比如多个治疗组与一个共同的对照组比较的研究和多个结局变量的研究就可能产生效量之间的相关^[2]。多变量Meta分析(multivariate meta-analysis)作为单变量Meta分析的一个拓展, 可合并估计多个研究的多个相关参数, 这些参数可以是多个结局或多组间的比较^[3]。当同一总体中的测量结局相关时, 分别对每个结局进行Meta分析, 测量结局之间的相关结构就可能被忽略^[4]。统计推断和模拟研究显示:研究内相关可能会通过“ 借力” 影响终点结局, 忽略它可能会使Meta分析的结果存在偏差, 如合并效应量均方和标准误的估计过高^[5]。

多变量Meta分析在随机对照研究中有多种应用, 最简单的是在临床试验中把每个组的结局分别处理, 其他的应用还有同时探索两个临床结局的治疗效应, 或同时探索成本效益的治疗效应, 比较多个治疗的联合试验, 以及在观察性研究中评估暴露量与疾病之间的相关性, 还有在诊断试验和网络干预中的应用^{[3, 4, 6]}。

现在多变量Meta分析变得越来越流行, 在许多统计软件中都含有正式的程序或者自动程序化的函数^[4]。目前国内Meta分析主要考虑单变量的形式, 本文的目的在于介绍多变量Meta分析的原理及其在R语言mvmeta包的实现。

1 原理和方法

假定有n个研究(i=1, 2, …, n), p个结局变量(j=1, 2, …, p), 每一个研究y_i=(y_i1, y_i2, …, y_ij, …, y_ip)’ 。为了简化, 先从两个结局变量介绍, 之后再拓展到p个结局变量。在固定效应的背景下, 假定y_i服从多元正态分布:

$\begin{matrix} (\begin{matrix} y_{i 1} \\ y_{i 2} \end{matrix}) ~ MVN ((\begin{matrix} μ_{1} \\ μ_{2} \end{matrix}), (\begin{matrix} σ_{i 1}^{2} & ρ_{i} σ_{i 1} σ_{i 2} \\ ρ_{i} σ_{i 1} σ_{i 2} & σ_{i 2}^{2} \end{matrix})) \end{matrix}$ (1)

式中, $\begin{matrix} ρ_{i} \end{matrix}$ 是第i个研究内两个结局的相关系数, $\begin{matrix} μ = \end{matrix}$ ( $\begin{matrix} μ_{1} \end{matrix}$ , $\begin{matrix} μ_{2} \end{matrix}$ )’ 是每个结局的均值向量, $\begin{matrix} S_{i} = (\begin{matrix} σ_{i 1}^{2} & ρ_{i} σ_{i 1} σ_{i 2} \\ ρ_{i} σ_{i 1} σ_{i 2} & σ_{i 2}^{2} \end{matrix}) \end{matrix}$ 是研究内的协方差矩阵。总变异可分解为研究内变异和研究间变异, 总相关系数也可分解为研究内相关系数 $\begin{matrix} ρ_{i} \end{matrix}$ 和研究间相关系数 $\begin{matrix} ρ_{τ} \end{matrix}$ 。

$\begin{matrix} (\begin{matrix} y_{i 1} \\ y_{i 2} \end{matrix}) | (\begin{matrix} θ_{i 1} \\ θ_{i 2} \end{matrix}) ~ MVN ((\begin{matrix} θ_{i 1} \\ θ_{i 2} \end{matrix}), (\begin{matrix} σ_{i 1}^{2} & ρ_{i} σ_{i 1} σ_{i 2} \\ ρ_{i} σ_{i 1} σ_{i 2} & σ_{i 2}^{2} \end{matrix})) \end{matrix}$ (2)

$\begin{matrix} θ = (θ_{i 1}, θ_{i 2}) ’ \end{matrix}$ 是研究的特定效应, 假定也服从正态分布:

$\begin{matrix} (\begin{matrix} θ_{i 1} \\ θ_{i 2} \end{matrix}) ~ MVN ((\begin{matrix} μ_{1} \\ μ_{2} \end{matrix}), (\begin{matrix} τ_{1}^{2} & ρ_{τ} τ_{1} τ_{2} \\ ρ_{τ} τ_{1} τ_{2} & τ_{2}^{2} \end{matrix})) \end{matrix}$ (3)

其中 $\begin{matrix} τ_{j} \end{matrix}$ 是效应量j的研究间变异(异质性), 研究间协方差矩阵为 $\begin{matrix} Δ = (\begin{matrix} τ_{1}^{2} & ρ_{τ} τ_{1} τ_{2} \\ ρ_{τ} τ_{1} τ_{2} & τ_{2}^{2} \end{matrix}) \end{matrix}$ , 它是未知的, $\begin{matrix} ρ_{τ} \end{matrix}$ 为不同研究间的两个效应量的相关系数。合并公式(2)和(3)可得:

$\begin{matrix} (\begin{matrix} y_{i 1} \\ y_{i 2} \end{matrix}) ~ MVN ((\begin{matrix} μ_{1} \\ μ_{2} \end{matrix}), (\begin{matrix} σ_{i 1}^{2} + τ_{1}^{2} & ρ_{i} σ_{i 1} σ_{i 2} + ρ_{τ} τ_{1} τ_{2} \\ ρ_{i} σ_{i 1} σ_{i 2} + ρ_{τ} τ_{1} τ_{2} & σ_{i 2}^{2} + τ_{2}^{2} \end{matrix})) \end{matrix}$ (4)

当有p个结局时, 可以拓展到:

进一步可表达为:y_i~ (μ , Δ + S)。 $\begin{matrix} Δ \end{matrix}$ 和 $\begin{matrix} S \end{matrix}$ 是p× p的(协)方差矩阵。应用广义最小二乘法(generalized least square, GLS)回归估计 $\begin{matrix} μ \end{matrix}$ 和 $\begin{matrix} Δ \end{matrix}$ 。固定效应的Meta分析模型, $\begin{matrix} Δ \end{matrix}$ 是不存在的, 研究间的变异仅仅由于研究内估计误差^[7]; 随机效应Meta分析难点在于估计研究间的方差^[3]。

2 实例分析

2.1 数据的选择和整理

数据来源于Lu和Ades^[8]的文献, 主要研究肝硬化的非手术治疗方式预防其出血的危险性, 以初次出血的例数为指标, 其中三个组分别是:β -受体阻滞剂(A), 硬化疗法(B), 对照组(C), 目的是评价这三种非手术治疗方式预防肝硬化出血的效果。

对原始数据进行处理, Bled表示初次出血的例数, Total表示干预组的总例数。Y_AC和Y_BC分别表示A、B两组相对于C组估计的ln(OR), 即干预组的肝硬化初次出血的危险性是对照组的倍数的自然对数; S_AA、S_BB和S_AB则表示其对应方差及两者之间的协方差。对于包含0的研究(研究10和研究20), 在每个组增加0.5个初次出血的例数。整理后见表1。

表1 肝硬化非手术治疗初次出血的随机试验

2.2 多变量Meta分析在R语言mvmeta包的实现

2.2.1 数据的读入

读入Excel数据:先把肝硬化初次出血整理后的数据集cirrhosis(至少包含Y_AC、Y_BC、S_AA、S_AB、S_BB变量)保存为csv格式, 然后利用下面命令将其导入R语言。先调用foreign包, 然后再把CSV格式的数据集读入。

> library(foreign)

> data< -read.csv (“ c:\cirrhosis.csv” )

2.2.2 R语言mvmeta程序包的加载和调用

R语言mvmeta包包含一系列函数可执行固定和随机效应的多变量和单变量的Meta分析和Meta回归。首先要对程序包进行加载:install.packages(‘ mvmeta’ )。然后调用mvmeta包:library(mvmeta)。

mvmeta的语句:mvmeta(formula, S, data, subset, method=“ reml” , bscov=“ unstr” , model=TRUE, contrasts=NULL, offset, na.action, control=list( ))

其中formula表示结局变量名称(即Y_AC、Y_BC); S表示研究内(协)方差(即S_AA、S_AB、S_BB); data表示数据集名称; method表示所用的估计方法:固定效应模型时选择FIXED; 随机效应模型时则选择限制性最大似然估计(REML)。最大似然估计(ML)、矩估计(MM)、方差成分法(VC)的其中之一, 默认为REML。由输出结果中Q检验的P值和I²统计量来判断异质性以及选择何种效应模型。

mvmeta包中主要提供了多变量Meta分析与多变量的Meta回归, 另外也提供了单变量的Meta分析和Meta回归。但对于后两者, 在R语言中的metafor、meta、rmeta及metalik等包提供了更多、更详尽和有效的功能、多变量Meta程序为library(mvmeta), 调用mvmeta软件包。

模型1:多变量Meta分析的程序

model< -mvmeta(cbind(Ya, Yb), S=S, data=cirrhosis)

模型2:多变量的Meta回归

model < - mvmeta(cbind(Ya, Yb)~X, S=S, data=cirrhosis), 此处X代表协变量。

模型3:单变量Meta分析

model< -mvmeta (Y, S=S, data=cirrhosis), 此处Y为单变量的效应量, S为效应量方差。

模型4:单变量Meta回归

model< -mvmeta (Y~X, S=S, data=cirrhosis), 此处X代表协变量。

运行以上程序后, 最后将结果输出。

2.3 结果及分析

单变量和多变量Meta分析都是采用ln(OR)值做分析、单变量Meta分析时Y_AC和Y_BC的Q检验P值均小于0.05, I²统计量分别为57.7%和77.8%。多变量Meta分析Q检验P< 0.05, I²统计量为73.9%。可知两种Meta分析均存在异质性, 都用随机效应模型、估计方法选择默认的REML法。

2.3.1 单变量Meta分析

表2是单变量Meta分析结果, 可得:AC与BC的OR值及95%可信区间分别为0.528 1(0.280 2, 0.995 5)、0.540 6(0.309 5, 0.944 3), 表明初次出血的危险性由于干预而降低, 即β -受体阻滞剂、硬化疗法可以预防肝硬化出血, 两者为保护因素。

表2 单变量Meta分析结果

2.3.2 多变量Meta分析

多变量Meta分析的结果:Y_AC 为-0.675 5(-1.307 3, -0.043 8), Y_BC为-0.593 8(-1.144 4, -0.043 2), 研究间相关系数为0.436 5(见表3), A组与B组的治疗效果呈正相关。OR值及95%可信区间分别为0.508 9(0.270 5, 0.957 1)、0.552 2(0.318 4, 0.957 7), 多变量Meta分析的结果说明β -受体阻滞剂预防肝硬化出血的效果是最好, 其次是硬化疗法。OR值的95%可信区间不包含1, 上下限均小于1, 说明两种疗法与对照组比较的初次出血危险性均小于1, 差异有统计学意义。

表3 多变量Meta分析结果

3 讨论

Lu和Ades对该数据在WinBUGS中采用贝叶斯分层模型进行混合治疗比较的分析, 结果d_AB的均值为-0.649, 它是A和B的ln(OR)的平均治疗效应^[8]。与Lu和Ades 的结果相比, 本研究中多变量Meta分析的结果(d_AB为-0.634 7)比单变量Meta分析结果的(d_AB为-0.626 7)更为接近。表1中较多的研究并不是完整地设置了3个组, 此时, 多变量Meta分析是先构造协方差矩阵, 然后根据协方差矩阵来估计参数。因此, 当部分研究不完整时, 它先是利用完整部分的数据计算相关系数或协方差矩阵, 然后再将这种相关性应用到整个研究中^[9]。缺失组别对结果的影响应考虑缺失比例及完整数据的代表性, 如果完整数据的代表性不好, 则可能会对结果造成影响。

多变量Meta分析的优势在于可以描述和利用多个效应的关系, 获得具有较好统计特征的参数估计并减少偏倚等; 不足之处是研究内相关未知, 在这种情况下, 可以用近似公式。替代模型、稳健的方差估计等来解决^[4]。尽管多变量Meta分析很早就有研究, 但由于研究者的习惯, 缺乏对研究效应相关结构的理解, 低估了忽略相关的影响等原因导致它还不能普及^{[4, 6]}、R、STATA、SAS均可以实现多变量Meta分析。STATA的mvmeta命令的原理和用法与R语言mvmeta包相同; SAS中Proc Mixed采用广义线性混合模型也可实现多变量Meta分析, 但程序比R和STATA复杂, 且不易理解、不同的多变量Meta分析实现软件的结果呈现的差异也需要做更深入的研究。近几年多变量Meta分析处于发展时期, 研究也越来越多, 但多变量Meta分析的普及仍然还有一段时间^[6]。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	孙振球. 医学统计学[M]. 第3版. 北京：人民卫生出版社， 2010： 674. [本文引用:1]
[2]	Gloudemans JA, Owens CM, Kromrey JD. MV_META: A SAS® macro for multivariate Meta-analysis. SAS Global Forum 2012. http://www.doc88.com/p-5681646124913.html [本文引用:1]
[3]	White IR. Multivariate rand om-effects Meta-analysis[J]. Stata J, 2009, 9(1): 40-56. [本文引用:3]
[4]	Mavridis D, Salanti G. A practical introduction to multivariate Meta-analysis[J]. Stat Methods Med Res, 2012, 22(2): 133-158. [本文引用:5]
[5]	R within-study correlation[J]. J Roy Stat, 2009, 172(4): 789-811. [本文引用:1]
[6]	otential and promise[J]. Stat Med, 2011, 30(20): 2481-2498. [本文引用:3]
[7]	Gasparrini A, Package ‘mvMeta’. Multivariate Meta-analysis and Meta⁃regression. http://cran.r-project.org/web/packages/mvmeta/mvmeta.pdf [本文引用:1]
[8]	Lu G, Ades AE. Combination of direct and indirect evidence in mixed treatment comparisons[J]. Stat Med, 2004, 23(20): 3105-3124. [本文引用:2]
[9]	Jackson D， White IR， Riley RD. A matrix⁃based method of moments for fitting the multivariate rand om effects model for Meta--analysis and Meta-regression[J]. Biom J, 2013, 55(2), 231-245. [本文引用:1]

2010

0.0

... 定量综合的一种研究方法^[1] ...

2012

0.0

... 然而,许多研究不能满足独立性的假设,比如多个治疗组与一个共同的对照组比较的研究和多个结局变量的研究就可能产生效量之间的相关^[2] ...

2009

0.0

... 多变量Meta分析(multivariate meta-analysis)作为单变量Meta分析的一个拓展,可合并估计多个研究的多个相关参数,这些参数可以是多个结局或多组间的比较^[3] ...

... 多变量Meta分析在随机对照研究中有多种应用,最简单的是在临床试验中把每个组的结局分别处理,其他的应用还有同时探索两个临床结局的治疗效应,或同时探索成本效益的治疗效应,比较多个治疗的联合试验,以及在观察性研究中评估暴露量与疾病之间的相关性,还有在诊断试验和网络干预中的应用^[3,4,6] ...

... 随机效应Meta分析难点在于估计研究间的方差^[3] ...

2012

0.0

... 当同一总体中的测量结局相关时,分别对每个结局进行Meta分析,测量结局之间的相关结构就可能被忽略^[4] ...

... 现在多变量Meta分析变得越来越流行,在许多统计软件中都含有正式的程序或者自动程序化的函数^[4] ...

... 替代模型、稳健的方差估计等来解决^[4] ...

... 尽管多变量Meta分析很早就有研究,但由于研究者的习惯,缺乏对研究效应相关结构的理解,低估了忽略相关的影响等原因导致它还不能普及^[4,6]、R、STATA、SAS均可以实现多变量Meta分析 ...

2009

0.0

... 影响终点结局,忽略它可能会使Meta分析的结果存在偏差,如合并效应量均方和标准误的估计过高^[5] ...

2011

0.0

... 近几年多变量Meta分析处于发展时期,研究也越来越多,但多变量Meta分析的普及仍然还有一段时间^[6] ...

0.0

... 固定效应的Meta分析模型, Δ是不存在的,研究间的变异仅仅由于研究内估计误差^[7] ...

2004

0.0

... 1 数据的选择和整理数据来源于Lu和Ades^[8]的文献,主要研究肝硬化的非手术治疗方式预防其出血的危险性,以初次出血的例数为指标,其中三个组分别是:#cod#x003b2 ...

... 649, 它是A和B的ln(OR)的平均治疗效应^[8] ...

2013

0.0

... 因此,当部分研究不完整时,它先是利用完整部分的数据计算相关系数或协方差矩阵,然后再将这种相关性应用到整个研究中^[9] ...