【r<-包】metafor——元分析技术实现（一）
最近在科研中使用到元分析，零零散散开始接触一些理论与实现。元分析的主要作用就是汇总某一个方面所有符合条件研究的统计结果，评估这个方向研究一个平均效力是什么（假设有很多研究探索喝可乐对男子精子能力的影响，不同的研究可能有不同的结果，元分析就是对这些研究进行一个汇总，并给出判断）。涉及一些研究异质性的评估与回归分析、固定效应模型与随机效应模型的使用。
我最近使用的元分析处理主要是依靠Stata的一些模块，但其实我对Stata知之甚少。Stata语法有点像命令而不是一门系统的编程语法，难懂又难写，所以这里学习R中Meta分析的实现。R里面进行元分析的包非常多，包与介绍可以参考任务列表。这里我选择metafor包，它涵盖了常用的元分析建模与分析、可视化。
里面的概念解释不一定到位，读者可以参考文档理解。如果有错误，请见谅并评论提出，谢谢。
实例
metafor包提供了数据集dat.bcg，它是13个不同研究BCG疫苗因tuberculosis状态效力差异的研究。
p_load(metafor)
data("dat.bcg", package = "metafor")
查看数据：
print(dat.bcg, row.names=FALSE)
## trial author year tpos tneg cpos cneg ablat alloc
## 1 Aronson 1948 4 119 11 128 44 random
## 2 Ferguson & Simes 1949 6 300 29 274 55 random
## 3 Rosenthal et al 1960 3 228 11 209 42 random
## 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
## 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
## 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate
## 7 Vandiviere et al 1973 8 2537 10 619 19 random
## 8 TPT Madras 1980 505 87886 499 87892 13 random
## 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random
## 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic
## 11 Comstock et al 1974 186 50448 141 27197 18 systematic
## 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic
## 13 Comstock et al 1976 27 16886 29 17825 33 systematic
tpos与tneg分别记录了治疗组中tuberculosis阳性和阴性的目标数；cpos与cneg记录控制组中tuberculosis阳性和阴性的目标数。除此之外，ablat记录研究的绝对维度，alloc介绍研究的病人分组方式。
这个结果可以表示为以下2x2的表格：
| TB+ | TB- | |
|---|---|---|
| Treated | tpos | tneg |
| Control | cpos | cneg |
下面我们以相关风险（对数化）作为结果的度量，根据数据集中的相关数据进行计算：
dat = escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg,
append = TRUE)
print(dat[, -c(4:7)], row.names = FALSE)
## trial author year ablat alloc yi vi
## 1 Aronson 1948 44 random -0.8893 0.3256
## 2 Ferguson & Simes 1949 55 random -1.5854 0.1946
## 3 Rosenthal et al 1960 42 random -1.3481 0.4154
## 4 Hart & Sutherland 1977 52 random -1.4416 0.0200
## 5 Frimodt-Moller et al 1973 13 alternate -0.2175 0.0512
## 6 Stein & Aronson 1953 44 alternate -0.7861 0.0069
## 7 Vandiviere et al 1973 19 random -1.6209 0.2230
## 8 TPT Madras 1980 13 random 0.0120 0.0040
## 9 Coetzee & Berjak 1968 27 random -0.4694 0.0564
## 10 Rosenthal et al 1961 42 systematic -1.3713 0.0730
## 11 Comstock et al 1974 18 systematic -0.3394 0.0124
## 12 Comstock & Webster 1969 33 systematic 0.4459 0.5325
## 13 Comstock et al 1976 33 systematic -0.0173 0.0714
在这个结果中，yi是对数化的相对风险，<0表示治疗组中风险要更小，除了两个研究，其他研究都有一致的趋势。
measure可以指定很多种不同的度量方式，常见相对风险（对数化）为RR，优势比（对数化）为OR，其他请参考函数文档。对数化的目的是使得度量尽量服从正态分布。
如果要使用escalc()函数的公式形式，需要提供长格式数据：
k = length(dat.bcg$trial)
dat.fm = data.frame(study=factor(rep(1:k, each=4)))
dat.fm$grp = factor(rep(c("T","T", "C", "C"), k), levels = c("T", "C"))
dat.fm$out = factor(rep(c("+","-","+","-"), k), levels = c("+", "-"))
dat.fm$freq = with(dat.bcg, c(rbind(tpos, tneg, cpos, cneg)))
dat.fm
## study grp out freq
## 1 1 T + 4
## 2 1 T - 119
## 3 1 C + 11
## 4 1 C - 128
## 5 2 T + 6
## 6 2 T - 300
## 7 2 C + 29
## 8 2 C - 274
## 9 3 T + 3
## 10 3 T - 228
## 11 3 C + 11
## 12 3 C - 209
## 13 4 T + 62
## 14 4 T - 13536
## 15 4 C + 248
## 16 4 C - 12619
## 17 5 T + 33
## 18 5 T - 5036
## 19 5 C + 47
## 20 5 C - 5761
## 21 6 T + 180
## 22 6 T - 1361
## 23 6 C + 372
## 24 6 C - 1079
## 25 7 T + 8
## 26 7 T - 2537
## 27 7 C + 10
## 28 7 C - 619
## 29 8 T + 505
## 30 8 T - 87886
## 31 8 C + 499
## 32 8 C - 87892
## 33 9 T + 29
## 34 9 T - 7470
## 35 9 C + 45
## 36 9 C - 7232
## 37 10 T + 17
## 38 10 T - 1699
## 39 10 C + 65
## 40 10 C - 1600
## 41 11 T + 186
## 42 11 T - 50448
## 43 11 C + 141
## 44 11 C - 27197
## 45 12 T + 5
## 46 12 T - 2493
## 47 12 C + 3
## 48 12 C - 2338
## 49 13 T + 27
## 50 13 T - 16886
## [到达getOption("max.print") -- 略过2行]]
使用：
escalc(out ~ grp | study, weights = freq, data = dat.fm, measure = "RR")
## yi vi
## 1 -0.8893 0.3256
## 2 -1.5854 0.1946
## 3 -1.3481 0.4154
## 4 -1.4416 0.0200
## 5 -0.2175 0.0512
## 6 -0.7861 0.0069
## 7 -1.6209 0.2230
## 8 0.0120 0.0040
## 9 -0.4694 0.0564
## 10 -1.3713 0.0730
## 11 -0.3394 0.0124
## 12 0.4459 0.5325
## 13 -0.0173 0.0714
一般元分析的输入数据都是宽格式，所以使用escalc的默认接口最好。但长格式使用的公式接口可以更好地处理更复杂地数据结构。
拟合模型
通过rma.uni()可以拟合众多元分析模型（函数别名rma()）。
参数：
args(rma)
## function (yi, vi, sei, weights, ai, bi, ci, di, n1i, n2i, x1i,
## x2i, t1i, t2i, m1i, m2i, sd1i, sd2i, xi, mi, ri, ti, sdi,
## r2i, ni, mods, measure = "GEN", intercept = TRUE, data, slab,
## subset, add = 1/2, to = "only0", drop00 = FALSE, vtype = "LS",
## method = "REML", weighted = TRUE, test = "z", level = 95,
## digits = 4, btt, tau2, verbose = FALSE, control, ...)
## NULL
指定数据
该函数可以与元分析常见的效应大小或者结果度量（例如，log优势比，标准化的均值差，相关系数）。我们只需要简单指定yi观察到的效应值以及指定vi对应的样本方差（或者通过sei指定标准误，样本方差的平方根）。如果这样输入数据，必须指定measure="GEN"（函数默认）。
另外，rma函数可以输入escalc函数一样的参数，然后会自动完成结果度量的计算，measure参数可以用来指定想要的输出度量，vtype参数与escalc一致。
也就是说，可以用rma完成原始输入数据到结果度量的计算，可想而知，rma函数中应该包含了escalc函数的实现。
指定模型
假设你已经通过参数yi与vi指定了效应值与对应的样本方差，随机效应模型就直接可以用rma(yi, vi, data=dat)进行拟合。默认使用严格最大似然评估方法（restricted maximum-likelihood estimation）用来评估真实效应的方差（REML评估器近似无偏并非常有效，参考Viechtbauer 2005）。多种（残差）异质性评估器可以通过method参数指定，有以下选项：
- “HS”: Hunter-Schmidt estimator
- “HE”: Hedges estimator
- “DL”: DerSimonian-Laird estimator
- “SJ”: Sidik-Jonkman estimator
- “ML”: Maximum-likelihood estimator
- “REML”: Restricted maximum-likelihood estimator
- “EB”: Empirical Bayes estimator
一个或多个协变量（moderators）可以通过参数mods指定。单个协变量可以以长度为k的行/列向量给出值。多个则需要创建一个k行p列的设计矩阵（比如，使用mods = cbind(mod1, mod2, mod3)），模型默认包含截距，如果不要，需要用intercept = FALSE指定。
许多R用户熟悉用公式形式来指定想要构建的模型，比如在使用lm()与glm()时（参见??formula）。mods支持指定单边的公式，比如mods = ~ mod1 + mod2 + mod3。交互项、多项式、因子等都可以用这种形式添加。注意一旦使用公式指定mods参数，intercept参数会被忽略，我们需要通过mods = mod1 + mod2 + mod3 - 1删除交互项。
固定效应模型可以通过method = "FE"进行构建（rma(yi, vi, data = dat, method = "FE），注意加权最小二乘法与不加权~的使用，默认weighted=TRUE，协变量也时通过mods参数指定。
参数的多类题测试（Omnibus test)
Omnibus tests are a kind of statistical test. They test whether the explained variance in a set of data is significantly greater than the unexplained variance, overall
对于包含协变量的模型，会对所有模型系数进行Omnibus测试（排除第一个系数β0β0如果它包含在模型中，也就是截距）。检验原假设为H9:β1=...βp=0H9:β1=...βp=0。如果模型中没有截距，那么检验时会包含所有的系数（包含截距）。备择，我们可以通过btt参数手动指定要检验的协变量。例如btt = c(3, 4)会在检验中包含第3、4个系数（如果模型包含截距，那么截距对应模型的第1个系数）。
分类协变量
分类协变量可以以类似的方式参入模型，但需要手动转换为线性数值形式，或者用factor转换为因子。
Knapp和Hartung校正
默认，模型单个系数（置信区间）的检验统计量都是基于正态分布，然而全类题测试是基于m个自由度的卡方χ2χ2。m是检验系数的数目。Knapp和Hartung方法（2013）knha = TRUE是对所评估的系数标准误差的校正，它可以帮助解释真实效应方差的不确定性以及生成不同的参考分布。然后单个系数和置信区间通过n-p个自由度的t分布进行评估，而全类题检验是用m和k-p个自由度的F分布。该校正只能用于随机或者混合效应模型。
例子
随机效应模型
我们将用BCG数据拟合随机效应模型开始学习。
这与下面的命令等效，不过下面的代码是用的提前算好的效应值大小与方差。
res = rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat, measure = "RR")
res
res = rma(yi, vi, data=dat)
res
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
## tau (square root of estimated tau^2 value): 0.5597
## I^2 (total heterogeneity / total variability): 92.22%
## H^2 (total variability / sampling variability): 12.86
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
结果显示了估计的评估log相对风险是μ=−0.7145μ=−0.7145（96% CI: -1.0669 到 -0.3622）。想要更好地解释该结果，我们进行指数转换：
exp(c(-0.7145, -1.0669, -0.3622))
## [1] 0.489 0.344 0.696
这显示的是相对风险（RR，另有很多人在比例风险模型中称HR）及95%置信区间，说明打了疫苗个人有tuberculosis感染的风险是没有打疫苗的一半左右。可以拒绝零假设H0:μ=0H0:μ=0 （z=−3.97,p<0.0001z=−3.97,p<0.0001）。
在log的相对风险中，异质性评估为τ2=0.3132τ2=0.3132。在Higgins和Thompson(2002)一文中建议了大量的异质性度量方式。I2I2统计值（百分比）评估了在效应大小的（由异质性和取样变异组成）度量中总的变异有多少可以归根于真实效应的异质性（如果τ2=0τ2=0则I2=0I2=0）。而H2H2统计量则是在所观察到结果中总变异的量与取样变异量的比值（如果τ2=0τ2=0则显示H2=1H2=1）。然而我们需要意识到，τ2τ2、I2I2和H2H2在研究数目比较少是估计不准确。
使用confint()我们可以获得这些统计量对应的置信区间（95%）。
confint(res)
##
## estimate ci.lb ci.ub
## tau^2 0.3132 0.1197 1.1115
## tau 0.5597 0.3460 1.0543
## I^2(%) 92.2214 81.9177 97.6781
## H^2 12.8558 5.5303 43.0680
可以看到置信区间都比较宽，说明我们不能太相信点估计的结果。而且，下置信区间非常大并且异质性检验结果表明真实效应有相当多的异质性（变异，可以理解为方差很大）。
为了使结果更易于观察，我们可以用forest函数创建forest图来可视化结果。虽然forest(res)足以给出必要的信息，添加额外的代码可以补充更多信息。默认，观察到的效应值会按比例绘制到图上，基于随机效应模型汇总的估计也自动添加到图上。以log尺度显示可以获得更好的解释。
运行下面代码：
forest(res, slab = paste(dat$author, dat$year, sep = ", "),
xlim = c(-16, 6), at = log(c(0.05, 0.25, 1, 4)), atransf = exp,
ilab = cbind(dat$tpos, dat$tneg, dat$cpos, dat$cneg),
ilab.xpos = c(-9.5, -8, -6, -4.5), cex = 0.75)
op = par(cex = 0.75, font = 2)
text(c(-9.5, -8, -6, -4.5), 15, c("TB+", "TB-", "TB+", "TB-"))
text(c(-8.75, -5.25), 16, c("Vaccinated", "Control"))
text(-16, 15, "Author(s) and Year", pos = 4)
text(6, 15, "Relative Risk [95% CI]", pos = 2)
img
par(op)
混合效应模型
部分异质性可能是由协变量（moderator）引起。例如，BCG疫苗的有效性取决于研究的地点，不同的纬度病菌环境不一样。另外，使用疫苗的时间也有可能改变疫苗的效力。我们这里使用一个混合效应模型进行拟合，探究绝对纬度与发表年份作为协变量对结果的影响。
res = rma(yi, vi, mods = cbind(ablat, year), data = dat)
res
与下面结果一致
res = rma(yi, vi, mods = ~ ablat + year, data = dat)
res
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1108 (SE = 0.0845)
## tau (square root of estimated tau^2 value): 0.3328
## I^2 (residual heterogeneity / unaccounted variability): 71.98%
## H^2 (unaccounted variability / sampling variability): 3.57
## R^2 (amount of heterogeneity accounted for): 64.63%
##
## Test for Residual Heterogeneity:
## QE(df = 10) = 28.3251, p-val = 0.0016
##
## Test of Moderators (coefficient(s) 2:3):
## QM(df = 2) = 12.2043, p-val = 0.0022
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -3.5455 29.0959 -0.1219 0.9030 -60.5724 53.4814
## ablat -0.0280 0.0102 -2.7371 0.0062 -0.0481 -0.0080 **
## year 0.0019 0.0147 0.1299 0.8966 -0.0269 0.0307
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
残余异质性（除了被协变量解释的）τ2=0.1108τ2=0.1108，表明(0.3132 - 0.1108) / 0.3132 = 65%的总异质性可以归根于两个协变量。基于Test of Moderator的结果我们可以拒绝原假设，但仅有协变量绝对纬度会显著性影响结果。残余异质性的检验是显著的，说明存在模型中没有包含的因素会影响疫苗的效力。
根据Log后的相对风险显示，绝对纬度1度的上升对应-0.03单位风险的变化。为了进一步促进解释，我们可以控制发表年份，生成纬度数据进行预测。
这里，我们使用predict()函数，在newmods参数中指定协变量新的数值结果，并用transf参数控制结果的转换，使用addx=TRUE将所有结果放在一起。
predict(res, newmods = cbind(seq(from=10, to=60, by=10), 1970),
transf = exp, addx = TRUE)
## pred ci.lb ci.ub cr.lb cr.ub X.intrcpt X.ablat X.year
## 1 0.9345 0.5833 1.4973 0.4179 2.0899 1 10 1970
## 2 0.7062 0.5149 0.9686 0.3421 1.4579 1 20 1970
## 3 0.5337 0.4196 0.6789 0.2663 1.0697 1 30 1970
## 4 0.4033 0.2956 0.5502 0.1958 0.8306 1 40 1970
## 5 0.3048 0.1916 0.4848 0.1369 0.6787 1 50 1970
## 6 0.2303 0.1209 0.4386 0.0921 0.5761 1 60 1970
结果显示了预测的平均相对风险pred和对应的95%置信区间。只有当不使用transf参数时才会给出标准误。
可以看到，纬度越高，相对风险越小。更为一般地，我们画一个图来看：
preds = predict(res, newmods = cbind(0:60, 1970), transf = exp)
wi = 1/sqrt(dat$vi)
size = 0.5 + 3 * (wi - min(wi)) / (max(wi) - min(wi))
plot(dat$ablat, exp(dat$yi), pch = 19, cex = size,
xlab = "Absolute Latitude", ylab = "Relative Risk",
las = 1, bty = "l", log = "y")
lines(0:60, preds$pred)
lines(0:60, preds$ci.lb, lty = "dashed")
lines(0:60, preds$ci.ub, lty = "dashed")
abline(h = 1, lty = "dotted")
img
分类协变量
协变量经常是分类的，因此基于分类协变量的水平对研究进行分组是元分析经常遇到的问题。一种方法是分别在每个水平拟合随机效应模型。例如，我们可以根据治疗的分配方式将协变量分为几组：
rma(yi, vi, data = dat, subset = (alloc == "random"))
rma(yi, vi, data = dat, subset = (alloc == "alternate"))
rma(yi, vi, data = dat, subset = (alloc == "systematic"))
这里使用了subset参数进行子集的提取，然而这并不是个好办法，除非感兴趣的异质性在每一个水平都存在着，而且一旦这样分组，我们将处理更小的研究数据。
与此相反，我们使用单一的混合效应模型来研究。
首先，我们创建必要的变量：
dat$a.random = ifelse(dat$alloc == "random", 1, 0)
dat$a.alternate = ifelse(dat$alloc == "alternate", 1, 0)
dat$a.systematic = ifelse(dat$alloc == "systematic", 1, 0)
然后为每一个水平计算单独的效应：
rma(yi, vi, mods = cbind(a.random, a.alternate, a.systematic), intercept = FALSE,
data = dat)
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
## tau (square root of estimated tau^2 value): 0.6013
## I^2 (residual heterogeneity / unaccounted variability): 88.77%
## H^2 (unaccounted variability / sampling variability): 8.91
##
## Test for Residual Heterogeneity:
## QE(df = 10) = 132.3676, p-val < .0001
##
## Test of Moderators (coefficient(s) 1:3):
## QM(df = 3) = 15.9842, p-val = 0.0011
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## a.random -0.9658 0.2672 -3.6138 0.0003 -1.4896 -0.4420 ***
## a.alternate -0.5180 0.4412 -1.1740 0.2404 -1.3827 0.3468
## a.systematic -0.4289 0.3449 -1.2434 0.2137 -1.1050 0.2472
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
相对于手动编码，我们可以使用factor()函数处理，结果一致
rma(yi, vi, mods = ~factor(alloc) - 1, data = dat)
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
## tau (square root of estimated tau^2 value): 0.6013
## I^2 (residual heterogeneity / unaccounted variability): 88.77%
## H^2 (unaccounted variability / sampling variability): 8.91
##
## Test for Residual Heterogeneity:
## QE(df = 10) = 132.3676, p-val < .0001
##
## Test of Moderators (coefficient(s) 1:3):
## QM(df = 3) = 15.9842, p-val = 0.0011
##
## Model Results:
##
## estimate se zval pval ci.lb
## factor(alloc)alternate -0.5180 0.4412 -1.1740 0.2404 -1.3827
## factor(alloc)random -0.9658 0.2672 -3.6138 0.0003 -1.4896
## factor(alloc)systematic -0.4289 0.3449 -1.2434 0.2137 -1.1050
## ci.ub
## factor(alloc)alternate 0.3468
## factor(alloc)random -0.4420 ***
## factor(alloc)systematic 0.2472
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
根据这个结果，仅有随机治疗分配获得了显著的治疗效果。然而，想要检验是否分配因子变量真的统计显著，我们需要使用对模型使用不同的参数。
rma(yi, vi, mods = cbind(a.alternate, a.systematic), data = dat)
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
## tau (square root of estimated tau^2 value): 0.6013
## I^2 (residual heterogeneity / unaccounted variability): 88.77%
## H^2 (unaccounted variability / sampling variability): 8.91
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 10) = 132.3676, p-val < .0001
##
## Test of Moderators (coefficient(s) 2:3):
## QM(df = 2) = 1.7675, p-val = 0.4132
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.9658 0.2672 -3.6138 0.0003 -1.4896 -0.4420 ***
## a.alternate 0.4478 0.5158 0.8682 0.3853 -0.5632 1.4588
## a.systematic 0.5369 0.4364 1.2303 0.2186 -0.3184 1.3921
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
或者：
rma(yi, vi, mods = ~ relevel(factor(alloc), ref = "random"), data = dat)
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.3615 (SE = 0.2111)
## tau (square root of estimated tau^2 value): 0.6013
## I^2 (residual heterogeneity / unaccounted variability): 88.77%
## H^2 (unaccounted variability / sampling variability): 8.91
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 10) = 132.3676, p-val < .0001
##
## Test of Moderators (coefficient(s) 2:3):
## QM(df = 2) = 1.7675, p-val = 0.4132
##
## Model Results:
##
## estimate se
## intrcpt -0.9658 0.2672
## relevel(factor(alloc), ref = "random")alternate 0.4478 0.5158
## relevel(factor(alloc), ref = "random")systematic 0.5369 0.4364
## zval pval ci.lb
## intrcpt -3.6138 0.0003 -1.4896
## relevel(factor(alloc), ref = "random")alternate 0.8682 0.3853 -0.5632
## relevel(factor(alloc), ref = "random")systematic 1.2303 0.2186 -0.3184
## ci.ub
## intrcpt -0.4420 ***
## relevel(factor(alloc), ref = "random")alternate 1.4588
## relevel(factor(alloc), ref = "random")systematic 1.3921
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
这里使用了random分配作为参考，结果β0=−0.9658β0=−0.9658是使用随机分配的对数化相对风险估计，而β1β1与β2β2分别是使用alternate与systematic分配会相对多大。比如β0+β1=−0.5180β0+β1=−0.5180是使用alternate分配的相对风险，而β0+β2=−0.4289β0+β2=−0.4289是使用systematic分配的相对风险。但协变量检验H0:β1=β2=0H0:β1=β2=0并不统计显著（QM=1.77,df=2,p=0.41QM=1.77,df=2,p=0.41），表明分配的方法实际没有影响疫苗的平均效力。
Knapp与Hartung校正
下面介绍校正方法的使用：
rma(yi, vi, mods = ~ factor(alloc) + ablat, data = dat, knha = TRUE)
##
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1446 (SE = 0.1124)
## tau (square root of estimated tau^2 value): 0.3803
## I^2 (residual heterogeneity / unaccounted variability): 70.11%
## H^2 (unaccounted variability / sampling variability): 3.35
## R^2 (amount of heterogeneity accounted for): 53.84%
##
## Test for Residual Heterogeneity:
## QE(df = 9) = 26.2034, p-val = 0.0019
##
## Test of Moderators (coefficient(s) 2:4):
## F(df1 = 3, df2 = 9) = 3.4471, p-val = 0.0650
##
## Model Results:
##
## estimate se tval pval ci.lb
## intrcpt 0.2932 0.4188 0.7000 0.5016 -0.6543
## factor(alloc)random -0.2675 0.3624 -0.7381 0.4793 -1.0873
## factor(alloc)systematic 0.0585 0.3925 0.1490 0.8849 -0.8294
## ablat -0.0273 0.0095 -2.8669 0.0186 -0.0488
## ci.ub
## intrcpt 1.2407
## factor(alloc)random 0.5523
## factor(alloc)systematic 0.9463
## ablat -0.0058 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
omnibus检验基于m=3和k-p=9的自由度，而基于k-p=9自由度的t分布用于单个系数和检验的参考分布。添加参数btt = c(2, 3)针对性检验H0:β1=β2=0H0:β1=β2=0。
通常，该校正会计算一个更为保守的p值。
其他函数与方法
下表提供了许多函数与方法的概览，其中一些我们这里加以详细讨论。
| 函数 | 描述 |
|---|---|
| print() | 标准输出方法 |
| summary() | 提供模型统计量的额外输出方法 |
| coef() | 抽取计算模型的系数、对应的标准误、检验统计量、p值以及置信区间边界 |
| vcov() | 抽取模型系数的方差-协方差矩阵 |
| fitstats() | 抽取严格的log似然度、偏差、AIC与BIC |
| fitted() | 拟合值 |
| predict() | 拟合/预测值（没有置信区间），可以用于新数据 |
| blup() | 真实结果的最好的线性无偏预测（BLUPs） |
| residuals() | 原始残差 |
| rstandard() | 内部标准化残差 |
| rstudent() | 外部标准化残差 |
| hatvalues() | 抽取hat矩阵对角元素 |
| weights() | 抽取用于模型拟合的权重 |
| influence() | 大量案例与删除诊断 |
| leave1out() | 固定/随机效应模型留一诊断 |
| forest() | 森林图 |
| funnel() | 漏斗图 |
| radial() | 雷达图 |
| qqnorm() | 正态分位数-分位数图 |
| plot() | 模型对象的通用绘图函数 |
| addpoly() | 为森林图添加多边形 |
| ranktest() | 漏斗图不对称等级相关检验 |
| regtest() | 漏斗图不对称回归检验 |
| trimfill() | 修剪和填充方法 |
| confint() | 随机/混合模型（残差）异质性置信区间（也会获取模型系数置信区间） |
| cumul() | 固定/随机效应模型累积元分析 |
| anova() | 根据模型统计量和似然度进行模型比较 |
| permutest() | 模型系数置换检验 |
拟合/预测值
fitted()函数可以用来获取k个研究的拟合值。predict()函数提供了拟合值加上标准误、置信区间边界。
例如，我们可以首先拟合一个随机效应模型获取对数化的相对风险，然后转换估计的评估log相对风险（即mumu）为指数形式。
res = rma(yi, vi, data = dat)
predict(res, transf = exp, digits = 2)
## pred ci.lb ci.ub cr.lb cr.ub
## 0.49 0.34 0.70 0.15 1.55
原始和标准化残差
许多元分析包含的一些研究观察效应值是异常乃至极端的，可视化观察是一种区分的方式，但这种方式在处理1个或多个协变量时存在问题。而且，元分析的研究通常一旦变化数量，方差也就跟着变了。
一种更正式的方式是检查残差与对应标准误直接的关系。
在线性回归中存在许多中残差形式，可以被采纳到元分析。最重要的，rstandard()与rstudent提供了内部与外部标准化残差（residual()提供的是原始残差）。如果一个研究与模型拟合得很好，它的标准化残差会服从标准正态分布。
res = rma(yi, vi, mods = cbind(ablat, year), data = dat)
rstudent(res)
## resid se z
## 1 0.2229 0.7486 0.2978
## 2 -0.2828 0.6573 -0.4303
## 3 -0.3826 0.7501 -0.5100
## 4 -1.0900 0.7768 -1.4032
## 5 -0.0774 0.5192 -0.1490
## 6 0.4519 0.4283 1.0551
## 7 -1.4061 0.5416 -2.5961
## 8 0.2321 0.4842 0.4793
## 9 0.0974 0.4803 0.2027
## 10 -0.4496 0.4554 -0.9872
## 11 -0.0557 0.4656 -0.1197
## 12 1.1864 0.8084 1.4677
## 13 0.7972 0.3742 2.1302
影响点诊断
跟一般地回归分析类似，影响点会影响整个模型的拟合效果，我们需要进行诊断和处理。influence()函数提供下面的一些诊断措施用于元分析：
- 外部标准化残差
- DFFITS值
- Cook距离
- 协方差比率
- DFBETAS值
- 轮流移除单个研究然后评估τ2τ2
- 轮流移除单个研究然后为（残差）异质性评估统计量
- hat矩阵对象元素
- 模型拟合时给观测结果的权重
例如，对于之前我们处理的混合效应模型：
res = rma(yi, vi, mods = cbind(ablat, year), data = dat)
inf = influence(res)
inf
## $inf
## rstudent dffits cook.d cov.r tau2.del QE.del hat weight inf
## 1 0.298 0.1785 0.0348 1.800 0.1317 28.3 0.1725 3.37
## 2 -0.430 -0.2368 0.0620 1.921 0.1308 27.6 0.2367 4.81
## 3 -0.510 -0.1094 0.0125 1.235 0.1191 27.8 0.0487 2.79
## 4 -1.403 -2.9415 7.3179 3.522 0.0906 23.2 0.8082 11.23 *
## 5 -0.149 -0.0263 0.0032 2.634 0.1497 27.3 0.2483 9.07
## 6 1.055 0.8926 0.7205 1.362 0.0994 21.3 0.4061 12.48
## 7 -2.596 -0.6815 0.4173 0.238 0.0544 19.1 0.0766 4.40
## 8 0.479 0.3703 0.1899 2.998 0.1498 24.1 0.3627 12.80
## 9 0.203 0.1305 0.0237 2.207 0.1501 28.3 0.1030 8.78
## 10 -0.987 -0.3870 0.1470 1.070 0.1072 24.8 0.1310 7.99
## 11 -0.120 -0.0030 0.0052 2.834 0.1583 25.5 0.2214 11.92
## 12 1.468 0.2171 0.0469 0.927 0.1059 26.1 0.0235 2.28
## 13 2.130 0.8150 0.4994 0.218 0.0498 21.5 0.1612 8.06
##
## $dfbs
## intrcpt ablat year
## 1 0.1492 -0.0622 -0.1491
## 2 -0.0949 -0.1001 0.0963
## 3 -0.0280 -0.0403 0.0283
## 4 2.2248 -2.5380 -2.2161
## 5 0.0350 0.0063 -0.0354
## 6 0.5935 -0.0008 -0.5954
## 7 -0.2659 0.5368 0.2591
## 8 -0.0394 -0.2110 0.0431
## 9 0.0728 -0.0865 -0.0718
## 10 -0.1188 -0.0973 0.1194
## 11 0.0631 -0.0343 -0.0631
## 12 -0.0413 0.0465 0.0419
## 13 -0.8276 0.6269 0.8279
相比于输出，我们可以绘制：
plot(inf, plotdfbs = TRUE)
img
img
上图结果中，发现研究7和13为模型引入了额外的异质性（即轮流移除研究7或13然后评估τ2τ2，发现减少），但这两者只有少量的影响。但研究4对模型异质性有极少的影响，但却极大影响了模型的拟合（Cook距离图非常明显，在hat图中也反应了出来）。
对于没有协变量的模型，我们可以使用leave1out()函数重复拟合模型，每次留下一个研究。例如：
res = rma(yi, vi, data = dat)
leave1out(res, transf = exp, digits = 3)
## estimate zval pval ci.lb ci.ub Q Qp tau2 I2 H2
## 1 0.493 -3.722 0.000 0.340 0.716 151.583 0.000 0.336 93.226 14.762
## 2 0.520 -3.620 0.000 0.365 0.741 145.318 0.000 0.293 92.254 12.910
## 3 0.504 -3.692 0.000 0.350 0.725 150.197 0.000 0.321 92.935 14.155
## 4 0.533 -3.558 0.000 0.377 0.754 96.563 0.000 0.263 90.412 10.430
## 5 0.466 -3.984 0.000 0.320 0.678 151.320 0.000 0.328 92.763 13.819
## 6 0.491 -3.550 0.000 0.332 0.727 128.187 0.000 0.360 90.912 11.003
## 7 0.519 -3.631 0.000 0.365 0.740 145.830 0.000 0.293 92.278 12.950
## 8 0.452 -4.418 0.000 0.317 0.643 67.986 0.000 0.273 87.031 7.711
## 9 0.477 -3.769 0.000 0.324 0.701 152.205 0.000 0.349 93.213 14.735
## 10 0.520 -3.544 0.000 0.363 0.747 139.827 0.000 0.299 92.232 12.874
## 11 0.469 -3.871 0.000 0.319 0.688 151.466 0.000 0.340 91.811 12.211
## 12 0.468 -4.173 0.000 0.327 0.668 150.787 0.000 0.308 92.678 13.658
## 13 0.460 -4.191 0.000 0.319 0.661 149.788 0.000 0.304 92.344 13.062
绘图函数（森林、漏斗、雷达与QQ图）
metafor包提供了一些频繁用于元分析的图形函数。
我们在之前就已经展示了森林图的绘制，下面用另一个例子展示如何绘制常见的效应大小与对象的取样方差，并解释addpoly()如何添加多边形。
forest(dat$yi, dat$vi, atransf = exp, ylim = c(-3.5, 16),
at = log(c(0.05, 0.25, 1, 4, 20)), xlim = c(-9, 7),
slab = paste(dat$author, dat$year, sep = ", "))
res = rma(yi, vi, mods = cbind(ablat), data = dat)
preds = predict(res, newmods = c(10, 30, 50))
addpoly(preds$pred, sei = preds$se, atransf = exp,
mlab = c("10 Degrees", "30 Degrees", "50 Degrees"))
text(-9, 15, "Author(s) and Year", pos = 4, font =2)
text(7, 15, "Relative Risk [95% CI]", pos = 2, font =2)
abline(h=0)
img
函数funnel()可以用来创建漏斗图，用于诊断异质性的出现以及确定的出版偏差。对于没有协变量的模型，该图为观察到的结果作为水平轴，对应的标准误作为垂直轴。一条垂直线显示了基于模型的结果。伪置信区间和它的1.96 SE边界绘制出来了，SE是垂直轴的标准误。对于涉及协变量的模型，该图显示了残差作为水平轴，标准误作为垂直轴。
下面分别展示没有和有协变量的漏斗图：
res = rma(yi, vi, data = dat)
funnel(res, main = "Random-Effects Model")
img
res = rma(yi, vi, mods = cbind(ablat), data = dat)
funnel(res, main = "Mixed-Effects Model")
img
雷达图用来评估有不同精度结果之间的一致性。固定效应模型与随机效应模型的评估方式有所不同，在图中也会显示。
res = rma(yi, vi, data = dat, method = "FE")
radial(res, main = "Fixed-Effects Model")
img
res = rma(yi, vi, data = dat, method = "REML")
radial(res, main = "Random-Effects Model")
img
qqnorm()函数创建的Q-Q正态图是元分析非常有用的诊断工具。
res = rma(yi, vi, data = dat)
qqnorm(res, main = "Random-Effects Model")
img
res = rma(yi, vi, mods = cbind(ablat), data = dat)
qqnorm(res, main = "Mixed-Effects Model")
img
漏斗图非对称检验
ranktest()与regtest()函数用来执行秩（等级）相关检验与回归检验。
对于回归检验，函数参数：
regtest(x, model = "rma", predictor = "sei", ni = NULL, ...)
x是拟合的模型对象。我们通过model参数选定模型，model= "lm"是加权回归，而model="rma"是标准的元分析模型（默认）。predictor指定检验的度量，比如标准误（默认）、取样方差（vi），样本大小的倒数（ninv）等。
res = rma(yi, vi, data = dat)
regtest(res, model = "lm")
##
## Regression Test for Funnel Plot Asymmetry
##
## model: weighted regression with multiplicative dispersion
## predictor: standard error
##
## test for funnel plot asymmetry: t = -1.4013, df = 11, p = 0.1887
res = rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat,
measure = "RR", mods = cbind(ablat, year))
regtest(res, predictor = "ni")
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: sample size
##
## test for funnel plot asymmetry: z = 0.7470, p = 0.4550
