Background
These R code chunks are related to the following research article which has been published recently. The article can be found via the following link: manuscript link. As the article is sufficiently detailed about the objective and results of the study, this will only serve as a tutorial to perform the analysis and draw the graphs included in the article. Authors believe in reproducible research and this is their way of contributing to the current knowledge. If this code helps you in any way we kindly ask you to cite the aforementioned article (that’s the best way of saying thanks). This tutorial will be focused on the realization of the forest plots (including subgroup analysis).
Prerequisites
We first load the corresponding data we will be using for the analysis, then we load required packages.
BWGF <- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/glutamine/excels/csvs/BWGF.csv")
GPX <- read.csv("D:/Chris/ETUDES1/PhD thesis/R files/DATABASE (excel)/glutamine/excels/csvs/GPX.csv")
head(BWGF)## Studies trt.n trt trt.sd con.n con con.sd conc dur hsper
## 1 Ayazi 2014a 50 2001.04 197.2 50 1975.25 197.2 0.25 42 42
## 2 Ayazi 2014b 50 2179.41 197.2 50 1975.25 197.2 0.50 42 42
## 3 Ayazi 2014c 50 2053.05 197.2 50 1975.25 197.2 1.00 42 42
## 4 Bai et al. 2019a 60 804.60 117.4 60 776.30 117.4 0.50 14 14
## 5 Bai et al. 2019b 60 852.30 117.4 60 776.30 117.4 1.00 14 14
## 6 Dai et al. 2011 18 1487.00 813.7 18 1087.00 813.7 0.50 42 21head(GPX)## Studies trt.n trt trt.sd con.n con con.sd conc dur hsper
## 1 Sifa et al. 2018a1 6 12.36 5.641175 6 6.18 5.641175 0.5 42 1
## 2 Sifa et al. 2018b1 6 13.61 5.641175 6 6.18 5.641175 1.5 42 1
## 3 Sifa et al. 2018c1 6 14.08 5.641175 6 6.18 5.641175 2.0 42 1
## 4 Hu et al. 2020 a1 18 16.19 3.568061 18 11.78 3.568061 0.5 21 21
## 5 Hu et al. 2020 a2 18 15.96 3.568061 18 11.78 3.568061 1.0 21 21
## 6 Hu et al. 2020 a3 18 18.46 3.568061 18 11.78 3.568061 1.5 21 21
## organ
## 1 Breast muscle
## 2 Breast muscle
## 3 Breast muscle
## 4 Breast muscle
## 5 Breast muscle
## 6 Breast muscleThere are two categories of data sets.The growth parameters (an example is BWGF), and antioxidant parameters (an example is GPX). we used the head function to show you the general shape of the data sets. Due to confidentiality reason, the original data sets will be made available upon request. Here is the description of the features:
Studies: name of each studies included;
trt.n: sample size of the treatment group;
trt: mean of the treatment group;
trt.sd: standard deviation of the treatment group;
con.n: sample size of the control group;
con: mean of the control group;
con.sd: standard deviation of the control group;
dur: duration of heat stress (used as covariate);
conc: dietary glutamine concentration;
organ: category of tissue;
As you might see, we will use the “metafor” library to pool effect sizes, perform subgroup analysis and draw forest plots.
Growth parameters forest plots
The growth parameters data was pooled without performing any subgroup analysis. Thus, the forest plots are drawn relatively easily. PS: prior analysis such as sensitivity and baujat plot were performed to determine outlying studies. We will have a separate tutorial discussing this approach. for the time being we assume that the reader is just interested in drawing forest plots after pooling effect sizes. Only the first code chunk will be detailed step by step.
TR<- escalc(m1i = trt, m2i= con,sd1i = trt.sd, sd2i = con.sd, n1i=trt.n,n2i = con.n,
data= BWGF,measure = "SMD") #we calculate effect size
TR %>%
as_tibble() %>%
head() ## then we print the results obtained## # A tibble: 6 x 12
## Studies trt.n trt trt.sd con.n con con.sd conc dur hsper yi vi
## <chr> <int> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int> <int> <dbl> <dbl>
## 1 Ayazi 20~ 50 2001. 197. 50 1975. 197. 0.25 42 42 0.130 0.0401
## 2 Ayazi 20~ 50 2179. 197. 50 1975. 197. 0.5 42 42 1.03 0.0453
## 3 Ayazi 20~ 50 2053. 197. 50 1975. 197. 1 42 42 0.391 0.0408
## 4 Bai et a~ 60 805. 117. 60 776. 117. 0.5 14 14 0.240 0.0336
## 5 Bai et a~ 60 852. 117. 60 776. 117. 1 14 14 0.643 0.0351
## 6 Dai et a~ 18 1487 814. 18 1087 814. 0.5 42 21 0.481 0.114New columns yi and vi were added to the original dataset.These are the effect sizes (yi) and their standard error (vi). We will use them to perform the meta-analysis (pool effect sizes) as well as to draw the forest plots. For details about technical terms such as random effects or fixed effects please refer to the manuscript.
res <-rma(yi,vi,method="REML",data = TR,slab = Studies) ## pooling effect sizes (ES) using random effects
res##
## Random-Effects Model (k = 17; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.1085 (SE = 0.0608)
## tau (square root of estimated tau^2 value): 0.3294
## I^2 (total heterogeneity / total variability): 67.25%
## H^2 (total variability / sampling variability): 3.05
##
## Test for Heterogeneity:
## Q(df = 16) = 50.0152, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub <U+200B>
## 0.7029 0.1023 6.8724 <.0001 0.5024 0.9033 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1You can now see the results we will use to draw the forest plot
forest(res,xlim=c(-6, 6), header = "Studies", )
### add text with Q-value, dfs, p-value, and I^2 statistic
text(-0.30, -1, pos=2, cex=0.75, bquote(paste(" for All Studies (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")), "%)")))
Antioxidant parameters forest plots
Since subgroup analysis was used in this case, we will have to perform the analysis and then include sections in the forest plot for each subgroups. While the code will get longer, the reasoning stays the same.
TR5<- escalc(m1i = trt, m2i= con,sd1i = trt.sd, sd2i = con.sd, n1i=trt.n,n2i = con.n,
data=GPX,measure = "SMD")
TR5 %>%
as_tibble() %>%
head()## # A tibble: 6 x 13
## Studies trt.n trt trt.sd con.n con con.sd conc dur hsper organ yi
## <chr> <int> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <int> <int> <chr> <dbl>
## 1 Sifa et a~ 6 12.4 5.64 6 6.18 5.64 0.5 42 1 Brea~ 1.01
## 2 Sifa et a~ 6 13.6 5.64 6 6.18 5.64 1.5 42 1 Brea~ 1.22
## 3 Sifa et a~ 6 14.1 5.64 6 6.18 5.64 2 42 1 Brea~ 1.29
## 4 Hu et al.~ 18 16.2 3.57 18 11.8 3.57 0.5 21 21 Brea~ 1.21
## 5 Hu et al.~ 18 16.0 3.57 18 11.8 3.57 1 21 21 Brea~ 1.15
## 6 Hu et al.~ 18 18.5 3.57 18 11.8 3.57 1.5 21 21 Brea~ 1.83
## # ... with 1 more variable: vi <dbl>res<-rma(yi,vi,method="REML",data = TR5,slab = Studies)
res##
## Random-Effects Model (k = 18; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.4386 (SE = 0.2395)
## tau (square root of estimated tau^2 value): 0.6623
## I^2 (total heterogeneity / total variability): 65.07%
## H^2 (total variability / sampling variability): 2.86
##
## Test for Heterogeneity:
## Q(df = 17) = 44.4697, p-val = 0.0003
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub <U+200B>
## 1.1219 0.1982 5.6620 <.0001 0.7336 1.5103 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Next we perform the subgroup analysis
res1 <- rma(yi, vi, mods = ~ organ, data=TR5, method = "REML") #test for significant difference in the ES between subgroup
##actually fitting subgroups##
res2 <- rma(yi, vi, subset=(organ=="Thigh muscle"), data=TR5, method = "REML")
res3 <- rma(yi, vi, subset=(organ=="Serum"), data=TR5, method = "REML")
res4 <- rma(yi, vi, subset=(organ=="Breast muscle"), data=TR5, method = "REML")Now comes the tricky part, we draw the forest plot step by step.
### separate subgroups in forest plot now ###
forest(res,xlim=c(-6, 6), ylim = c(-2,34), order = organ, rows=c(3:8,13:18,23:28), header = "Studies", )
##add text on the graph##
### add text with Q-value, dfs, p-value, and I^2 statistic for the whole model##
text(-1.4, -1, pos=2, cex=0.75, bquote(paste(" for All Studies (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")), "%)")))
### add text for the test of subgroup differences
text(-1.4, -2, pos=2, cex=0.75, bquote(paste("Test for Subgroup Differences: ",
Q[M], " = ", .(formatC(res1$QM, digits=2, format="f")), ", df = ", .(res$p - 1),
", p = ", .(formatC(res1$QMp, digits=2, format="f")))))
### add text for the headers of subgroups####
op <- par(cex=1.2, font=4)
text(-6, c(10,20,30), cex=1, pos=4, c("Breast Muscle",
"Serum", "Thigh muscle"))
par(op)
### add text with Q-value, dfs, p-value, and I^2 statistic for the subgroups ###
text(-6, 21.5, pos=4, cex=1, bquote(paste("RE Model for Subgroup (Q = ",
.(formatC(res2$QE, digits=2, format="f")), ", df = ", .(res2$k - res2$p),
", p = ", .(formatC(res2$QEp, digits=2, format="f")), "; ", I^2, " = ",
.(formatC(res2$I2, digits=1, format="f")), "%)")))
text(-6, 11.5, pos=4, cex=1, bquote(paste("RE Model for Subgroup (Q = ",
.(formatC(res3$QE, digits=2, format="f")), ", df = ", .(res3$k - res3$p),
", p = ", .(formatC(res3$QEp, digits=2, format="f")), "; ", I^2, " = ",
.(formatC(res3$I2, digits=1, format="f")), "%)")))
text(-6, 1.5, pos=4, cex=1, bquote(paste("RE Model for Subgroup (Q = ",
.(formatC(res4$QE, digits=2, format="f")), ", df = ", .(res4$k - res4$p),
", p = ", .(formatC(res4$QEp, digits=2, format="f")), "; ", I^2, " = ",
.(formatC(res4$I2, digits=1, format="f")), "%)")))
##add polygons to the subgroups##
addpoly(res2, row=21.5, cex=1, mlab="")
addpoly(res3, row= 11.5, cex=1, mlab="")
addpoly(res4, row= 1.5, cex=1, mlab="")
By using the same approaches all remaining forest plots of the manuscript can be plotted.