#エビデンス計量評価学 Rを使って統合を行う
#必要なパッケージの読み込み
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.5
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(metafor)
## 要求されたパッケージ Matrix をロード中です
##
## 次のパッケージを付け加えます: 'Matrix'
##
## 以下のオブジェクトは 'package:tidyr' からマスクされています:
##
## expand, pack, unpack
##
## 要求されたパッケージ metadat をロード中です
##
## Loading the 'metafor' package (version 3.8-1). For an
## introduction to the package please type: help(metafor)
library(meta)
## Loading 'meta' package (version 6.0-0).
## Type 'help(meta)' for a brief overview.
## Readers of 'Meta-Analysis with R (Use R!)' should install
## older version of 'meta' package: https://tinyurl.com/dt4y5drs
#データセットの読み込み
ds1 <- read_csv("homework.csv")
## Rows: 10 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Source
## dbl (5): N1, Event1, N2, Event2, N
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ds2 <- ds1 %>%
mutate(pos.i = Event1,
pos.c = Event2,
total.i = N1,
total.c = N2,
neg.i = N1 - Event1,
neg.c = N2 - Event2)
2.必要な統計量(平均値差や標準誤差)を計算する。
ds3 <-
escalc(measure="RR",
ai=pos.i,
bi=neg.i,
ci=pos.c,
di=neg.c,
data = ds2,
append=TRUE)
ds3 |> as.data.frame()
## Source N1 Event1 N2 Event2 N pos.i pos.c total.i total.c neg.i
## 1 Cook 1990 25 12 24 11 49 12 11 25 24 13
## 2 Conan McCaul 2003 35 6 36 4 71 6 4 35 36 29
## 3 Shin 2007 27 9 36 10 63 9 10 27 36 18
## 4 Shin 2007 26 5 35 4 61 5 4 26 35 21
## 5 Dabu-Bondoc 2012 30 14 32 20 62 14 20 30 32 16
## 6 Patel P 2013 87 9 75 7 162 9 7 87 75 78
## 7 Jain 2013 75 23 75 43 150 23 43 75 75 52
## 8 Pin On 2016 42 2 44 2 86 2 2 42 44 40
## 9 Mishra 2017 50 6 50 14 100 6 14 50 50 44
## 10 Rao 2017 56 17 59 32 115 17 32 56 59 39
## neg.c yi vi
## 1 13 0.04618938 0.09257576
## 2 32 0.43363599 0.36031746
## 3 26 0.18232156 0.14629630
## 4 31 0.52039507 0.38296703
## 5 12 -0.29213642 0.05684524
## 6 68 0.10289442 0.22914067
## 7 32 -0.62570590 0.04006741
## 8 42 0.04652002 0.95346320
## 9 36 -0.84729786 0.19809524
## 10 27 -0.58033681 0.05526723
統合する統計量を指定する 平均値差: “MD” リスク比: “RR” オッズ比: “OR” リスク差: “RD” 異質性(τ)を見積もる様々な方法: 良く使われるのはDLとREM L DL : DerSimonian-Laird REML : Restricted Maximum-Likelihood
ds4 <- metafor::rma(measure="RR", yi=yi, vi=vi,
data = ds3,
method="DL")
ds4
##
## Random-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0445 (SE = 0.0770)
## tau (square root of estimated tau^2 value): 0.2109
## I^2 (total heterogeneity / total variability): 27.93%
## H^2 (total variability / sampling variability): 1.39
##
## Test for Heterogeneity:
## Q(df = 9) = 12.4872, p-val = 0.1872
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2695 0.1315 -2.0496 0.0404 -0.5272 -0.0118 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
metafor::forest(ds4, atransf = exp, digits=2)
##metaパッケージをつかった統合と可視化 metaforパッケージにおいてForest plotをお手軽に綺麗にするのは難しいので、ここはあきらめてmetaパッケージを使ったほうが無難。 metaパッケージでは「2.必要な統計量(平均値差や標準誤差)」と「3.上記で算出した統計量を統合する」は同時に行う
dic.HDS <- meta::metabin(event.e = pos.i,
n.e = total.i,
event.c = pos.c,
n.c = total.c,
data = ds2,
studlab = paste(Source),
fixed = FALSE,
random = TRUE,
method.tau = "DL",
sm = "RR")
dic.HDS
## Number of studies combined: k = 10
## Number of observations: o = 919
## Number of events: e = 250
##
## RR 95%-CI z p-value
## Random effects model 0.7639 [0.5902; 0.9886] -2.05 0.0406
##
## Quantifying heterogeneity:
## tau^2 = 0.0446; tau = 0.2113; I^2 = 28.0% [0.0%; 65.4%]; H = 1.18 [1.00; 1.70]
##
## Test of heterogeneity:
## Q d.f. p-value
## 12.50 9 0.1866
##
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - DerSimonian-Laird estimator for tau^2
## - Mantel-Haenszel estimator used in calculation of Q and tau^2 (like RevMan 5)
まずはシンプルな作図(これだけでも十分綺麗)
meta::forest(dic.HDS) #meta::でmetaのforest関数を使うってこと。forest関数はmetaにも、metaforにもあるため混乱を避ける。指定しないと直近使用になる。
ラベルや小数点以下の記載を修正する(適宜)
forest.dic <-
meta::forest(dic.HDS,
sortvar=TE, #効果純にソート あまりすすめられない。
xlim = c(0.1,8),
rightlabs = c("RR","95% CI","weight"),
leftlabs = c("Author", "Event","N","Event","N"),
lab.e = "Intervention",
pooled.totals = TRUE,
smlab = "",
text.random = "Random effect",
print.tau2 = FALSE,
col.diamond = "blue",
col.diamond.lines = "black",
print.I2.ci = TRUE,
digits = 2)
meta::forest(dic.HDS,
layout = "RevMan5")
meta::forest(dic.HDS,
layout = "JAMA")
#Funnel Plot フォレストプロットはmetaが綺麗だが、それ以外の機能はmetaforが使いやすい。
等高線強調Funnel plot
metafor::funnel(ds4,
level=c(90, 95, 99),
shade=c("white", "gray", "darkgray"),
refline=0)
legend("topright",
c("p > 0.1", "0.1 > p > 0.05", "0.05 > p > 0.01"),
fill=c("white", "gray", "darkgray"))
###左右対称性の検定
Beggの方法(ランクテスト)
metafor::ranktest(ds4)
##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.4222, p = 0.1083
###Eggerの方法(回帰)
metafor::regtest(ds4, model="lm") #左右非対称となる。閾値がひっかかる
##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 2.2715, df = 8, p = 0.0528
## Limit Estimate (as sei -> 0): b = -0.8342 (CI: -1.3992, -0.2691)