-
Notifications
You must be signed in to change notification settings - Fork 1
/
magnesium-meta.R
211 lines (173 loc) · 6.21 KB
/
magnesium-meta.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
# project 2
rm(list=ls())
##### set up data ##############################################################
setwd("C:/Users/hw1220/Desktop/project")
mag <- read.csv("magnesium.csv")
# deleting unnecessary column
mag$X <- NULL
# rearrange to fit table 2 in Higgins paper
colnames(mag) <- c("trial", "name", "year", "nm", "rm", "nc", "rc")
mag <- mag[c("trial", "name", "year", "rm", "nm", "rc", "nc")]
mag[13, 1] <- 8
mag[8:12, 1] <- c(9:13)
mag <- mag[order(mag$trial),]
# examine the studies
mag$total <- mag$nm + mag$nc
sum(mag[1:14, 8]) # 4388 patients in other 14 studies
58050/4388 # how much bigger ISIS-4 is than the others
##### frequentist analysis #####################################################
# install.packages("metafor")
# install.packages("rmeta")
# install.packages("formattable")
library(metafor)
library(rmeta)
library(formattable)
### reproduce table 2
# OR in Peto FE for first 8 trials
or.p8 <- rma.peto(ai=rm, n1i=nm, ci=rc, n2i=nc,
data=mag[1:8,], slab=name,
add=1/2, to="only0", drop00=TRUE,
level=95, digits=2, verbose=FALSE)
summary(or.p8) # OR 0.65, CI(0.51, 0.82)
# OR in Peto FE for first 14 trials
or.p14 <- rma.peto(ai=rm, n1i=nm, ci=rc, n2i=nc,
data=mag[1:14,], slab=name,
add=1/2, to="only0", drop00=TRUE,
level=95, digits=2, verbose=FALSE)
summary(or.p14) # OR 0.57, CI(0.46, 0.71)
# OR in Peto FE for all trials incl ISIS-4
or.p15 <- rma.peto(ai=rm, n1i=nm, ci=rc, n2i=nc,
data=mag[1:15,], slab=name,
add=1/2, to="only0", drop00=TRUE,
level=95, digits=2, verbose=FALSE)
summary(or.p15) # OR 1.01, CI(0.95, 1.07)
# OR in DSL RE for first 8 trials
or.dl8 <- meta.DSL(ntrt=nm, nctrl=nc, ptrt=rm, pctrl=rc,
conf.level=0.95, names=name,
data=mag[1:8,], na.action=na.fail,statistic="OR")
summary(or.dl8) # OR 0.55, CI(0.34, 0.89)
# OR in DSL RE for first 14 trials
or.dl14 <- meta.DSL(ntrt=nm, nctrl=nc, ptrt=rm, pctrl=rc,
conf.level=0.95, names=name,
data=mag[1:14,], na.action=na.fail,statistic="OR")
summary(or.dl14) # OR 0.47, CI(0.32, 0.68)
# OR in DSL RE for all trials incl ISIS-4
or.dl15 <- meta.DSL(ntrt=nm, nctrl=nc, ptrt=rm, pctrl=rc,
conf.level=0.95, names=name,
data=mag[1:15,], na.action=na.fail,statistic="OR")
summary(or.dl15) # OR 0.53, CI(0.36, 0.77)
# forest plots for all 3 FE models
forest(or.p8, main="Odds Ratio for the First 8 Trials")
forest(or.p14, main="Odds Ratio for the First 14 Trials")
forest(or.p15, main="Odds Ratio for All Trials")
##### bayesian #################################################################
# install.packages("rstan", repos = "http://cran.rstudio.com", dependencies = TRUE)
library(rstan)
k <- length(mag$trial)
nc <- mag$nc
nm <- mag$nm
rc <- mag$rc
rm <- mag$rm
# Bayesian model with reference prior
model_ref <- "
data{
int<lower=0> k;
int<lower=0> nc[k];
int<lower=0> nm[k];
int<lower=0> rc[k];
int<lower=0> rm[k];
}
parameters {
vector[k] delta;
real<lower=0, upper=1> pc[k];
real<lower=0> sigma;
real deltanew;
real mu;
}
transformed parameters{
real<lower=0, upper=1> pm[k];
for (i in 1:k) {
pm[i] <- exp(log(pc[i]/(1-pc[i]))+delta[i])/(1+exp(log(pc[i]/(1-pc[i]))+delta[i]));
}
}
model{
for (i in 1:k) {
rc[i] ~ binomial(nc[i], pc[i]);
rm[i] ~ binomial(nm[i], pm[i]);
delta[i] ~ normal(mu, sigma);
pc[i] ~ uniform(0, 1);
}
deltanew ~ normal(mu, sigma);
mu ~ normal(0, 100);
sigma ~ uniform(0, 100);
}"
# Bayesian model with skeptical prior
model_skp <- "
data{
int<lower=0> k;
int<lower=0> nc[k];
int<lower=0> nm[k];
int<lower=0> rc[k];
int<lower=0> rm[k];
}
parameters {
vector[k] delta;
real<lower=0, upper=1> pc[k];
real<lower=0> sigma;
real deltanew;
real mu;
}
transformed parameters{
real<lower=0, upper=1> pm[k];
for (i in 1:k) {
pm[i] <- exp(log(pc[i]/(1-pc[i]))+delta[i])/(1+exp(log(pc[i]/(1-pc[i]))+delta[i]));
}
}
model{
for (i in 1:k) {
rc[i] ~ binomial(nc[i], pc[i]);
rm[i] ~ binomial(nm[i], pm[i]);
delta[i] ~ normal(mu, sigma);
pc[i] ~ uniform(0, 1);
}
deltanew ~ normal(mu, sigma);
mu ~ normal(0, 1/sqrt(32.69));
sigma ~ uniform(0, 100);
}"
# fit reference model in stan
mod_fit_ref <- stan(model_code=model_ref, data=c("k", "nc", "nm", "rc", "rm"),
pars=c("delta", "mu","sigma","deltanew"),
iter=500000, chains=3, warmup=500, verbose=FALSE)
print(mod_fit_ref)
# fit skeptical model in stan
mod_fit_skp <- stan(model_code=model_skp, data=c("k", "nc", "nm", "rc", "rm"),
pars=c("delta", "mu","sigma","deltanew"),
iter=500000, chains=3, warmup=500, verbose=FALSE)
print(mod_fit_skp)
# trace plots for deltanew
traceplot(mod_fit_ref, pars= "deltanew")
# plot(mod_fit_ref)
traceplot(mod_fit_skp, pars= "deltanew")
# plot(mod_fit_skp)
# histogram for reference posterior distribution
mag.sim.ref <- extract(mod_fit_ref, permuted=TRUE)
# hist(mag.sim.ref$deltanew, breaks=200, xlim=range(-4, 2),
# xlab="Posterior distribution of log of odds ratio",
# main="Reference prior histogram")
hist(exp(mag.sim.ref$deltanew), breaks=20000, xlim = range(0, 2),
xlab="Posterior distribution of odds ratio",
main="Reference prior histogram")
# histogram for skeptical posterior distribution
mag.sim.skp <- extract(mod_fit_skp, permuted=TRUE)
# hist(mag.sim.skp$deltanew, breaks=200, xlim=range(-3, 3),
# xlab="Posterior distribution of log of odds ratio",
# main="Skeptical prior histogram")
hist(exp(mag.sim.skp$deltanew), breaks=200000, xlim = range(0, 4),
xlab="Posterior distribution of odds ratio",
main="Skeptical prior histogram")
# statistical superiority
sum(exp(mag.sim.ref$deltanew) < 1) / length(mag.sim.ref$deltanew) # 0.90
sum(exp(mag.sim.skp$deltanew) < 1) / length(mag.sim.skp$deltanew) # 0.68
# clinical superiority
sum(exp(mag.sim.ref$deltanew) < 0.9) / length(mag.sim.ref$deltanew) # 0.87
sum(exp(mag.sim.skp$deltanew) < 0.9) / length(mag.sim.skp$deltanew) # 0.62