library(metafor) ; library(brms); library(bayesmeta)
## Warning: package 'metafor' was built under R version 4.0.2
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.0.2
## Loading 'metafor' package (version 2.4-0). For an overview
## and introduction to the package please type: help(metafor).
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 4.0.2
## Loading 'brms' package (version 2.13.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:stats':
##
## ar
## Warning: package 'bayesmeta' was built under R version 4.0.2
## Loading required package: forestplot
## Warning: package 'forestplot' was built under R version 4.0.2
## Loading required package: grid
## Loading required package: magrittr
## Warning: package 'magrittr' was built under R version 4.0.2
## Loading required package: checkmate
# Full data from Andy Hill paper
Trial = c("Mahmud", "Niaee", "Hashim", "Elgazzar", "Okumus", "Kirti", "Rezai", "Abd-Elsalam", "Gonzalez", "Lopez-Medina", "Fonseca")
x1 = c(0, 4, 2, 2, 6, 0, 1, 3, 5, 0, 12)
n1 = c(183, 120, 70, 200, 30, 55, 35, 82, 36, 200, 53)
y1 = n1-x1
x0 = c(3, 11, 6, 24, 9, 4, 0.5, 4, 6, 1, 25)
n0 = c(180, 60, 70, 200, 30, 57, 34, 82, 37, 98, 115)
y0 = n0-x0
df1 = data.frame(Trial, x1, y1, n1, x0, y0, n0)
# Data without Elgazzar study
Trial = c("Mahmud", "Niaee", "Hashim", "Okumus", "Kirti", "Rezai", "Abd-Elsalam", "Gonzalez", "Lopez-Medina", "Fonseca")
x1 = c(0, 4, 2, 6, 0, 1, 3, 5, 0, 12)
n1 = c(183, 120, 70, 30, 55, 35, 82, 36, 200, 53)
y1 = n1-x1
x0 = c(3, 11, 6, 9, 4, 0, 4, 6, 1, 25)
n0 = c(180, 60, 70, 30, 57, 34, 82, 37, 98, 115)
y0 = n0-x0
df2= data.frame(Trial, x1, y1, n1, x0, y0, n0)
Analysis
dat = escalc(measure="RR", ai=x1, bi=y1, ci=x0, di=y0, data=df1, append=T)
res = rma(yi, vi, data = dat)
forest(res, slab = paste(dat$Trial), atransf = exp)
# Bayesian analysis
bayes = bayesmeta(y=dat[,"yi"], sigma=dat[,"vi"], labels=dat[,"Trial"], mu.prior.mean=0, mu.prior.sd=4, tau.prior=function(t){dhalfnormal(t,scale=0.5)} )
forestplot(bayes)
bayes$summary
## tau mu theta
## mode 0.6597662 -0.7895313 -0.7749340
## median 0.7031287 -0.8091607 -0.8053067
## mean 0.7249432 -0.8187005 -0.8187005
## sd 0.2013534 0.3119465 0.8154821
## 95% lower 0.3621767 -1.4491730 -2.4750790
## 95% upper 1.1293846 -0.2119467 0.8026270
exp(bayes$summary[,"mu"])
## mode median mean sd 95% lower 95% upper
## 0.4540576 0.4452316 0.4410044 1.3660816 0.2347644 0.8090078