大家应该很熟悉meta分析,所谓meta分析就是一个全面收集所有相关研究并逐个进行严格评价和分析,再用定量合成的方法对资料进行统计学处理得出综合结论的整个过程。今天我们给大家介绍一个在R语言中进行meta分析的工具metafor包。我们通过这个包把相应的meta分析的常规的一些图为大家介绍下。
1. 森林图,主要是对研究的一致性进行评估的可视化展示形式,以竖线为界,总结结果在线左认为是研究的因素降低,或者对研究的因素不利。此处我们使用此包自带的卡介苗抵抗肺结核(TB)的研究数据进行森林他绘制。
首先看下数据构成:
具体的每一个变量的解释如下:
然后,我们利用meta分析模型对数据进行评估,此处我们需要用到一个函数rma,主要是利用线性模型对数据进行一致性评估。
其中主要的参数有meature,一般我们都选用‘RR’,‘OR‘(具体解释可以自行百度或者参见我们前面的教程);
method,具体的方法如下:
还有一个参数slab他主要是列出你想要在图中显示的每个研究的标签。
接下来我们看下实例:
library(metafor)
res <-rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, measure="RR",
slab=paste(author, year, sep=","), method="REML")
接下来我们看下森林图的绘制函数forest:
其中主要的参数有at用于标记X轴的值;ilab,ilab.pos用于显示研究的相关数据及相应的位置。具体实例如下:
par(mar=c(4,4,1,2))
forest(res,xlim=c(-16, 6), at=log(c(0.05, 0.25, 1, 4)), atransf=exp,
ilab=cbind(dat.bcg$tpos, dat.bcg$tneg,dat.bcg$cpos, dat.bcg$cneg),
ilab.xpos=c(-9.5,-8,-6,-4.5), cex=0.75,ylim=c(-1, 27),
order=order(dat.bcg$alloc),rows=c(3:4,9:15,20:23),
xlab="Risk Ratio",mlab="", psize=1)
接下来就是对森林图的进一步数据的完善,美化,代码如下:
text(-16, -1, pos=4, cex=0.75,bquote(paste("RE Model for All Studies (Q = ",
.(formatC(res$QE, digits=2, format="f")), ", df = ",.(res$k - res$p),
", p = ", .(formatC(res$QEp, digits=2, format="f")),"; ", I^2, " = ",
.(formatC(res$I2, digits=1, format="f")), "%)")))
### set font expansion factor (as in forest()above) and use bold italic
### font and save original settings in object'op'
op <- par(cex=0.75, font=4)
### add text for the subgroups
text(-16, c(24,16,5), pos=4, c("SystematicAllocation",
"RandomAllocation",
"AlternateAllocation"))
### switch to bold font
par(font=2)
### add column headings to the plot
text(c(-9.5,-8,-6,-4.5), 26, c("TB+","TB-", "TB+", "TB-"))
text(c(-8.75,-5.25), 27, c("Vaccinated","Control"))
text(-16, 26, "Author(s) andYear", pos=4)
text(6, 26, "Risk Ratio [95%CI]", pos=2)
### set par back to the original settings
par(op)
### fit random-effects model in the threesubgroups
res.s <- rma(ai=tpos, bi=tneg, ci=cpos,di=cneg, data=dat.bcg, measure="RR",
subset=(alloc=="systematic"), method="REML")
res.r <- rma(ai=tpos, bi=tneg, ci=cpos,di=cneg, data=dat.bcg, measure="RR",
subset=(alloc=="random"), method="REML")
res.a <- rma(ai=tpos, bi=tneg, ci=cpos,di=cneg, data=dat.bcg, measure="RR",
subset=(alloc=="alternate"), method="REML")
### add summary polygons for the threesubgroups
addpoly(res.s, row=18.5, cex=0.75, atransf=exp,mlab="")
addpoly(res.r, row= 7.5, cex=0.75, atransf=exp,mlab="")
addpoly(res.a, row= 1.5, cex=0.75, atransf=exp,mlab="")
### add text with Q-value, dfs, p-value, andI^2 statistic for subgroups
text(-16, 18.5, pos=4, cex=0.75,bquote(paste("RE Model for Subgroup (Q = ",
.(formatC(res.s$QE, digits=2, format="f")), ", df =", .(res.s$k - res.s$p),
", p = ", .(formatC(res.s$QEp, digits=2,format="f")), "; ", I^2, " = ",
.(formatC(res.s$I2, digits=1, format="f")), "%)")))
text(-16, 7.5, pos=4, cex=0.75,bquote(paste("RE Model for Subgroup (Q = ",
.(formatC(res.r$QE, digits=2, format="f")), ", df =", .(res.r$k - res.r$p),
", p = ", .(formatC(res.r$QEp, digits=2,format="f")), "; ", I^2, " = ",
.(formatC(res.r$I2, digits=1,format="f")), "%)")))
text(-16, 1.5, pos=4, cex=0.75,bquote(paste("RE Model for Subgroup (Q = ",
.(formatC(res.a$QE, digits=2, format="f")), ", df =", .(res.a$k - res.a$p),
", p = ", .(formatC(res.a$QEp, digits=2,format="f")), "; ", I^2, " = ",
.(formatC(res.a$I2, digits=1, format="f")), "%)")))
以上的森林图对子集合照样做了相关的分析。如果一个整体的进行累积的结果展示。我们可以看下下面的代码:
par(mar=c(4,4,1,2))
### calculate(log) risk ratios and corresponding sampling variances
dat <-escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg,data=dat.bcg)
### fitrandom-effects models
res <-rma(yi, vi, data=dat, slab=paste(author, year, sep=", "))
### cumulativemeta-analysis (in the order of publication year)
tmp <-cumul(res, order=order(dat$year))
### cumulativeforest plot
forest(tmp,xlim=c(-4,2), at=log(c(0.125, 0.25, 0.5, 1, 2)),
atransf=exp, digits=c(2,3), cex=0.75)
### switch tobold font
par(cex=0.75,font=2)
### add columnheadings to the plot
text(-4, 15,"Author(s) and Year", pos=4)
text( 2, 15,"Risk Ratio [95% CI]", pos=2)
2. CaterpillarPlot,这个中文实在是没发现是啥意思。就叫他履带图吧,这个图其实就是森林图的一种,去掉了所有的标签,但是如果在研究的数量比较大时,此图就发挥了它的作用。接下来我们看下具体的应用:
set.seed(5132)
k <- 250
vi <-rchisq(k, df=1) * .03
yi <-rnorm(k, rnorm(k, 0.5, 0.4), sqrt(vi))
### fit RE model
res <-rma(yi, vi)
### create plot
forest(yi, vi,
xlim=c(-2.5,3.5), ### adjust horizontal plot regionlimits
subset=order(yi), ### order by size of yi
slab=NA, annotate=FALSE, ### removestudy labels and annotations
efac=0, ### remove vertical bars at end of CIs
pch=19, ### changing point symbol tofilled circle
col="gray40", ### change color of points/CIs
psize=2, ### increase point size
cex.lab=1, cex.axis=1, ### increase size of x-axis title/labels
lty=c("solid","blank")) ### remove horizontal line at top of plot
### draw pointsone more time to make them easier to see
points(sort(yi),k:1, pch=19, cex=0.5)
### add summarypolygon at bottom and text
addpoly(res,mlab="", annotate=FALSE, cex=1)
text(-2, -2,"RE Model", pos=4, offset=0, cex=1)
3. 漏斗图,主要是反应的样本效应的波动性,一般横轴为效应量(RR,RO等),纵轴为样本量或者标准误。那么随着样本的增加标准误和效应量相对会越来越集中,也就形成了倒置的漏斗;反之,那么说明数据的结果存在一定的问题。其中Y轴的数据选择:
具体的实现实例如下:
library(metafor)
### fitfixed-effects model
res <-rma(yi, vi, data=dat.hackshaw1998, measure="OR",method="FE")
### set up 2x2array for plotting
par(mfrow=c(2,2))
### draw funnelplots
funnel(res,main="Standard Error")
funnel(res,yaxis="vi", main="Sampling Variance")
funnel(res,yaxis="seinv", main="Inverse Standard Error")
funnel(res,yaxis="vinv", main="Inverse Sampling Variance")
4. 散点图,主要是绘制预测结果,观察是否在置信区间内部。由于其运算的简单,我们就不把程序一一展开了。具体实例如下:
par(mar=c(5,5,1,2))
### calculate(log) risk ratios and corresponding sampling variances
dat <-escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg,data=dat.bcg)
### fitmixed-effects model with absolute latitude as predictor
res <-rma(yi, vi, mods = ~ ablat, data=dat)
### calculatepredicted risk ratios for 0 to 60 degrees absolute latitude
preds <-predict(res, newmods=c(0:60), transf=exp)
### calculatepoint sizes by rescaling the standard errors
wi <- 1/sqrt(dat$vi)
size <- 0.5 + 3.0 * (wi - min(wi))/(max(wi) -min(wi))
### plot therisk ratios against absolute latitude
plot(dat$ablat,exp(dat$yi), pch=19, cex=size,
xlab="Absolute Latitude",ylab="Risk Ratio",
las=1, bty="l",log="y")
### addpredicted values (and corresponding CI bounds)
lines(0:60,preds$pred)
lines(0:60,preds$ci.lb, lty="dashed")
lines(0:60,preds$ci.ub, lty="dashed")
### dotted lineat RR=1 (no difference between groups)
abline(h=1,lty="dotted")
### labels somepoints in the plot
ids <-c(4,7,12,13)
pos <-c(3,3,1,1)
text(dat$ablat[ids],exp(dat$yi)[ids], ids, cex=0.9, pos=pos)
5. Plot ofInfluence Diagnostics 主要是评估模型的中研究质量,从而发现对分析主要的影响的研究以及偏差很大的研究。其中residuals主要越大说明研究偏差越大;diffits主要评估偏离均值的大小,越大偏离均值越大;cook.d也可以称马氏距离,其值越大则对研究影响越大;cov.r协方差比率,其绝对值越大则相互作用越强;tau2.del,QE.del主要是发现存在异质性的研究;hat,weight主要体现标准差系数,反应数据的离散程度。
### fitrandom-effects model with r-to-z transformed correlations
res <-rma(ri=ri, ni=ni, measure="ZCOR", data=dat.mcdaniel1994)
### calculateinfluence diagnostics
inf <-influence(res)
### plot theinfluence diagnostics
plot(inf,layout=c(8,1))
6. Radial(Galbraith) Plot 主要反应各研究的异质性,从而发现异质性的点。此处以固定效应模型为例。
### adjustmargins so the space is better used
par(mar=c(5,4,0,2))
### fitfixed-effects model
res <-rma(yi, vi, data=dat.hackshaw1998, method="FE")
### draw radialplot
radial(res)
图中横轴是标准差/标准差倒数,纵轴是效应评估结果的对数值。弧线对应的效应评估大小分布。图中右侧的直线指示了1附近的效应值。
7. Baujat图主要通过Q检验进行研究异质性分析。横轴是每个样本相对于总体的异质性大小;纵轴是表示有或者没有这个研究之后对总体效应的标准化的平方差。
par(mar=c(5,4,2,2))
### get datafrom Pignon et al. (2000)
dat <-dat.pignon2000
### computeestimated log hazard ratios and sampling variances
dat$yi <-with(dat, OmE/V)
dat$vi <-with(dat, 1/V)
###meta-analysis based on all 65 trials
res <-rma(yi, vi, data=dat, method="FE", slab=id)
### createBaujat plot
baujat(res,xlim=c(0,20), ylim=c(0,0.2))
8. L’Abbe图,同样是异质性可视化定性图。图中圈的大小代表了研究样本的大小。横轴是对照组效应值对数,纵轴是实验组效应值对数。分布在对角实线代表没有差异的研究,分布在线下代表风险比对照组小。
### decreasemargins so the full space is used
par(mar=c(5,4,1,2))
### fitrandom-effects model
res <-rma(ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg, measure="RR")
### draw L'Abbéplot
labbe(res)