Topic 1. SCI 文章中 Meta 分析之 metafor_桓峰基因的博客-程序员宅基地_r语言metafor包

技术标签: r语言  开发语言  

最近搞了一下 Meta 分析,蛮有意思的,掌握了这些技能未来发文章不用愁,Meta 分析也是医学上很重要的方法之一,这期基于现有数据,分析几个经典的图形,包括森林图(forest)、漏洞图(funnel)、星状图(radial)、拉贝图(labbe)、以及Q-Q正态分位图(Q-Q normal)。

01 Meta 分析的概念

1976年,Glass首次将这一概念命名为 Meta-analysis (荟萃分析)。并定义为一种对不同研究结果进行收集、合并及统计分析的方法。这种方法逐渐发展为一门新兴学科——“循证医学”的主要内容和研究手段。荟萃分析的主要目的是将以往的研究成果更为客观地综合反映出来。研究者并不进行原始的研究,而是将研究以获得的结果进行综合分析。

图片

02 输入数据模式

图片 我们使用 R 包 metafor,由Viechtbauer开发,除可完成二分类和连续性变量的 Meta 分析外,还可进行回归分析,这个软件包是唯一可以进行混合效应模型拟合运算的程序包,其包括:

  1. 单个变量

  2. 多分类变量

  3. 连续性变量

还可以检测模型系数并获得可信区间,以及对参数进行精准检验如置换检验(permutation)。

程序包安装,如下:

install.packages("metafor")
library("metafor")

读取软件包数据,来自13项检验卡介苗抗结核效果的研究的结果,meta分析的目的是检查卡介苗预防结核病的总体有效性,并检查可能影响效果大小的缓减剂。该数据集已经在一些出版物中被用来说明元分析方法,如下:

data("dat.bcg")
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

#######数据说明
?dat.bcg
  Format
The data frame contains the following columns:

trial  numeric  trial number
author  character  author(s)
year  numeric  publication year
tpos  numeric  number of TB positive cases in the treated (vaccinated) group
tneg  numeric  number of TB negative cases in the treated (vaccinated) group
cpos  numeric  number of TB positive cases in the control (non-vaccinated) group
cneg  numeric  number of TB negative cases in the control (non-vaccinated) group
ablat  numeric  absolute latitude of the study location (in degrees)
alloc  character  method of treatment allocation (random, alternate, or systematic assignment)  

这13项研究以下列表格提供数据细节,如下:


TB positive TB negative
vaccinated group tpos tneg
control group cpos cneg

03 实例解析

计算log相对风险和相应的抽样方差,如下:

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

演示公式界面,如下:

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
51    13   C   +    29
52    13   C   - 17825
  • 随机效应模型,如下:
#默认模型为REML模型Random-effect model using Restricted maximum likelihood estimator
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

res <- rma(ai = tpos, bi = tneg, ci = cpos, 
           di = cneg, data = dat, measure = "RR")
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

异质性的置信区间可用metafor中confint()可以求异质性系数(tau2,tau,I2(%),H^2)的置信区间,如下:

### confidence intervals for tau^2, I^2, and H^2

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 
# 估算效应量的预测/可信度区间
predict(res, digits = 3)
 pred    se  ci.lb  ci.ub  pi.lb pi.ub 
 -0.715 0.180 -1.067 -0.362 -1.867 0.438
  • forest 森林图

森林图是以统计指标和统计分析方法为基础,用数值运算结果绘制出的图型。它在平面直角坐标系中,以一条垂直的无效线(横坐标刻度为1或0)为中心,用平行于横轴的多条线段描述了每个被纳入研究的效应量和可信区间,用一个棱形(或其它图形)描述了多个研究合并的效应量及可信区间。它非常简单和直观地描述了Meta分析的统计结果,是Meta分析中最常用的结果表达形式。如下:

### forst plot
forest(res, slab = paste(dat$author, dat$year, sep = ", "), 
       xlim = c(-16, 6), at = log(c(.05, .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 = .75)
op <- par(cex = .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)
par(op)

图片

混合效应模型,如下:

### mixed-effects model with absolute latitude and publication year as moderators

res <- rma(yi, vi, mods = cbind(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 (coefficients 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

### predicted average relative risks for 10, 20, 30, 40, and 50 degrees 
### absolute latitude holding publication year constant at 1970

predict(res, newmods = cbind(seq(from = 10, to = 60, by = 10), 1970), 
        transf = exp, addx = TRUE)

### plot of the predicted average relative risk as a function of absolute latitude
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")

图片

  • funnel 漏斗图

漏斗图是由Light等在1984年提出,一般以单个研究的效应量为横坐标,样本含量为纵坐标做的散点图。效应量可以为RR、OR和死亡比或者其对数值等。理论上讲,被纳入Meta分析的各独立研究效应的点估计,在平面坐标系中的集合应为一个倒置的漏斗形,因此称为漏斗图。

  1. 样本量小,研究精度低,分布在漏斗图的底部,向周围分散;

  2. 样本量大,研究精度高,分布在漏斗图的顶部,向中间集中。

实际使用时,做Meta分析的研究个数较少时不宜做漏斗图,一般推荐Meta分析的研究个数在10个及以上才需做漏斗图。漏斗图主要用于观察某个系统评价或Meta分析结果是否存在偏倚,如发表偏倚或其他偏倚,如下图:

### funnel plots
res <- rma(yi, vi, data = dat)
funnel(res, main = "Random-Effects Model")
res <- rma(yi, vi, mods = cbind(ablat), data = dat)
funnel(res, main = "Mixed-Effects Model")
#######Begger's检验
ranktest(res)
Rank Correlation Test for Funnel Plot Asymmetry

Kendall's tau = 0.0256, p = 0.9524
##########Egger's检验
regtest(res)
Regression Test for Funnel Plot Asymmetry

Model:     mixed-effects meta-regression model
Predictor: standard error

Test for Funnel Plot Asymmetry: z = -0.4623, p = 0.6439

图片

  • radial 星状图

主要反映各研究的异质性,从而发现异质性的点。弧线对应的效应评估大小分布。如下:

### radial plots
res <- rma(yi, vi, data = dat, method = "REML")
radial(res, main = "Random-Effects Model")
res <- rma(yi, vi, data = dat, method = "ML")
radial(res, main = "Mixed-Effects Model")

图片

  • Q-Q分位图

绘制预测结果,观察是否在置信区间内部,如下:

### qq-normal plots
res <- rma(yi, vi, data = dat)
qqnorm(res, main = "Random-Effects Model")
res <- rma(yi, vi, mods = cbind(ablat), data = dat)
qqnorm(res, main = "Mixed-Effects Model")

图片

  • labbe 拉贝图

该函数计算两组的臂水平结果(例如,治疗和控制),并将它们相互比较。特别是功能块的原始比例两组互相在分析风险差异,比例在分析的日志(日志)风险比率,日志赔率在分析优势比(日志),平方根反正弦转换比例在分析平方根反正弦转换风险差异,分析发病率差异时用原始发病率,分析发病率比时用发病率的对数(log),分析发病率差异时用发病率的平方根转换,如下:

res <- rma(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

### default plot
labbe(res)

### funnel plot with risk values on the x- and y-axis and add grid
labbe(res, transf=exp, grid=TRUE)
# }

图片

04 分析结果解读

可以看到I^2为71.98%,属于高度异质性,可采用随机效应模型。异质性低的时候可以采用固定效应模型和随机效应模型,结果差别不大,但高异质性只能选择随机效应模型,否则会使结果外推性受到约束。此处选择随机效应模型是出于保守情况考虑。

random-effect model是基于跨研究间存在异质性的假设,该合并模型承认研究间异质性的存在,但是不对异质性加以处理;如果纳入合并的研究间存在异质性,尽管未达到我们常规设定的I^2>50%,但是在用 fixed-effect model合并时,默认运算直接忽略这一部分异质性的存在,这样合并的结果会造成假阳性误差,而选用random-effect model合并时,尽管不处理异质性,但是其默认运算承认异质性的存在,合并结果更可信!

从漏斗图以及Begg’s检验及Egger’s 检验的结果可以看出,P值都是大于0.05的,也就意味着没有发表偏倚。

Meta 分析同样是发文章的一种选择,搞清这些图解以及内部的方法,至于SCI能发到多少分,还需要看我们研究的课题的新颖程度,而方法也就是这些,表现形式基本都有,只要您有好的数据,顺利发表个SCI高分还是指日可待的,下期在聊聊其他软件!

关注公众号,扫码进群,免费解答,后期会有免费直播教程,敬请期待!
在这里插入图片描述

Reference:

  1. Viechtbauer, W. (2010). Conducting Meta-Analyses in R with the metafor Package. Journal of Statistical Software, 36(3), 1–48.

  2. Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., & Mosteller, F. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Meta-analysis of the published literature. Journal of the American Medical Association, 271(9), 698–702.

  3. Berkey, C. S., Hoaglin, D. C., Mosteller, F., & Colditz, G. A. (1995). A random-effects regression model for meta-analysis. Statistics in Medicine, 14(4), 395–411.

  4. van Houwelingen, H. C., Arends, L. R., & Stijnen, T. (2002). Advanced methods in meta-analysis: Multivariate approach and meta-regression. Statistics in Medicine, 21(4), 589–624.

  5. Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48.

桓峰基因

生物信息分析,SCI文章撰写及生物信息基础知识学习:R语言学习,perl基础编程,linux系统命令,Python遇见更好的你

34篇原创内容

公众号

图片

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/weixin_41368414/article/details/122830991

智能推荐

mysql集群-主从同步(自用笔记)_the_pure的博客-程序员宅基地_mysql集群修改master_port

主从同步原理:TODO数据库准备:至少需要两台服务器,以后肯定是独立两台电脑,当然也可以使用虚拟机。本次采用安装多个服务(使用不同的端口)来代替服务器.一、环境准备:准备主库master:1.解压mysql安装包。2.修改my.ini文件:在[mysqld]下添加如下配置[mysqld]#二进制日志文件,主从重点文件log-bin=mysql-bin#服务ID随便取但是要唯一server-id=1#端口号 在单台机器上通过不同端口来模拟集群...

Typora添加LaTeX内联公式_粽子小黑的博客-程序员宅基地_内联公式

在Typora中添加LaTeX内联公式所谓LaTeX内联公式:O(nlogn)O(nlogn)O(nlogn),公式和文字可以显示在一行,在数学公式两端各敲一个$。如果想要把公式显示为单独一行a+b=ca+b=ca+b=c如上所示,需要在公式两端各加上两个$$。Typora 是一个离线的Markdown编辑软件。如想写如H2O等Markdown拓展语法,需要进行偏好设置。如下图:根据自己...

HDU-1533 Going Home(二分图匹配)_AC_Arthur的博客-程序员宅基地_hdu1533时间复杂度

最近开始做最小费用流的题目,该题是二分图完美匹配下的最小权匹配,所谓完美匹配就是说从源点流入的总流量等于从汇点流出的总流量,在这种状态下的最小费用 。 那么显然是要套用最小费用流模板,另外二分图匹配的第一步就是要划分集合,划分两个集合,集合A与源点相连,集合B与汇点相连,至于容量和权值就要依据题目而定 。比如该题,因为每个小人恰好能对应一个房子,所以每个小人与汇点的容量为1,房子与汇点的容

loadrunner Web_类函数之web_add_auto_header()_testingstar的博客-程序员宅基地_web_add_auto_header

web_add_auto_header()--常用函数将指定的标头添加到所有后续HTTP请求。该函数与Run-time SettingsàBrowser/Browser Emuletion/User-Agent里面的设置功能相似,如下图; 在以下示例中,web_add_auto_header函数通知服务器重播接受压缩数据:C语言写法:web_add_auto_header(“...

简述js中的面向对象编程_天心天地生的博客-程序员宅基地_介绍面向对象编在js中的运用

1.object类型到目前为止,我们看到的大多数引用类型值都是 Object 类型的实例;而且, Object 也是 ECMAScript 中使用最多的一个类型。虽然 Object 的实例不具备多少功能,但对于在应用程序中存储 和传输数据而言,它们确实是非常理想的选择。2.继承继承的基本思想:利用原型让一个引用类型继承另一个引用类型的属性和方法.3.构造函数(constructor...

随便推点

sublime 自动触发补全_信行合一的博客-程序员宅基地

Preferences- settings user -Preferences.sublime-settings中设置自动触发:{"auto_complete": true,"auto_mathch_enabled": true,"color_scheme": "Packages/User/SublimeLinter/Monokai (SL).tmTheme","f

安装laravel缺少autoload.php_不是皮卡chiu的博客-程序员宅基地

亲测按照 这篇文章 安装composer和laravel是否可行到了这一步问题就出现了:一旦安装完成后,就可以使用 laravel new 命令在你指定的目录中建立一份全新安装的 Laravel 应用。例如: laravel new blog 命令会在当前目录下建立一个名为 blog 的目录, 此目录里面存放着全新安装的 Laravel ,并且所有依赖包也已经安装好了。此方法的安装速度会比通...

定时任务框架Quartz_灰灰2016的博客-程序员宅基地_c# quartz demo

深入解读Quartz的原理https://blog.csdn.net/scgyus/article/details/79360316Quartz入门与Demo搭建https://blog.csdn.net/noaman_wgs/article/details/80984873

Spring系列之集成Druid连接池及监控配置_程序员阿牛啊的博客-程序员宅基地_datasource.druid.url

前言前一篇文章我们熟悉了HikariCP连接池,也了解到它的性能很高,今天我们讲一下另一款比较受欢迎的连接池:Druid,这是阿里开源的一款数据库连接池,它官网上声称:为监控而生!他可以实现页面监控,看到SQL的执行次数、时间和慢SQL信息,也可以对数据库密码信息进行加密,也可以对监控结果进行日志的记录,以及可以实现对敏感操作实现开关,杜绝SQL注入,下面我们详细讲一下它如何与Spring集成,并且顺便了解一下它的监控的配置。 文章要点: Spring集成Druid 监控Filters配置(

使用JiaoZiVideoPlayer播放视频方向不对问题_tenda的博客-程序员宅基地_jzvdstd怎么设置横屏或者竖屏播放

JzvdStd--JiaoZiVideoPlayer 播放视频横竖不对问题先上两张图片看问题:两个视频播放始终都是一个方向,找了好久都没找到解决方案,查看源码也没找到啥提示。最终还是幸亏公司里有个大佬给解决了,所以在此记录下,以后用这个播放器的肯定多少会遇到这个问题,废话不多说,上代码。public class JZplayer extends JzvdStd { public...

Maven 各命令执行流程解析和说明_风舞松林涧的博客-程序员宅基地

1: 本机安装Maven,可参考其它网站或者我的博客,修改MAVEN( setting.xml文件) 存放本机资源库的位置:D:\Repositories\Maven       命令说明,对应的命令执行流程可见流程如图3、图4 console台打印信息,能知一二。  mvn compile:编译源代码,生成对应的CLASS文件,执行流程可见流程如图...

推荐文章

热门文章

相关标签