HDS_nakamura

0-1. パッケージの起動

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ 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

library(svglite)
library(gridExtra)

## 
##  次のパッケージを付け加えます: 'gridExtra' 
## 
##  以下のオブジェクトは 'package:dplyr' からマスクされています:
## 
##     combine

library(readr)

0-2.データのインポート

0-3. リスク比やリスク差の追加

ds2<-ds1 %>% mutate(upper_right=(N1-Event1), lower_right=(N2-Event2),
                    incidence1=(Event1/N1), incidence2=(Event2/N2))
ds3<-ds2 %>% mutate(Risk_Ratio=(incidence1/incidence2), Risk_Difference=incidence1-incidence2)
ds3

##               Source N1 Event1 N2 Event2   N upper_right lower_right incidence1
## 1          Cook 1990 25     12 24     11  49          13          13 0.48000000
## 2  Conan McCaul 2003 35      6 36      4  71          29          32 0.17142857
## 3          Shin 2007 27      9 36     10  63          18          26 0.33333333
## 4          Shin 2007 26      5 35      4  61          21          31 0.19230769
## 5   Dabu-Bondoc 2012 30     14 32     20  62          16          12 0.46666667
## 6       Patel P 2013 87      9 75      7 162          78          68 0.10344828
## 7          Jain 2013 75     23 75     43 150          52          32 0.30666667
## 8        Pin On 2016 42      2 44      2  86          40          42 0.04761905
## 9        Mishra 2017 50      6 50     14 100          44          36 0.12000000
## 10          Rao 2017 56     17 59     32 115          39          27 0.30357143
##    incidence2 Risk_Ratio Risk_Difference
## 1  0.45833333  1.0472727     0.021666667
## 2  0.11111111  1.5428571     0.060317460
## 3  0.27777778  1.2000000     0.055555556
## 4  0.11428571  1.6826923     0.078021978
## 5  0.62500000  0.7466667    -0.158333333
## 6  0.09333333  1.1083744     0.010114943
## 7  0.57333333  0.5348837    -0.266666667
## 8  0.04545455  1.0476190     0.002164502
## 9  0.28000000  0.4285714    -0.160000000
## 10 0.54237288  0.5597098    -0.238801453

0-4.リスク比(RR)の効果量の計算

res1<- escalc(measure="RR",
              n1i=N1, n2i=N2,
              ai=Event1, bi=upper_right, ci=Event2, di=lower_right,
              data=ds2)
res1 |> as.data.frame()

##               Source N1 Event1 N2 Event2   N upper_right lower_right incidence1
## 1          Cook 1990 25     12 24     11  49          13          13 0.48000000
## 2  Conan McCaul 2003 35      6 36      4  71          29          32 0.17142857
## 3          Shin 2007 27      9 36     10  63          18          26 0.33333333
## 4          Shin 2007 26      5 35      4  61          21          31 0.19230769
## 5   Dabu-Bondoc 2012 30     14 32     20  62          16          12 0.46666667
## 6       Patel P 2013 87      9 75      7 162          78          68 0.10344828
## 7          Jain 2013 75     23 75     43 150          52          32 0.30666667
## 8        Pin On 2016 42      2 44      2  86          40          42 0.04761905
## 9        Mishra 2017 50      6 50     14 100          44          36 0.12000000
## 10          Rao 2017 56     17 59     32 115          39          27 0.30357143
##    incidence2          yi         vi
## 1  0.45833333  0.04618938 0.09257576
## 2  0.11111111  0.43363599 0.36031746
## 3  0.27777778  0.18232156 0.14629630
## 4  0.11428571  0.52039507 0.38296703
## 5  0.62500000 -0.29213642 0.05684524
## 6  0.09333333  0.10289442 0.22914067
## 7  0.57333333 -0.62570590 0.04006741
## 8  0.04545455  0.04652002 0.95346320
## 9  0.28000000 -0.84729786 0.19809524
## 10 0.54237288 -0.58033681 0.05526723

1-1. 統計量の統合

res2 <- rma(yi, vi, data=res1, method="DL")
res2

## 
## 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-2. metaforでのフォレストプロットの作成

metafor::forest(res2)

1-3. metaでのフォレストプロットの作成

metacont: ２群の連続変数

metabin: ２群の２値変数

metainc: ２群の発症 (incidence) 率

metacor: 相関係数 (１群)

metamean: 平均値 (１群)

metaprop: proportion (１群) 引数 event(イベント数) と n(観察数)

metarate: 発症 (incidence) 率 (１群) 引数 event(イベント数)と person time

RR<-metabin(event.e=Event1, n.e=N1, event.c=Event2, n.c=N2, 
        data=ds3, method.tau = "DL",sm="RR")
RR

## Number of studies combined: k = 10
## Number of observations: o = 919
## Number of events: e = 250
## 
##                          RR           95%-CI     z p-value
## Common effect model  0.7294 [0.5958; 0.8931] -3.06  0.0022
## 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(RR)

1-4. 既存のレイアウトでのフォレストプロット(RevMan, JAMA)

meta::forest(RR,layout="RevMan5")

meta::forest(RR, layout="JAMA")

2-1. ファンネルプロット

metafor::funnel(RR)
meta::funnel(RR)

2-2. 等高線強調ファンネルプロット

metafor::funnel(res2, 
                level=c(90, 95, 99), 
                shade=c("white", "gray55", "gray75"), 
                refline=0)

legend("topright",
       c("p > 0.1", "0.1 > p > 0.05", "0.05 > p > 0.01"),
       fill=c("white", "gray55", "gray75"))

2-3. 左右対称性の検定

Beggの方法

metafor::ranktest(res2)

## 
## Rank Correlation Test for Funnel Plot Asymmetry
## 
## Kendall's tau = 0.4222, p = 0.1083

Eggerの方法（回帰)

metafor::regtest(res2, 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)