基于中位数为效应量的Meta分析及R实现

黄 桥，邝心颖，童铁军，毛 智，郑 亮，陈 昊，贺小峰，董圣杰，张 超，仇成凤，翁 鸿，桂裕亮，任学群＊，曾宪涛,,＊

1.武汉大学中南医院循证与转化医学中心（武汉 430071）

2.武汉大学第二临床学院循证医学与临床流行病学教研室（武汉 430071）

3.日本圣路加国际大学看护学研究科全球健康看护系（日本东京 1040045）

4.香港浸会大学数学系（香港 999077）

5.中国人民解放军总医院重症医学科（北京 100853）

6.同济大学附属东方医院转化医学中心（上海 200120）

7.南京中医药大学针灸推拿学院针灸临床教研室（南京 210029）

8.长治医学院附属和平医院科教科（山西长治 046000）

9.长治医学院第一临床学院教学科研科（山西长治 046000）

10.烟台市烟台山医院骨关节科（山东烟台 264001）

11.十堰市太和医院（湖北医药学院附属）循证医学中心（湖北十堰 442000）

12.南华大学附属怀化市第一人民医院循证医学与临床研究中心（湖南怀化 418000）

13.河南大学循证医学与临床转化研究院（河南开封 475000）

14.武汉大学中南医院人事处（武汉 430071）

15.武汉大学全球健康研究中心（武汉 430072）

连续变量是临床研究中非常常见的资料类型，如血压和空腹血糖测量值。对于连续变量的集中趋势和变异趋势进行描述统计学时，正态分布的数据常采用均数±标准差进行描述，偏态分布数据常采用中位数（四分位间距）或中位数（第一四分位数～第三四分位数）进行描述[1]。一般情况下，当文献中汇报效应量为中位数时，表明文献的作者知晓正态分布对连续变量描述统计的重要意义，根据数据的分布形态采用合适的方法。

Meta分析是循证医学中系统评价的定量合成的重要部分，而数据提取是Meta分析中的重要步骤。Meta分析的制作者需要首先根据“PICO原则”和数据的类型明确收集效应量指标的类型，并制定详细的提取流程。效应量的数据主要来源于已发表的相关文献、书籍、会议摘要等。目前在Meta实践中，对于连续变量信息，制作者常基于均数或标准化均数及其置信区间的数据信息进行数据合并，当原始研究未提供均数（标准差）而仅提供中位数（四分位间距）时，常采用的策略：（1）尝试联系作者索取详细数据后重新计算；（2）Meta分析时排除该文献；（3）根据中位数（四分位间距等）估计均数和标准差。其中，简单转化的策略为使用中位数代替均值，根据四分位间距计算标准差（标准差=四分位间距/1.35），但是明确获取的数据为偏态分布，均数和中位数相差较大，四分位间距计算标准差时要求样本量足够大同时满足正态分布时才能计算。Hozo SP及罗德惠等基于结局变量为正态分布的假设实现从样本量、中位数、极值或四分位数到均数与标准差的估算[2-3]，估算基于复杂的数学公式进行。此外，当样本数据已经确定为偏态分布，对均值和标准进行近似估计时受数据分布的偏度影响较大，不可避免存在一定的偏差，在Meta分析中进行多个中位数的转换后进行Meta合并，对结果的大小和方向可能产生无法预估的影响。

均数和中位数均是连续性资料的效应量的估计和表达，均可计算相应的置信区间[4]，为Meta分析进行的效应量合并提供可能。2018年，McGrath S等提出一系列方法实现单组的中位数效应合并和两组中位数差值的Meta合并（对应于两组均数差的Meta合并）[5-6]，并开发相应metamedian程序包，本文将为实现的原理以及实现方法进行解读，并介绍metamedian包的使用方法。

1 中位数作为效应量的统计基础

中位数是对连续变量集中趋势的统计量，可进行统计描述，相应的置信区间可用于统计推断。中位数与均数应用场景不同，但具有相同的统计地位。

1.1 非正态分布的常见描述统计

研究者根据研究目的设计并完成数据收集后，首先需要区分数据中各个变量的类型。连续变量是常见的变量类型，在进行描述统计前，需要先正态分布检验，其常见的方法有图示法、Shapiro-Wilk检验和Kolmogorov-Smirnov检验。若证据显示数据呈现为偏态分布，研究者常见的汇报显示为以下三种：S1=｛中位数、第一四分位数、第三四分位数、样本量｝，S2=｛中位数、最小值、最大值、样本量｝，S3=｛中位数、第一四分位数、第三四分位数、最小值、最大值、样本量｝。

1.2 中位数的点估计及其置信区间估计（推断统计）

当进行统计推断时，若连续变量的分布不满足参数分布的形态，常采用非参数检验方法，推断总体分布位置和已知分布位置、两个和多个总体分布位置是否有差异，对效应量进行点估计时，不再采用均值，而是中位数，例如配对样本的符号秩和检验、两独立样本的秩和检验。单个样本（配对样本）的中位数可通过简单排序方法获得，2个独立样本差值的中位数常采用Hodges-Lehmann（HL）法进行估计，需要注意的是进行2个独立样本的分布位置差异的估计时，不是先计算各组中位数后计算差异，而是估计各组间差异后再计算中位数。

均值的置信区间估计基于明确的概率分布和大样本，但中位数的分布未知，置信区间估计方法有所不同。目前常用基于样本数据模拟抽样的Bootstrap法，根据bootstrap样本统计量的分布可以估计总参参数（中位数），然后根据百分位数，如2.5%～97.5%对应95%置信区间，或者BCa（bias corrected,accelerated）法等构建置信区间[7]。R软件中的“rcompanion”和“boot”等程序包可进行单组和成组数据的中位数的置信区间的估计。需要注意的是，对于等级资料，特别是小样本且等级水平较少时，基于Boostrap法估计的置信区间可能不可信。

2 连续变量的效应量的Meta合并

2.1 基于均数的Meta合并

当原始研究采用均数和标准差进行汇报时，常采用的效应量为未标准化均数差或标准化均数差，效应量估计的精确性采用方差进行衡量，可进一步计算出对应的置信区间。根据研究设计类型常分为两种情况：

2.1.1 成组设计

原始研究常报告各组的均数、标准差、样本量。

其中n1和n2是各组样本量，S1和S2是各组的标准差。

2.1.2 配对设计

Sd为前后测量差值的标准差，r为前后测量的相关系数。

因其方差的计算较为复杂，在此不作呈现，具体可参考相关书籍。计算每个纳入的原始研究的均数差和方差后，再基于倒方差法采用固定效应模型或随机相应模型进行均值差的合并，得出合并效应的点估计和置信区间估计。

2.2 基于中位数的Meta合并

当原始研究采用S1、S2或S3进行汇报时，常采用的效应量为中位数。合并效应量的点估计和置信区间根据研究设计类型可归纳为以下2大类[5]：

2.2.1 成组设计

原始研究常报告各组的中位数、分位数、极值等，有3种估计方法。

第一，基于HL法和符号检验的方法，该方法未采用倒方差法进行效应合并。合并效应计算中位数差值的中位数（Median of the Difference of Medians，MDM）进行估计，主要采用HL法，具体实现如下：共纳入k个研究，每个研究有2个组，从原始研究中观测值主要采用中位数M1i，M2i（i=1，2，3，…k），计算M1i和M2i之间所有的差值，共k×k个，计算k2个差值的中位数，即为合并效应的点估计值，然后采用符号检验（sign test）构建相应的置信区间。因为符号检验基于离散数据的二项式分布（i=1，2，3，…k）计算累积概率，寻找第k个数满足条件，如5%累积概率，所以采用该方法估计的95%置信区间可能不完全等于95%，受纳入Meta研究个数的影响较大。

第二，基于线性分位数混合效应模型（Linear quantile mixed model，LQMM） 进行估计，相比于分位数回归允许纳入随机效应，相比于线性混合效应模型的非参数的方法，实现基于中位数的效应合并。在进行Meta合并时，将每个研究的中位数的差值作为因变量，设置随机截距（random effect），这意味着默认每个研究的效应之间存在差异，由固定效应部分和随机效应部分构成，随机部分满足均数为0，方差为σ2的正态分布，估计截距的固定效应部分和95%置信区间即为Meta合并效应及其95%置信区间。

第三，基于分位数估计（Quantile estimation,QE）的方差估计。原始研究中的2组（一般为试验组和对照组）均采用S1或S2或S3的方式进行统计描述，基于现有信息拟合正态分布、对数正态分布、Weibull分布和Gamma分布4种参数分布，其概率密度函数为f(*)，计算基于参数分布的分位数和已观察样本中分位数之间的差异的平方和，根据差异最小法确定最合适的参数和分布类型。中位数的差值满足近似正态分布（公式1），则可以计算每个研究中位数差值的方差，最后采用倒方差的方法进行效应量的合并和置信区间的估计。

mi1、ni1、f(mi1)为第i个研究的组1的中位数、样本量、估计的最合适概率密度函数。

2.2.2 配对设计

原始研究常报告单组中位数、分位数、极值等，有3种估计方法。第一、基于中位数的中位数（Median of Medians,MM）估计。合并效应的点估计采用纳入研究的中位数效应的中位数，置信区间的估计采用符号检验（sign test），因为符合检验基于P=0.5的二项式分布，当实验次数较大时，二项分布趋近正态分布（Normal approximation），当采用渐近正态分布时，将纳入的效应从小到大排序（m1≤m2≤m3…mk），然后计算以下百分位数（公式2）作为95%置信区间。当其估计其 他置信区间范围时，替换其中的Z0.025即可。采用该方法估计的95%置信区间为渐近区间。

k为纳入Meta分析的研究个数，Z0.025为标准正态分布的97.5%对应的Z分数。

第二、基于中位数的加权中位数（Weighted Median of Medians,WMM）估计，相比前一种方法考虑权重对合并效应估计的影响，常用的权重wi为各研究样本量所占的比例（公式3）。按照第一种方法计算所需的分位数，根据加权分位数法计算合并效应的点估计值和置信区间。其中，加权分位数法首先将纳入的效应从小到大排序（m1≤m2≤m3…..mk），根据排序后的中位数对应的权重用阶梯函数绘制经验累积概率分布，x轴为中位数效应，y轴为累积概率，最后按y轴对比前述的分位数确定合并效应的点估计和置信区间。

第三、基于分位数估计（Quantile estimation,QE）的方差估计，实现思路与成组设计中的QE法类似，不在进行阐述。

2.3 基于中位数的Meta合并的异质性评价

对纳入研究的异质性进行评价是Meta分析中的重要步骤，它决定我们收集的原始资料能否进行合并。在基于中位数的Meta分析中，QE法计算出每个研究的中位数效应量及对应的方差，根据DerSimonian-Laird法进行异质性评价，可计算常见异质性评价指标τ2和I2等。其他方法，如MM、WMM和MDM暂未提供合适的异质性评价。

3 metamedian程序包在R软件实现

3.1 R软件和metamedian程序包安装

为保证程序的正常运行，建议从https://www.r-project.org/下载R软件[8]，选择国内镜像下载安装文件（https://mirrors.tuna.tsinghua.edu.cn/CRAN/）。本文选择R 3.5.3版本，在64位的win10系统上进行安装，建议读者选择安装RStudio。它是针对R开始的IDE，能有效编写和存储R脚本，管理数据、R程序包和绘制的图片等。在菜单栏找到R x64 3.5.3并启动，在交互式窗口（R-GUI）中的提示符“＞”后输入命令install.packages（"metamedian"），安装时会自动安装“estmeansd”，“Hmisc”，“metafor”和“stats”程序包，等待数分钟后完成[9]。

3.2 数据分析

metamedian程序包中主要有pool.med和qe两个函数，经过对源代码分析和阅读帮助文档，pool.med函数可实现MM、WMM和MDM三种方法，用于合并来自配对设计和成组设计中获取的单组的中位数效应。qe函数采用QE法进行单组数据和2组独立数据的中位数效应量合并，其先调用qe.study.level函数计算各研究效应的方差，然后调用metafor程序包中的rma.uni函数实现效应量合并。成组设计中的MDM法和LQMM法未在该程序包实现，但MDM法可结合pool.med函数中的代码和HL法修改完成，LQMM法则可在lqmm程序包中实现。

3.2.1 单组数据的合并

数据来源于metamedian程序包的帮助文件，提供单组数据表示效应量，采用MM、WMM和qe三种方法实现效应量的合并。代码如下：

分析结果详见表1。需要注意的是：在pool.med函数中设置norm.approx中设置为FALSE时，表明置信区间的估计不采用渐近正态分布估计，而采用符号检验。此时，如果纳入的研究小于6，程序会提示“Not enough studies for exact coverage ＞= 95% CI”，该警告未在帮助文件中说明。此外，当使用加权分位数估计法时，文献中说明权重是与纳入研究的样本量成比例且之和为1，但在实际的R实践中因pool.med函数调用的是“Hmisc”程序包中的wtd.quantile函数，需要公式3中的权重乘以100后作为参数传入。

3.2.2 成组数据的合并

数据来源于metamedian程序包的帮助文件，常见S1、S2或S3的格式进行汇报时，本例以S1格式的数据为例，其他形式仅修改和添加对应参数即可，代码如下：

分析结果见表2。其中，qe函数中可通过single.family参数指定两组的参数分布是否一致，默认为FALSE。loc.shift参数指定两组的参数分布的差异是否仅限定为位置参数偏倚，而不考虑参数分数的其他特征，默认为FALSE。loc.shift参数指定两组的参数分布的差异是否仅限定为位置参数偏倚，而不考虑参数分数的其他特征，默认为FALSE。

表1 三种方法实现单组中位数数据的Meta合并Table 1.Three methods to achieve meta-analysis of single group median data

表2 QE法实现成组数据的中位数Meta合并Table 2.Performing meta-analysis of group median data using the QE method

4 讨论

随着生物统计分析在研究报告规范中的不断强化，强调针对数据的分布形态采用合适的方法进行统计描述和统计推断。循证医学的快速发展，Meta分析的需求越来越多，针对连续变量的效应量合并时，目前大部分研究采用均数（标准差）的效应量进行Meta合并。而当纳入的文献因数据为非正态分布时，汇报的中位数（四分位间距等）为效应量，对Meta分析造成一定困扰。均数和中位数均为连续变量集中趋势的体现，均可计算相应的置信区间，从而衡量效应量的大小、准确度和精度。目前对中位数作为效应量进行Meta研究不多，本文基于前人研究对中位数合并的原理进行总结，并对技术细节进行阐述，以帮助读者理解实现的前提、原理和常用方法。成组设计的数据中的Meta合并主要有MDM、LQMM和QE三种方法，单组设计的数据中主要有MM、WMM和QE三种方法[5]。QE法适用广泛，基于4种参数分布进行尝试，选择最佳的参数分布来估计每个研究效应的方差，可实现异质性的定性定量评价。

R软件作为一款免费的统计分析及绘图软件，其开放性使得研究者能快速实现不同算法并分享，越来越受到研究者的青睐。R软件中有metafor和metaplus等程序包可以实现均数（标准差）效应的Meta合并[10-11]。本文对R软件metamedian程序包实现中位数效应量的分析过程和结果进行了展示。读者在实际操作过程中根据纳入数据的不同形态可进行亚组分析，效应量为均数及95% CI和效应量为中位数及95% CI分别进行合并。目前该程序包中还未加入森林图和异质性的图形评价的绘制，读者可根据forest等程序包实现。基于R的开源特征，metamedian程序包的功能会在未来逐步完善。

利益冲突：所有作者声明不存在利益冲突。