BMJ Booksの“Systematic Reviews in Health Care"の中にある
"Chapter 18: Meta-analysis in Stata"をRを用いて再現。
http://www.blackwellpublishing.com/medicine/bmj/systreviews/samplechap.asp
http://www.blackwellpublishing.com/medicine/bmj/systreviews/pdfs/chapter18.pdf
使い方
コードをコピーして手元のRのエディタに貼付ける。
Windows版では左上 [ファイル] → [新しいスクリプト]でエディタ。
Mac版では真っ白い紙の形のアイコンでエディタ。
実行する行、あるいは範囲を選択して
WindowsであればCtrl＋R、MacであればCommand＋Returnで実行。
## 存在しないパッケージのみインストールするための関数を定義
inst <- function(PKG) {
if(!(PKG %in% rownames(installed.packages()))) {
install.packages(pkgs=PKG, dependencies=TRUE)
}
}
## Stataのデータを読むための library(foreign)を使用。
library(foreign)
## Webから直にデータを読み込みます。
strepto <- read.dta("http://www.blackwellpublishing.com/medicine/bmj/systreviews/datasets/strepto.dta")
magnes <- read.dta("http://www.blackwellpublishing.com/medicine/bmj/systreviews/datasets/magnes.dta")
bcgtrial <- read.dta("http://www.blackwellpublishing.com/medicine/bmj/systreviews/datasets/bcgtrial.dta")
## rmeta, metafor, metaパッケージを使用するので、
## なければインストールして読み込みます。
## いろいろありますがmetaが一番気が利いてるかも。
inst("rmeta") ; library(rmeta)
inst("metafor"); library(metafor)
inst("meta"); library(meta)
## MIに対するStreptokinase使用/不使用と死亡数のデータ。p 349のものです。
strepto
trial trialnam year pop1 deaths1 pop0 deaths0
1 1 Fletcher 1959 12 1 11 4
2 2 Dewar 1963 21 4 21 7
3 3 1st European 1969 83 20 84 15
4 4 Heikinheimo 1971 219 22 207 17
5 5 Italian 1971 164 19 157 18
6 6 2nd European 1971 373 69 357 94
7 7 2nd Frankfurt 1973 102 13 104 29
8 8 1st Australian 1973 264 26 253 32
9 9 NHLBI SMIT 1974 53 7 54 3
10 10 Valere 1975 49 11 42 9
11 11 Frank 1975 55 6 53 6
12 12 UK Collab 1976 302 48 293 52
13 13 Klein 1976 14 4 9 1
14 14 Austrian 1977 352 37 376 65
15 15 Lasierra 1977 13 1 11 3
16 16 N German 1977 249 63 234 51
17 17 Witchitz 1977 32 5 26 5
18 18 2nd Australian 1977 112 25 118 31
19 19 3rd European 1977 156 25 159 50
20 20 ISAM 1986 859 54 882 63
21 21 GISSI-1 1986 5860 628 5852 758
22 22 ISIS-2 1988 8592 791 8595 1029
## データ構造の表示はSTRucture()を使います。
str(strepto)
'data.frame': 22 obs. of 7 variables:
$ trial : int 1 2 3 4 5 6 7 8 9 10 ...
$ trialnam: chr "Fletcher" "Dewar" "1st European" "Heikinheimo" ...
$ year : int 1959 1963 1969 1971 1971 1971 1973 1973 1974 1975 ...
$ pop1 : int 12 21 83 219 164 373 102 264 53 49 ...
$ deaths1 : int 1 4 20 22 19 69 13 26 7 11 ...
$ pop0 : int 11 21 84 207 157 357 104 253 54 42 ...
$ deaths0 : int 4 7 15 17 18 94 29 32 3 9 ...
- attr(*, "datalabel")= chr "Streptokinase after MI"
- attr(*, "time.stamp")= chr "15 Feb 2001 10:07"
- attr(*, "formats")= chr "%8.0g" "%14s" "%8.0g" "%12.0g" ...
- attr(*, "types")= int 98 141 105 105 105 105 105
- attr(*, "val.labels")= chr "" "" "" "" ...
- attr(*, "var.labels")= chr "Trial number" "Trial name" "Year of publication" "Treated population" ...
- attr(*, "version")= int 6
## rmetaパッケージのmeta.Mantel-Haenszelを使用します。以下は1行でもかけます。
## p 351の例にあるStataの解析を再現。
## trial名にスペースがあるとうまくいかないので、sub()で" "を"_"に入れ替えておく。
## パケージ名::コマンド名 の表記にしていますがコマンド名だけでもok。
strepto_rmeta <-
rmeta::meta.MH(
data=strepto, # streptoデータを使用。
ntrt=pop1, nctrl=pop0, # Number in TReaTmentとNumber in ConTRoLを指定します。
ptrt=deaths1, pctrl=deaths0, # TReaTment群とConTRoL群のイベントの数を指定します。
names=sub(pattern=" ", replacement="_", strepto$trialnam),
statistic="RR" # 本にあわせてRelative Riskを使います。
)
## summary()すると結果を見ることができます。
## 有効数字が違いますが、同じ結果が再現できているかと思います。
summary(strepto_rmeta)
Fixed effects ( Mantel-Haenszel ) meta-analysis
Call: rmeta::meta.MH(ntrt = pop1, nctrl = pop0, ptrt = deaths1, pctrl = deaths0,
names = sub(pattern = " ", replacement = "_", strepto$trialnam),
data = strepto, statistic = "RR")
------------------------------------
RR (lower 95% upper)
Fletcher 0.23 0.03 1.75
Dewar 0.57 0.20 1.66
1st_European 1.35 0.74 2.45
Heikinheimo 1.22 0.67 2.24
Italian 1.01 0.55 1.85
2nd_European 0.70 0.53 0.92
2nd_Frankfurt 0.46 0.25 0.83
1st_Australian 0.78 0.48 1.27
NHLBI_SMIT 2.38 0.65 8.71
Valere 1.05 0.48 2.28
Frank 0.96 0.33 2.80
UK_Collab 0.90 0.63 1.28
Klein 2.57 0.34 19.48
Austrian 0.61 0.42 0.89
Lasierra 0.28 0.03 2.34
N_German 1.16 0.84 1.60
Witchitz 0.81 0.26 2.51
2nd_Australian 0.85 0.54 1.34
3rd_European 0.51 0.33 0.78
ISAM 0.88 0.62 1.25
GISSI-1 0.83 0.75 0.91
ISIS-2 0.77 0.70 0.84
------------------------------------
Mantel-Haenszel RR =0.8 95% CI ( 0.75,0.85 )
Test for heterogeneity: X^2( 21 ) = 30.41 ( p-value 0.084 )
## おなじものをEZRでも使用されているmetaパッケージで施行。
strepto_meta <-
meta::metabin(data=strepto,
event.e=deaths1, n.e=pop1, # exposure群のeventとn
event.c=deaths0, n.c=pop0, # control群のeventとn
sm="RR", studlab=trialnam, # summary measureとstudy label
comb.fixed=TRUE, comb.random=FALSE)
## summaryは不要でそのまま表示。
strepto_meta
RR 95%-CI %W(fixed)
Fletcher 0.2292 [0.0300; 1.7499] 0.18
Dewar 0.5714 [0.1962; 1.6647] 0.30
1st European 1.3494 [0.7429; 2.4509] 0.64
Heikinheimo 1.2232 [0.6688; 2.2371] 0.75
Italian 1.0105 [0.5510; 1.8531] 0.78
2nd European 0.7026 [0.5338; 0.9247] 4.10
2nd Frankfurt 0.4571 [0.2522; 0.8282] 1.22
1st Australian 0.7786 [0.4780; 1.2683] 1.39
NHLBI SMIT 2.3774 [0.6490; 8.7086] 0.13
Valere 1.0476 [0.4809; 2.2821] 0.41
Frank 0.9636 [0.3316; 2.8005] 0.26
UK Collab 0.8956 [0.6261; 1.2809] 2.25
Klein 2.5714 [0.3394; 19.4813] 0.05
Austrian 0.6080 [0.4173; 0.8861] 2.68
Lasierra 0.2821 [0.0340; 2.3403] 0.14
N German 1.1609 [0.8403; 1.6038] 2.24
Witchitz 0.8125 [0.2634; 2.5062] 0.24
2nd Australian 0.8497 [0.5369; 1.3446] 1.29
3rd European 0.5096 [0.3327; 0.7805] 2.11
ISAM 0.8801 [0.6195; 1.2503] 2.65
GISSI-1 0.8274 [0.7491; 0.9138] 32.34
ISIS-2 0.7690 [0.7044; 0.8395] 43.86
Number of studies combined: k=22
RR 95%-CI z p.value
Fixed effect model 0.7988 [0.7546; 0.8455] -7.747 < 0.0001
Quantifying heterogeneity:
tau^2 = 0.0112; H = 1.2 [1; 1.56]; I^2 = 30.9% [0%; 58.9%]
Test of heterogeneity:
Q d.f. p.value
30.41 21 0.084
Details on meta-analytical method:
- Mantel-Haenszel method
## metaパッケージでは、p 352のようなforest plotもいきなりきれいに描いてくれます。
meta::forest.meta(strepto_meta)
### Funnel plot
meta::funnel(strepto_meta)
### Cumulative meta-analysis
strepto_meta.cuml <- meta::metacum(strepto_meta)
strepto_meta.cuml
Cumulative meta-analysis (Fixed effect model)
RR 95%-CI p.value tau^2 I^2
Adding Fletcher (k=1) 0.2292 [0.0300; 1.7499] 0.1555 0
Adding Dewar (k=2) 0.4436 [0.1740; 1.1309] 0.0887 0 0.0%
Adding 1st European (k=3) 0.9614 [0.5905; 1.5651] 0.8741 0.3145 50.7%
Adding Heikinheimo (k=4) 1.0664 [0.7296; 1.5587] 0.7398 0.08 30.3%
Adding Italian (k=5) 1.0498 [0.7610; 1.4483] 0.7671 0.0128 7.8%
Adding 2nd European (k=6) 0.8387 [0.6812; 1.0327] 0.0976 0.0557 37.2%
Adding 2nd Frankfurt (k=7) 0.7800 [0.6415; 0.9485] 0.0128 0.082 47.6%
Adding 1st Australian (k=8) 0.7798 [0.6503; 0.9351] 0.0073 0.0525 38.9%
Adding NHLBI SMIT (k=9) 0.8012 [0.6697; 0.9585] 0.0154 0.0711 43.8%
Adding Valere (k=10) 0.8115 [0.6814; 0.9663] 0.019 0.0597 38.7%
Adding Frank (k=11) 0.8154 [0.6863; 0.9687] 0.0203 0.0482 32.3%
Adding UK Collab (k=12) 0.8299 [0.7106; 0.9693] 0.0186 0.0312 26.7%
Adding Klein (k=13) 0.8372 [0.7172; 0.9772] 0.0243 0.0325 26.0%
Adding Austrian (k=14) 0.7966 [0.6906; 0.9189] 0.0018 0.0351 29.7%
Adding Lasierra (k=15) 0.7919 [0.6868; 0.9132] 0.0013 0.0344 27.9%
Adding N German (k=16) 0.8392 [0.7368; 0.9558] 0.0082 0.0462 37.3%
Adding Witchitz (k=17) 0.8388 [0.7371; 0.9545] 0.0077 0.0404 33.1%
Adding 2nd Australian (k=18) 0.8395 [0.7413; 0.9508] 0.0059 0.0322 28.9%
Adding 3rd European (k=19) 0.8066 [0.7160; 0.9087] 0.0004 0.0459 37.6%
Adding ISAM (k=20) 0.8148 [0.7278; 0.9122] 0.0004 0.0376 34.6%
Adding GISSI-1 (k=21) 0.8220 [0.7629; 0.8857] < 0.0001 0.0201 31.2%
Adding ISIS-2 (k=22) 0.7988 [0.7546; 0.8455] < 0.0001 0.0112 30.9%
Pooled estimate 0.7988 [0.7546; 0.8455] < 0.0001 0.0112 30.9%
Details on meta-analytical method:
- Mantel-Haenszel method
### Cumulative forest plot
meta::forest(strepto_meta.cuml, xlim = c(0.5, 2))
### Funnel plotの検定?
meta::metabias(strepto_meta)
Linear regression test of funnel plot asymmetry
data: strepto_meta
t = 0.4681, df = 20, p-value = 0.6448
alternative hypothesis: asymmetry in funnel plot
sample estimates:
bias se.bias slope
0.1670 0.3567 -0.2389
## Smoking cessationデータでのDerSimonian-Laird random effect。
## A Handbook of Statistical Analyses Using R — 2nd Edより
## http://cran.r-project.org/web/packages/HSAUR2/vignettes/Ch_meta_analysis.pdf
inst("HSAUR2"); library(HSAUR2)
data(smoking)
## Quit Treated, Total Treated, Quit Control, Total Controlの略です。
smoking
qt tt qc tc
Blondal89 37 92 24 90
Campbell91 21 107 21 105
Fagerstrom82 30 50 23 50
Fee82 23 180 15 172
Garcia89 21 68 5 38
Garvey00 75 405 17 203
Gross95 37 131 6 46
Hall85 18 41 10 36
Hall87 30 71 14 68
Hall96 24 98 28 103
Hjalmarson84 31 106 16 100
Huber88 31 54 11 60
Jarvis82 22 58 9 58
Jensen91 90 211 28 82
Killen84 16 44 6 20
Killen90 129 600 112 617
Malcolm80 6 73 3 121
McGovern92 51 146 40 127
Nakamura90 13 30 5 30
Niaura94 5 84 4 89
Pirie92 75 206 50 211
Puska79 29 116 21 113
Schneider85 9 30 6 30
Tonnesen88 23 60 12 53
Villa99 11 21 10 26
Zelman92 23 58 18 58
## metaパッケージでやってみます。
smoking_meta <-
meta::metabin(data=smoking,
event.e=qt, n.e=tt,
event.c=qc, n.c=tc,
sm="OR", studlab=rownames(smoking))
## そのまま使うとMantel-Haenszel (Fixed effect)とDerSimonian-Laird (Random effect)両方施行。
## それぞれHSAUR2のPDFのp 4とp 6の結果に一致していると思います。
smoking_meta
OR 95%-CI %W(fixed) %W(random)
Blondal89 1.8500 [0.9892; 3.460] 3.96 4.82
Campbell91 0.9767 [0.4971; 1.919] 4.66 4.33
Fagerstrom82 1.7609 [0.7965; 3.893] 2.51 3.41
Fee82 1.5333 [0.7713; 3.048] 3.66 4.23
Garcia89 2.9489 [1.0094; 8.615] 1.21 2.08
Garvey00 2.4866 [1.4256; 4.337] 5.04 5.63
Gross95 2.6241 [1.0265; 6.708] 1.74 2.60
Hall85 2.0348 [0.7829; 5.289] 1.63 2.53
Hall87 2.8223 [1.3289; 5.994] 2.26 3.69
Hall96 0.8687 [0.4614; 1.636] 5.63 4.75
Hjalmarson84 2.1700 [1.1005; 4.279] 3.18 4.30
Huber88 6.0040 [2.5721; 14.015] 1.21 3.07
Jarvis82 3.3272 [1.3706; 8.077] 1.53 2.86
Jensen91 1.4345 [0.8429; 2.441] 6.32 5.95
Killen84 1.3333 [0.4279; 4.155] 1.43 1.88
Killen90 1.2349 [0.9310; 1.638] 23.69 10.55
Malcolm80 3.5224 [0.8531; 14.543] 0.57 1.26
McGovern92 1.1676 [0.7040; 1.937] 7.61 6.31
Nakamura90 3.8235 [1.1500; 12.713] 0.77 1.70
Niaura94 1.3449 [0.3487; 5.188] 1.00 1.38
Pirie92 1.8435 [1.2044; 2.822] 8.59 7.61
Puska79 1.4603 [0.7750; 2.752] 4.36 4.74
Schneider85 1.7143 [0.5228; 5.621] 1.15 1.74
Tonnesen88 2.1239 [0.9285; 4.858] 2.15 3.19
Villa99 1.7600 [0.5489; 5.643] 1.16 1.80
Zelman92 1.4603 [0.6791; 3.140] 2.97 3.60
Number of studies combined: k=26
OR 95%-CI z p.value
Fixed effect model 1.670 [1.468; 1.899] 7.819 < 0.0001
Random effects model 1.752 [1.483; 2.069] 6.604 < 0.0001
Quantifying heterogeneity:
tau^2 = 0.0475; H = 1.18 [1; 1.5]; I^2 = 28.4% [0%; 55.8%]
Test of heterogeneity:
Q d.f. p.value
34.9 25 0.09
Details on meta-analytical method:
- Mantel-Haenszel method
meta::forest.meta(smoking_meta)
meta::funnel(smoking_meta)
meta::metabias(smoking_meta)
Linear regression test of funnel plot asymmetry (efficient score)
data: smoking_meta
t = 2.326, df = 24, p-value = 0.02878
alternative hypothesis: asymmetry in funnel plot
sample estimates:
bias se.bias slope
1.2912 0.5551 0.1127
smoking_meta.cuml <- meta::metacum(smoking_meta)
smoking_meta.cuml
Cumulative meta-analysis (Fixed effect model)
OR 95%-CI p.value tau^2 I^2
Adding Blondal89 (k=1) 1.850 [0.9892; 3.460] 0.0541 0
Adding Campbell91 (k=2) 1.378 [0.8736; 2.175] 0.1678 0.0936 45.9%
Adding Fagerstrom82 (k=3) 1.465 [0.9866; 2.175] 0.0583 0.0078 5.8%
Adding Fee82 (k=4) 1.482 [1.0520; 2.087] 0.0244 0 0.0%
Adding Garcia89 (k=5) 1.593 [1.1515; 2.203] 0.0049 0 0.0%
Adding Garvey00 (k=6) 1.807 [1.3682; 2.386] < 0.0001 0.0123 8.8%
Adding Gross95 (k=7) 1.869 [1.4325; 2.440] < 0.0001 0.0022 1.6%
Adding Hall85 (k=8) 1.881 [1.4550; 2.430] < 0.0001 0 0.0%
Adding Hall87 (k=9) 1.960 [1.5376; 2.499] < 0.0001 0 0.0%
Adding Hall96 (k=10) 1.770 [1.4134; 2.216] < 0.0001 0.0549 28.3%
Adding Hjalmarson84 (k=11) 1.806 [1.4586; 2.236] < 0.0001 0.0399 22.6%
Adding Huber88 (k=12) 1.944 [1.5820; 2.390] < 0.0001 0.1186 45.9%
Adding Jarvis82 (k=13) 2.000 [1.6359; 2.444] < 0.0001 0.1162 44.8%
Adding Jensen91 (k=14) 1.919 [1.5910; 2.316] < 0.0001 0.103 43.3%
Adding Killen84 (k=15) 1.901 [1.5798; 2.288] < 0.0001 0.0933 39.9%
Adding Killen90 (k=16) 1.675 [1.4354; 1.954] < 0.0001 0.1045 48.6%
Adding Malcolm80 (k=17) 1.690 [1.4496; 1.969] < 0.0001 0.1034 47.1%
Adding McGovern92 (k=18) 1.639 [1.4153; 1.897] < 0.0001 0.0978 46.9%
Adding Nakamura90 (k=19) 1.660 [1.4355; 1.920] < 0.0001 0.102 47.0%
Adding Niaura94 (k=20) 1.656 [1.4333; 1.914] < 0.0001 0.0947 44.1%
Adding Pirie92 (k=21) 1.674 [1.4602; 1.920] < 0.0001 0.0796 41.7%
Adding Puska79 (k=22) 1.664 [1.4560; 1.902] < 0.0001 0.0708 39.0%
Adding Schneider85 (k=23) 1.665 [1.4577; 1.901] < 0.0001 0.0646 36.1%
Adding Tonnesen88 (k=24) 1.675 [1.4692; 1.910] < 0.0001 0.0594 33.9%
Adding Villa99 (k=25) 1.676 [1.4713; 1.910] < 0.0001 0.0536 31.0%
Adding Zelman92 (k=26) 1.670 [1.4684; 1.899] < 0.0001 0.0475 28.4%
Pooled estimate 1.670 [1.4684; 1.899] < 0.0001 0.0475 28.4%
Details on meta-analytical method:
- Mantel-Haenszel method
meta::forest(smoking_meta.cuml)
