In martinlaw/dcv: Meta-Analysis on Postoperative Pneumonia

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(metafor)
library(xtable)

initial.data <- read.csv(file="../data/secondary-outcomes.csv")
dat <- escalc(measure="MD", data=initial.data, m1i=m1i, sd1i=sd1i, n1i=n1i, m2i=m2i, sd2i=sd2i, n2i=n2i, slab=study)

Create models

# Fig 2: Respiratory compliance (Cdyn)
cdyn.olv <- rma(yi, vi, data=dat, subset = (outcome1=="Cdyn" & time=="OLV"), method="PM")

# Fig 3: Plateau pressure (Plat)
plat.tlv <- rma(yi, vi, data=dat, subset = (outcome1=="Plat" & time=="TLV"), method="PM")
plat.olv <- rma(yi, vi, data=dat, subset = (outcome1=="Plat" & time=="OLV"), method="PM")

# Fig 4: Peak inspiratory pressure (Peak)
peak.tlv <- rma(yi, vi, data=dat, subset = (outcome1=="Peak" & time=="TLV"), method="PM")
peak.olv <- rma(yi, vi, data=dat, subset = (outcome1=="Peak" & time=="OLV"), method="PM")

# Fig 5: Arterial oxygen pressure (pO2)
po2.tlv <- rma(yi, vi, data=dat, subset = (outcome1=="pO2" & time=="TLV"), method="PM")
# po2.olv <- rma(yi, vi, data=dat, subset = (outcome1=="pO2" & time=="OLV"), method="PM")
# Error when using default options. Fix:
po2.olv <- rma(yi, vi, data=dat, subset = outcome1=="pO2" & time=="OLV", method="PM", control=list(tau2.max=1000))

Collated results

est.ci <- function(mod){
  est.with.ci <- paste(round(mod$b, 1), " (", round(mod$ci.lb, 1), ", ", round(mod$ci.ub, 1), ")", sep="")
  est.with.ci
}
olv.est <- c(est.ci(cdyn.olv), est.ci(plat.olv), est.ci(peak.olv), est.ci(po2.olv))
tlv.est <- c(est.ci(plat.tlv), est.ci(peak.tlv), est.ci(po2.tlv))

olv.pval <- round(c(cdyn.olv$pval, plat.olv$pval, peak.olv$pval, po2.olv$pval), 3)
tlv.pval <- round(c(plat.tlv$pval, peak.tlv$pval, po2.tlv$pval), 3)

olv.I2 <- round(c(cdyn.olv$I2, plat.olv$I2, peak.olv$I2, po2.olv$I2), 1)
tlv.I2 <- round(c(plat.tlv$I2, peak.tlv$I2, po2.tlv$I2), 1)

olv.results <- data.frame(est=olv.est, p.value=olv.pval, I2=olv.I2)
rownames(olv.results) <- c("Cdyn", "Plateau", "Peak", "PaO2")

tlv.results <- data.frame(est=tlv.est, p.value=tlv.pval, I2=tlv.I2)
rownames(tlv.results) <- c("Plateau", "Peak", "PaO2")

olv.tab <- xtable(olv.results, align="rrrr", digits=c(0,1,3,1), caption="OLV results")
tlv.tab <- xtable(tlv.results, align="rrrr", digits=c(0,1,3,1), caption="TLV results")

print(olv.tab, type="html")
print(tlv.tab, type="html")

Evidence of effect in all secondary outcomes except Plateau+TLV. These effects persist when the HKSJ modification is used, with the exception of Plateau+OLV (results not shown).

Notes on secondary outcome figures in draft

Figure 3 in paper has different headings compared to Figs 4 and 5.
Fig 4 has TLV first, while Figs 3 and 5 have OLV first.
TLV and OLV data have been mixed up in the document (assuming the Excel sheet is correct)

Forest plots

Fig 2: Respiratory compliance (Cdyn)

forest(cdyn.olv, header="Cdyn (OLV)", mlab="")
text(x=-16, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(cdyn.olv$I2, digits=1, format="f")))))

Fig 3: Plateau pressure (Plat)

forest(plat.tlv, header="Plateau (TLV)", mlab="")
text(x=-20, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(plat.tlv$I2, digits=1, format="f")))))
forest(plat.olv, header="Plateau (OLV)", mlab="")
text(x=-20, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(plat.olv$I2, digits=1, format="f")))))

Fig 4: Peak inspiratory pressure (Peak)

forest(peak.tlv, header="Peak (TLV)", mlab="")
text(x=-18, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(peak.tlv$I2, digits=1, format="f")))))
forest(peak.olv, header="Peak (OLV)", mlab="")
text(x=-22, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(peak.olv$I2, digits=1, format="f")))))

Fig 5: Arterial oxygen pressure (pO2)

forest(po2.tlv, header="pO2 (TLV)", mlab="")
text(x=-200, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(po2.tlv$I2, digits=1, format="f")))))
forest(po2.olv, header="pO2 (OLV)", mlab="")
text(x=-200, y=-1, pos=4, labels=bquote(paste(I^2, "=", .(formatC(po2.olv$I2, digits=1, format="f")))))

Attempt to recreate paper results

# Fig 2: Respiratory compliance (Cdyn)
cdyn.olv <- rma(yi, vi, data=dat, subset = (outcome1=="Cdyn" & time=="OLV"), method="REML")

# Fig 3: Plateau pressure (Plat)
plat.tlv <- rma(yi, vi, data=dat, subset = (outcome1=="Plat" & time=="TLV"), method="REML")
plat.olv <- rma(yi, vi, data=dat, subset = (outcome1=="Plat" & time=="OLV"), method="REML")

# Fig 4: Peak inspiratory pressure (Peak)
peak.tlv <- rma(yi, vi, data=dat, subset = (outcome1=="Peak" & time=="TLV"), method="REML")
peak.olv <- rma(yi, vi, data=dat, subset = (outcome1=="Peak" & time=="OLV"), method="REML")

# Fig 5: Arterial oxygen pressure (pO2)
po2.tlv <- rma(yi, vi, data=dat, subset = (outcome1=="pO2" & time=="TLV"), method="REML")
# po2.olv <- rma(yi, vi, data=dat, subset = (outcome1=="pO2" & time=="OLV"), method="PM")
# Error when using default options. Fix:
po2.olv <- rma(yi, vi, data=dat, subset = outcome1=="pO2" & time=="OLV", control=list(tau2.max=1000), method="REML")

Collated results

est.ci <- function(mod){
  est.with.ci <- paste(round(mod$b, 1), " (", round(mod$ci.lb, 1), ", ", round(mod$ci.ub, 1), ")", sep="")
  est.with.ci
}
olv.est <- c(est.ci(cdyn.olv), est.ci(plat.olv), est.ci(peak.olv), est.ci(po2.olv))
tlv.est <- c(est.ci(plat.tlv), est.ci(peak.tlv), est.ci(po2.tlv))

olv.pval <- round(c(cdyn.olv$pval, plat.olv$pval, peak.olv$pval, po2.olv$pval), 3)
tlv.pval <- round(c(plat.tlv$pval, peak.tlv$pval, po2.tlv$pval), 3)

olv.I2 <- round(c(cdyn.olv$I2, plat.olv$I2, peak.olv$I2, po2.olv$I2), 1)
tlv.I2 <- round(c(plat.tlv$I2, peak.tlv$I2, po2.tlv$I2), 1)

olv.results <- data.frame(est=olv.est, p.value=olv.pval, I2=olv.I2)
rownames(olv.results) <- c("Cdyn", "Plateau", "Peak", "PaO2")

tlv.results <- data.frame(est=tlv.est, p.value=tlv.pval, I2=tlv.I2)
rownames(tlv.results) <- c("Plateau", "Peak", "PaO2")

olv.tab <- xtable(olv.results, align="rrrr", digits=c(0,1,3,1), caption="OLV results")
tlv.tab <- xtable(tlv.results, align="rrrr", digits=c(0,1,3,1), caption="TLV results")

print(olv.tab, type="html")
print(tlv.tab, type="html")