エビデンス計量評価論講義:メタ解析演習
options(repos="https://cran.rstudio.com" )必要パッケージのインポート
pacman::p_load(tidyverse, metafor, meta, svglite)データのインポート
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.データの確認と、名前の変更
RStudio内のEmvironmentからデータの確認をする。
ヘッダー項目の意味
| Source | N1 | Dic1 | Mean1 | Sd1 | N2 | Dic2 | Mean2 | Sd2 | RoBsummary | RoB_dic | High_dose |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 著者 | 介入群 | 対照群 | RoB | RoB2値 | 高用量 | ||||||
| 症例数 | 2値 | 平均 | 標準偏差 | 症例数 | 2値 | 平均 | 標準偏差 | 3値 | High=1 | High=1 |
実データ
Dataset |>
head(5)データの統合と可視化（連続値）
metaforパッケージを使った統合と可視化
1. 平均値データのある行に限定したデータ作成
dmean <-
Dataset %>%
filter(!is.na(Mean1))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()3. 上記で算出した統計量を統合する
- 統合する統計量を指定する
- 平均値差: “MD”
- リスク比: “RR”
- オッズ比: “OR”
- リスク差: “RD”
- 平均値差: “MD”
- 異質性（τ）を見積もる様々な方法:
良く使われるのはDLとREML
- DL : DerSimonian-Laird
- REML : Restricted Maximum-Likelihood
- DL : DerSimonian-Laird
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 ' ' 14. 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")## Warning: Use argument 'common' instead of 'fixed' (deprecated).mean.HDS## Number of studies: k = 7
## Number of observations: o = 661 (o.e = 326, o.c = 335)
##
## MD 95%-CI z p-value
## Random effects model 72.0708 [14.5869; 129.5547] 2.46 0.0140
##
## Quantifying heterogeneity (with 95%-CIs):
## 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 of meta-analysis methods:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Jackson method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q4. フォレストプロット作成
まずはシンプルな作図（これだけでも十分綺麗）
meta::forest(mean.HDS)
ラベルや小数点以下の記載を修正する（適宜）
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
)## Warning: Use argument 'label.e' instead of 'lab.e' (deprecated).
既存レイアウトを使う方法
1. RevMan（コクラン）方式
meta::forest(mean.HDS,
layout = "RevMan5")
###### 2. JAMA方式
meta::forest(mean.HDS,
layout = "JAMA")
図の保存・保存した図の処理
pdf, png, svgなど好きな形式で保存できる （図の微修正はInkscape等のソフトで行
# png("forest01.png", width = 10, height = 3, units = "in", res = 600)
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"),
label.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)
# dev.off()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.5619Eggerの方法（回帰）
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)# metabias(mean.HDS, method.bias = "Egger") #少数のため計算不可サブグループ解析
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")## Warning: Use argument 'common' instead of 'fixed' (deprecated).mean.HDS.sub## Number of studies: k = 7
## Number of observations: o = 661 (o.e = 326, o.c = 335)
##
## MD 95%-CI z p-value
## Random effects model 72.0708 [14.5869; 129.5547] 2.46 0.0140
##
## Quantifying heterogeneity (with 95%-CIs):
## 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 of meta-analysis methods:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Jackson method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on QForest 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)## Warning: Use argument 'label.e' instead of 'lab.e' (deprecated).
データの統合と可視化（２値）
metaforパッケージを使った統合と可視化
0. （項目名の変更は適宜）
1. 2値データのある行に限定したデータ作成
Dataset |> head(6)Dataset <-
Dataset %>%
mutate(pos.i = Dic1,
pos.c = Dic2,
total.i = N1,
total.c = N2,
neg.i = N1 - Dic1,
neg.c = N2 - Dic2)
ddic <- Dataset %>%
filter(!is.na(Dic1))2. 統計量算出
陽性症例数と陰性症例数を使う場合はai, bi, ci, diを使う。 陽性症例数と総症例数の場合は、ai, n1i, ci, n2iを使う。
# res3 <-
# escalc(measure="RR",
# ai=pos.i,
# bi=neg.i,
# ci=pos.c,
# di=neg.c,
# data = ddic,
# append=TRUE)
res3 <-
escalc(measure = "RR",
ai = pos.i,
n1i = total.i,
ci = pos.c,
n2i = total.c,
data = ddic,
append = T)
res3 |> as.data.frame()3 and 4. 統合と可視化
res4 <- metafor::rma(measure="RR", yi=yi, vi=vi,
data = res3,
method="DL")
res4##
## 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 ' ' 1metafor::forest(res4, atransf = exp, digits=2) #digits=2とすると小数点以下第二位まで
metaパッケージを使った統合と可視化
dic.HDS <- meta::metabin(event.e = pos.i,
n.e = total.i,
event.c = pos.c,
n.c = total.c,
data = ddic,
studlab = paste(Source),
fixed = FALSE,
random = TRUE,
method.tau = "DL",
sm = "RR")## Warning: Use argument 'common' instead of 'fixed' (deprecated).dic.HDS## Number of studies: k = 10
## Number of observations: o = 919 (o.e = 453, o.c = 466)
## 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 (with 95%-CIs):
## 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 of meta-analysis methods:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Mantel-Haenszel estimator used in calculation of Q and tau^2 (like RevMan 5)
## - Calculation of I^2 based on Qmeta::forest(dic.HDS)
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)## Warning: Use argument 'label.e' instead of 'lab.e' (deprecated).
Funnel plot
metafor::funnel(res4,
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(res4)##
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.4222, p = 0.1083Eggerの方法（回帰）
metafor::regtest(res4, 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)# meta::metabias(dic.HDS, method.bias = "Egger") #同じ結果になるサブグループ解析
subgroup analysis: RoB_dicを用いて
dic.HDS.sub <-
meta::metabin(event.e = pos.i,
n.e = total.i,
event.c = pos.c,
n.c = total.c,
data = ddic,
studlab = paste(Source),
subgroup = RoB_dic , ###
subgroup.name = "RoB_dic", ###
fixed = FALSE,
random = TRUE,
method.tau = "DL",
sm = "RR")## Warning: Use argument 'common' instead of 'fixed' (deprecated).dic.HDS.sub## Number of studies: k = 10
## Number of observations: o = 919 (o.e = 453, o.c = 466)
## 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 (with 95%-CIs):
## 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
##
## Results for subgroups (random effects model):
## k RR 95%-CI tau^2 tau Q I^2
## RoB_dic = 0 5 0.8401 [0.6027; 1.1709] 0.0145 0.1203 4.42 9.4%
## RoB_dic = 1 5 0.7195 [0.4832; 1.0713] 0.0738 0.2716 6.57 39.1%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 0.34 1 0.5579
##
## Details of meta-analysis methods:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Mantel-Haenszel estimator used in calculation of Q and tau^2 (like RevMan 5)
## - Calculation of I^2 based on QForest plot
# pdf("forest_sub_rob.pdf", height = 6, width = 10)
# svglite::svglite("forest_sub_rob.svg")
meta::forest(dic.HDS.sub,
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",
test.subgroup.random = T,
print.I2.ci = FALSE,
digits = 2)## Warning: Use argument 'label.e' instead of 'lab.e' (deprecated).