#### R CODE FOR "A BAYESIAN APROACH TO SEQUENTIAL META-ANALYSIS" #### #### G. T. Spence, D. Steinsaltz, T. R. Fanshawe #### ### 1. R code to simulate meta-analysis data # Input: # n = number of studies # theta, tsq = true effect size and heterogeneity # t, Imax = parameters determining average study precision (/size) # Output: # y.i = observed study effects # sigma.sq.i = observed study variances theta.i <- rnorm(n, theta, sqrt(tsq)) sigma.sq <- t/Imax sigma.sq.i <- runif(n, 0.25 * sigma.sq, 1.75 * sigma.sq) y.i <- rnorm(n, theta.i, sqrt(sigma.sq.i)) ### 2. R code to perform fully Bayesian sequential meta-analysis on ### 5,000 simulated meta-analyses applying different stopping rule thresholds # Priors: N(0, var) for theta, IG(eta, lambda) for tau.sq # Input: # n, theta, tsq, t, Imax = parameters for simulating meta-analysis data # delta = minimally relevant clinical effect size (e.g. logs odds ratio) # epsilon.e, epsilon.f = vectors of posterior decision thresholds # var, eta, lambda = prior parameters # iterations, chains, burn, thin = MCMC settings # Requires JAGS download: http://mcmc-jags.sourceforge.net # and `rjags' R package: https://cran.r-project.org/web/packages/rjags/ ## Define model and priors ## ================================================= write(" model { for (i in 1:n) { y.i[i] ~ dnorm(theta.i[i], 1 / sigma.sq.i[i]) theta.i[i] ~ dnorm(theta, inv.tsq)} # JAGS uses precision not variance theta ~ dnorm(mu, prec) inv.tsq ~ dgamma(eta, lambda) # Appropriate prior for precision tsq <- 1 / inv.tsq } ", "Bsim.bug") ## Bayesian sequential meta-analysis function with different thresholds ## ==== Bayes.sim.thres <- function(runs=5000, theta=0, tsq=0.0625, n=50, t=10, Imax=44.316, # Imax only defines sigma.sq epsilon.e=0.006, epsilon.f=0.01, delta=0.5, eta=1.5, lambda=0.08, var=5, chains=3, iterations=1E5, burn=100, thin=10){ epsilons <- data.frame(e=rep(epsilon.e, each=length(epsilon.f)), f=rep(epsilon.f, length(epsilon.e))) n.eps <- dim(epsilons)[1] # Settings settings <- data.frame(runs=runs, t=t, Imax=Imax, theta=theta, tsq=tsq, n=n, delta=delta, eta=eta, lambda=lambda, var=var) effect <- futility <- benefit <- concl <- stopped <- studies <- matrix(NA, runs, n.eps) diag <- matrix(NA, runs, n) # Simulation loop for (j in 1:runs){ if (j %% 100 == 0) cat(j, ", ", sep="") # Simulate meta-analysis data theta.i <- rnorm(n, theta, sqrt(tsq)) sigma.sq <- t/Imax sigma.sq.i <- runif(n, 0.25 * sigma.sq, 1.75 * sigma.sq) y.i <- rnorm(n, theta.i, sqrt(sigma.sq.i)) i <- 1; concl[j, ] <- 0 # Perform Bayesian inference, sequentially adding trials within loop while(sum(concl[j, ] == 0) > 0 & i <= n){ Y <- y.i[1:i] v <- sigma.sq.i[1:i] # Set-up and run simulation Bsim.jags <- jags.model("Bsim.bug", n.chains=chains, quiet=TRUE, list(y.i=Y, sigma.sq.i=v, n=i, eta=eta, lambda=lambda, mu=0, prec=1 / var), inits=list(theta=runif(1, -2, 2), inv.tsq=1 / runif(1, 0, 0.5))) update(Bsim.jags, burn, progress.bar="none") # Burn-in Bsim.sims <- coda.samples(Bsim.jags, c("theta", "tsq"), n.iter=iterations, thin=thin, progress.bar="none") all.chains <- as.mcmc(do.call(rbind, Bsim.sims)) # Check for convergence diag[j,i] <- gelman.diag(Bsim.sims)$mpsrf # At each sequential analysis, apply each set of threshold pairs for (k in 1:n.eps) { if (concl[j, k] == 0) { # HPD CI theta.HPD.eps <- HPDinterval(all.chains, 1 - epsilons$e[k])[1, ] # Status effect[j, k] <- ifelse(theta.HPD.eps[[1]] < 0 & theta.HPD.eps[[2]] > 0, 0, 1) benefit[j, k] <- ifelse(theta.HPD.eps[[1]] > 0, 1, 0) futility[j, k] <- ifelse(quantile(all.chains[, 1], 1 - epsilons$f[k]) < delta & quantile(all.chains[, 1], epsilons$f[k]) > -delta, 1, 0) stopped[j, k] <- ifelse(effect[j, k] == 0 & futility[j, k] == 0 & i == n, 1, 0) concl[j, k] <- ifelse(effect[j, k] == 1 | futility[j, k] == 1, 1, 0) studies[j, k] <- i } } # Warning if run truncated (i.e. n not big enough) if (sum(stopped[j, ] == 1) > 0) cat("Warning: meta-analysis truncated , ") i <- i + 1 } # Warning if possible problem with convergence if (sum(diag > 1.2, na.rm=TRUE) > 0) { cat("Warning: problem with convergence") } } # Function output x <- list(Settings=settings, Epsilon=epsilons, Concl=concl, Effect=effect, Benefit=benefit, Futility=futility, Stopped=stopped, Studies=studies, Diag=diag) class(x) <- "Bsimthres" x } ## Summary function ## ======================================================== summary.Bsimthres <- function(a) { list(a$Settings, data.frame("Epsilon.e"=a$Epsilon$e, "Epsilon.f"=a$Epsilon$f, "Studies"=apply(a$Studies, 2, mean), "Benefit"=apply(a$Benefit, 2, mean), "Effect"=apply(a$Effect, 2, mean))) } ## Code to replicate data in paper ## ========================================= t <- rep(c(5, 10, 20), each=6) theta <- rep(rep(c(0, 0.5), each=3), 3) tsq <- rep(c(0, 0.0625, 0.25), 6) n <- rep(c(50, 60, 75), each=6) set <- data.frame(t, theta, tsq, n) for (i in 1:18){ sim <- Bayes.sim.thres(epsilon.e=seq(0.002, 0.01, 0.002), epsilon.f=seq(0.002, 0.01, 0.002), t=set$t[i], theta=set$theta[i], tsq=set$tsq[i], n=set$n[i]) assign(paste0("sim.thres", i), sim) print(summary(sim)) } ### 3. R code to perform semi-Bayes sequential meta-analysis on 5,000 ### simulated meta-analyses # Prior: IG(eta, lambda) for tau.sq # Input: # n, theta, tsq, t, Imax = parameters for simulating meta-analysis data # delta = minimally relevant clinical effect size (e.g. logs odds ratio) # alpha, beta = target error levels for construction of Whitehead boundaries # eta, lambda = prior parameters # Requires JAGS download: http://mcmc-jags.sourceforge.net # and `rjags' R package: https://cran.r-project.org/web/packages/rjags/ ## Functions for semi-Bayes integration ## ==================================== num <- function(t2, eta, lambda, mu, y, var) { t2^(-eta) * exp(-lambda / t2) * prod(1 / sqrt(2 * pi * (t2 + var))) * exp(-0.5 * sum((y - mu)^2 / (var + t2))) } denom <- function(t2, eta, lambda, mu, m, y, var) { t2^(-eta - 1) * exp(-lambda / t2) * prod(1 / sqrt(2 * pi * (t2 + var))) * exp(-0.5 * sum((y - mu)^2 / (var + t2))) } ## Calculate restricted Whitehead boundaries with discrete corrections ## ===== Whitehead.boundaries <- function(H, Ik, Imax) { dk <- H - 0.583 * sqrt(pmax(0, diff(c(0, Ik)))) dk[(which(Ik >= Imax)[1] + 1):length(Ik)] <- NA dk } ## Record simulation measures ## ============================================== Measures <- function(stat, Ik, dk, tau.sq.method) { stopnumber <- ifelse(sum(abs(stat) > dk, na.rm=TRUE) > 0, which(abs(stat) > dk)[1], sum(!is.na(dk))) benefit <- ifelse(sum(stat > dk, na.rm=TRUE) > 0, 1, 0) effect <- ifelse(sum(abs(stat) > abs(dk), na.rm=TRUE) > 0, 1, 0) tau.sq <- tau.sq.method[stopnumber] info <- Ik[stopnumber] list(studies=stopnumber, benefit=benefit, effect=effect, tau.sq=tau.sq, info=info) } ## Semi-Bayes simulation function ## ========================================== SB.sim <- function(runs=5000, t=5, theta=0, tsq=0, n=50, alpha=0.05, beta=0.1, delta=0.5, eta=1.5, lambda=0.08) { # Restricted Whitehead boundary values val <- data.frame(alpha=c(0.001, 0.01, 0.05), beta=rep(c(0.2, 0.1, 0.05), each=3), H=c(14.576, 9.779, 6.457, 16.12, 11.029, 7.461, 17.394, 12.061, 8.288), Imax=c(17.535, 12.138, 8.299, 21.447, 15.438, 11.079, 24.972, 18.461, 13.673)) if (sum(alpha == c(0.001, 0.01, 0.05)) == 0 | sum(beta == c(0.2, 0.1, 0.05)) == 0) { stop("Boundaries not able to be calclulated from non-standard error values") } H <- val$H[val$alpha == alpha & val$beta == beta] / delta Imax <- val$Imax[val$alpha == alpha & val$beta == beta] / delta^2 # Output structures settings <- data.frame(runs=runs, t=t, theta=theta, tsq=tsq, n=n, alpha=alpha, beta=beta, delta=delta, eta=eta, lambda=lambda) studies <- benefit <- effect <- tau.sq <- info <- c() # Simulation loop for (j in 1:runs) { if (j %% 100 == 0) cat(j, ", ") # Simulate meta-analysis data theta.i <- rnorm(n, theta, sqrt(tsq)) sigma.sq <- t / Imax sigma.sq.i <- runif(n, 0.25 * sigma.sq, 1.75 * sigma.sq) y.i <- rnorm(n, theta.i, sqrt(sigma.sq.i)) ## Semi-Bayes random effects model # Heterogeneity, Information and Score tau.sq <- c(); mu <- c(); Ik <- c(); Sk <- c() for (i in 1:n) { if (i == 1) { num.val <- integrate(Vectorize(num, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=0, y=y.i[1:i], var=sigma.sq.i[1:i])$value denom.val <- integrate(Vectorize(denom, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=0, y=y.i[1:i], var=sigma.sq.i[1:i])$value } else { num.val <- integrate(Vectorize(num, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=mu[i - 1], y=y.i[1:i], var=sigma.sq.i[1:i])$value denom.val <- integrate(Vectorize(denom, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=mu[i - 1], y=y.i[1:i], var=sigma.sq.i[1:i])$value } tau.sq[i] <- num.val / denom.val Ik[i] <- sum(1 / (sigma.sq.i[1:i] + tau.sq[i])) Sk[i] <- sum(y.i[1:i] / (sigma.sq.i[1:i] + tau.sq[i])) mu[i] <- Sk[i] / Ik[i] } Zk <- Sk / sqrt(Ik) # Boundaries and Measures dk <- Whitehead.boundaries(H=H, Ik=Ik, Imax=Imax) B <- Measures(stat=Sk, Ik=Ik, dk=dk, tau.sq.method=tau.sq) ## Outputs studies[j] <- B$studies benefit[j] <- B$benefit effect[j] <- B$effect tau.sq[j] <- B$tau.sq info[j] <- B$info } x <- list(Settings=settings, Studies=studies, Benefit=benefit, Effect=effect, Tau.sq=tau.sq, Information=info) class(x) <- "SBsim" x } ## Summary function ## ======================================================== summary.SBsim <- function(a) { x <- data.frame(SemiBayes=c(mean(a$Studies), mean(a$Info), mean(a$Tau.sq), mean(a$Benefit), mean(a$Effect)), row.names=c("Studies at stopping", "Info at stopping", "Tau.sq mean", "Probability of benefit", "Probability of effect")) list(a$Settings, x) } ## Code to replicate data in paper ## ========================================= for (i in 1:18){ sim <- SB.sim(t=set$t[i], theta=set$theta[i], tsq=set$tsq[i], n=set$n[i]) assign(paste0("SBsim", i), sim) print(summary(sim)) } ### 4. R function for fully Bayesian sequential meta-analysis on simulated ### meta-analyses, using one set of thresholds, and getting parameter estimates # Priors: N(0, var) for theta, IG(eta, lambda) for tau.sq # Input: # n, theta, tsq, t, Imax = parameters for simulating meta-analysis data # delta = minimally relevant clinical effect size (e.g. logs odds ratio) # epsilon.e, epsilon.f = posterior decision thresholds # var, eta, lambda = prior parameters # iterations, chains, burn, thin = MCMC settings # Requires JAGS download: http://mcmc-jags.sourceforge.net # and `rjags' R package: https://cran.r-project.org/web/packages/rjags/ Bayes.sim <- function(runs=5000, theta=0, tsq=0.0625, n=50, t=10, Imax=44.316, # Imax only defines sigma.sq epsilon.e=0.006, epsilon.f=0.01, delta=0.5, eta=1.5, lambda=0.08, var=5, chains=3, iterations=1E5, burn=100, thin=10){ # Output structures settings <- data.frame(runs=runs, t=t, Imax=Imax, theta=theta, tsq=tsq, n=n, delta=delta, eta=eta, lambda=lambda, var=var, epsilon.e=epsilon.e, epsilon.f=epsilon.f) theta.cover <- tsq.cover <- tsq.mode <- tsq.median <- theta.mode <- theta.median <- c() effect <- futility <- benefit <- concl <- stopped <- studies <- c() diag <- matrix(NA, runs, n) # Simulation loop for (j in 1:runs) { if (j %% 100 == 0) cat(j, ", ") # Simulate meta-analysis data theta.i <- rnorm(n, theta, sqrt(tsq)) sigma.sq <- t / Imax sigma.sq.i <- runif(n, 0.25 * sigma.sq, 1.75 * sigma.sq) y.i <- rnorm(n, theta.i, sqrt(sigma.sq.i)) # Perform Bayesian inference, sequentially adding trials within loop i <- 1; concl[j] <- 0 while(concl[j] == 0 & i <= n){ Y <- y.i[1:i] v <- sigma.sq.i[1:i] # Set-up and run simulation Bsim.jags <- jags.model("Bsim.bug", n.chains=chains, quiet=TRUE, list(y.i=Y, sigma.sq.i=v, n=i, eta=eta, lambda=lambda, mu=0, prec=1 / var), inits=list(theta=runif(1, -2, 2), inv.tsq=1 / runif(1, 0, 0.5))) update(Bsim.jags, burn, progress.bar="none") # Burn-in Bsim.sims <- coda.samples(Bsim.jags, c("theta", "tsq"), n.iter=iterations, thin=thin, progress.bar="none") all.chains <- as.mcmc(do.call(rbind, Bsim.sims)) # Check for convergence diag[j,i] <- gelman.diag(Bsim.sims)$mpsrf # HPD CIs all.chains <- as.mcmc(do.call(rbind, Bsim.sims)) theta.HPD.eps <- HPDinterval(all.chains, 1 - epsilon.e)[1, ] theta.HPD.95 <- HPDinterval(all.chains)[1, ] tsq.HPD <- HPDinterval(all.chains)[2, ] # Posterior estimates for tsq tsq.post <- density(all.chains[, 2]) tsq.median[j] <- median(all.chains[, 2]) tsq.mode[j] <- tsq.post$x[tsq.post$y == max(tsq.post$y)] # Posterior estimates for theta theta.post <- density(all.chains[, 1]) theta.median[j] <- median(all.chains[, 1]) theta.mode[j] <- theta.post$x[theta.post$y == max(theta.post$y)] # Status effect[j] <- ifelse(theta.HPD.eps[[1]] < 0 & theta.HPD.eps[[2]] > 0, 0, 1) benefit[j] <- ifelse(theta.HPD.eps[[1]] > 0, 1, 0) futility[j] <- ifelse(quantile(all.chains[, 1], 1 - epsilon.f) < delta & quantile(all.chains[, 1], epsilon.f) > -delta, 1, 0) stopped[j] <- ifelse(effect[j] == 0 & futility[j] == 0 & i == n, 1, 0) concl[j] <- ifelse(effect[j] == 1 | futility[j] == 1, 1, 0) studies[j] <- i i <- i + 1 } # Warning if run truncated (i.e. n not big enough) if (stopped[j] == 1) cat("Warning: meta-analysis truncated , ") # Coverage at end point theta.cover[j] <- ifelse(theta.HPD.95[[1]] < theta & theta.HPD.95[[2]] > theta, 1, 0) tsq.cover[j] <- ifelse(tsq.HPD[[1]] < tsq & tsq.HPD[[2]] > tsq, 1, 0) } # Warning if possible problem with convergence if (sum(diag > 1.2, na.rm=TRUE) > 0) { cat("Warning: problem with convergence") } x <- list(Settings=settings, Studies=studies, Benefit=benefit, Effect=effect, Futility=futility, Stopped=stopped, Concl=concl, Tsq.mode=tsq.mode, Tsq.median=tsq.median, Theta.mode=theta.mode, Theta.median=theta.median, Theta.cover=theta.cover, Tsq.cover=tsq.cover, Diag=diag) class(x) <- "Bsim" x } ## Summary function ## ======================================================== summary.Bsim <- function(a) { x <- data.frame("Bayes"=c(mean(a$Studies), mean(a$Theta.median), mean(a$Tsq.median), mean(a$Theta.cover), mean(a$Tsq.cover), mean(a$Benefit)), row.names=c("Studies at stopping", "Theta median", "Tau.sq median", "Theta coverage", "Tau.sq coverage", "Probability of benefit")) list(a$Settings, x) } ## Code to replicate data in paper ## ========================================= for (i in 1:18){ sim <- Bayes.sim(epsilon.e=0.006, epsilon.f=0.01, t=set$t[i], theta=set$theta[i], tsq=set$tsq[i], n=set$n[i]) assign(paste0("Bsim", i), sim) print(summary(sim)) } ### 5. R code to perform frequentist, semi-Bayes and fully Bayesian sequential ### meta-analysis methods on real data from Cochrane Library # Priors: N(0, var) for theta, IG(eta, lambda) for tau.sq # Input: # data with 'yearofstudy', 'events1', 'total1', 'events2' & 'total2' columns # delta.OR = minimally relevant clinical effect size on odds ratio scale # epsilon.e, epsilon.f = vectors of posterior decision thresholds # alpha, beta = target error levels for construction of Whitehead boundaries # var, eta, lambda = prior parameters # iterations, chains, burn, thin = MCMC settings # Requires JAGS download: http://mcmc-jags.sourceforge.net # `rjags' R package: https://cran.r-project.org/web/packages/rjags/ # and `rmeta' R package: https://cran.r-project.org/web/packages/rmeta/ ## Adjusted measures for Whitehead methods (frequentist and semi-Bayes) ## ==== Measures.W <- function(stat, Ik, Imax, dk, tau.sq.method, mu) { result <- "Uncertain" if(sum(stat > dk, na.rm=TRUE) > 0) result <- "Harm" if(sum(stat < (-dk), na.rm=TRUE) > 0) result <- "Benefit" if(result == "Uncertain" & tail(Ik, 1) > Imax) result <- "Futility" if (result == "Benefit" | result == "Harm") { stopnumber <- which(abs(stat) > dk)[1] } if (result == "Futility") stopnumber <- sum(!is.na(dk)) if (result == "Uncertain") stopnumber <- length(stat) tau.sq <- tau.sq.method[stopnumber] theta <- mu[stopnumber] list(result=result, studies=stopnumber, theta=theta, tau.sq=tau.sq) } ## Apply sequential meta-analysis methods to data ## ========================== seq.MA <- function(data, delta.OR=1.65, alpha=0.05, beta=0.1, epsilon.e=0.006, epsilon.f=0.01, eta=1.5, lambda=0.08, prior.var=5, iterations=1E5, chains=3) { #Order by year d <- data[order(data$yearofstudy), ] #Continuity correction for 0/full counts, drop if full counts in both groups if(sum(d$events1==0)!=0) d[d$events1==0,][c(7, 10, 11, 14)] <- d[d$events1==0,][c(7, 10, 11, 14)] + c(0.5, 1) if(sum(d$events2==0)!=0) d[d$events2==0,][c(7, 10, 11, 14)] <- d[d$events2==0,][c(7, 10, 11, 14)] + c(0.5, 1) if(sum(d$events1==d$total1)!=0) d[d$events1==d$total1,][c(7, 10, 11, 14)] <- d[d$events1==d$total1,][c(7, 10, 11, 14)] + c(rep(0.5, sum(d$events1==d$total1)), rep(1, sum(d$events1==d$total1))) if(sum(d$events2==d$total2)!=0) d[d$events2==d$total2,][c(7, 10, 11, 14)] <- d[d$events2==d$total2,][c(7, 10, 11, 14)] + c(rep(0.5, sum(d$events2==d$total2)), rep(1, sum(d$events2==d$total2))) #Convert into log(OR) y.i <- log(d$events1*(d$total2-d$events2)/(d$events2*(d$total1-d$events1))) sigma.sq.i <- 1/d$events1 + 1/d$events2 + 1/(d$total1-d$events1) + 1/(d$total2-d$events2) ## Minimum clinically relevant effect size delta <- log(delta.OR) # Restricted Whitehead boundary values val <- data.frame(alpha=c(0.001, 0.01, 0.05), beta=rep(c(0.2, 0.1, 0.05), each=3), H=c(14.576, 9.779, 6.457, 16.12, 11.029, 7.461, 17.394, 12.061, 8.288), Imax=c(17.535, 12.138, 8.299, 21.447, 15.438, 11.079, 24.972, 18.461, 13.673)) if (sum(alpha == c(0.001, 0.01, 0.05)) == 0 | sum(beta == c(0.2, 0.1, 0.05)) == 0) { stop("Boundaries not able to be calclulated from non-standard error values") } H <- val$H[val$alpha == alpha & val$beta == beta] / delta Imax <- val$Imax[val$alpha == alpha & val$beta == beta] / delta^2 ## Frequentist (Non-Bayes) Method # Heterogeneity, Information and Score tau.sq.DL <- 0; Ik.NB <- 1 / sigma.sq.i[1] Sk.NB <- y.i[1] / sigma.sq.i[1] for (i in 2:length(y.i)) { tau.sq.DL[i] <- meta.summaries(y.i[1:i], sqrt(sigma.sq.i[1:i]), method="random")$tau2 Ik.NB[i] <- sum(1 / (sigma.sq.i[1:i] + tau.sq.DL[i])) Sk.NB[i] <- sum(y.i[1:i] / (sigma.sq.i[1:i] + tau.sq.DL[i])) } Zk.NB <- Sk.NB / sqrt(Ik.NB) mu.NB <- Sk.NB / Ik.NB # Boundaries and Measures dk.NB <- Whitehead.boundaries(H=H, Ik=Ik.NB, Imax=Imax) NB <- Measures.W(stat=Sk.NB, Ik=Ik.NB, Imax=Imax, dk=dk.NB, tau.sq.method=tau.sq.DL, mu=mu.NB) ## Semi-Bayes model # Heterogeneity, Information and Score tau.sq.SB <- c(); mu.SB <- c(); Ik.SB <- c(); Sk.SB <- c() for (i in 1:length(y.i)) { if (i == 1) { num.val <- integrate(Vectorize(num, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=0, y=y.i[1:i], var=sigma.sq.i[1:i])$value denom.val <- integrate(Vectorize(denom, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=0, y=y.i[1:i], var=sigma.sq.i[1:i])$value } else { num.val <- integrate(Vectorize(num, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=mu.SB[i - 1], y=y.i[1:i], var=sigma.sq.i[1:i])$value denom.val <- integrate(Vectorize(denom, vectorize.args='t2'), 0, Inf, eta=eta, lambda=lambda, mu=mu.SB[i - 1], y=y.i[1:i], var=sigma.sq.i[1:i])$value } tau.sq.SB[i] <- num.val / denom.val Ik.SB[i] <- sum(1 / (sigma.sq.i[1:i] + tau.sq.SB[i])) Sk.SB[i] <- sum(y.i[1:i] / (sigma.sq.i[1:i] + tau.sq.SB[i])) mu.SB[i] <- Sk.SB[i] / Ik.SB[i] } Zk.SB <- Sk.SB / sqrt(Ik.SB) # Boundaries and Measures dk.SB <- Whitehead.boundaries(H=H, Ik=Ik.SB, Imax=Imax) SB <- Measures.W(stat=Sk.SB, Ik=Ik.SB, Imax=Imax, dk=dk.SB, tau.sq.method=tau.sq.SB, mu=mu.SB) ## Fully Bayesian Method i <- 1; result.FB <- "Uncertain"; diag <- 0 while(result.FB == "Uncertain" & i <= length(y.i)) { Y <- y.i[1:i] v <- sigma.sq.i[1:i] # Set-up and run simulation Bsim.jags <- jags.model("Bsim.bug", n.chains = chains, quiet=TRUE, list(y.i=Y, sigma.sq.i=v, n=i, eta=eta, lambda=lambda, mu=0, prec=1 / prior.var), inits=list(theta=runif(1, -2, 2), inv.tsq=1 / runif(1, 0, 0.5))) update(Bsim.jags, 100, progress.bar="none") # Burn-in Bsim.sims <- coda.samples(Bsim.jags, c("theta", "tsq"), n.iter=iterations, thin=10, progress.bar="none") # Check for convergence diag[i] <- gelman.diag(Bsim.sims)$mpsrf # HPD CIs all.chains <- as.mcmc(do.call(rbind, Bsim.sims)) theta.HPD.eps <- HPDinterval(all.chains, 1 - epsilon.e)[1, ] theta.HPD.95 <- HPDinterval(all.chains)[1, ] tsq.HPD <- HPDinterval(all.chains)[2, ] # Status if (theta.HPD.eps[[1]] > 0) result.FB <- "Harm" if (theta.HPD.eps[[2]] < 0) result.FB <- "Benefit" if (quantile(all.chains[, 1], 1 - epsilon.f)[[1]] < abs(delta) & quantile(all.chains[, 1], epsilon.f)[[1]] > -abs(delta)) { result.FB <- "Futility" } studies.FB <- i i <- i + 1 } # Warning if possible problem with convergence if (sum(diag > 1.2, na.rm=TRUE) > 0) { cat("\n Warning: problem with convergence") } # Fully Bayesian Summary measures tau.sq.FB <- median(all.chains[, 2]) theta.FB <- median(all.chains[, 1]) FB <- list(result=result.FB, studies=studies.FB, theta=theta.FB, tau.sq=tau.sq.FB) ## Output data.frame(Result=c(NB$result, SB$result, FB$result), Studies=c(NB$studies, SB$studies, FB$studies), Theta=c(NB$theta, SB$theta, FB$theta), Tau.sq=c(NB$tau.sq, SB$tau.sq, FB$tau.sq), row.names=c("Non-Bayes", "Semi-Bayes", "Fully Bayes")) } ### 6. R code to perform fully Bayesian sequential meta-analysis with a ### log-normal heterogeneity prior on real data from Cochrane Library # Priors: N(0, var) for theta, LN(log.mean, log.var) for tau.sq # Input: # data with 'yearofstudy', 'events1', 'total1', 'events2' & 'total2' columns # delta.OR = minimally relevant clinical effect size on odds ratio scale # epsilon.e, epsilon.f = vectors of posterior decision thresholds # var, log.mean, log.var = prior parameters # iterations, chains, burn, thin = MCMC settings # Requires JAGS download: http://mcmc-jags.sourceforge.net # and `rjags' R package: https://cran.r-project.org/web/packages/rjags/ ### Define Bayesian model with log-normal prior for tsq ### =================== write(" model { for (i in 1:n) { y.i[i] ~ dnorm(theta.i[i], 1 / sigma.sq.i[i]) theta.i[i] ~ dnorm(theta, 1 / tsq)} theta ~ dnorm(mu, prec) tsq ~ dlnorm(log.mean, log.prec) } ", "BsimLN.bug") ## Apply Bayesian sequential meta-analysis using log-normal tsq prior## ======= Bayes.seq.MA.LN <- function(data, delta.OR=1.65, epsilon.e=0.006, epsilon.f=0.01, log.mean=-2.56, log.var=1.74^2, prior.var=5, iterations=1E5, chains=3) { #Order by year d <- data[order(data$yearofstudy), ] #Continuity correction for 0/full counts, drop if full counts in both groups if(sum(d$events1==0)!=0) d[d$events1==0,][c(7, 10, 11, 14)] <- d[d$events1==0,][c(7, 10, 11, 14)] + c(0.5, 1) if(sum(d$events2==0)!=0) d[d$events2==0,][c(7, 10, 11, 14)] <- d[d$events2==0,][c(7, 10, 11, 14)] + c(0.5, 1) if(sum(d$events1==d$total1)!=0) d[d$events1==d$total1,][c(7, 10, 11, 14)] <- d[d$events1==d$total1,][c(7, 10, 11, 14)] + c(rep(0.5, sum(d$events1==d$total1)), rep(1, sum(d$events1==d$total1))) if(sum(d$events2==d$total2)!=0) d[d$events2==d$total2,][c(7, 10, 11, 14)] <- d[d$events2==d$total2,][c(7, 10, 11, 14)] + c(rep(0.5, sum(d$events2==d$total2)), rep(1, sum(d$events2==d$total2))) #Convert into log(OR) y.i <- log(d$events1*(d$total2-d$events2)/(d$events2*(d$total1-d$events1))) sigma.sq.i <- 1/d$events1 + 1/d$events2 + 1/(d$total1-d$events1) + 1/(d$total2-d$events2) ## Minimum clinically relevant effect size delta <- log(delta.OR) ## Fully Bayesian Method i <- 1; result.FB <- "Uncertain"; diag <- 0 while(result.FB == "Uncertain" & i <= length(y.i)) { Y <- y.i[1:i] v <- sigma.sq.i[1:i] # Set-up and run simulation Bsim.jags <- jags.model("BsimLN.bug", n.chains = chains, quiet=TRUE, list(y.i=Y, sigma.sq.i=v, n=i, log.mean=log.mean, log.prec=1/log.var, mu=0, prec=1 / prior.var), inits=list(theta=runif(1, -2, 2), tsq=runif(1, 0, 0.5))) update(Bsim.jags, 100, progress.bar="none") # Burn-in Bsim.sims <- coda.samples(Bsim.jags, c("theta", "tsq"), n.iter=iterations, thin=10, progress.bar="none") # Check for convergence diag[i] <- gelman.diag(Bsim.sims)$mpsrf # HPD CIs all.chains <- as.mcmc(do.call(rbind, Bsim.sims)) theta.HPD.eps <- HPDinterval(all.chains, 1 - epsilon.e)[1, ] theta.HPD.95 <- HPDinterval(all.chains)[1, ] tsq.HPD <- HPDinterval(all.chains)[2, ] # Status if (theta.HPD.eps[[1]] > 0) result.FB <- "Harm" if (theta.HPD.eps[[2]] < 0) result.FB <- "Benefit" if (quantile(all.chains[, 1], 1 - epsilon.f)[[1]] < abs(delta) & quantile(all.chains[, 1], epsilon.f)[[1]] > -abs(delta)) { result.FB <- "Futility" } studies.FB <- i i <- i + 1 } # Warning if possible problem with convergence if (sum(diag > 1.2, na.rm=TRUE) > 0) { cat("\n Warning: problem with convergence") } # Fully Bayesian Summary measures tau.sq.FB <- median(all.chains[, 2]) theta.FB <- median(all.chains[, 1]) FB <- list(result=result.FB, studies=studies.FB, theta=theta.FB, tau.sq=tau.sq.FB) ## Output data.frame(Result=FB$result, Studies=FB$studies, Theta=FB$theta, Tau.sq=FB$tau.sq, row.names="Fully Bayes") }