As published in Singh S, Loke YK, Spangler JG, Furberg CD, Risk of serious adverse cardiovascular events associated with varenicline: a systematic review and meta-analysis, Canadian Medical Association Journal, 2011;183;1359-1366.

(This is a screen shot, so it may not be perfectly sharp.)

PDF produced by R/metafor/forest:

After editing the PDF file usingInkscape:

This code, and only this code, produced the raw graph above.

require("metafor"

options(show.signif.stars=FALSE, stringsAsFactors=FALSE)

# Meta analysis of 14 double-blind RCTs reported in:

# Singh S, Loke YK, Spangler JG, Furberg CD, Risk of serious adverse

# cardiovascular events associated with varenicline: a systematic review

# and meta-analysis, Canadian Medical Association Journal, 2011;183;1359-1366.

#

# Varenicline is "Chantix, the best-selling prescription drug for

# smoking cessation."

Study <- c("A3051080 2010", "A3051095 2010", "Fagerstrom 2010",

"Gonzalez 2006", "Jorenby 2006", "Nakamura 2007", "Niaura 2008",

"Nides 2006", "Oncken 2006", "Rigotti 2010", "Tashkin 2010",

"Tonstad 2006", "Tsai 2007", "Williams 2007")

# V: Varenicline P: Placebo SACE: Severe Adverse Cardiac Event

n.V <- c(394, 493, 214, 352, 344, 465, 160, 383, 518, 355, 250, 603, 126, 251)

SACE.V <- c(1, 1, 0, 2, 1, 1, 2, 1, 2, 25, 5, 4, 1, 6)

e_n.V <- paste(SACE.V,"/",n.V, sep="")

n.P <- c(199, 166, 218, 344, 341, 154, 160, 127, 129, 359, 254, 607, 124, 126)

SACE.P <- c(0, 0, 1, 2, 1, 0, 0, 0, 0, 20, 2, 0, 0, 1)

e_n.P <- paste(SACE.P,"/",n.P, sep="")

(dframe <- data.frame(Study, n.V, SACE.V, n.P, SACE.P, e_n.V, e_n.P))

meta2 <- escalc(measure="PETO", data=dframe,

ai=SACE.V, n1i=n.V, ci=SACE.P, n2i=n.P, add=0)

dframe <- data.frame(dframe, logOR=meta2$yi, Var.logOR=meta2$vi)

orderES <- order(dframe$logOR)

dframe <- dframe[orderES,]

# check against article

OR.est <- round(exp(dframe$logOR),2)

CLim95Lo <- round(exp(dframe$logOR - 1.96*sqrt(dframe$Var.logOR)),2)

CLim95Hi <- round(exp(dframe$logOR + 1.96*sqrt(dframe$Var.logOR)),2)

data.frame(dframe, OR.est, CLim95Lo, CLim95Hi)

#This agrees with values in Figure 2.

(overall.2 <- rma(dframe$logOR, dframe$Var.logOR))

data.frame(estimate=round(exp(overall.2$b),2),

CLim95Lo=round(exp(overall.2$ci.lb),2),

CLim95Hi=round(exp(overall.2$ci.ub),2), row.names="Peto OR")

#This agrees with values in Figure 2.

graphics.off()

try(windows(width=7, height=4.5), silent=TRUE)

try(quartz(width=7, height=4.5), silent=TRUE)

forest(overall.2, slab=dframe$Study, efac=4, lwd=2,

xlim=c(-13,12), alim=c(-4.3,6.0), digits=1,

at=log(c(1/64, 1/32, 1/16, 1/8, 1/4, 1/2, 1, 2, 4, 8, 16, 32, 64, 128, 256)), atransf=exp,

ilab=cbind(dframe$e_n.V, dframe$e_n.P), xlab="Odds Ratio (Peto method, log-scaled)",

ilab.xpos=c(-8,-5.1), cex=0.70, )

text(1.7,14.7, "increased risk\nwithvarenicline", pos=3, cex=0.65)

arrows(x0=3.2, y0=16, x1=4, length=0.05)

op <- par(cex=0.70,font=2)

text(c(-8,-5.1),16,c("events/n","events/n"))

text(c(-8,-5.1),17,c("Varenicline","Placebo"))

text(-12,16,"Study",pos=4)

text(12,16,"Odds Ratio [95%CI]",pos=2)

par(op)

graphics.off()

try(windows(width=6, height=4), silent=TRUE)

try(quartz(width=6, height=4), silent=TRUE)

funnel(overall.2, main="Funnel Plot for Simple Random Effects Analysis")

abline(v=0, col="red")

# How could this be improved?