It's interesting to do meta analysis recently. Once you master these skills, you don't have to worry about publishing articles in the future. Meta analysis is also one of the very important methods in medicine. This issue analyzes several classic graphics based on existing data, including forest, vulnerability, radial, labe And QQ normal.
01 concept of meta # analysis
In 1976, Glass named this concept meta analysis for the first time. It is defined as a method of collecting, merging and statistical analysis of different research results. This method has gradually developed into a new discipline  "evidencebased medicine", the main content and research means. The main purpose of Metaanalysis is to comprehensively reflect the previous research results more objectively. Researchers do not conduct the original research, but comprehensively analyze the results obtained from the research.
02 input data mode
We use R package metafor, developed by Viechtbauer, which can not only complete the Metaanalysis of binary classification and continuity variables, but also carry out regression analysis. This software package is the only package that can carry out the fitting operation of mixed effect model, including:

Single variable

Multi classification variable

Continuity variable
It can also detect the model coefficients and obtain the confidence interval, and accurately test the parameters, such as permutation test.
Package installation, as follows:
install.packages("metafor") library("metafor")
The software package data were obtained from 13 studies that examined BCG's anti tuberculosis effects. The purpose of meta analysis was to examine the overall effectiveness of BCG in preventing TB and to examine possible mitigation agents. This data set has been used in some publications to illustrate the metaanalysis method, as follows:
data("dat.bcg") print(dat.bcg, row.names = FALSE) trial author year tpos tneg cpos cneg ablat alloc 1 Aronson 1948 4 119 11 128 44 random 2 Ferguson & Simes 1949 6 300 29 274 55 random 3 Rosenthal et al 1960 3 228 11 209 42 random 4 Hart & Sutherland 1977 62 13536 248 12619 52 random 5 FrimodtMoller et al 1973 33 5036 47 5761 13 alternate 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate 7 Vandiviere et al 1973 8 2537 10 619 19 random 8 TPT Madras 1980 505 87886 499 87892 13 random 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic 11 Comstock et al 1974 186 50448 141 27197 18 systematic 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic 13 Comstock et al 1976 27 16886 29 17825 33 systematic #######Data description ?dat.bcg Format The data frame contains the following columns: trial numeric trial number author character author(s) year numeric publication year tpos numeric number of TB positive cases in the treated (vaccinated) group tneg numeric number of TB negative cases in the treated (vaccinated) group cpos numeric number of TB positive cases in the control (nonvaccinated) group cneg numeric number of TB negative cases in the control (nonvaccinated) group ablat numeric absolute latitude of the study location (in degrees) alloc character method of treatment allocation (random, alternate, or systematic assignment)
The 13 studies provide data details in the following tables:
TB positive  TB negative  
vaccinated group  tpos  tneg 
control group  cpos  cneg 
03 example analysis
Calculate the log relative risk and the corresponding sampling variance as follows:
dat < escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, append = TRUE) print(dat[,c(4:7)], row.names = FALSE trial author year ablat alloc yi vi 1 Aronson 1948 44 random 0.8893 0.3256 2 Ferguson & Simes 1949 55 random 1.5854 0.1946 3 Rosenthal et al 1960 42 random 1.3481 0.4154 4 Hart & Sutherland 1977 52 random 1.4416 0.0200 5 FrimodtMoller et al 1973 13 alternate 0.2175 0.0512 6 Stein & Aronson 1953 44 alternate 0.7861 0.0069 7 Vandiviere et al 1973 19 random 1.6209 0.2230 8 TPT Madras 1980 13 random 0.0120 0.0040 9 Coetzee & Berjak 1968 27 random 0.4694 0.0564 10 Rosenthal et al 1961 42 systematic 1.3713 0.0730 11 Comstock et al 1974 18 systematic 0.3394 0.0124 12 Comstock & Webster 1969 33 systematic 0.4459 0.5325 13 Comstock et al 1976 33 systematic 0.0173 0.0714
The demo formula interface is as follows:
k < length(dat.bcg$trial) dat.fm < data.frame(study = factor(rep(1:k, each = 4))) dat.fm$grp < factor(rep(c("T", "T", "C", "C"), k), levels = c("T", "C")) dat.fm$out < factor(rep(c("+", "", "+", ""), k), levels = c("+", "")) dat.fm$freq < with(dat.bcg, c(rbind(tpos, tneg, cpos, cneg))) dat.fm study grp out freq 1 1 T + 4 2 1 T  119 3 1 C + 11 4 1 C  128 5 2 T + 6 6 2 T  300 7 2 C + 29 8 2 C  274 9 3 T + 3 10 3 T  228 11 3 C + 11 12 3 C  209 13 4 T + 62 14 4 T  13536 15 4 C + 248 16 4 C  12619 17 5 T + 33 18 5 T  5036 19 5 C + 47 20 5 C  5761 21 6 T + 180 22 6 T  1361 23 6 C + 372 24 6 C  1079 25 7 T + 8 26 7 T  2537 27 7 C + 10 28 7 C  619 29 8 T + 505 30 8 T  87886 31 8 C + 499 32 8 C  87892 33 9 T + 29 34 9 T  7470 35 9 C + 45 36 9 C  7232 37 10 T + 17 38 10 T  1699 39 10 C + 65 40 10 C  1600 41 11 T + 186 42 11 T  50448 43 11 C + 141 44 11 C  27197 45 12 T + 5 46 12 T  2493 47 12 C + 3 48 12 C  2338 49 13 T + 27 50 13 T  16886 51 13 C + 29 52 13 C  17825
 The random effect model is as follows:
#The default model is REML model random effect model using restricted maximum likelihood estimator res < rma(yi, vi, data = dat) res RandomEffects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664) tau (square root of estimated tau^2 value): 0.5597 I^2 (total heterogeneity / total variability): 92.22% H^2 (total variability / sampling variability): 12.86 Test for Heterogeneity: Q(df = 12) = 152.2330, pval < .0001 Model Results: estimate se zval pval ci.lb ci.ub 0.7145 0.1798 3.9744 <.0001 1.0669 0.3622 ***  Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 res < rma(ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat, measure = "RR") res RandomEffects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664) tau (square root of estimated tau^2 value): 0.5597 I^2 (total heterogeneity / total variability): 92.22% H^2 (total variability / sampling variability): 12.86 Test for Heterogeneity: Q(df = 12) = 152.2330, pval < .0001 Model Results: estimate se zval pval ci.lb ci.ub 0.7145 0.1798 3.9744 <.0001 1.0669 0.3622 ***  Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For the confidence interval of heterogeneity, the confidence interval of heterogeneity coefficient (tau2,tau,I2(%),H^2) can be calculated by confint () in metafor, as follows:
### confidence intervals for tau^2, I^2, and H^2 confint(res) estimate ci.lb ci.ub tau^2 0.3132 0.1197 1.1115 tau 0.5597 0.3460 1.0543 I^2(%) 92.2214 81.9177 97.6781 H^2 12.8558 5.5303 43.0680 # Prediction / confidence interval of estimated effect quantity predict(res, digits = 3) pred se ci.lb ci.ub pi.lb pi.ub 0.715 0.180 1.067 0.362 1.867 0.438
 forest map
Forest map is based on statistical indicators and statistical analysis methods and drawn with numerical calculation results. In the plane rectangular coordinate system, it takes a vertical invalid line (the abscissa scale is 1 or 0) as the center, uses multiple line segments parallel to the horizontal axis to describe the effect quantity and confidence interval of each research, and uses a prism (or other graphics) to describe the effect quantity and confidence interval of multiple research combinations. It is very simple and intuitive to describe the statistical results of Metaanalysis. It is the most commonly used result expression in Metaanalysis. As follows:
### forst plot forest(res, slab = paste(dat$author, dat$year, sep = ", "), xlim = c(16, 6), at = log(c(.05, .25, 1, 4)), atransf = exp, ilab = cbind(dat$tpos, dat$tneg, dat$cpos, dat$cneg), ilab.xpos = c(9.5, 8, 6, 4.5), cex = .75) op < par(cex = .75, font = 2) text(c(9.5, 8, 6, 4.5), 15, c("TB+", "TB", "TB+", "TB")) text(c(8.75, 5.25), 16, c("Vaccinated", "Control")) text(16, 15, "Author(s) and Year", pos = 4) text(6, 15, "Relative Risk [95% CI]", pos = 2) par(op)
The mixed effect model is as follows:
### mixedeffects model with absolute latitude and publication year as moderators res < rma(yi, vi, mods = cbind(ablat, year), data = dat) res MixedEffects Model (k = 13; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.1108 (SE = 0.0845) tau (square root of estimated tau^2 value): 0.3328 I^2 (residual heterogeneity / unaccounted variability): 71.98% H^2 (unaccounted variability / sampling variability): 3.57 R^2 (amount of heterogeneity accounted for): 64.63% Test for Residual Heterogeneity: QE(df = 10) = 28.3251, pval = 0.0016 Test of Moderators (coefficients 2:3): QM(df = 2) = 12.2043, pval = 0.0022 Model Results: estimate se zval pval ci.lb ci.ub intrcpt 3.5455 29.0959 0.1219 0.9030 60.5724 53.4814 ablat 0.0280 0.0102 2.7371 0.0062 0.0481 0.0080 ** year 0.0019 0.0147 0.1299 0.8966 0.0269 0.0307  Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ### predicted average relative risks for 10, 20, 30, 40, and 50 degrees ### absolute latitude holding publication year constant at 1970 predict(res, newmods = cbind(seq(from = 10, to = 60, by = 10), 1970), transf = exp, addx = TRUE) ### plot of the predicted average relative risk as a function of absolute latitude preds < predict(res, newmods = cbind(0:60, 1970), transf = exp) wi < 1/sqrt(dat$vi) size < 0.5 + 3 * (wi  min(wi))/(max(wi)  min(wi)) plot(dat$ablat, exp(dat$yi), pch = 19, cex = size, xlab = "Absolute Latitude", ylab = "Relative Risk", las = 1, bty = "l", log = "y") lines(0:60, preds$pred) lines(0:60, preds$ci.lb, lty = "dashed") lines(0:60, preds$ci.ub, lty = "dashed") abline(h = 1, lty = "dotted")
 funnel diagram
Funnel chart is a scatter chart proposed by Light et al. In 1984, which generally takes the effect quantity of a single study as the abscissa and the sample content as the ordinate. The effect quantity can be RR, OR and death ratio OR its logarithm. Theoretically, the point estimation of each independent research effect included in Metaanalysis should be an inverted funnel shape in the plane coordinate system, so it is called funnel diagram.

The sample size is small and the research accuracy is low. It is distributed at the bottom of the funnel diagram and scattered around;

The sample size is large and the research accuracy is high. It is distributed at the top of the funnel diagram and concentrated in the middle.
In actual use, when the number of studies for Metaanalysis is small, it is not suitable to make a funnel chart. Generally, it is recommended that the funnel chart should be made only when the number of studies for Metaanalysis is 10 or more. Funnel chart is mainly used to observe whether there is bias in the results of a system evaluation or Metaanalysis, such as publication bias or other bias, as shown in the following figure:
### funnel plots res < rma(yi, vi, data = dat) funnel(res, main = "RandomEffects Model") res < rma(yi, vi, mods = cbind(ablat), data = dat) funnel(res, main = "MixedEffects Model") #######Begger's test ranktest(res) Rank Correlation Test for Funnel Plot Asymmetry Kendall's tau = 0.0256, p = 0.9524 ##########Egger's test regtest(res) Regression Test for Funnel Plot Asymmetry Model: mixedeffects metaregression model Predictor: standard error Test for Funnel Plot Asymmetry: z = 0.4623, p = 0.6439
 radial star
It mainly reflects the heterogeneity of various studies, so as to find the points of heterogeneity. The effect corresponding to the arc is evaluated by size distribution. As follows:
### radial plots res < rma(yi, vi, data = dat, method = "REML") radial(res, main = "RandomEffects Model") res < rma(yi, vi, data = dat, method = "ML") radial(res, main = "MixedEffects Model")
 QQ bitmap
Draw the prediction results and observe whether they are within the confidence interval, as follows:
### qqnormal plots res < rma(yi, vi, data = dat) qqnorm(res, main = "RandomEffects Model") res < rma(yi, vi, mods = cbind(ablat), data = dat) qqnorm(res, main = "MixedEffects Model")
 labbe labetu
This function calculates the arm level results of the two groups (e.g., treatment and control) and compares them with each other. Incidence rate of incidence rate is incidence of the two groups of functional blocks, especially in the analysis of risk differences. The log odds ratio is analyzed in the log log log ratio. The log odds ratio is analyzed in the analysis of the odds ratio (log), the square root inverse sine transformation ratio is used to analyze the risk difference of the square root sine transform, and the incidence rate is analyzed with the original incidence rate, and the incidence ratio is analyzed by the logarithm of incidence (log). The square root conversion of incidence rate is analyzed as follows: incidence rate difference.
res < rma(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg) ### default plot labbe(res) ### funnel plot with risk values on the x and yaxis and add grid labbe(res, transf=exp, grid=TRUE) # }
04 interpretation of analysis results
It can be seen that I^2 is 71.98%, which belongs to high heterogeneity, and the random effect model can be used. When heterogeneity is low, fixed effect model and random effect model can be used, and the results are not different, but high heterogeneity can only choose random effect model, otherwise the extrapolation of results will be constrained. The random effect model is selected here for conservative consideration.
The random effect model is based on the assumption that there is heterogeneity across studies. The merger model recognizes the existence of heterogeneity among studies, but does not deal with heterogeneity; If there is heterogeneity among the studies included in the merger, although it does not reach the I ^ 2 > 50% we normally set, when merging with fixed effect model, the default operation directly ignores the existence of this part of heterogeneity, so the merger result will cause false positive error. When merging with random effect model, although the heterogeneity is not handled, However, its default operation recognizes the existence of heterogeneity, and the merging result is more reliable!
From the funnel chart and the results of Begg's test and Egger's test, we can see that the P value is greater than 0.05, which means there is no publication bias.
Meta analysis is also an option for publishing articles. To find out these diagrams and internal methods, as for the score of SCI, we still need to see the novelty of the subject we study, and the methods are basically these. As long as you have good data, it's just around the corner to publish a high score of SCI. In the next issue, we'll talk about other software!
Pay attention to the official account, scan the code into groups, free answers, and later there will be free live tutorials. Please look forward to it!
Reference:

Viechtbauer, W. (2010). Conducting MetaAnalyses in R with the metafor Package. Journal of Statistical Software, 36(3), 1–48.

Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., & Mosteller, F. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Metaanalysis of the published literature. Journal of the American Medical Association, 271(9), 698–702.

Berkey, C. S., Hoaglin, D. C., Mosteller, F., & Colditz, G. A. (1995). A randomeffects regression model for metaanalysis. Statistics in Medicine, 14(4), 395–411.

van Houwelingen, H. C., Arends, L. R., & Stijnen, T. (2002). Advanced methods in metaanalysis: Multivariate approach and metaregression. Statistics in Medicine, 21(4), 589–624.

Viechtbauer, W. (2010). Conducting metaanalyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48.
Huanfeng gene
Bioinformatics analysis, SCI article writing and bioinformatics basic knowledge learning: R language learning, perl basic programming, linux system commands, Python meet you better
34 original content
official account