if (FALSE) {
#################################################################
# perform a meta-analysis using binary ("indicator") covariables:
# 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 (binary indicator variables):
X <- cbind("basiliximab"=as.numeric(crins.es$IL2RA=="basiliximab"),
"daclizumab" =as.numeric(crins.es$IL2RA=="daclizumab"))
print(X)
# perform meta-analysis:
bmr01 <- bmr(crins.es, X=X,
tau.prior=function(t){dhalfnormal(t, scale=0.5)})
# show forest plot:
forestplot(bmr01)
# show forest plot including contrast
# (difference between the two groups):
forestplot(bmr01,
X.mean=rbind("basiliximab" = c(1, 0),
"daclizumab" = c(0, 1),
"group difference" = c(-1, 1)))
##############################################
# perform the meta-analysis using a different
# ("intercept / slope") regressor setup:
X <- cbind("intercept"=1,
"offset.dac"=as.numeric(crins.es$IL2RA=="daclizumab"))
print(X)
# perform meta-analysis:
bmr02 <- bmr(crins.es, X=X,
tau.prior=function(t){dhalfnormal(t, scale=0.5)})
# show default forest plot:
forestplot(bmr02)
# show forest plot including both group means and their difference:
forestplot(bmr02,
X.mean=rbind("basiliximab" = c(1, 0),
"daclizumab" = c(1, 1),
"group difference" = c(0, 1)))
###############################################################
# a meta analysis using a continuous regressor
# and including prediction:
help("NicholasEtAl2019")
# load data:
data("NicholasEtAl2019")
# compute effect sizes (logarithic odds) from count data:
es <- escalc(measure="PLO",
xi=patients*(prog.percent/100), ni=patients,
slab=study, data=NicholasEtAl2019)
# set up regressor matrix:
X <- cbind("intercept2000" = 1, "year" = (es$year-2000))
# perform analysis:
bmr03 <- bmr(es, X=X)
# show forest plot including some mean estimates for the
# years from 1990 to 2018, and a prediction for 2019:
forestplot(bmr03,
X.mean=rbind("intercept (2000)" = c(1, 0),
"annual change" = c(0, 1),
"change per decade" = c(0, 10),
"mean 1990" = c(1, -10),
"mean 2000" = c(1, 0),
"mean 2010" = c(1, 10),
"mean 2018" = c(1, 18)),
X.predict=rbind("prediction 2019" = c(1, 19)),
xlab="log-odds",
txt_gp = fpTxtGp(ticks = gpar(cex=1), xlab = gpar(cex=1)))
# the shown summaries and predictions may also be computed "manually";
# mean effect (year 2018), posterior median and 95 percent CI:
bmr03$qpredict(p=0.5, x=c(1, 18))
bmr03$pred.interval(level=0.95, x=c(1, 18))
# prediction (year 2019), posterior median and 95 percent CI:
bmr03$qpredict(p=0.5, x=c(1, 19), mean=FALSE)
bmr03$pred.interval(level=0.95, x=c(1, 19), mean=FALSE)
# means and predictions may also be derived
# using the "summary()" function:
summary(bmr03,
X.mean=rbind("intercept (2000)" = c(1, 0),
"annual change" = c(0, 1),
"change per decade" = c(0, 10),
"mean 1990" = c(1, -10),
"mean 2000" = c(1, 0),
"mean 2010" = c(1, 10),
"mean 2018" = c(1, 18)),
X.predict=rbind("prediction 2019" = c(1, 19)))
##########################################################
# the tabular part of the forest plot may also be changed;
# draw a default plot:
forestplot(bmr03)
# don't plot, only extract the tabular bits:
fp <- forestplot(bmr03, plot=FALSE)
labtxt <- fp$labeltext
head(labtxt)
# drop two columns:
labtxt <- labtxt[,-c(2,3)]
# add two new columns:
labtxt <- cbind(labtxt[,1],
c("year", es$year, "",""),
c("events / total",
paste(round(es$patients*(es$prog.percent/100)),
"/", es$patients), "",""),
labtxt[,2:3])
head(labtxt)
# draw new forest plot:
forestplot(bmr03, labeltext=labtxt, xlab="log-odds")
}
Run the code above in your browser using DataLab