首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >R语言实现meta分析过程中的可视化展示

R语言实现meta分析过程中的可视化展示

作者头像
一粒沙
发布2019-07-31 14:19:17
3.5K0
发布2019-07-31 14:19:17
举报
文章被收录于专栏:R语言交流中心R语言交流中心

大家应该很熟悉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)

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-01-07,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 R语言交流中心 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档