### Reproduce ### Prior densities of Figure 1. library("bayesmeta") data("CrinsEtAl2014", package = "bayesmeta") CrinsEtAl2014 library("metafor") crins.es <- escalc(measure = "OR", ai = exp.AR.events, n1i = exp.total, ci = cont.AR.events, n2i = cont.total, slab = publication, data = CrinsEtAl2014) crins.es # compute analyses # (eventually only the corresponding priors are required): ma01.shrinkage <- bayesmeta(crins.es, tau.prior="shrinkage") ma01.dumouchel <- bayesmeta(crins.es, tau.prior="dumouchel") ma01.conventional <- bayesmeta(crins.es, tau.prior="conventional") ma01.jeffreys <- bayesmeta(crins.es, tau.prior="jeffreys") ma01.bergerdeely <- bayesmeta(crins.es, tau.prior="bergerdeely") # compute harmonic mean of squared standard errors: harmonic <- nrow(crins.es)/sum(1/crins.es$vi) # derive corresponding "Turner et al." prior: tp <- TurnerEtAlPrior(outcome = "surgical", comparator1 = "pharmacological", comparator2 = "placebo / control") # specify colours and line styles: colvec <- c(rep(c("blue","green3","red3","black"),c(3,3,2,1))) ltyvec <- c("solid","dashed","12", "solid","dashed","12", "solid","dashed", "solid") # "tau" argument for plotting: tau <- seq(0, 2.2, le=200) # create empty plot: plot(c(0,2.1), c(0,2), type = "n", axes = FALSE, xlab = "", ylab = "") # add axis labels: mtext(expression(tau), side = 1, line = par("mgp")[1]) mtext(expression(p(tau)), side = 2, line = par("mgp")[2]) # draw axes: abline(h = 0, v = 0, col = "darkgray") # show s_0 value: abline(v = sqrt(harmonic), col = "gray", lty = "dotted") # add density lines: lines(tau, dhalfnormal(tau,scale = 0.5), col = colvec[1], lty = ltyvec[1]) lines(tau, dhalfcauchy(tau,scale = 0.5), col = colvec[2], lty = ltyvec[2]) lines(tau, tp$dprior(tau), col = colvec[3], lty = ltyvec[3]) lines(tau, ma01.shrinkage$dprior(tau), col = colvec[4], lty = ltyvec[4]) lines(tau, ma01.dumouchel$dprior(tau), col = colvec[5], lty = ltyvec[5]) lines(tau, ma01.conventional$dprior(tau), col = colvec[6], lty = ltyvec[6]) lines(tau, 0.5*ma01.jeffreys$dprior(tau), col = colvec[7], lty = ltyvec[7]) # (note the arbitrary scaling) lines(tau, ma01.bergerdeely$dprior(tau), col = colvec[8], lty = ltyvec[8]) # add tick marks and box: axis(3, at = sqrt(harmonic), lab = expression(italic(s)[0])) axis(1); box() # add legend: legend("topright", col = colvec, lty = ltyvec, bg = "white", c("half-normal(0.5)", "half-Cauchy(0.5)", "log-normal(-1.07, 0.87)", "uniform shrinkage", "DuMouchel", "conventional", "Jeffreys", "Berger-Deely")) ############## # Section 3.2 # Example data: a systematic review in immunosuppression # ------------------------------------------------------------------------- data("CrinsEtAl2014", package = "bayesmeta") CrinsEtAl2014 crins.es <- escalc(measure = "OR", ai = exp.AR.events, n1i = exp.total, ci = cont.AR.events, n2i = cont.total, slab = publication, data = CrinsEtAl2014) crins.es ############## # Section 3.3 # Performing a Bayesian random-effects meta-analysis # The bayesmeta() function ma01 <- bayesmeta(y = crins.es[, "yi"], sigma = sqrt(crins.es[, "vi"]), labels = crins.es[, "publication"], mu.prior.mean = 0, mu.prior.sd = 4, tau.prior = function(t) dhalfnormal(t, scale = 0.5)) ma01 <- bayesmeta(crins.es, mu.prior.mean = 0, mu.prior.sd = 4, tau.prior = function(t) dhalfnormal(t, scale = 0.5)) ma01 # The forestplot() function forestplot(ma01) # The plot() function plot(ma01) # Elements of the bayesmeta() output ma01$summary ma01$y ma01$sigma ma01$labels ma01$k x <- seq(-3, 0.5, length = 200) plot(x, ma01$dposterior(mu = x), type = "l", xlab = "effect", ylab = "posterior density") abline(h = 0, v = 0, col = "gray") 1 - ma01$pposterior(mu = 0) x <- seq(-3, 0.5, length = 200) plot(x, ma01$pposterior(mu = x), type = "l", xlab = "effect", ylab = "posterior CDF") abline(h = 0:1, v = 0, col = "gray") 1 - ma01$pposterior(tau = 1) ma01$qposterior(tau.p = 0.99) set.seed(123) ma01$rposterior(n = 5) ma01$rposterior(n = 5, tau.sample = FALSE) prob.control <- 0.5 logodds.control <- log(prob.control / (1 - prob.control)) logodds.treat <- (logodds.control + ma01$rposterior(n = 10000, tau.sample = FALSE)) prob.treat <- exp(logodds.treat) / (1 + exp(logodds.treat)) riskdiff <- (prob.treat - prob.control) median(riskdiff) quantile(riskdiff, c(0.025, 0.975)) # Credible intervals ma01$post.interval(tau.level = 0.99) ma01$post.interval(tau.level = 0.99, method = "central") ma01$qposterior(tau.p = c(0.005, 0.995)) # Prediction x <- seq(-3.5, 0.5, length = 200) plot(x, ma01$dposterior(mu = x), type = "n", xlab = "effect", ylab = "probability density") abline(h = 0, v = 0, col = "gray") lines(x, ma01$dposterior(mu = x), col = "red") lines(x, ma01$dposterior(mu = x, predict = TRUE), col = "blue") # Shrinkage ma01$theta[,1:2] x <- seq(-4, 0.5, length = 200) plot(x, ma01$dposterior(theta = x, individual = 1), type = "n", xlab = "effect", ylab = "probability density") abline(h = 0, v = 0, col = "gray") lines(x, dnorm(x, mean = ma01$y[1], sd = ma01$sigma[1]), col = "green", lty = "dashed") lines(x, ma01$dposterior(theta = x, individual = 1), col = "green") ############## # Section 3.4 # Investigating prior variations # Prior predictive distributions hn05 <- normalmixture(cdf = function(t) phalfnormal(t, scale = 0.5)) hn10 <- normalmixture(cdf = function(t) phalfnormal(t, scale = 1.0)) hc05 <- normalmixture(cdf = function(t) phalfcauchy(t, scale = 0.5)) x <- seq(-1, 3, by = 0.01) plot(x, hn05$cdf(x), type = "l", col = "blue", ylim = 0:1, xlab = expression(theta[i]), ylab = "prior predictive CDF") lines(x, hn10$cdf(x), col = "green") lines(x, hc05$cdf(x), col = "red") abline(h = 0:1, v = 0, col = "gray") axis(3, at = log(c(0.5, 1, 2, 5, 10, 20)), lab = c(0.5, 1, 2, 5, 10, 20)) mtext(expression(exp(theta[i])), side = 3, line = 2.5) q975 <- c("half-normal(0.5)" = hn05$quantile(0.975), "half-normal(1.0)" = hn10$quantile(0.975), "half-Cauchy(0.5)" = hc05$quantile(0.975)) print(cbind("theta" = q975, "exp(theta)" = exp(q975))) # Informative heterogeneity priors tp <- TurnerEtAlPrior(outcome = "surgical", comparator1 = "pharmacological", comparator2 = "placebo / control") tp$qprior(c(0.025, 0.5, 0.975)) ma02 <- bayesmeta(crins.es, mu.prior.mean = 0, mu.prior.sd = 4, tau.prior = tp$dprior) # Non-informative priors ma03 <- bayesmeta(crins.es, tau.prior = "Jeffreys") ############## # Section 3.5 # Making the connection with frequentist results ma04 <- rma(crins.es) sqrt(ma04$tau2) ma04$b ma04$se ma03$cond.moment(tau = sqrt(ma04$tau2)) ############## # Section 3.6 # Posterior predictive checks # A meta-analysis of two studies ma05 <- bayesmeta(crins.es[crins.es[, "randomized"] == "yes",], mu.prior.mean = 0, mu.prior.sd = 4, tau.prior = function(t) dhalfnormal(t, scale = 0.5)) # The forestplot() function forestplot(ma05) # Posterior predictive p-values for the effect (mu) ma05$pposterior(mu = 0) p1 <- pppvalue(ma05, parameter = "mu", value = 0, alternative = "less", statistic = "cdf", n = 1000) p1 plot(ma05, which = 2, mulim = c(-3.5, 1), taulim = c(0, 2)) abline(h = p1$null.value) # (the null-hypothesized mu value) points(p1$replicates$tau, p1$replicates$mu, col = "cyan") # (the samples) plot(ecdf(p1$replicates$statistic[, 1])) abline(v = p1$statistic, col = "red") abline(h = 1 - p1$p.value, col = "green") # Posterior predictive p-values for the heterogeneity (tau) p2 <- pppvalue(ma05, parameter = "tau", value = 0, alternative = "greater", statistic = "q", n = 1000) BF <- function(y, sigma) { bm <- bayesmeta(y = y, sigma = sigma, mu.prior.mean = 0, mu.prior.sd = 4, tau.prior = function(t) dhalfnormal(t, scale = 0.5), interval.type = "central") return(bm$bayesfactor[1, "tau=0"]) } p3 <- pppvalue(ma05, parameter = "tau", value = 0, alternative = "greater", statistic = BF, rejection.region = "lower.tail", n = 1000, sigma = ma05$sigma) # Posterior predictive p-values for individual effects (theta[i]) p4 <- pppvalue(ma05, parameter = "Spada", value = 0, alternative = "less", statistic = "cdf", n = 1000) ############### # Appendix A.5 # Calibration check mupriormean <- 0.0 mupriorsd <- 4.0 taupriorscale <- 0.5 Nsim <- 1000 pit <- matrix(NA_real_, nrow = Nsim, ncol = 2, dimnames = list(NULL, c("mu", "tau"))) for (i in 1:Nsim) { # generate data: mu <- rnorm(n = 1, mean = mupriormean, sd = mupriorsd) # effect tau <- rhalfnormal(n = 1, scale = taupriorscale) # heterogeneity k <- sample(c(2,3,5,10,20), size = 1) # number of studies sigma <- runif(n = k, min = 0.2, max = 1.0) # standard errors y <- rnorm(n = k, mean = mu, sd = sqrt(sigma^2+tau^2)) # estimates # perform analysis: bma <- try(bayesmeta(y = y, sigma = sigma, tau.prior = function(t) dhalfnormal(t, scale = taupriorscale), mu.prior = c(mupriormean, mupriorsd))) # log probability integral transform (PIT) values: if (!is(bma, "try-error")) { pit[i, "mu"] <- bma$pposterior(mu = mu) pit[i, "tau"] <- bma$pposterior(tau = tau) } } plot(ecdf(pit[, "mu"]), col = "blue", main = "effect (mu)", xlab = "PIT value", ylab = "CDF") lines(0:1, 0:1, lty = "dashed", lwd = 2) plot(ecdf(pit[, "tau"]), col = "red", main = "heterogeneity (tau)", xlab = "PIT value", ylab = "CDF") lines(0:1, 0:1, lty = "dashed", lwd = 2)