if (FALSE) {
######################################################################
# (1) A simple example with two groups of studies
# load data:
data("CrinsEtAl2014")
# compute effect measures (log-OR):
crins.es <- escalc(measure="OR",
ai=exp.AR.events, n1i=exp.total,
ci=cont.AR.events, n2i=cont.total,
slab=publication, data=CrinsEtAl2014)
# show data:
crins.es[,c("publication", "IL2RA", "exp.AR.events", "exp.total",
"cont.AR.events", "cont.total", "yi", "vi")]
# specify regressor matrix:
X <- cbind("bas"=as.numeric(crins.es$IL2RA=="basiliximab"),
"dac"=as.numeric(crins.es$IL2RA=="daclizumab"))
print(X)
print(cbind(crins.es[,c("publication", "IL2RA")], X))
# NB: regressor matrix specifies individual indicator covariates
# for studies with "basiliximab" and "daclizumab" treatment.
# perform regression:
bmr01 <- bmr(y=crins.es$yi, sigma=sqrt(crins.es$vi),
labels=crins.es$publication, X=X)
# alternatively, one may simply supply the "escalc" object
# (yields identical results):
bmr01 <- bmr(crins.es, X=X)
# show results:
bmr01
bmr01$summary
plot(bmr01)
pairs(bmr01)
# NOTE: there are many ways to set up the regressor matrix "X"
# (also affecting the interpretation of the involved parameters).
# See the above specification and check out the following alternatives:
X <- cbind("bas"=1, "offset.dac"=c(1,0,1,0,0,0))
X <- cbind("intercept"=1, "offset"=0.5*c(1,-1,1,-1,-1,-1))
# One may also use the "model.matrix()" function
# to specify regressor matrices via the "formula" interface; e.g.:
X <- model.matrix( ~ IL2RA, data=crins.es)
X <- model.matrix( ~ IL2RA - 1, data=crins.es)
######################################################################
# (2) A simple example reproducing a "bayesmeta" analysis:
data("CrinsEtAl2014")
crins.es <- escalc(measure="OR",
ai=exp.AR.events, n1i=exp.total,
ci=cont.AR.events, n2i=cont.total,
slab=publication, data=CrinsEtAl2014)
# a "simple" meta-analysis:
bma02 <- bayesmeta(crins.es,
tau.prior=function(t){dhalfnormal(t, scale=0.5)},
mu.prior.mean=0, mu.prior.sd=4)
# the equivalent "intercept-only" meta-regression:
bmr02 <- bmr(crins.es,
tau.prior=function(t){dhalfnormal(t, scale=0.5)},
beta.prior.mean=0, beta.prior.sd=4)
# the corresponding (default) regressor matrix:
bmr02$X
# compare computation time and numbers of bins used internally:
cbind("seconds" = c("bayesmeta" = unname(bma02$init.time),
"bmr" = unname(bmr02$init.time)),
"bins" = c("bayesmeta" = nrow(bma02$support),
"bmr" = nrow(bmr02$support$tau)))
# compare heterogeneity estimates:
rbind("bayesmeta"=bma02$summary[,1], "bmr"=bmr02$summary[,1])
# compare effect estimates:
rbind("bayesmeta"=bma02$summary[,2], "bmr"=bmr02$summary[,2])
######################################################################
# (3) An example with binary as well as continuous covariables:
# load data:
data("RobergeEtAl2017")
str(RobergeEtAl2017)
head(RobergeEtAl2017)
?RobergeEtAl2017
# compute effect sizes (log odds ratios) from count data:
es.pe <- escalc(measure="OR",
ai=asp.PE.events, n1i=asp.PE.total,
ci=cont.PE.events, n2i=cont.PE.total,
slab=study, data=RobergeEtAl2017,
subset=complete.cases(RobergeEtAl2017[,7:10]))
# show "bubble plot" (bubble sizes are
# inversely proportional to standard errors):
plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
col=c("blue","red")[as.numeric(es.pe$onset)],
xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
# set up regressor matrix:
# (individual intercepts and slopes for two subgroups):
X <- model.matrix(~ -1 + onset + onset:dose, data=es.pe)
colnames(X) <- c("intEarly", "intLate", "slopeEarly", "slopeLate")
# check out regressor matrix (and compare to original data):
print(X)
# perform regression:
bmr03 <- bmr(es.pe, X=X)
bmr03$summary
# derive predictions from the model;
# specify corresponding "regressor matrices":
newx.early <- cbind(1, 0, seq(50, 150, by=5), 0)
newx.late <- cbind(0, 1, 0, seq(50, 150, by=5))
# (note: columns correspond to "beta" parameters)
# compute predicted medians and 95 percent intervals:
pred.early <- cbind("median"=bmr03$qpred(0.5, x=newx.early),
bmr03$pred.interval(x=newx.early))
pred.late <- cbind("median"=bmr03$qpred(0.5, x=newx.late),
bmr03$pred.interval(x=newx.late))
# draw "bubble plot":
plot(es.pe$dose, es.pe$yi, cex=1/sqrt(es.pe$vi),
col=c("blue","red")[as.numeric(es.pe$onset)],
xlab="dose (mg)", ylab="log-OR (PE)", main="Roberge et al. (2017)")
legend("topright", col=c("blue","red"), c("early onset", "late onset"), pch=1)
# add predictions to bubble plot:
matlines(newx.early[,3], pred.early, col="blue", lty=c(1,2,2))
matlines(newx.late[,4], pred.late, col="red", lty=c(1,2,2))
}
Run the code above in your browser using DataLab