PUFA Meta-Analysis Reanalysis

First, we load the CVD events data from Excel.

Data Frame Legend:

ā€œPā€ refers to the polyunsaturated fat diet group

ā€œCā€ refers to the control diet group

pufaMETA<-read.csv("PUFA.csv", header=TRUE)
pufaMETA
##         Study P_Events P_Total C_Events C_Total
## 1       DARTS      132    1018      144    1015
## 2  LAVeterans       53     424       71     422
## 3 MinnesotaCS      131    4541      121    4516
## 4     MRC Soy       45     199       51     194
## 5        Oslo       61     206       81     206
## 6       STARS        2      27        5      28
### Here's our dataframe. Now we load "metafor." ###
library(metafor)
## Loading required package: Matrix
## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

Next, we produce an effect estimate for each study.

Along with the corresponding variance.

dat <- escalc(measure="RR", 
              ai=P_Events, n1i=P_Total, 
              ci=C_Events, n2i=C_Total, 
              data=pufaMETA)
dat
##         Study P_Events P_Total C_Events C_Total      yi     vi
## 1       DARTS      132    1018      144    1015 -0.0900 0.0126
## 2  LAVeterans       53     424       71     422 -0.2971 0.0282
## 3 MinnesotaCS      131    4541      121    4516  0.0739 0.0155
## 4     MRC Soy       45     199       51     194 -0.1506 0.0317
## 5        Oslo       61     206       81     206 -0.2836 0.0190
## 6       STARS        2      27        5      28 -0.8799 0.6272
res <- rma(yi, vi, data=dat)
confint(res)
## 
##        estimate  ci.lb   ci.ub
## tau^2    0.0066 0.0000  0.3010
## tau      0.0812 0.0000  0.5486
## I^2(%)  21.4317 0.0000 92.5694
## H^2      1.2728 1.0000 13.4578
par(mar=c(4,4,1,2))

Then, we pool the effects to produce a summary effect.

We use a random-effects model, as done originally.

res <- rma(ai=P_Events, n1i=P_Total, 
           ci=C_Events, n2i=C_Total, 
           data=dat, measure="RR",
           slab=paste(Study, sep=", "), 
           method="REML")

We plot the results and their corresponding confidence intervals.

forest(res, xlim=c(-16, 6), 
       at=log(c(0.05, 0.25, 1, 2)), atransf=exp,
       ilab=cbind(dat$P_Events, dat$P_Total, 
                  dat$C_Events, dat$C_Total),
       ilab.xpos=c(-9.5,-8,-6,-4.5), 
       cex=0.85, ylim=c(-1, 8.5),
       xlab="Risk Ratio", mlab="", psize=1.7, pch=16)

##We also paste the heterogeneity statistics.## 

text(-16, -1, pos=4, cex=0.85, 
     bquote(paste("Random Effects (Q = ",
      .(formatC(res$QE, digits=2, format="f")),
      ", df = ", .(res$k - res$p),
      ", p = ", .(formatC(res$QEp, digits=2, 
      format="f")), "; ", I^2, " = ",
     .(formatC(res$I2, digits=1, format="f")),
     "%)")))

op <- par(cex=0.85, font=4)

par(font=2)

##Then we label the columns in the forest plot.## 

text(c(-9.5,-8,-6,-4.5), 7.5, c("Events ", " Total", 
                                "Events ", " Total"))
text(c(-8.75,-5.25),     8.5, c("PUFA Diet", 
                                "Control Diet"))
text(-16,                7.5, "Study Name",  pos=4)
text(6,                  7.5, "Risk Ratio [95% CI]", 
     pos=2)

These are the results when excluding the Finnish studies.