#エビデンス計量評価学  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.

2値データの統合

metaforパッケージを使った統合と可視化

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
  1. 上記で算出した統計量を統合する

統合する統計量を指定する   平均値差: “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
  1. metaforでフォレストプロット
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)

4. フォレストプロット作成

まずはシンプルな作図(これだけでも十分綺麗)

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)

既存レイアウトを使う方法

  1. RevMan(コクラン)方式
meta::forest(dic.HDS,
       layout = "RevMan5")

  1. JAMA方式
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)