##################################################
#
#  R code for presentation
#  "Meta-analysis using the bayesmeta package"
#  by Christian Röver
#  2023-08-24
#

library("bayesmeta")


################################################################################
# check out and load example data set
# (Karner et al. (2014); https://doi.org/10.1002/14651858.CD009285.pub3):

?KarnerEtAl2014
data("KarnerEtAl2014")
str(KarnerEtAl2014)
head(KarnerEtAl2014)


################################################################################
# exacerbation endpoint, long studies only:
# (Analysis 1.10.2, page 73) 

# effect measure (log-OR, 1st study, Bateman 2010a):
KarnerEtAl2014[1,c("study","tiotropium.exa","tiotropium.total","placebo.exa","placebo.total")]

# compute OR:
(685/1304) / (842/1160)

# compute log-OR:
log((685/1304) / (842/1160))

# standard error:
sqrt(sum(1 / c(685, 1304, 842, 1160)))

# same calculation using the metafor package's "escalc()" function:
escalc(measure="OR", ai=685, bi=1304, ci=842, di=1160)
escalc(measure="OR", ai=685, n1i=1989, ci=842, n2i=2002)

# effect sizes and standard errors for all eight "long" studies:
exa.long <- escalc(measure="OR",
                   ai=tiotropium.exa, n1i=tiotropium.total,
                   ci=placebo.exa,    n2i=placebo.total,
                   slab=study, subset = (duration=="1 year or longer"),
                   data=KarnerEtAl2014)

# illustrate in forest plot:
forestplot(exa.long, xlab="log-OR")


################################################################################
# analyze:
bma01 <- bayesmeta(exa.long,
                   mu.prior=c(0,4),
                   tau.prior=function(tau){dhalfnormal(tau,scale=0.5)})

# show plain output:
bma01
bma01$summary

# plots:
plot(bma01)
plot(bma01, prior=TRUE) # (including prior densities in marginal posteriors)

forestplot(bma01, xlab="log-OR")

# posterior density, CDF, quantiles, etc.:

bma01$pposterior(mu=0)

bma01$qposterior(mu.p=0.90)

x <- seq(-1, 0.5, le=200)
plot(x, bma01$dposterior(mu=x), type="l")
abline(h=0, v=0, col="grey")

# trace plot:
traceplot(bma01)

# shrinkage estimates:
bma01$theta


################################################################################
# sensitivity analyses;
# analysis with uniform effect prior:
bma02 <- bayesmeta(exa.long,
                   tau.prior=function(tau){dhalfnormal(tau,scale=0.5)})

# analysis with uniform effect + heterogeneity priors:
bma03 <- bayesmeta(exa.long)

# compare resulting effect estimates:
rbind("original"      =bma01$summary[,"mu"],
      "uniform effect"=bma02$summary[,"mu"],
      "both uniform"  =bma03$summary[,"mu"])

# compare resulting heterogeneity estimates:
rbind(bma01$summary[,"tau"],
      bma02$summary[,"tau"],
      bma03$summary[,"tau"])


################################################################################
#
#  Meta-regression
#
################################################################################

# compute effect sizes (log-ORs) for ALL 22 studies:
exa.all <- escalc(measure="OR",
                  ai=tiotropium.exa, n1i=tiotropium.total,
                  ci=placebo.exa,    n2i=placebo.total,
                  slab=study, # subset = (duration=="1 year or longer"),
                  data=KarnerEtAl2014)

################################
# show data, including duration
# (as a binary covariable):
exa.all[,c("study", "duration", "yi", "vi")]

# assemble regressor matrix:
X1 <- cbind("short"=as.numeric(exa.all$duration=="up to 1 year"),
            "long" =as.numeric(exa.all$duration=="1 year or longer"))
head(X1)

# perform analysis:
bmr01 <- bmr(exa.all, X=X1,
             tau.prior=function(tau){dhalfnormal(tau, scale=0.5)})

# show forest plot:
forestplot(bmr01, xlab="log-OR")

# show forest plot with additional linear combination:
forestplot(bmr01, xlab="log-OR",
           X.mean=rbind("short"=c(1,0),
                        "long"=c(0,1),
                        "difference"=c(-1,1)))

# show estimates only:
summary(bmr01, X.mean=rbind("short"=c(1,0),
                            "long"=c(0,1),
                            "difference"=c(-1,1)))

summary(bmr01, X.mean=rbind("short"=c(1,0),
                            "long"=c(0,1),
                            "difference"=c(-1,1)))$mean


##########################################
# show data, including FEV1 (% predicted)
# (a continuous covariable):
exa.all[,c("study", "baseline.fev1pp", "yi", "vi")]
# (note: missing data for 13th study)

# assemble regressor matrix (intercept/slope):
X2 <- cbind("intercept" = 1,
            "FEV1"      = exa.all$baseline.fev1pp)
head(X2)

# perform analysis:
bmr02 <- bmr(exa.all[-13,], X=X2[-13,],
             tau.prior=function(tau){dhalfnormal(tau, scale=0.5)})

# show forestplot (default layout):
forestplot(bmr02, xlab="log-OR")

# show forestplot (with sensible covariable settings):
forestplot(bmr02, xlab="log-OR",
           X.mean=rbind("intercept" = c(1,0),
                        "slope"     = c(0,1),
                        "FEV1=40%"  = c(1,40),
                        "FEV1=50%"  = c(1,50),
                        "FEV1=60%"  = c(1,60),
                        "FEV1=70%"  = c(1,70)))

# estimates for pre-specified covariable values:
summary(bmr02,
        X.mean=rbind("intercept" = c(1,0),
                     "slope"     = c(0,1),
                     "FEV1=40%"  = c(1,40),
                     "FEV1=50%"  = c(1,50),
                     "FEV1=60%"  = c(1,60),
                     "FEV1=70%"  = c(1,70)))

summary(bmr02,
        X.mean=rbind("intercept" = c(1,0),
                     "slope"     = c(0,1),
                     "FEV1=40%"  = c(1,40),
                     "FEV1=50%"  = c(1,50),
                     "FEV1=60%"  = c(1,60),
                     "FEV1=70%"  = c(1,70)))$mean