metafor
Package
metafor
package (Viechtbauer,
2010)
# install the metafor package
install.packages("metafor")
# load the metafor package
library(metafor)
# for the example data, we also need the bayesmeta package
install.packages("bayesmeta")
escalc()
allows for the calculation of 70+ different types of effect size
measures / outcome measures# copy the data to 'dat'
dat <- get(data(KarnerEtAl2014, package="bayesmeta"))
# examine the dataset (for a selection of variables)
dat[c(1,10:11,17:18)]
# study tiotropium.total tiotropium.exa placebo.total placebo.exa
# 1 Bateman2010a 1989 685 2002 842
# 2 Bateman2010b 1337 495 653 288
# 3 Beeh2006 1236 180 403 80
# 4 Brusasco2003 402 129 400 156
# 5 Casaburi2002 550 198 371 156
# 6 Chan2007 608 268 305 125
# 7 Cooper2010 260 112 259 102
# 8 Covelli2005 100 9 96 12
# 9 Dusser2006 500 213 510 272
# 10 Freeman2007 200 19 195 35
# 11 Johansson2008 107 2 117 4
# 12 Magnussen2008 228 13 244 26
# 13 Moita2008 147 6 164 6
# 14 NCT00144326 123 11 127 12
# 15 Niewoehner2005 914 255 915 296
# 16 Powrie2007 69 30 73 47
# 17 Sun2007 30 0 30 2
# 18 Tashkin2008 2987 2001 3006 2049
# 19 Tonnel2008 266 101 288 130
# 20 Trooster2011 238 11 219 24
# 21 Verkindre2006 46 10 54 8
# 22 Voshaar2008 360 43 181 21
# tiotropium.total - total number of patients in the treatment group
# tiotropium.exa - patients with ≥1 exacerbations in the treatment group
# placebo.total - total number of patients in the placebo group
# placebo.exa - patients with ≥1 exacerbations in the placebo group
# calculate log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=tiotropium.exa, n1i=tiotropium.total,
ci=placebo.exa, n2i=placebo.total,
slab=study, data=dat)
# examine the dataset
cbind(dat[1], "..."="...", dat[28:29])
# study ... yi vi
# 1 Bateman2010a ... -0.3234 0.0043
# 2 Bateman2010b ... -0.2943 0.0094
# 3 Beeh2006 ... -0.3737 0.0221
# 4 Brusasco2003 ... -0.3023 0.0219
# 5 Casaburi2002 ... -0.2546 0.0190
# 6 Chan2007 ... 0.1267 0.0202
# 7 Cooper2010 ... 0.1526 0.0319
# 8 Covelli2005 ... -0.3677 0.2173
# 9 Dusser2006 ... -0.4317 0.0161
# 10 Freeman2007 ... -0.7342 0.0930
# 11 Johansson2008 ... -0.6197 0.7684
# 12 Magnussen2008 ... -0.6793 0.1246
# 13 Moita2008 ... 0.1138 0.3468
# 14 NCT00144326 ... -0.0606 0.1919
# 15 Niewoehner2005 ... -0.2117 0.0104
# 16 Powrie2007 ... -0.8544 0.1187
# 17 Sun2007 ... -1.6773 2.4679
# 18 Tashkin2008 ... -0.0536 0.0030
# 19 Tonnel2008 ... -0.2958 0.0300
# 20 Trooster2011 ... -0.9321 0.1421
# 21 Verkindre2006 ... 0.4683 0.2745
# 22 Voshaar2008 ... 0.0329 0.0803
yi
– log odds ratios (negative values indicate lower
odds of exacerbations in the treatment group)vi
– corresponding sampling variances
rma()
is one of the main modeling functions for fitting equal-, fixed-,
random-, and mixed-effects meta-regression models# Random-Effects Model (k = 22; tau^2 estimator: REML)
#
# tau^2 (estimated amount of total heterogeneity): 0.0204 (SE = 0.0150)
# tau (square root of estimated tau^2 value): 0.1427
# I^2 (total heterogeneity / total variability): 48.61%
# H^2 (total variability / sampling variability): 1.95
#
# Test for Heterogeneity:
# Q(df = 21) = 42.4592, p-val = 0.0037
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# -0.2451 0.0544 -4.5076 21 0.0002 -0.3582 -0.1320 ***
#
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# pred ci.lb ci.ub pi.lb pi.ub
# 0.78 0.70 0.88 0.57 1.08
# forest plot
par(mar=c(5,4,2,2))
forest(res, xlim=c(-5,3.8), atransf=exp, at=log(2^(-4:3)),
addpred=TRUE, digits=c(2L,4L), header=TRUE, top=2, shade=TRUE)
# risk-of-bias ratings extracted from Figure 2 from Karner et al (2014)
rob <- read.table(header=TRUE, text = "
rb.a rb.b rb.c rb.d rb.e rb.f
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '-' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '-' '+'
'+' '+' '+' '+' '?' '?'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '?' '+' '?' '+' '+'
'+' '+' '+' '+' '-' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '+' '+'")
# turn the risk of bias items into factors with levels +, -, and ?
rob <- lapply(rob, factor, levels=c("+", "-", "?"))
# estimated average odds ratio (and 95% CI/PI)
pred <- predict(res, transf=exp, digits=2)
# need the rounded estimate and CI bounds further below
pred <- fmtx(c(pred$pred, pred$ci.lb, pred$ci.ub), 2)
# total number of studies
k <- nrow(dat)
# get the weights and format them as will be used in the forest plot
weights <- paste0(fmtx(weights(res), 1), "%")
# adjust the margins
par(mar=c(10.6,0,2.3,1.3), mgp=c(3,0.4,0), tcl=-0.3)
# forest plot with extra annotations
sav <- forest(res, atransf=exp, at=log(c(0.1, 0.2, 0.5, 1, 2, 5, 10)), xlim=c(-15,5),
xlab="", efac=c(0,1,1), textpos=c(-15,-2.5), lty=c(1,1,0), refline=NA,
ilab=cbind(tiotropium.exa, tiotropium.total, placebo.exa, placebo.total, weights),
ilab.xpos=c(-10.8,-9.6,-8.1,-6.9,-5.2), ilab.pos=2,
cex=0.79, header=c("Study","95% CI"), mlab="")
# add the horizontal line at the top
segments(sav$xlim[1], k+1, sav$xlim[2], k+1)
# add the vertical reference line at 0
segments(0, -3, 0, k+1)
# now we add a bunch of text; since some of the text falls outside of the
# plot region, we set xpd=NA so nothing gets clipped
par(xpd=NA)
# adjust cex as used in the forest plot and use a bold font
par(cex=sav$cex, font=2)
text(sav$ilab.xpos, k+2, pos=2, c("Events","Total","Events","Total","Weight"))
text(c(mean(sav$ilab.xpos[1:2])+0.2,mean(sav$ilab.xpos[3:4])), k+3, pos=2, c("tiotropium","placebo"))
text(sav$textpos[2], k+3, "Odds ratio", pos=2)
text(0, k+3, "Odds ratio")
text(sav$xlim[2]-0.35, k+3, "Risk of Bias", pos=2)
text(0, k+2, "IV, Random, 95% CI")
text(c(sav$xlim[1],sav$ilab.xpos[c(2,4,5)]), -1, pos=c(4,2,2,2,2),
c("Total (95% CI)", sum(dat$tiotropium.total), sum(dat$placebo.total), "100.0%"))
text(sav$xlim[1], -7, pos=4, "Risk of bias legend")
# first hide the non-bold summary estimate text and then add it back in bold font
rect(sav$textpos[2], -1.5, sav$ilab.xpos[5], -0.5, col="white", border=NA)
text(sav$textpos[2], -1, paste0(pred[1], " [", pred[2], ", ", pred[3], "]"), pos=2)
# use a non-bold font for the rest of the text
par(cex=sav$cex, font=1)
# add favours text below the x-axis
text(log(c(1,1)), -4.5, c("Favours tiotropium","Favours placebo"), pos=c(2,4), offset=0.5)
# add text for total events
text(sav$xlim[1], -2, pos=4, "Total events:")
text(sav$ilab.xpos[c(1,3)], -2, c(sum(dat$tiotropium.exa),sum(dat$placebo.exa)), pos=2)
# add text with heterogeneity statistics
text(sav$xlim[1], -3, pos=4, bquote(paste("Heterogeneity: ",
"Tau"^2, " = ", .(fmtx(res$tau2, 2)), "; ",
"Chi"^2, " = ", .(fmtx(res$QE, 2)),
", df = ", .(res$k - res$p),
" (", .(fmtp(res$QEp, digits=3, pname="P", add0=TRUE, sep=TRUE, equal=TRUE)), "); ",
I^2, " = ", .(round(res$I2)), "%")))
# add text for test of overall effect
text(sav$xlim[1], -4, pos=4, bquote(paste("Test for overall effect: Z = ",
.(fmtx(res$zval, 2)), " (", .(fmtp(res$pval, digits=5, pname="P", add0=TRUE, sep=TRUE, equal=TRUE)), ")")))
# add risk of bias points and symbols
cols <- c("#00cc00", "#cc0000", "#eeee00")
syms <- levels(rob$rb.a)
pos <- seq(sav$xlim[2]-2.4,sav$xlim[2]-0.5,length=6)
for (i in 1:6) {
points(rep(pos[i],k), k:1, pch=19, col=cols[rob[[i]]], cex=2.2)
text(pos[i], k:1, syms[rob[[i]]], font=2)
}
text(pos, k+2, c("A","B","C","D","E","F"), font=2)
# add risk of bias legend
text(sav$xlim[1], -8:-13, pos=4, c(
"(A) Random sequence generation (selection bias)",
"(B) Allocation concealment (selection bias)",
"(C) Blinding of participants and personnel (performance bias)",
"(D) Blinding of outcome assessment (detection bias)",
"(E) Incomplete outcome data (attrition bias)",
"(F) Selective reporting (reporting bias)"))
# funnel plot
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(2^(-5:5)), ylim=c(0,2),
las=1, digits=list(5L,1))
# Regression Test for Funnel Plot Asymmetry
#
# Model: mixed-effects meta-regression model
# Predictor: standard error
#
# Test for Funnel Plot Asymmetry: t = -1.1076, df = 20, p = 0.2812
# Limit Estimate (as sei -> 0): b = -0.1675 (CI: -0.3509, 0.0159)
# Rank Correlation Test for Funnel Plot Asymmetry
#
# Kendall's tau = -0.0130, p = 0.9556
# Estimated number of missing studies on the right side: 3 (SE = 3.1378)
#
# Random-Effects Model (k = 25; tau^2 estimator: REML)
#
# tau^2 (estimated amount of total heterogeneity): 0.0239 (SE = 0.0165)
# tau (square root of estimated tau^2 value): 0.1546
# I^2 (total heterogeneity / total variability): 49.75%
# H^2 (total variability / sampling variability): 1.99
#
# Test for Heterogeneity:
# Q(df = 24) = 50.0371, p-val = 0.0014
#
# Model Results:
#
# estimate se zval pval ci.lb ci.ub
# -0.2194 0.0534 -4.1089 <.0001 -0.3241 -0.1148 ***
#
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# funnel plot
par(mar=c(5,4,2,2))
funnel(taf, atransf=exp, at=log(2^(-5:5)), ylim=c(0,2),
las=1, digits=list(5L,1))
legend("topright", inset=.01, pch=c(19,21), pt.bg="gray90",
legend=c("Observed Studies", "Filled Studies"))
# contour-enhanced funnel plot (zoomed in at the top)
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(2^(-2:2)), xlim=c(-1,1), ylim=c(0,0.6),
las=1, digits=list(2L,1), refline=0,
level=c(90,95), shade=c("white","gray75"), legend=TRUE)
# Test of Excess Significance
#
# Observed Number of Significant Findings: 9 (out of 22)
# Expected Number of Significant Findings: 6.9080
# Observed Number / Expected Number: 1.3028
#
# Estimated Power of Tests (based on theta = -0.2451)
#
# min q1 median q3 max
# 0.0537 0.1042 0.2489 0.4509 0.8251
#
# Test of Excess Significance: p = 0.1683 (X^2 = 0.9235, df = 1)
# Limit Estimate (theta_lim): -0.2187 (where p = 0.1)
# fit step function selection model
sel <- selmodel(res, type="stepfun", alternative="less", steps=c(.025,.05,.5,1))
sel
# Random-Effects Model (k = 22; tau^2 estimator: ML)
#
# tau^2 (estimated amount of total heterogeneity): 0.0144 (SE = 0.0124)
# tau (square root of estimated tau^2 value): 0.1200
#
# Test for Heterogeneity:
# LRT(df = 1) = 6.3227, p-val = 0.0119
#
# Model Results:
#
# estimate se zval pval ci.lb ci.ub
# -0.1264 0.0887 -1.4260 0.1539 -0.3002 0.0473
#
# Test for Selection Model Parameters:
# LRT(df = 3) = 7.2109, p-val = 0.0655
#
# Selection Model Results:
#
# k estimate se zval pval ci.lb ci.ub
# 0 < p <= 0.025 9 1.0000 --- --- --- --- ---
# 0.025 < p <= 0.05 3 0.8065 0.5991 -0.3230 0.7467 0.0000 1.9807
# 0.05 < p <= 0.5 5 0.1468 0.1161 -7.3478 <.0001 0.0000 0.3744 ***
# 0.5 < p <= 1 5 0.2096 0.2216 -3.5675 0.0004 0.0000 0.6438 ***
#
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# meta-regression model using the mean FEV1 at baseline as predictor
res <- rma(yi, vi, mods = ~ baseline.fev1, data=dat, test="knha")
res
# Mixed-Effects Model (k = 22; tau^2 estimator: REML)
#
# tau^2 (estimated amount of residual heterogeneity): 0.0109 (SE = 0.0107)
# tau (square root of estimated tau^2 value): 0.1044
# I^2 (residual heterogeneity / unaccounted variability): 33.77%
# H^2 (unaccounted variability / sampling variability): 1.51
# R^2 (amount of heterogeneity accounted for): 46.40%
#
# Test for Residual Heterogeneity:
# QE(df = 20) = 29.6562, p-val = 0.0756
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 20) = 11.4396, p-val = 0.0030
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.6740 0.2728 2.4708 20 0.0226 0.1050 1.2430 *
# baseline.fev1 -0.7935 0.2346 -3.3822 20 0.0030 -1.2829 -0.3041 **
#
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# scatter/bubbeplot of the results
par(mar=c(5,4,2,2))
tmp <- regplot(res, xlab="Baseline Mean FEV1", bty="l", las=1, digits=1, legend=TRUE)
# label some points
lab <- c(11,17,18,20)
text(tmp$xi[lab], tmp$yi[lab], tmp$ids[lab], pos=c(1,3,4,1), offset=0.6)
# predicted pooled odds ratio at FEV1 equal to 1, 1.5, and 2
predict(res, newmods=c(1, 1.5, 2), transf=exp, digits=2)
# pred ci.lb ci.ub pi.lb pi.ub
# 1 0.89 0.79 0.99 0.69 1.13
# 2 0.60 0.49 0.72 0.45 0.80
# 3 0.40 0.26 0.61 0.25 0.65
rma.mh()
- meta-analysis via the Mantel-Haenszel methodrma.peto()
- meta-analysis via Peto’s one-step methodrma.glmm()
- meta-analysis via generalized linear models# tiotropium.total tiotropium.deaths placebo.total placebo.deaths
# 1 1989 52 2002 38
# 2 1337 34 653 10
# 3 1236 2 403 2
# 4 402 1 400 5
# 5 550 7 371 7
# 6 608 13 305 2
# 7 260 6 259 6
# 8 100 0 96 0
# 9 500 7 510 8
# 10 200 1 195 4
# 11 107 0 117 1
# 12 228 0 244 2
# 13 147 2 164 0
# 14 123 0 127 0
# 15 914 22 915 19
# 16 69 1 73 2
# 17 30 0 30 0
# 18 2987 381 3006 411
# 19 266 3 288 6
# 20 238 0 219 0
# 21 46 0 54 0
# 22 360 2 181 0
# fit a random-effects model using a 'binomial-normal model'
res <- rma.glmm(measure="OR", data=dat,
ai=tiotropium.deaths, n1i=tiotropium.total,
ci=placebo.deaths, n2i=placebo.total)
# Warning: 5 studies with NAs omitted from model fitting.
# Warning: Some yi/vi values are NA.
# Random-Effects Model (k = 17; tau^2 estimator: ML)
# Model Type: Unconditional Model with Fixed Study Effects
#
# tau^2 (estimated amount of total heterogeneity): 0
# tau (square root of estimated tau^2 value): 0
# I^2 (total heterogeneity / total variability): 0.00%
# H^2 (total variability / sampling variability): 1.00
#
# Tests for Heterogeneity:
# Wld(df = 16) = 14.6596, p-val = 0.5497
# LRT(df = 16) = 25.4335, p-val = 0.0625
#
# Model Results:
#
# estimate se zval pval ci.lb ci.ub
# -0.0234 0.0653 -0.3584 0.7200 -0.1513 0.1045
#
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# pred ci.lb ci.ub pi.lb pi.ub
# 0.98 0.86 1.11 0.86 1.11
rma.mv()
can be used to fit multilevel and multivariate modelsand \(\left[ \begin{array}{c} \varepsilon_{i1} \\ \varepsilon_{i2} \\ \vdots \end{array} \right] \sim N\left( \left[ \begin{array}{c} 0 \\ 0 \\ \vdots \end{array} \right], \left[ \begin{array}{ccc} v_{i1} & cov_{i12} & \ldots \\ cov_{i12} & v_{i2} & \ldots \\ \vdots & \vdots & \ddots \end{array} \right] \right)\)
# copy the data to 'dat' again
dat <- get(data(KarnerEtAl2014, package="bayesmeta"))
# select only relevant variables (also remove two studies)
dat <- dat[-c(17,20),c(1,10:12,17:19)]
# abbreviate the variable names
names(dat) <- abbreviate(names(dat))
# examine the dataset
dat
# stdy ttrpm.t ttrpm.x ttrpm.s plcb.t plcb.x plcb.s
# 1 Bateman2010a 1989 685 161 2002 842 198
# 2 Bateman2010b 1337 495 78 653 288 44
# 3 Beeh2006 1236 180 29 403 80 7
# 4 Brusasco2003 402 129 48 400 156 90
# 5 Casaburi2002 550 198 30 371 156 35
# 6 Chan2007 608 268 51 305 125 19
# 7 Cooper2010 260 112 21 259 102 16
# 8 Covelli2005 100 9 2 96 12 1
# 9 Dusser2006 500 213 28 510 272 33
# 10 Freeman2007 200 19 2 195 35 1
# 11 Johansson2008 107 2 0 117 4 0
# 12 Magnussen2008 228 13 4 244 26 0
# 13 Moita2008 147 6 1 164 6 1
# 14 NCT00144326 123 11 2 127 12 0
# 15 Niewoehner2005 914 255 64 915 296 87
# 16 Powrie2007 69 30 2 73 47 3
# 18 Tashkin2008 2987 2001 759 3006 2049 811
# 19 Tonnel2008 266 101 11 288 130 8
# 21 Verkindre2006 46 10 0 54 8 3
# 22 Voshaar2008 360 43 2 181 21 0
# ttrpm.t - total number of patients in the treatment group
# ttrpm.x - patients with ≥1 exacerbations in the treatment group
# ttrpm.s - patients with ≥1 severe exacerbations in the treatment group
# plcb.t - total number of patients in the placebo group
# plcb.x - patients with ≥1 exacerbations in the placebo group
# plcb.s - patients with ≥1 severe exacerbations in the placebo group
# calculate log odds ratios for exacerbations
dat <- escalc(measure="OR", ai=ttrpm.x, n1i=ttrpm.t,
ci=plcb.x, n2i=plcb.t,
data=dat, var.names=c("y1i","v1i"))
# calculate log odds ratios for severe exacerbations
dat <- escalc(measure="OR", ai=ttrpm.s, n1i=ttrpm.t,
ci=plcb.s, n2i=plcb.t,
data=dat, var.names=c("y2i","v2i"))
y1i
and y2i
are calculate based on
the same sample of subjects, they are not independentttrpm.s
is a subset of ttrpm.x
and
plcb.s
is a subset of plcb.x
, we can calculate
their covariances (cov12i
) without further information
(Trikalinos &
Olkin, 2012, see Appendix)# calculate the covariances
dat$cov12i <- with(dat, ttrpm.t / (ttrpm.x * (ttrpm.t-ttrpm.s)) +
plcb.t / (plcb.x * (plcb.t-plcb.s)))
# plot of the log odds ratios against each other (with 50% confidence ellipses)
# set up the plot
par(mar=c(5,4,2,2))
plot(NA, xlim=c(-1.7,1), ylim=c(-3.6,4), bty="l", xlab="Log(OR) for Exacerbations",
ylab="Log(OR) for Severe Exacerbations")
# add the confidence ellipses (requires that the 'ellipse' package is loaded)
invisible(sapply(1:nrow(dat), function(i) {
x <- dat[i,]
xy <- ellipse(matrix(c(x$v1i,x$cov12i,x$cov12i,x$v2i), nrow=2),
centre=c(x$y1i,x$y2i), level=0.5)
lines(xy[,1],xy[,2], col="gray80")
}))
# add the points
points(dat$y1i, dat$y2i, pch=21, bg="gray60", cex=1.5)
# reshape dataset into long format
dat.long <- data.frame(study = rep(dat$stdy, each=2),
outcome = c("exa","sev"),
yi = c(t(dat[c("y1i","y2i")])),
vi = c(t(dat[c("v1i","v2i")])))
# examine the first six rows of the reshaped dataset
head(dat.long)
# study outcome yi vi
# 1 Bateman2010a exa -0.3233776 0.004276443
# 2 Bateman2010a sev -0.2200787 0.012363055
# 3 Bateman2010b exa -0.2942854 0.009419799
# 4 Bateman2010b sev -0.1537356 0.037984103
# 5 Beeh2006 exa -0.3736609 0.022098500
# 6 Beeh2006 sev 0.3069067 0.180693654
# construct the 2x2 var-cov matrices of the studies
V <- tapply(dat, dat$stdy, function(x) matrix(c(x$v1i, x$cov12i,
x$cov12i, x$v2i), nrow=2))
# examine the first three elements
V[1:3]
# $Bateman2010a
# [,1] [,2]
# [1,] 0.004276443 0.00290643
# [2,] 0.002906430 0.01236305
#
# $Bateman2010b
# [,1] [,2]
# [1,] 0.009419799 0.00586845
# [2,] 0.005868450 0.03798410
#
# $Beeh2006
# [,1] [,2]
# [1,] 0.0220985 0.0184100
# [2,] 0.0184100 0.1806937
# fit a multivariate random-effects model
res <- rma.mv(yi, V, mods = ~ factor(outcome) - 1,
random = ~ outcome | study, struct="UN",
data=dat.long)
res
# Multivariate Meta-Analysis Model (k = 40; method: REML)
#
# Variance Components:
#
# outer factor: study (nlvls = 20)
# inner factor: outcome (nlvls = 2)
#
# estim sqrt k.lvl fixed level
# tau^2.1 0.0187 0.1369 20 no exa
# tau^2.2 0.0494 0.2222 20 no sev
#
# rho.exa rho.sev exa sev
# exa 1 - 20
# sev 0.4931 1 no -
#
# Test for Residual Heterogeneity:
# QE(df = 38) = 68.8813, p-val = 0.0016
#
# Test of Moderators (coefficients 1:2):
# QM(df = 2) = 20.2765, p-val < .0001
#
# Model Results:
#
# estimate se zval pval ci.lb ci.ub
# factor(outcome)exa -0.2311 0.0513 -4.5018 <.0001 -0.3317 -0.1305 ***
# factor(outcome)sev -0.1825 0.0923 -1.9782 0.0479 -0.3634 -0.0017 *
#
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit the model with cvvc="varcov" to get the var-cov matrix of the elements in G
res <- update(res, cvvc="varcov")
# regression of the log(OR) for severe exacerbations on the log(OR) for exacerbations
sav <- matreg(R=res$G, cov=TRUE, y=2, x=1, V=res$vvc, means=coef(res))
# set up the plot
par(mar=c(5,4,2,2))
plot(NA, xlim=c(-1.7,1), ylim=c(-3.6,4), bty="l", xlab="Log(OR) for Exacerbations",
ylab="Log(OR) for Severe Exacerbations")
# add the confidence ellipses for the studies
invisible(sapply(1:nrow(dat), function(i) {
x <- dat[i,]
xy <- ellipse(matrix(c(x$v1i,x$cov12i,x$cov12i,x$v2i), nrow=2),
centre=c(x$y1i,x$y2i), level=0.5)
lines(xy[,1],xy[,2], col="gray90")
}))
# add the points for the studies
points(dat$y1i, dat$y2i, pch=21, col="gray80", bg="gray90", cex=1)
# add the regression line
abline(a=coef(sav)[1], b=coef(sav)[2], lwd=3, col="gray60", lty="dashed")
# add the 95% PI ellipsis based on the model
xy <- ellipse(res$G, centre=coef(res), level=0.95)
lines(xy[,1],xy[,2], col="gray30", lwd=3, lty="dotted")
# add the 95% CI ellipsis based on the model
xy <- ellipse(vcov(res), centre=coef(res), level=0.95)
lines(xy[,1],xy[,2], col="gray30", lwd=3)
# add the point for the pooled effects
points(coef(res)[1], coef(res)[2], pch=21, bg="gray40", cex=2)
# add the legend
legend("topright", inset=0.01, lty=c("solid","dotted","dashed"), lwd=3,
col=c("gray30","gray30","gray60"),
legend=c("95% Confidence Ellipse", "95% Prediction Interval Ellipse", "Latent Regression Model"))
# Hypothesis:
# 1: factor(outcome)exa - factor(outcome)sev = 0
#
# Results:
# estimate se zval pval
# 1: -0.0486 0.0847 -0.5731 0.5666
#
# Test of Hypothesis:
# QM(df = 1) = 0.3285, p-val = 0.5666
metafor
provides a comprehensive collection of
functions for conducting meta-analyses in R