随机对照研究(二分类变量)的Meta分析——R软件实现(metafor包)

关键词：R; Meta分析; 随机对照试验; 二分类数据; metafor包

一、案例介绍

本文以模拟数据为例，为研究补钙对预防骨质疏松发生的效果，检索纳入9篇RCT研究文献，具体数据见表1，干预组为钙补充剂组，对照组为安慰剂组，具体指标含义见表1，针对这9篇文献的相对危险度进行Meta分析。案例数据可从“附件下载”处下载。

二、软件操作及结果解读

(一) 数据整理及读取

二分类数据主要是符合2×2四格表形式，在metafor包中，变量分别为试验组和对照组中发生事件人数和未发生例数。表1中t1、n1、c1、n2分别表示实验组结局事件发生数、试验组总人数(注意不是实验组未发生例数)、对照组结局事件发生数、对照组总人数。故需要生成试验组未发生例数t0,对照组未发生例数c0。

mydata <- read.csv("随机对照研究(二分类数据).csv")
mydata$t0 <- (mydata$n1-mydata$t1)  #生成试验组未发生例数
mydata$c0 <- (mydata$n2-mydata$c1)  #生成对照组未发生例数
View (mydata)

(二) 异质性检验

选择RR为效应量，任意选用固定或随机效应模型进行异质性检验，可以直接得到Q统计量和I²统计量，再根据异质性检验的结果选择相应的效应模型。

代码如下：

library (metafor) 
mydata1 <- rma (ai = t1, bi = t0, ci = c1, di = c0, data = mydata, measure = "RR", method = "REML")
summary (mydata1)

从结果图2可知，异质性检验Q统计量=34.2016，P <0.001，I² = 77.67%，均提示研究间异质性较大，故选用随机效应模型。

(三) 合并效应量，绘制森林图

根据具体要求，合并数据，进行Meta分析，采用倒方差法拟合随机效应模型，合并RR，采用REML法估计研究间异质性方差。以print()函数打印Meta分析结果，predict()函数将结果返回RR，以forest()函数绘制森林图。结果见图3—图5。

代码如下：

mydata.meta <- escalc(measure="RR", ai=t1, bi=t0, ci=c1, di=c0, data=mydata)
mydata.meta2<-rma.uni (yi, vi, data=mydata.meta, slab=paste (mydata.meta$Author, mydata.meta$Year, sep= ","), digits=3)
print (summary(mydata.meta2), digits=3)
predict (mydata.meta2, transf=exp, digits=2)

该命令参数介绍如下：

• measure：表示选择效应值，本例选择RR 为效应值；
• ai，bi，ci，di：分别代表试验组事件发生人数和未发生人数，对照组事件发生人数和未发生人数；
• data：为所要分析的数据对象；
• yi：表示选择效应值，如“RR”，“OR”等；
• vi：表示抽样方差的大小；
• slab：表示标签；
• digits：表示选择结果的小数位数；

forest (mydata.meta2, slab = paste(mydata.meta$Author, mydata.meta$Year, sep = ", "), xlim = c(-16, 6), at = log(c(0.05, 0.25, 1, 4,20)), atransf = exp,ilab = cbind(mydata.meta$t1,mydata.meta$n1,mydata.meta$c1,mydata.meta$n2), ilab.xpos = c(-9.5, -8, -6, -4.5), cex = 0.75)
text (-16, 11, "Author(Year)", pos = 4, font =2)  #在森林图上添加表头
text (6, 11, "Relative Risk [95% CI]", pos = 2, font =2)
text (c(-8.5, -7, -4.6, -3.1), 11, c("Events", "Total", "Events", "Total"), pos = 2, font =2)
text (c(-7.3, -4), 12, c("Experiental", "Control"), pos = 2, font =2)

森林图中水平横线表示该效应量的95%置信区间范围，点大小表示该研究的权重，即该项研究对Meta分析的贡献度；图中的菱形为合并效应值；图中的垂直实线用于判定结果差异有无统计学意义，若单个研究或合并效应量的95%置信区间与该直线相交，则代表两组的差异没有统计学意义。

(四) 结果解读

无论是数字化结果，还是森林图，均给出了每一研究效应量的点估计及95%可信区间。随机效应模型得到的合并效应量点估计及95%可信区间为0.69 (95%CI: 0.57~0.84)，相应z = -3.705，P<0.001。表明钙补充剂组发生骨质疏松的风险是对照组的0.69倍，即补充钙剂可降低31%的骨质疏松发生风险。

三、知识小贴士

(一) metafor包中进行Meta分析的部分函数用法

metafor包中进行Meta分析主要有两种使用途径，第一种途径为用ram.uni()函数联用escalc()函数，具体用法为先运行ram.uni()函数会得到效应值对数及方差的结果，分别以yi，vi表示。然后再采用ram.uni函数拟合经典的固定效应模型和随机效应模型。第二种途径为直接利用ram.uni()函数。

escalc()函数和ram.uni()函数设置选项及说明如下，更多设置可查看meta包中escalc()函数和ram.uni()函数自带帮助文件：

escalc()函数的主要用法为：escalc (measure, ai, bi, ci, di, nli, n2i, data=NULL, add=1/2,…)

• measure：表示选择效应值，如“RR”，“OR”等；
• ai，bi，ci，di：分别代表试验组事件发生人数和未发生人数，对照组事件发生人数和未发生人数；
• n1i，n2i：表示试验组和对照组的总人数；
• data：为所要分析的数据对象；
• add：表示四格表中有“0”值时函数的赋值，运行该函数后，会得到效应值对数及方差的结果。

ram.uni()函数的主要用法为：ram.uni(yi, vi, ai, bi, ci, di, n1i, n2i, mods, measure, data, intercept=TRUE, data, slab, subset, add=1/2, to=“only0”, drop00=FALSE, vtype=“LS”, method=“REML”, weighted=TRUE, test=“z”, level=95, digits, btt, tau2, verbose=FALSE, control, …)

• ai，bi，ci，di：分别代表试验组事件发生人数和未发生人数，对照组事件发生人数和未发生人数；
• n1i，n2i：表示试验组和对照组的总人数；
• yi：表示选择效应值，如“RR”，“OR”等；
• vi：表示抽样方差的大小；
• mods：用于指定将调节因子纳入模型，用于Meta回归模型；
• measure：表示选择效应值，如“RR”，“OR”等；
• data：表示所要分析的数据对象；
• intercept：表示该模型中是否包含截距，默认是TRUE；
• slab：表示标签；
• subset：可选的(逻辑或数字)向量，用于指定将包括在该函数返回的数据框中的研究子集；
• add：一个非负数，用于指定添加到零单元格、计数或频率的数量；
• to：一个字符串，用于指定何时添加添加项下的值(“全部”、“仅0”、“if0all”或“无”)；
• drop00：逻辑参数，指定在计算观察到的效应大小或结果时，是否应将两组中没有病例/事件(或只有病例)的研究剔除；
• vtype：字符串，用于指定要计算的抽样变异的类型；
• method：代表选择方差分量tau2的估计方法，默认为限制极大似然估计(REML)；
• digits：表示选择结果的小数位数；
• btt：在含有解释变量的模型中，模型中除了截距之外的其它参数都是通过多类题测试(omnibus test)进行估计的，btt是用来指定要测试的系数的索引；
• level：表示置信区间；
• tau2：用于指定随机效应模型和混合效应模型中的残留异质性的大小，如果不特别指定，tau2会从数据中测得；
• control：迭代估计算法的控制值的可选列表。如果没有指定，默认值将在函数内部定义。

(二) 森林图代码参数

森林图的部分设置选项及说明如下，更多设置可查看meta包中forest()函数自带帮助文件：

forest(x, annotate=TRUE, addfit=TRUE, addpred=FALSE,showweights=FALSE, header=FALSE, xlim, alim, olim, ylim, at, steps=5, level=x$level, refline=0, digits=2L, width,xlab, slab, mlab, ilab, ilab.xpos, ilab.pos,order, plim=c(0.5, 1.5), colout,col, cex, cex.lab, cex.axis, ...)

• x：表示目标对象；
• annotate：逻辑参数，用于指定是否应将注释添加到绘图中(默认为TRUE)；
• addfit：逻辑上指定是否应该将摘要估计(对于没有调节器的模型)或拟合值(对于有调节器的模型)添加到图中(默认为TRUE)；
• addpred：逻辑上指定预测区间的边界是否应该被添加到图中(默认为FALSE)；
• showweights：逻辑上指定注释是否应该包括模型拟合过程中给予观察结果的权重(默认为FALSE)；
• header：指定是否应将列标题添加到绘图中(默认为FALSE)。也可以是一个字符向量来指定左右标题(或者只指定左边的)；
• xlim：绘图区域的水平限制。如果没有指定，该函数会尝试将绘图的水平限制设置为一些合理的值；
• alim：X轴的限制。如果没有指定，该函数会尝试将X轴的极限值设置为一些合理的值；
• olim：可选参数，用于指定观察/结果限制。如果没有指定，则不使用限制；
• ylim：绘图的Y轴极限。如果没有指定，函数会尝试将Y轴限制设置为一些合理的值；
• at：X轴刻度线的位置和相应的标签。如果没有指定，该函数会尝试将刻度线的位置/标签设置为一些合理的值；
• steps：x轴的刻度线的数量(默认为5)。当通过at参数指定位置时被忽略；
• level：0到100之间的数值，用于指定置信区间的水平(默认是取自对象的值)；
• refline：数值，用于指定垂直"参考"线的位置(默认为0)。可以通过将此参数设置为NA来抑制该线；
• digits：整数，用于指定X轴刻度线标签和注释的小数位数；
• width：用于手动调整注释列的宽度；
• xlab：X轴的标题。如果没有指定，该函数会尝试设置一个合适的轴的标题；
• slab：可选的向量，包含k个研究报告的标签。如果没有指定，函数将尝试从x中提取研究标签，或者在函数中创建简单的标签。要抑制标签，请将此参数设置为NA；
• mlab：可选的字符串，为来自等效或随机效应模型的汇总估计值提供一个标签。如果没有指定，标签将在函数中创建；
• ilab：可选的向量、矩阵或数据框，提供有关研究的额外信息，这些信息应被添加到图表中；
• ilab.xpos：数字向量，用于指定通过ilab给出的变量的水平位置(如果指定了ilab，必须指定)；
• ilab.pos：整数(1、2、3或4)，用于指定通过ilab给出的向量的对齐方式(2表示右对齐，4表示左对齐)。如果没有指定，默认是将标签居中；
• order：可选的字符串，用于指定研究报告的排序方式。也可以是一个变量，将根据该变量对研究进行排序；
• plim：长度为2的数值向量，用于缩放点的大小(当指定psize时被忽略)；
• colout：可选的字符串，用于指定绘制观察结果的颜色(如果没有指定，默认使用 "黑色")。也可以是一个矢量；
• col：可选的字符串，用于指定用于总结多边形或拟合值的颜色。如果没有指定，函数会设置一个默认的颜色；
• cex：可选的字符和符号扩展因子。如果没有指定，该函数将尝试将其设置为一个合理的值；
• cex.lab：可选的X轴标题扩展因子。如果没有指定，函数会尝试将其设置为一个合理的值；
• cex.axis：X轴标签的可选扩展因子。如果没有指定，函数会尝试将其设置为一个合理的值。