エビデンス計量評価 統合
chikamasa
2022-10-19
#エビデンス計量評価学 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#データセットの読み込み
Dataset <- read_csv("HDS_metaanalysis.csv")## Rows: 12 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Source, RoBsummary
## dbl (10): N1, Dic1, Mean1, Sd1, N2, Dic2, Mean2, Sd2, RoB_dic, High_dose
##
## ℹ 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.#連続量データの統合と可視化 ##metaforパッケージを使った統合と可視化 1.平均値データのある行に限定したデータの作成。
dmean <- Dataset %>%
filter(!is.na(Mean1)) #is.na(特定の列)関数を利用すれば欠測値を簡単に除去・抽出できる。2.必要な統計量(平均値差や標準誤差)を計算する。
res1 <-
escalc(measure="MD",
n1i=N1,
n2i=N2,
m1i=Mean1,
m2i=Mean2,
sd1i=Sd1,
sd2i=Sd2,
data = dmean,
append=TRUE)
res1 |> as.data.frame()## Source N1 Dic1 Mean1 Sd1 N2 Dic2 Mean2 Sd2 RoBsummary RoB_dic
## 1 Takeshima 2013 87 9 137.00 54.70 75 7 107.0 33.1 Some 0
## 2 Takahashi 2017 60 NA 175.33 36.47 61 NA 124.4 10.5 High 1
## 3 Nakayama 2007 27 9 123.90 32.70 36 10 89.1 16.1 High 1
## 4 Kawano 2007 26 5 228.60 59.40 35 4 88.5 14.5 High 1
## 5 Nitta 2017 56 17 282.00 32.00 59 32 92.0 8.0 Some 1
## 6 Yao 1990 25 12 88.20 37.80 24 11 79.2 23.2 Some 0
## 7 Yagome 2016 45 NA 150.20 35.50 45 NA 100.3 30.2 Some 0
## High_dose yi vi
## 1 0 30.00 48.99997
## 2 1 50.93 23.97506
## 3 0 34.80 46.80361
## 4 1 140.10 141.71330
## 5 1 190.00 19.37046
## 6 1 9.00 79.58027
## 7 1 49.90 48.27311- 上記で算出した統計量を統合する
統合する統計量を指定する 平均値差: “MD” リスク比: “RR” オッズ比: “OR” リスク差: “RD” 異質性（τ）を見積もる様々な方法: 良く使われるのはDLとREM L DL : DerSimonian-Laird REML : Restricted Maximum-Likelihood
res2<-rma(measure="MD", yi=yi, vi=vi,
data = res1,
method="DL") #fixed effect model: "FE"
res2##
## Random-Effects Model (k = 7; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 5963.2132 (SE = 3932.0797)
## tau (square root of estimated tau^2 value): 77.2218
## I^2 (total heterogeneity / total variability): 99.30%
## H^2 (total variability / sampling variability): 141.87
##
## Test for Heterogeneity:
## Q(df = 6) = 851.1951, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 72.0708 29.3291 2.4573 0.0140 14.5869 129.5547 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- metaforでフォレストプロット
metafor::forest(res2)
##metaパッケージをつかった統合と可視化 metaforパッケージにおいてForest plotをお手軽に綺麗にするのは難しいので、ここはあきらめてmetaパッケージを使ったほうが無難。 metaパッケージでは「2.必要な統計量（平均値差や標準誤差）」と「3.上記で算出した統計量を統合する」は同時に行う
mean.HDS <- metacont(n.e = N1,
mean.e = Mean1,
sd.e = Sd1,
n.c = N2,
mean.c = Mean2,
sd.c = Sd2,
data = dmean,
studlab = paste(Source),
fixed = FALSE,
random = TRUE,
method.tau = "DL",
sm = "MD")
mean.HDS## Number of studies combined: k = 7
## Number of observations: o = 661
##
## MD 95%-CI z p-value
## Random effects model 72.0708 [14.5869; 129.5547] 2.46 0.0140
##
## Quantifying heterogeneity:
## tau^2 = 5963.2132 [2015.1128; 25562.0042]; tau = 77.2218 [44.8900; 159.8812]
## I^2 = 99.3% [99.1%; 99.5%]; H = 11.91 [10.48; 13.54]
##
## Test of heterogeneity:
## Q d.f. p-value
## 851.20 6 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Jackson method for confidence interval of tau^2 and tau4. フォレストプロット作成
まずはシンプルな作図（これだけでも十分綺麗）
meta::forest(mean.HDS) #meta::でmetaのforest関数を使うってこと。forest関数はmetaにも、metaforにもあるため混乱を避ける。指定しないと直近使用になる。
ラベルや小数点以下の記載を修正する（適宜）
meta::forest(mean.HDS,
sortvar=TE, #arrangement by effect size エフェクトサイズにてソートする。
xlim = c(-100,300),
rightlabs = c("MD","95% CI","weight"),
leftlabs = c("Author", "N","Mean","SD","N","Mean","SD"),
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.mean=1, #小数点第一位まで記載
digits.sd = 1, #小数点第一位まで記載
digits = 1
)
既存レイアウトを使う方法
- RevMan（コクラン）方式
meta::forest(mean.HDS,
layout = "RevMan5")
- JAMA方式
meta::forest(mean.HDS,
layout = "JAMA")
#Funnel Plot フォレストプロットはmetaが綺麗だが、それ以外の機能はmetaforが使いやすい。
まずは普通のFunnel plot
metafor::funnel(res2)
等高線強調Funnel plot
metafor::funnel(res2, 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(res2)##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = -0.2381, p = 0.5619###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 = -1.1616, df = 5, p = 0.2978
## Limit Estimate (as sei -> 0): b = 193.7771 (CI: -51.0285, 438.5826)#サブグループ解析 subgroup analysis: High_dose (or not)
mean.HDS.sub <-
metacont(n.e = N1,
mean.e = Mean1,
sd.e = Sd1,
n.c = N2,
mean.c = Mean2,
sd.c = Sd2,
data = dmean,
studlab = paste(Source),
subgroup = High_dose , ###
subgroup.name = "High_dose", ###
fixed = FALSE,
random = TRUE,
method.tau = "DL",
sm = "MD")
mean.HDS.sub## Number of studies combined: k = 7
## Number of observations: o = 661
##
## MD 95%-CI z p-value
## Random effects model 72.0708 [14.5869; 129.5547] 2.46 0.0140
##
## Quantifying heterogeneity:
## tau^2 = 5963.2132 [2015.1128; 25562.0042]; tau = 77.2218 [44.8900; 159.8812]
## I^2 = 99.3% [99.1%; 99.5%]; H = 11.91 [10.48; 13.54]
##
## Test of heterogeneity:
## Q d.f. p-value
## 851.20 6 < 0.0001
##
## Results for subgroups (random effects model):
## k MD 95%-CI tau^2 tau Q I^2
## High_dose = 0 2 32.4550 [22.8655; 42.0445] 0 0 0.24 0.0%
## High_dose = 1 5 87.9779 [13.9580; 161.9978] 7069.0343 84.0775 682.32 99.4%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 2.13 1 0.1448
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Jackson method for confidence interval of tau^2 and tau###Forest plot
# pdf(filename = "forest_sub_dose.pdf", height = 6, width = 10)
meta::forest(mean.HDS.sub,
sortvar=TE,
xlim = c(-100,300),
rightlabs = c("MD","95% CI","weight"),
leftlabs = c("Author", "N","Mean","SD","N","Mean","SD"),
lab.e = "Intervention",
pooled.totals = TRUE,
smlab = "",
text.random = "Random effect",
print.tau2 = FALSE,
col.diamond = "blue",
col.diamond.lines = "black",
print.I2.ci = FALSE,
digits.mean=1,
digits.sd = 1,
digits = 2)
# dev.off()