用R做meta分析（附效应量计算神器）

Hello，这里是行上行下，我是喵君姐姐~众所周知，R具有免费、源代码开放，以及出色的统计计算和绘图表现能力等一系列优点，颇受科研人员的喜爱。当然，这些优势在meta分析上也表现得极为突出。本教程以标准平均差(cohen’s d)作为效应量的meta分析为例进行演示，适用于比较两种实验条件之间因变量的差异。例如“有、无干预条件下社交焦虑程度的差异”、“集中注意、分散注意条件下n-back任务的正确率差异”等等。简单来讲，实证研究中，可以进行t检验的结果，都可以参考以下教程，对数据进行meta分析。接下来，我们邀请到FarAway将逐步介绍用R做meta分析的基本步骤。本文以软件实操讲解为主，如果读者对meta分析的原理感兴趣，推荐以下两本书自行学习。

（扫码即可购买，满100减50，还可叠加优惠券哟~）一、 准备工欲善其事必先利其器，先来说说准备工作：1. 下载并安装R studio，在往期推文中我们已经详细介绍了R和Rstudio的安装教程，按步骤操作即可！2. 安装三个工具包：meta、metafor和xlsx3. 准备好上述需要进行meta分析的数据，存放在电脑中，以备调用。

为方便讲解，我们引用了《meta分析导论》里的一组数据，如下：从左到右，依次是纳入的研究效应量的编号(No)、第一作者的姓名(author)、发表年份(year)、实验组因变量均值(exp_mean)、实验组均值的标准差(exp_sd)、实验组样本量(exp_n)、控制组因变量均值(con_mean)、控制组均值的标准差(con_sd)、控制组样本量(con_n)、效应量(cohen’s d)、效应量95%置信区间下限(lower)、效应量95%置信区间上限(upper)、施测的地区(area)。从上到下，每一行数据都代表一个原始研究(single study)的一个效应量，如果一个原始研究中包括多个效应量（如：有好几个实验（或成对比较）及其结果），那么每一个效应量应该作为独立的一行列在表格中。我们注意到，很多实证研究中只报告均值、标准差/标准误、样本量等描述统计结果，并未提供效应量，因此，通常在meta分析前，我们需要自己计算效应量及其置信区间。文末有效应量计算工具的推荐。二、 配置环境和工具包相信用过R软件的小伙伴一定对工具包的安装不陌生（可参见：R语言安装及入门；以及RStudio介绍及入门）。如图输入代码-用鼠标选中代码内容-点击Run运行代码，R会自动运行下载安装步骤。

**如果xlsx安装不上，可能是电脑java环境配置的问题，可以参考这篇文章《解决在R Studio中无法加载xlsx包的教程及方法配置java环境》。**如果仍然无法安装xlsx package，也可以尝试其他调用数据的方法：例如：采用Rstudio导入数据的工具

三、导入数据以代码调用数据为例，如图输入数据存放的路径代码。

此条代码中，路径后的蓝色1代表调用xlsx里sheet1的数据，如果需要调用sheet2里的数据，则将1改成2即可。这里对需要分析的数据赋名为“metagen”，此后需要用metagen()函数执行meta分析。选中8-10行代码，点击Run，则完成数据导入。**查询路径的方法是找到数据文档，右键属性，查看位置。

数据导入成功后，可以看到Console栏里出现数据内容。

在Environment栏里也能看到metagen数据集，点击红框处，可以看到编列如表格的数据。

四、计算综合效应量终于到了小伙伴们最关心的综合效应量计算环节了！对cohen’s d进行综合之前，我们还需要完成一步数据转换的工作。Metagen()函数实现cohen’s d的meta分析是基于d值及其标准误，以下是效应量标准误计算的代码：

运行后，可以看到Console窗口里增加了一列se。

接着，就需要载入meta工具，运行如下代码：

分析之后可以看到console报告的结果：

红框显示的是meta分析里纳入的每一个研究（或效应量）所对应的作者名、年份、Cohen’s d值，置信区间，以及固定效应模型(fixed)和随机效应模型(random)计算出的权重。本次纳入meta分析的研究（或效应量）有6个，因此k=6。黄框显示的是综合效应量的结果。如图可见：1. 固定效应模型（fixed effect model）计算的综合效应量为0.4166，95%CI为[0.2906; 0.5426]，z检验值为6.48，p<0.0001。【结果显著】2. 随机效应模型（random effects model）计算的综合效应量为0.3602，95%CI为[0.1534; 0.5671]，z检验值为3.41，p=0.0006。【结果显著】**可以看到这里采用两个模型计算出的显著结果略有不同，研究中通常会采用更加保守的随机效应模型。感兴趣的小伙伴可以进一步了解两个模型的区别。蓝框显示的是异质性检验的结果。主要看Q检验的结果，以及I^2(I2)的值，一般来讲，I2 值0%-25%为低异质性, 50%-75%为中等异质性，>75% 为高异质性。异质性越高，代表纳入的研究（或效应量）间差异越大，在实证研究中，这种差异可能是由于实验范式、被试特征、施策时间、施策地点差异等造成，因此大部分meta分析会进一步假设可能的影响因素，再做亚组分析（详见第七步）。五、绘制森林图在meta分析的论文里森林图是最为常见的结果可视化图形，它可以直观表现每个效应量以及综合效应量的结果。利用R绘制森林图十分简单，运行如下代码即可。

在Plots窗口就可以看到自动生成的森林图，可以通过export导出图片，简单的处理后就可以插入论文啦！

六、发表偏倚检验Meta分析的发表偏倚通常会通过：肉眼看漏斗图的对称分布主观判断和基于egger’s test的统计结果来判断，以下就是两种方法的代码1. 漏斗图代码：

运行后在Plots可以查看并导出漏斗图。

2. egger’s test代码

运行后可以在console窗口可以查看分析结果，p>0.05，说明不存在发表偏倚，对应的检验图在plots窗口可以查看并导出。

假设上述分析发现，存在发表偏倚，可以再做trim-and-fill method剪补法来检验这一发表偏倚是否会影响到综合效应量的结果。3. 剪补法运行如下代码：

剪补法分析后，结果可以在Console窗口查看：

从分析结果看，原数据还是存在一定的偏差，漏斗右侧少了两个效应量（此时k=9，剪补法分析后的漏斗图右侧多了3个空心圆，代表补上的效应量）。通过算法自动补齐后，得到了一个新的综合效应量0.56，显著性跟剪补法分析前的没有变化，因此，在一定程度上可以说明综合效应量不会受发表偏倚的影响。七、亚组分析并绘制森林图在meta分析中，我们通常需要对单个研究进行分类讨论，类似于将单个研究按照变量分成不同的水平。在本范例中，我们假设纳入的部分单个效应量数据采集自A区域，另一部分采集自B区域，我们想要考察两类问卷测量出的结果是否存在差异，可以对问卷进行亚组分析，并绘制森林图。运行如下代码：

Console窗口显示两个亚组各自的综合效应量，可以看到：固定效应模型和随机效应模型的亚组分析都没有发现组间(between groups)显著的差异，表明数据采集的地区并不调节综合效应的结果。

八、元回归分析在meta分析里，我们对分类变量的调节作用的考察可以通过亚组分析来完成，对连续变量的调节作用，我们可以通过元分析来考察。例如：我们想知道文章发表年份(year)是否调节综合效应量的结果，我们可以采用如下代码进行元回归分析：

可以看到元回归并没有发现显著的结果，即year并不调节综合效应量，回归系数为0.1035，p值为0.6395，截距为0.2993，R方为0.00%。九、敏感性分析我们也看到很多的meta分析论文会做额外的敏感性分析，尤其在发现综合效应量的异质性较大的情况下，通过敏感性分析可以查看是否存在某些极端的纳入研究，产生对综合效应量的影响。

运行代码后，我们看到结果显示依次逐个剔除纳入的研究，综合效应量的显著性并未发生明显的改变，表明纳入的研究不存在极端的现象。

**meta工具包还可以自动生成一些符合期刊要求的图型，假设你想投医学期刊JAMA，可以在绘图前设置jama参数，代码如下：

接着再绘制图形，就会生成符合JAMA要求的图形啦！例如，在设置后再次运行绘制森林图，就可以获得如下图形：

想要恢复默认设置，输入以下代码即可：

效应量计算工具有很多，在这里推荐一款专业的在线计算效应量的神器：Practical Meta-Analysis Effect Size Calculator，基本上可以满足大多数效应量计算的需要。网址如下：http://www.campbellcollaboration.org/escalc/html/EffectSizeCalculator-SMD-main.php这款工具最近也开发了R的程序包，感兴趣的朋友可以关注，它可以实现一些比网页版更高级的功能。esc package网址如下：https://strengejacke.github.io/esc作者：FarAway排版：喵君姐姐*图片与书籍来源网络，其目的在于促进信息分享交流，不做商业用途。若有侵权，请在后台留言，管理人员删侵。