如何用 R 语言进行元分析？(Part1)
文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。
元分析在许多研究学科领域的循证证据中发挥着关键作用，如医学、社会科学和经济学等等。之前我们讲过用 CMA 软件进行元分析，今天主要介绍如何用 R 语言进行元分析。R 语言的魅力在于其快速、强大、完善且开源的算法基础，获取网址：http://www.r-project.org/。
mada (version 0.5.10)
https://cran.r-project.org/web/packages/mada/index.html
meta (version 5.1-1)
https://cran.r-project.org/web/packages/meta/index.html
metafor (version 3.0-2)
https://cran.r-project.org/web/packages/metafor/index.html
metasens (version 1.0-0)
https://cran.r-project.org/web/packages/metasens/index.html
mvmeta (version 1.0.3)
https://cran.r-project.org/web/packages/mvmeta/index.html
netmeta (version 2.0-1)
https://cran.r-project.org/web/packages/netmeta/index.html
rmeta (version 3.0)
https://cran.r-project.org/web/packages/rmeta/index.html
ellipse (version 0.4.2)
https://cran.r-project.org/web/packages/ellipse/index.html
这些 R 包是目前最新的版本，由于 R 包的更新速度比较快，所以这些版本号也会不断更新，所以，在用 R 语言进行元分析时，要注意你输入的 R 命令是否正确)。
使用以下 R 命令来安装这些 R 包的最新版本：
pkgs <- c("mada", "meta", "metafor", "metasens", "mvmeta", "netmeta", "rmeta", "ellipse") install.packages(pkgs, repos="http://cran.rstudio.com/")
R 不是一个菜单或工具栏所驱动的程序，必须在命令行输入命令，或者执行一个预先写好的脚本。在 R中需要注意区分大小写，大多数函数名都是小写的。要执行一个函数命令，必须打上括号。
加载、保存和恢复数据。
R 中存储数据的方式与其他包或软件略有不同，在 R 中可以被认为是一个虚拟库“或虚拟目录”，而不是工作区(也称全局环境或环境)，R 可以包含许多数据集，以及其他 R 对象，如函数。存储文件名的后缀有所不同，在 Mac/Linux/Unix 下默认名称为 .rdata，在 Windows 下默认名称为 _RData。
在这里，以 .rdata 文件名为例，将其存储在当前的工作目录中。在 R 命令行中查找当前工作目录需要输入 getwd()。通过键入对象的名称来访问对象。典型的对象包括数字、变量、向量、矩阵、数据帧(近似于其他包中的数据集)和函数。查看附加的 .rdata 工作区中的对象，可用 ls 函数。
> ls()
[1] "g" "lev" "n" "opar" "pie.sales" "pin"
[7] "scale" "usr" "x" "xadd" "xdelta" "xscale"
[13] "xx" "y" "yadd" "ydelta" "yscale" "yy"
删除 R 中的对象 g，可以使用 R 函数 rm：
> rm(g)
> ls()
[1] "lev" "n" "opar" "pie.sales" "pin" "scale"
[7] "usr" "x" "xadd" "xdelta" "xscale" "xx"
[13] "y" "yadd" "ydelta" "yscale" "yy"
删除工作区中的所有 R 对象，可使用以下命令：
> rm(list=ls())
> ls()
character(0)
结果字符(0)表示在工作区中找不到R对象。注意，只要不将其保存在 .rdata 文件中，目前对工作区的更改就不是永久的。如果R由于内部错误而停止(基本上没发生过)，或者终止了 R 进程，那么对工作区的所有更改都将丢失。所以，你可以随时在 R 中使用命令 save.image() 保存工作空间。
要更改当前的工作目录，输入 setwd("my_directory")，将 my_directory 替换为已经存在的目录的路径，或者从 Mac OS 和 Windows 下的工具栏中更改当前的工作目录。
接下来，通过一系列命令来加载、查看和保存元分析中的数据。以 dataset1.csv 为例(示例数据来自 Spooner等人(2002) 研究中的数据)，并将其保存在当前的工作目录中。要检查该数据是否存在，使用以下命令：
> list.files(pattern="dataset1")
[1] "dataset1.csv"
这个命令能列出当前工作目录中包含“dataset1.csv”的所有文件，这里只需要列出文件dataset1.csv。用下方命令将这些数据加载到 R 中。str 命令显示 R 对象 data1 是一个数据帧，即一个数据集，包括 17 项研究。
最后，只需输入 save.image() 来保存数据。如果退出 R 不保存数据，输入 q(“no”)。注意，不是保存整个工作区，而是 R 对象的子集，例如数据集或函数。例如，为了将 R 对象 data1 保存到当前工作目录下的文件 data1.rda 中，可以使用 save(data1, file="data1.rda") 命令。对于包含 R 对象的文件，推荐使用 .rda 扩展名。可以使用命令 load(“data1.rda”)，将文件 data1 加载到R中。使用 read.csv 函数读取该元分析数据。使用 str 命令可以显示导入的数据集的结构。输入数据集的名称 (data1)，可以显示整个数据集信息。
> # 1\. Read in the data
> data1 <- read.csv("dataset1.csv", as.is=TRUE)
> # 2\. Print structure of R object data1
> str(data1)
’data.frame’: 17 obs. of 8 variables:
$ author: chr "Boner" "Boner" "Chudry" "Comis" ...
$ year : chr "1988" "1989" "1987" "1993" ...
$ Ne : int 13 20 12 12 17 8 13 12 12 12 ...
$ Me : num 13.5 15.7 21.3 14.5 14.4 ...
$ Se : num 13.8 13.1 13.1 12.2 11.1 ...
$ Nc : int 13 20 12 12 17 8 13 12 12 12 ...
$ Mc : num 20.8 22.7 39.7 31.3 27.4 ...
$ Sc : num 21.5 16.5 12.9 15.1 17.3 ...
> # 3\. To view an R object, just type its name:
> data1
author year Ne Me Se Nc Mc Sc
1 Boner 1988 13 13.54 13.85 13 20.77 21.46
2 Boner 1989 20 15.70 13.10 20 22.70 16.47
3 Chudry 1987 12 21.30 13.10 12 39.70 12.90
4 Comis 1993 12 14.50 12.20 12 31.30 15.10
5 DeBenedictis 1994a 17 14.40 11.10 17 27.40 17.30
6 DeBenedictis 1994b 8 14.80 18.60 8 31.40 20.60
7 DeBenedictis 1995 13 15.70 16.80 13 29.60 18.90
8 Debelic 1986 12 29.83 15.95 12 48.08 15.08
9 Henriksen 1988 12 17.50 13.10 12 47.20 16.47
10 Konig 1987 12 12.00 14.60 12 26.20 12.30
11 Morton 1992 16 15.83 13.43 16 38.36 18.01
12 Novembre 1994f 24 15.42 8.35 24 28.46 13.84
13 Novembre 1994s 19 11.00 12.40 19 26.10 14.90
14 Oseid 1995 20 14.10 9.50 20 28.90 18.00
15 Roberts 1985 9 18.90 17.70 9 38.90 18.90
16 Shaw 1985 8 10.27 7.02 8 34.43 10.96
17 Todaro 1993 13 10.10 8.90 13 23.50 4.00
从R数据集中选择变量。
如果不首先引用对象 data1 的名称，R 就不知道变量的名称 author、year 等。使用只存在于变量名中的 R 命令将会出现错误。
> author
Error: object ’author’ not found
显示 author 变量的数据有不同的方法。首先，可以使用美元符号“$”来选择数据集的单个变量。下面的 R 代码表示从数据集 data1 中选择变量 author：
> data1$author
[1] "Boner" "Boner" "Chudry" "Comis"
[5] "DeBenedictis" "DeBenedictis" "DeBenedictis" "Debelic"
[9] "Henriksen" "Konig" "Morton" "Novembre"
[13] "Novembre" "Oseid" "Roberts" "Shaw"
[17] "Todaro"
其次，可以使用方括号来选择数据集的行(观测值)和列(变量)，从数据集 data1 中提取变量 author，即一个列，可以通过以下方式提取：
> data1[, "author"]
[1] "Boner" "Boner" "Chudry" "Comis"
*** Output truncated ***
此外，数据集 data1 的前四个作者 / 第一行可以用：
> data1[1:4, "author"]
[1] "Boner" "Boner" "Chudry" "Comis"
或者，
> data1$author[1:4]
[1] "Boner" "Boner" "Chudry" "Comis"
第三，可以使用 with 函数来选择作者列表。
> # List the first four authors in data frame data1
> with(data1, author[1:4])
[1] "Boner" "Boner" "Chudry" "Comis"
> with(data1[1:4,], author)
[1] "Boner" "Boner" "Chudry" "Comis"
with 函数的第一个参数是一个数据集，这里分别是“data1”和“data1[1:4,]”。with 函数的第二个参数是在指定数据集内求值的表达式。
第四，附加一个数据集，即使用 R 命令 attach，使一个数据集的所有变量在搜索路径中是可用的：
> # 1\. Make variables of data frame data1 directly available:
> attach(data1)
> # 2\. List the first four authors in data frame data1
> author[1:4]
[1] "Boner" "Boner" "Chudry" "Comis"
> # 3\. Detach data frame data1
> detach(data1)
这种方法的一个缺点是工作空间中的 R 对象可能会在附加的数据集 data1 中隐藏一个变量，接下来会在示例 R 代码中展示这一点。
> # 1\. Create a new R object author (numeric variable)
> author <- 123
> # 2\. Attaching data frame data1 results in a warning
> attach(data1)
The following object is masked _by_ .GlobalEnv:
author
> # 3\. The following command prints the numeric variable author
> author[1:4]
[1] 123 NA NA NA
> # 4\. Search for R objects called "author"
> find("author")
[1] ".GlobalEnv" "data1"
> # 5\. Detach data frame data1
> detach(data1)
> # 6\. Remove R object author from global environment
> rm(author)
find 命令显示 R 对象 author 同时存在于工作区“.GlobalEnv”和 R 对象“data1”中。当首先搜索工作空间，这个 R 对象可以在 print 命令中使用。
在这里，使用前三种方法中的一种来从数据集中选择变量。如果对单个变量的信息感兴趣，则首选美元赋值。方括号可用于从数据集中选择多个变量。例如，可以使用下面的 R 代码显示出版年份、样本量以及作者姓名：
> data1[1:4, c("author", "year", "Ne", "Nc")]
author year Ne Nc
1 Boner 1988 13 13
2 Boner 1989 20 20
3 Chudry 1987 12 12
4 Comis 1993 12 12
另一方面，函数 with 在数据集的多个变量的计算中更简洁。
> # 1\. Calculate and display the total sample sizes:
> data1$Ne + data1$Nc
[1] 26 40 24 24 34 16 26 24 24 24 32 48 38 40 18 16 26
> with(data1, Ne + Nc)
[1] 26 40 24 24 34 16 26 24 24 24 32 48 38 40 18 16 26
> # 2\. Calculate and display the total sample size
> # for the Chudry study
> data1$Ne[data1$author=="Chudry"] +
+ data1$Nc[data1$author=="Chudry"]
[1] 24
> with(data1[data1$author=="Chudry",], Ne + Nc)
[1] 24
运行脚本。
这里，将演示如何创建和运行一个简单的源文件。从工具栏上单击File--New Script (用于Windows) 或者File--New Document (用于Mac OS)。将出现一个新的脚本文件，在脚本文件中输入以下命令：
getwd()
dir(pattern="example1")
data1 <- read.csv("dataset01.csv")
summary(data1)
print(summary(data1))
保存脚本文件。尽管它是一个文本文件，但它的后缀是“.R”。要运行第一行，将光标移动到这一行，按住 Ctrl 键并按 R(用于Windows)或 cmd 键和<RETURN>(用于Mac OS)。要运行部分代码或整个脚本，首先必须高亮显示相应的区域。
或者，如果文件作为example1.R保存在当前工作目录中，可以输入source("example1.R")来执行脚本文件中的所有命令。
安装和使用附加函数库。
在进行元分析之前的最后一步是安装元分析需要用到的 R 包。计算机需要连接到互联网，并且运行 R。在R命令提示符下输入 install.packages ("meta")或 install.packages()。运行第二种会打开一个对话框，按字母顺序显示可用的 R 包，向下滑动找到 meta 并单OK。要找到关于 R 包类型用 help(package=meta)。
需要注意的是，meta 包并不是唯一一个具有元分析功能的 R 包。其他具有一般元分析方法的 R 包是 metafor 和 rmeta。特别是，metafor 是大神们都建议安装的一个 R 包，可用于 meta 的某些方法，如研究间的方差和元回归方法中的极大似然估计方法，只有安装了 metafor 才可用。
关于 Spooner等人(2002) 元分析研究中所用到的森林图代码(前文提到过，本示例数据来自于Spooner)，生成如下所示的森林图：
> # 1\. Load add-on package meta
> library(meta)
Loading ’meta’ package (version 4.0-2).
> # 2\. Do meta-analysis
> m <- metacont(Ne, Me, Se, Nc, Mc, Sc,
+ studlab=paste(author, year),
+ data=data1)
> # 3\. Produce forest plot
> forest(m, xlab="Maximum % fall in FEV1")
详情请参见该研究的原文
使用 xlab 选项来标记 x 轴。在森林图的左侧可以看到每一项标有第一作者和日期的研究，以及用于计算的数据。在森林图的右侧，描述了均值差(MD) 及其 95% 的置信区间。默认情况下，固定效应和随机效应元分析的结果都会给出，标为 W(fixed) 和 W(random) 的列表示在各自的元分析中，这些研究所占的权重(百分比)。如果要从森林图中剔除固定效应或随机效应模型，可以用参数 comb.fixed=FALSE 或者 comb.random=FALSE。
标准分析方法
固定效应和随机效应元分析。
对于研究1 (Boner 1988)，实验组和控制组的均值分别是 13.54、20.77，均值差MD为 -7.23。对于研究2 (Boner 1989)，实验组和控制组的均值是 15.70、22.70，均值差MD为 -7.00。因此，(Boner 1988)和(Boner 1989)的 FEV1 最大平均下降约 7%。对于研究1 (Boner 1988)，95% 的置信区间为：
这里可以使用基础 R 来计算(Boner 1988)的均值差和 95% 的置信区间(假设文件 dataset1.csv 在当前工作目录中)：
> # 1\. Read in the data
> data1 <- read.csv("dataset01.csv", as.is=TRUE)
> # 2\. Calculate mean difference and its standard error for
> # study 1 (Boner 1988) of dataset data1:
> MD <- with(data1[1,], Me - Mc)
> seMD <- with(data1[1,], sqrt(Seˆ2/Ne + Scˆ2/Nc))
> # 3\. Print mean difference and limits of 95% confidence
> # interval using round function to show only two digits:
> round(c(MD, MD + c(-1,1) * qnorm(1-(0.05/2)) * seMD), 2)
[1] -7.23 -21.11 6.65
这里可以看到，95% 置信区间的均值差、下限和上限与手动公式计算出来的是一致的。也可以使用 R 中 meta 包的 metacont 函数来计算均值差和置信区间：
> with(data1[1, ],
+ print(metacont(Ne, Me, Se, Nc, Mc, Sc),
+ digits=2))
MD 95%-CI z p-value
-7.23 [-21.11; 6.65] -1.02 0.3074
Details:
- Inverse variance method
使用 metacont 函数的参数 sm="MD" 得到了相同的结果。与使用 with 函数相比，更方便的方法是使用带有参数 data 和 subset 的 metacont 函数。
> print(metacont(Ne, Me, Se, Nc, Mc, Sc,
+ data=data1, subset=1), digits=2)
MD 95%-CI z p-value
-7.23 [-21.11; 6.65] -1.02 0.3074
*** Output truncated ***
除了均值差及其 95% 置信区间外，metacont 函数还报告了用于检验整体效果的 z 分数和 p 值。这些值可用如下基础 R 函数 pnorm 和 abs 进行计算：
> zscore <- MD/seMD
> round(c(zscore, 2*pnorm(abs(zscore), lower.tail=FALSE
[1] -1.0206 0.3074
固定效应模型。
固定效应模型假设元分析中成分研究的估计效应来自单一同质群体。固定效应估计及其方差可以用如下 R 代码计算：
> # 1\. Calculate mean difference, variance and weights
> MD <- with(data1, Me - Mc)
> varMD <- with(data1, Seˆ2/Ne + Scˆ2/Nc)
> weight <- 1/varMD
> # 2\. Calculate the inverse variance estimator
> round(weighted.mean(MD, weight), 4)
[1] -15.514
> # 3\. Calculate the variance
> round(1/sum(weight), 4)
[1] 1.4126
使用 metacont 函数可以更容易地进行元分析，并且它产生的结果是相同的：
> mc1 <- metacont(Ne, Me, Se, Nc, Mc, Sc,
+ data=data1,
+ studlab=paste(author, year))
> round(c(mc1$TE.fixed, mc1$seTE.fixedˆ2), 4)
[1] -15.5140 1.4126
前面提到 MD 列，是实验组和控制组反应的均值差，以及固定效应 (%W(fixed)) 和随机效应 (%W(random)) 模型下，每个研究的权重。
研究1(Boner 1988)在固定效应元分析中所占的权重可用以下公式进行计算：
也可以用R来计算这些值：
> mc1$w.fixed[1]
[1] 0.01992783
> sum(mc1$w.fixed)
[1] 0.7079028
> round(100*mc1$w.fixed[1] / sum(mc1$w.fixed), 2)
[1] 2.82
由 R 命令生成的森林图和代码如下所示：
> forest(mc1, comb.random=FALSE, xlab=
+ "Difference in mean response (intervention - control)
+ units: maximum % fall in FEV1",
+ xlim=c(-50,10), xlab.pos=-20, smlab.pos=-20)
用 xlab 选项标记 x 轴，选项 xlim=c(-50,10) 用于指定 x 轴的范围为 -50 到 10 之间。选项 xlab.pos 和 smlab.pos 用于指定 x 轴上的标签中心和图形顶部的汇总测量区间，否则运行以上代码将以 0 为中心。元分析也可以使用 metagen 函数来进行基于一般逆方差法的元分析：
> # 1\. Apply generic inverse variance method
> mc1.gen <- metagen(mc1$TE, mc1$seTE, sm="MD")
> # 2\. Same result
> mc1.gen <- metagen(TE, seTE, data=mc1, sm="MD")
> # 3\. Print results for fixed effect and random effects method
> c(mc1$TE.fixed, mc1$TE.random)
[1] -15.51403 -15.64357
> c(mc1.gen$TE.fixed, mc1.gen$TE.random)
[1] -15.51403 -15.64357
随机效应模型。
随机效应模型试图解释这样一个事实，即研究效应估计
往往比固定效应模型中假设的变量更多。
①研究间方差估计
在 metagen 和其它 R 包 meta 中，可使用以下方法估计研究间方差τ2(参数为 method.tau)：
DerSimonian-Laird估计(method.tau="DL") (default)；
Paule-Mandel估计(method.tau="PM")；
Restricted maximum-likelihood估计(method.tau="REML")；
Maximum-likelihood估计(method.tau="ML")
Hunter-Schmidt估计(method.tau="HS")
Sidik-Jonkman估计(method.tau="SJ")
Hedges估计(method.tau="HE")
Empirical Bayes估计(method.tau="EB")
DerSimonian-Laird 估计法是目前比较流行的方法，尤其是在医学研究中。这个方法也是 R 包 meta 中的默认方法。下图显示了数据集中τ2的各种估计结果的森林图。DerSimonian-Laird、Restricted maximum-likelihood 和 Empirical Bayes 估计的结果相似。Sidik-Jonkman 估计值非常之大，而其他的估计值(如 Paule-Mandel，Maximum-likelihood，Hunter-Schmidt 和 Hedges)却又相当小。来自 Sidik-Jonkman 方法的非常大的估计值告诫我们不要完全依赖这种方法。
②Hartung-Knapp校正
可以使用 metacont 函数来进行 Hartung-Knap 校正：
> mc2.hk <- metacont(Ne, Me, Se, Nc, Mc, Sc, sm="SMD",
+ data=data2, comb.fixed=FALSE,
+ hakn=TRUE)
或者使用 metagen 函数：
> mc2.hk <- metagen(TE, seTE, data=mc2, comb.fixed=FALSE,
+ hakn=TRUE)
③预测区间
在 R 包 meta 中，预测期间的导入方式有几种。使用 metacont 函数创建元分析对象时的参数 prediction=TRUE，或者，在 summary，forest 或 print 命令中指定 prediction 参数。在下面的 R 代码中，使用的是 summary.meta 命令中的 prediction 参数。
> print(summary(mc1, prediction=TRUE), digits=2)
Number of studies combined: k=17
MD 95%-CI z p-value
Fixed effect model -15.51 [-17.84; -13.18] -13.05 < 0.0001
Random effects model -15.64 [-18.14; -13.15] -12.30 < 0.0001
Prediction interval [-19.94; -11.35]
*** Output truncated ***
使用以下命令可以很容易地生成显示预测区间的森林图：
> forest(mc1, prediction=TRUE, col.predict="black")
（未完待续~）
