#########################################################################
# Example: Geeganage2010 data
# Preparing covariances for a multivariate meta-analysis
# Choose variable SBP, DBP, DD, D with effect sizes "MD","MD","RD","lgOR"
#########################################################################
data(Geeganage2010)
## set the correlation coefficients list r
r12 <- 0.71
r13 <- 0.5
r14 <- 0.25
r23 <- 0.6
r24 <- 0.16
r34 <- 0.16
r <- vecTosm(c(r12, r13, r14, r23, r24, r34))
diag(r) <- 1
mix.r <- lapply(1:nrow(Geeganage2010), function(i){r})
attach(Geeganage2010)
## compute variance co-variances
computvcov <- mix.vcov(type = c("MD", "MD", "RD", "lgOR"),
d = cbind(MD_SBP, MD_DBP, NA, NA),
sdt = cbind(sdt_SBP, sdt_DBP, NA, NA),
sdc = cbind(sdc_SBP, sdc_DBP, NA, NA),
nt = cbind(nt_SBP, nt_DBP, nt_DD, nt_D),
nc = cbind(nc_SBP, nc_DBP, nc_DD, nc_D),
st = cbind(NA, NA, st_DD, st_D),
sc = cbind(NA, NA, sc_DD, sc_D),
r = mix.r,
name = c("MD.SBP", "MD.DBP", "RD.DD", "lgOR.D"))
# save different effect sizes in y
y <- computvcov$ef
head(y)
# save variances and covariances of all the effect sizes in a matrix S
S <- computvcov$matrix.vcov
S[1, ]
## fixed-effect model
MMA_FE <- summary(metafixed(y = y, Slist = computvcov$list.vcov))
#######################################################################
# Running random-effects model using package "mixmeta" or "metaSEM"
#######################################################################
# Restricted maximum likelihood (REML) estimator from the mixmeta package
# library(mixmeta)
# mvmeta_RE <- summary(mixmeta(cbind(MD.SBP, MD.DBP, RD.DD, lgOR.D) ~ 1,
# S = S, data = y, method = "reml"))
# mvmeta_RE
# maximum likelihood estimators from the metaSEM package
# library(metaSEM)
# metaSEM_RE <- summary(meta(y = y, v = S))
# metaSEM_RE
##############################################################
# Plotting the result:
##############################################################
obj <- MMA_FE
# obj <- mvmeta_RE
# obj <- metaSEM_RE
# pdf("CI.pdf", width = 4, height = 7)
plotCI(y = computvcov$ef, v = computvcov$list.vcov,
name.y = NULL, name.study = Geeganage2010$studyID,
y.all = obj$coefficients[,1],
y.all.se = obj$coefficients[,2],
hline = c(0, 0, 0, 1))
# dev.off()
Run the code above in your browser using DataLab