#' Influence Diagnostics #' #' Conducts an influence analysis of a meta-analysis generated by \code{\link[meta]{meta}} functions, #' and allows to produce influence diagnostic plots. #' #' @usage InfluenceAnalysis(x, random = FALSE, subplot.heights = c(30,18), #' subplot.widths = c(30,30), forest.lims = 'default', #' return.separate.plots = FALSE, text.scale = 1) #' #' @param x An object of class \code{meta}, generated by the \code{metabin}, \code{metagen}, #' \code{metacont}, \code{metacor}, \code{metainc}, \code{metarate} or \code{metaprop} function. #' @param random Logical. Should the random-effects model be used to generate the influence diagnostics? #' Uses the \code{method.tau} specified in the \code{meta} object if one #' of "\code{DL}", "\code{HE}", "\code{SJ}", "\code{ML}", "\code{REML}", "\code{EB}", "\code{PM}", "\code{HS}" or "\code{GENQ}" (to ensure compatibility with #' the \code{\link[metafor]{metafor}} package). Otherwise, the DerSimonian-Laird #' (\code{"DL"}; DerSimonian & Laird, 1986) estimator is used. \code{FALSE} by default. #' @param subplot.heights Concatenated array of two numerics. Specifies the heights of the #' first (first number) and second (second number) row of the overall plot generated when plotting the results. #' Default is \code{c(30,18)}. #' @param subplot.widths Concatenated array of two numerics. Specifies the widths of the #' first (first number) and second (second number) column of the overall results plot generated when plotting the results. #' Default is \code{c(30,30)}. #' @param forest.lims Concatenated array of two numerics. Specifies the x-axis limits of the forest plots #' generated when plotting the results. Use \code{"default"} if standard settings should be used (this is the default). #' @param return.separate.plots Logical. When plotted, should the influence plots be shown as separate plots in lieu #' of returning them in one overall plot? #' @param text.scale Positive numeric. Scaling factor for the text geoms used when plotting the results. Values <1 shrink the #' text, while values >1 increase the text size. Default is \code{1}. #' #' @details #' The function conducts an influence analysis using the "Leave-One-Out" paradigm internally #' and produces data for four influence diagnostics. Diagnostic plots can be produced by saving the output of the #' function to an object and plugging it into the \code{plot} function. #' These diagnostics may be used to determine which study or effect size #' may have an excessive influence on the overall results of a meta-analysis and/or contribute substantially to #' the between-study heterogeneity in an analysis. This may be used for outlier detection and to test #' the robustness of the overall results found in an analysis. Results for four diagnostics are calculated: #' \itemize{ #' \item \strong{Baujat Plot}: Baujat et al. (2002) proposed a plot to evaluate heterogeneity patterns in #' a meta-analysis. The x-axis of the Baujat plot shows the overall heterogeneity contribution of each effect size #' while the y-axis shows the influence of each effect size on the pooled result. The \code{\link[meta]{baujat}} #' function is called internally to produce the results. Effect sizes or studies with high values #' on both the x and y-axis may be considered to be influential cases; effect sizes or studies #' with high heterogeneity contribution (x-axis) and low influence on the overall results can be outliers #' which might be deleted to reduce the amount of between-study heterogeneity. #' \item \strong{Influence Characteristics}: Several influence analysis diagnostics #' proposed by Viechtbauer & Cheung (2010). Results are calculated by an internal call #' to \code{\link[metafor]{influence.rma.uni}}. In the console output, potentially influential studies are marked #' with an asterisk (\code{*}). When plotted, effect sizes/studies determined to be influential cases #' using the "rules of thumb" described in Viechtbauer & Cheung (2010) are shown in red. For further #' details, see the documentation of the \code{\link[metafor]{influence.rma.uni}} function. #' \item \strong{Forest Plot for the Leave-One-Out Analysis, sorted by Effect Size}: This #' displays the effect size and \eqn{I^2}-heterogeneity when omitting one of the \eqn{k} studies each time. #' The plot is sorted by effect size to determine which studies or effect sizes particularly #' affect the overall effect size. Results are generated by an internal call to \code{\link[meta]{metainf}}. #' \item \strong{Forest Plot for the Leave-One-Out Analysis, sorted by \eqn{I^2}}: see above; results are sorted #' by \eqn{I^2} to determine the study for which exclusion results in the greatest reduction of heterogeneity. #'} #' #' @references Harrer, M., Cuijpers, P., Furukawa, T.A, & Ebert, D. D. (2019). #' \emph{Doing Meta-Analysis in R: A Hands-on Guide}. DOI: 10.5281/zenodo.2551803. \href{https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/influenceanalyses.html}{Chapter 6.3} #' #' DerSimonian R. & Laird N. (1986), Meta-analysis in clinical trials. \emph{Controlled Clinical Trials, 7}, 177–188. #' #' Viechtbauer, W., & Cheung, M. W.-L. (2010). Outlier and influence diagnostics for meta-analysis. \emph{Research Synthesis Methods, 1}, 112–125. #' #' @author Mathias Harrer & David Daniel Ebert #' #' @return A \code{list} object of class \code{influence.analysis} containing the #' following objects is returned (if results are saved to a variable): #' \itemize{ #' \item \code{BaujatPlot}: The Baujat plot #' \item \code{InfluenceCharacteristics}: The Viechtbauer-Cheung influence characteristics plot #' \item \code{ForestEffectSize}: The forest plot sorted by effect size #' \item \code{ForestI2}: The forest plot sorted by between-study heterogeneity #' \item \code{Data}: A \code{data.frame} containing the data used for plotting. #'} #' Otherwise, the function prints out (1) the results of the Leave-One-Out Analysis (sorted by \eqn{I^2}), #' (2) the Viechtbauer-Cheung Influence Diagnostics and (3) Baujat Plot data (sorted by heterogeneity contribution), #' in this order. Plots can be produced manually by plugging a saved object of class \code{InfluenceAnalysis} generated by #' the function into the \code{plot} function. It is also possible to only produce one specific plot by #' specifying the name of the plot as a \code{character} in the second argument of the \code{plot} call (see Examples). #' #' @import ggplot2 ggrepel grid #' @importFrom gridExtra grid.arrange arrangeGrob #' @importFrom metafor rma.uni influence.rma.uni #' @import meta #' @import utils #' @importFrom graphics abline axis lines mtext par plot points rect segments text #' @importFrom stats as.formula hat influence ks.test optimize pbinom pchisq pf pnorm pt punif qchisq qf qnorm qt reformulate reorder setNames uniroot #' #' @export InfluenceAnalysis #' #' @seealso \code{\link[metafor]{influence.rma.uni}}, \code{\link[meta]{metainf}}, \code{\link[meta]{baujat}} #' #' @examples #' #' \dontrun{ #' # Load 'ThirdWave' data #' data(ThirdWave) #' #' # Create 'meta' meta-analysis object #' suppressPackageStartupMessages(library(meta)) #' meta = metagen(TE, seTE, studlab = paste(ThirdWave$Author), data=ThirdWave) #' #' # Run influence analysis; specify to return separate plots when plotted #' inf.an = InfluenceAnalysis(meta, random = TRUE, #' return.separate.plots = TRUE) #' #' # Show results in console #' inf.an #' #' # Generate all plots #' plot(inf.an) #' #' # For baujat plot #' plot(inf.an, "baujat") #' #' # For influence diagnostics plot #' plot(inf.an, "influence") #' #' # For forest plot sorted by effect size #' plot(inf.an, "ES") #' #' # For forest plot sorted by I-squared #' plot(inf.an, "I2")} ### Influence Analysis function for fixed-effect-model meta-analyses InfluenceAnalysis = function(x, random = FALSE, subplot.heights = c(30, 18), subplot.widths = c(30, 30), forest.lims = "default", return.separate.plots = FALSE, text.scale = 1) { update.meta = getFromNamespace("update.meta", "meta") # Validate x = x if (class(x)[1] %in% c("meta", "metabin", "metagen", "metacont", "metacor", "metainc", "metaprop", "metarate")) { x <- update.meta(x, subset = !(is.na(x$TE) | is.na(x$seTE))) } else { stop("Object 'x' must be of class 'meta', 'metabin', 'metagen', 'metacont', 'metacor', 'metainc', or 'metaprop'") } n.studies = x$k TE = x$TE seTE = x$seTE random = random # Make unique studlabs x$studlab = make.unique(as.character(x$studlab)) if (random %in% c(TRUE, FALSE)) { } else { stop("'random' must be set to either TRUE or FALSE.") } forest.lims = forest.lims if (forest.lims[1] == "default" | (class(forest.lims[1]) == "numeric" & class(forest.lims[2]) == "numeric")) { } else { stop("'forest.lims' must either be 'default' or two concatenated numerics for ymin and ymax.") } return.seperate.plots = return.separate.plots if (return.seperate.plots %in% c(TRUE, FALSE)) { } else { stop("'return.separate.plots' must be set to either TRUE or FALSE.") } heights = subplot.heights if (class(heights[1]) == "numeric" & class(heights[2]) == "numeric") { } else { stop("'subplot.heights' must be two concatenated numerics.") } widths = subplot.widths if (class(widths[1]) == "numeric" & class(widths[2]) == "numeric") { } else { stop("'widths' must be two concatenated numerics.") } text.scale = text.scale if (text.scale > 0) { } else { stop("'text.scale' must be a single number greater 0.") } if (length(unique(x$studlab)) != length(x$studlab)) { stop("'Study labels in the 'meta' object must be unique.") } ######################################## cat("[===============") ######################################## if (random == FALSE) { method.rma = "FE" method.meta = "fixed" } else { if (x$method.tau %in% c("DL", "HE", "SJ", "ML", "REML", "EB", "HS", "GENQ", "PM")) { method.rma = x$method.tau } else { method.rma = "DL" cat("Tau estimator is unkown to metafor::rma; DerSimonian-Laird ('DL') estimator used.") } method.meta = "random" } # Perform Meta-Analysis using metafor, get influence results res = metafor::rma.uni(yi = TE, sei = seTE, measure = "GEN", data = x, method = method.rma, slab = studlab) metafor.inf = influence(res) # Recode inf metafor.inf$is.infl = ifelse(metafor.inf$is.infl == TRUE, "yes", "no") cheungviechtdata = cbind(study = substr(rownames(as.data.frame(metafor.inf$inf)), 1, 3), as.data.frame(metafor.inf$inf), is.infl = metafor.inf$is.infl) rownames(cheungviechtdata) = NULL if (length(unique(cheungviechtdata$study)) < length(cheungviechtdata$study)) { i = 3 while (length(unique(cheungviechtdata$study)) < length(cheungviechtdata$study)) { i = i + 1 cheungviechtdata$study = substr(rownames(as.data.frame(metafor.inf$inf)), 1, i) } } ######################################## cat("===============") ######################################## cheungviechtdata = as.data.frame(cheungviechtdata) within(cheungviechtdata, {study = factor(study, levels = study)}) -> cheungviechtdata # Generate plots rstudent.plot = ggplot2::ggplot(cheungviechtdata, aes(y = rstudent, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("Stand. Residual") dffits.thresh = 3 * sqrt(metafor.inf$p/(metafor.inf$k - metafor.inf$p)) dffits.plot = ggplot2::ggplot(cheungviechtdata, aes(y = dffits, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("DFFITS") cook.d.plot = ggplot2::ggplot(cheungviechtdata, aes(y = cook.d, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("Cook's Distance") cov.r.plot = ggplot2::ggplot(cheungviechtdata, aes(y = cov.r, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("Covariance Ratio") tau2.del.plot = ggplot2::ggplot(cheungviechtdata, aes(y = tau2.del, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("tau-squared (L-0-0)") QE.del.plot = ggplot2::ggplot(cheungviechtdata, aes(y = QE.del, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("Q (L-0-0)") hat.thresh = 3 * (metafor.inf$p/metafor.inf$k) hat.plot = ggplot2::ggplot(cheungviechtdata, aes(y = hat, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("hat") weight.plot = ggplot2::ggplot(cheungviechtdata, aes(y = weight, x = study, color = is.infl, group = 1)) + geom_line(color = "black") + geom_point(size = 2) + scale_color_manual(values = c("blue", "red")) + theme_minimal() + theme(axis.title.x = element_blank(), legend.position = "none", axis.text.x = element_text(angle = 45, size = 5), axis.title.y = element_text(size = 7), axis.text.y = element_text(size = 5)) + ylab("weight") rma.influence.plot = arrangeGrob(rstudent.plot, dffits.plot, cook.d.plot, cov.r.plot, tau2.del.plot, QE.del.plot, hat.plot, weight.plot, ncol = 2) # Perform Influence Analysis on meta object, generate forests meta.inf = meta::metainf(x, pooled = method.meta) if (x$sm %in% c("RR", "OR", "IRR")) { effect = x$sm n.studies = n.studies # Create Sortdat data set for sorting sortdat = data.frame(studlab = meta.inf$studlab, mean = exp(meta.inf$TE), lower = exp(meta.inf$lower), upper = exp(meta.inf$upper), i2 = meta.inf$I2) sortdat2 = sortdat[1:(nrow(sortdat) - 2), ] lastline = sortdat[nrow(sortdat), ] # Change summary label if (random == TRUE) { lastline[1] = "Random-Effects Model" } else { lastline[1] = "Fixed-Effect Model" } for (i in 2:4) { lastline[i] = format(round(lastline[i], 2), nsmall = 2) } # Sort sortdat.es = sortdat2[order(sortdat2$mean), ] sortdat.es$studlab = factor(sortdat.es$studlab, levels = sortdat.es$studlab[order(-sortdat.es$mean)]) sortdat.i2 = sortdat2[order(sortdat2$i2), ] sortdat.i2$studlab = factor(sortdat.i2$studlab, levels = sortdat.i2$studlab[order(-sortdat.i2$i2)]) # Generate Forest Plots if (forest.lims[1] == "default") { if (min(sortdat.es$lower) > 0.5){ min = 0.5 } else { min = NA } if (max(sortdat.es$upper) <= 1){ max = 1.2 } else { max = round(max(sortdat.es$upper) + 6, 0) } } else { if (forest.lims[1] <= 0){ min = NA } else { min = forest.lims[1] } max = forest.lims[2] + 4 } if (method.meta == "fixed"){ plot.sum.effect = exp(x$TE.fixed) plot.sum.lower = exp(x$lower.fixed) plot.sum.upper = exp(x$upper.fixed) } else { plot.sum.effect = exp(x$TE.random) plot.sum.lower = exp(x$lower.random) plot.sum.upper = exp(x$upper.random) } ######################################## cat("===============") ######################################## title.es = with(sortdat.es, { paste0("hat(theta)['*']~'='~", paste0("'", format(round(mean, 2)), "'"), "~'['*", paste0("'", format(round(lower, 2), nsmall = 2), "'"), "*'-'*", paste0("'", format(round(upper, 2), nsmall = 2), "'"), "*']'*';'~italic(I)^2~'='~", paste0("'", format(round(i2, 2)*100,nsmall = 0), "'"), "*'%'")}) title.i2 = with(sortdat.i2,{ paste0("italic(I)^2~'='~", paste0("'", format(round(i2, 2)*100,nsmall = 0), "'"), "*'%'", "*';'~", "hat(theta)['*']~'='~", paste0("'", format(round(mean, 2)), "'"), "~'['*", paste0("'", format(round(lower, 2), nsmall = 2), "'"), "*'-'*", paste0("'", format(round(upper, 2), nsmall = 2), "'"), "*']'")}) forest.es = ggplot(sortdat.es, aes(x = studlab, y = mean, ymin = lower, ymax = upper)) + geom_text(aes(label = title.es, y = Inf), parse = T, hjust = "inward", size = 3 * text.scale) + geom_hline(yintercept = 1, color = "blue") + ylab(paste(effect, " (", as.character(lastline$studlab), ")", sep = "")) + ggtitle("Sorted by Effect Size") + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_text(color = "black", size = 12, face = "bold"), axis.text.y = element_text(color = "black", size = 9 * text.scale), plot.title = element_text(face = "bold", hjust = 0.5), axis.line.x = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.x = element_text(color = "black", size = 9 * text.scale)) + scale_y_continuous(trans = "log2", limits = c(min, max)) + geom_rect(aes(ymin=plot.sum.lower, ymax=plot.sum.upper, xmin=0, xmax=Inf), alpha=0.08, fill="lightgreen", color=NA) + geom_hline(yintercept = plot.sum.effect, color = "darkgreen", linetype="dotted", size=0.5) + geom_point(shape = 15, size = 4.5, color = "grey") + geom_linerange(size = 0.9) + geom_pointrange(shape = 3, size = 0.3) forest.i2 = ggplot(sortdat.i2, aes(x = studlab, y = mean, ymin = lower, ymax = upper)) + geom_text(aes(label = title.i2, y = Inf), parse = T, hjust = "inward", size = 3 * text.scale) + geom_hline(yintercept = 1, color = "blue") + ylab(paste(effect, " (", as.character(lastline$studlab), ")", sep = "")) + ggtitle(expression(Sorted~by~italic(I)^2)) + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_text(color = "black", size = 12, face = "bold"), axis.text.y = element_text(color = "black", size = 9 * text.scale), plot.title = element_text(face = "bold", hjust = 0.5), axis.line.x = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.x = element_text(color = "black", size = 9 * text.scale)) + scale_y_continuous(trans = "log2", limits = c(min, max)) + geom_rect(aes(ymin=plot.sum.lower, ymax=plot.sum.upper, xmin=0, xmax=Inf), alpha=0.08, fill="lightgreen", color=NA) + geom_hline(yintercept = plot.sum.effect, color = "darkgreen", linetype="dotted", size=0.5) + geom_point(shape = 15, size = 4.5, color = "grey") + geom_linerange(size = 0.9) + geom_pointrange(shape = 3, size = 0.3) } else if (class(x)[1] %in% c("metacor", "metaprop", "metarate")) { effect = x$sm n.studies = n.studies # Create Sortdat data set for sorting sortdat = data.frame(studlab = meta.inf$studlab, mean = meta.inf$TE, lower = meta.inf$lower, upper = meta.inf$upper, i2 = meta.inf$I2) sortdat2 = sortdat[1:(nrow(sortdat) - 2), ] lastline = sortdat[nrow(sortdat), ] # Change summary label if (random == TRUE) { lastline[1] = "Random-Effects Model" } else { lastline[1] = "Fixed-Effect Model" } for (i in 2:4) { lastline[i] = format(round(lastline[i], 2), nsmall = 2) } # Sort sortdat.es = sortdat2[order(sortdat2$mean), ] sortdat.es$studlab = factor(sortdat.es$studlab, levels = sortdat.es$studlab[order(-sortdat.es$mean)]) sortdat.i2 = sortdat2[order(sortdat2$i2), ] sortdat.i2$studlab = factor(sortdat.i2$studlab, levels = sortdat.i2$studlab[order(-sortdat.i2$i2)]) # Backtransform backtransformer = function(x, sm, n){ # Define functions z2cor = function(x) { res <- (exp(2 * x) - 1)/(exp(2 * x) + 1) res } logit2p = function(x) { res <- 1/(1 + exp(-x)) res } asin2p = function (x, n = NULL, value = "mean", warn = TRUE) { if (all(is.na(x))) return(x) if (is.null(n)) { minimum <- asin(sqrt(0)) maximum <- asin(sqrt(1)) } else { minimum <- 0.5 * (asin(sqrt(0/(n + 1))) + asin(sqrt((0 + 1)/(n + 1)))) maximum <- 0.5 * (asin(sqrt(n/(n + 1))) + asin(sqrt((n + 1)/(n + 1)))) } sel0 <- x < minimum sel1 <- x > maximum if (any(sel0, na.rm = TRUE)) { if (is.null(n)) { if (warn) warning("Negative value for ", if (length(x) > 1) "at least one ", if (value == "mean") "transformed proportion using arcsine transformation.\n Proportion set to 0.", if (value == "lower") "lower confidence limit using arcsine transformation.\n Lower confidence limit set to 0.", if (value == "upper") "upper confidence limit using arcsine transformation.\n Upper confidence limit set to 0.", sep = "") } else { if (warn) warning("Too small value for ", if (length(x) > 1) "at least one ", if (value == "mean") "transformed proportion using Freeman-Tukey double arcsine transformation.\n Proportion set to 0.", if (value == "lower") "lower confidence limit using Freeman-Tukey double arcsine transformation.\n Lower confidence limit set to 0.", if (value == "upper") "upper confidence limit using Freeman-Tukey double arcsine transformation.\n Upper confidence limit set to 0.", sep = "") } } if (any(sel1, na.rm = TRUE)) { if (is.null(n)) { if (warn) warning("Too large value for ", if (length(x) > 1) "at least one ", if (value == "mean") "transformed proportion using arcsine transformation.\n Proportion set to 1.", if (value == "lower") "lower confidence limit using arcsine transformation.\n Lower confidence limit set to 1.", if (value == "upper") "upper confidence limit using arcsine transformation.\n Upper confidence limit set to 1.", sep = "") } else { if (warn) warning("Too large value for ", if (length(x) > 1) "at least one ", if (value == "mean") "transformed proportion using Freeman-Tukey double arcsine transformation.\n Proportion set to 1.", if (value == "lower") "lower confidence limit using Freeman-Tukey double arcsine transformation.\n Lower confidence limit set to 1.", if (value == "upper") "upper confidence limit using Freeman-Tukey double arcsine transformation.\n Upper confidence limit set to 1.", sep = "") } } res <- rep(NA, length(x)) sel <- !(sel0 | sel1) sel <- !is.na(sel) & sel res[sel0] <- 0 res[sel1] <- 1 if (is.null(n)) { res[sel] <- sin(x[sel])^2 } else { res[sel] <- 0.5 * (1 - sign(cos(2 * x[sel])) * sqrt(1 - (sin(2 * x[sel]) + (sin(2 * x[sel]) - 1/sin(2 * x[sel]))/n[sel])^2)) } res } asin2ir = function (x, time = NULL, value = "mean", warn = TRUE) { if (all(is.na(x))) return(x) minimum <- 0.5 * (sqrt(0/time) + sqrt((0 + 1)/time)) sel0 <- x < minimum if (any(sel0, na.rm = TRUE)) { if (warn) warning("Too small value for ", if (length(x) > 1) "at least one ", if (value == "mean") "transformed proportion using Freeman-Tukey double arcsine transformation.\n Rate set to 0.", if (value == "lower") "lower confidence limit using Freeman-Tukey double arcsine transformation.\n Lower confidence limit set to 0.", if (value == "upper") "upper confidence limit using Freeman-Tukey double arcsine transformation.\n Upper confidence limit set to 0.", sep = "") } res <- rep(NA, length(x)) sel <- !sel0 sel <- !is.na(sel) & sel res[sel0] <- 0 res[sel] <- (1/time[sel] - 8 * x[sel]^2 + 16 * time[sel] * x[sel]^4)/(16 * x[sel]^2 * time[sel]) res[res < 0] <- 0 res } if(sm == "COR"){ res = x } if(sm == "IR"){ res = x } if(sm == "PRAW"){ res = x } if(sm == "ZCOR"){ res = z2cor(x) } if(sm == "PLOGIT"){ res = logit2p(x) } if (sm == "PAS"){ res <- asin2p(x, value = value, warn = FALSE) } if (sm == "PFT"){ res = asin2p(x, n, value = value, warn = FALSE) } if (sm == "IRS"){ res = x^2 } if (sm == "IRFT"){ res = asin2ir(x, time=n, value = value, warn = FALSE) } if (sm == "IRLN"){ res = exp(x) } if (sm == "PLN"){ res = exp(x) } res } if (class(x)[1] %in% c("metaprop", "metacor")){ n.h.m = meta.inf$n.harmonic.mean[1:(length(meta.inf$n.harmonic.mean)-2)] n.h.m.tot = meta.inf$n.harmonic.mean[length(meta.inf$n.harmonic.mean)] n.h.m.es = n.h.m[order(sortdat.es$mean)] n.h.m.i2 = n.h.m[order(sortdat.i2$mean)] sortdat.es$mean = backtransformer(sortdat.es$mean, sm=effect, n=n.h.m.es) sortdat.es$lower = backtransformer(sortdat.es$lower, sm=effect, n=n.h.m.es) sortdat.es$upper = backtransformer(sortdat.es$upper, sm=effect, n=n.h.m.es) sortdat.i2$mean = backtransformer(sortdat.i2$mean, sm=effect, n=n.h.m.i2) sortdat.i2$lower = backtransformer(sortdat.i2$lower, sm=effect, n=n.h.m.i2) sortdat.i2$upper = backtransformer(sortdat.i2$upper, sm=effect, n=n.h.m.i2) if (method.meta == "fixed"){ plot.sum.effect = backtransformer(x$TE.fixed, sm=effect, n=n.h.m.tot) plot.sum.lower = backtransformer(x$lower.fixed, sm=effect, n=n.h.m.tot) plot.sum.upper = backtransformer(x$upper.fixed, sm=effect, n=n.h.m.tot) } else { plot.sum.effect = backtransformer(x$TE.random, sm=effect, n=n.h.m.tot) plot.sum.lower = backtransformer(x$lower.random, sm=effect, n=n.h.m.tot) plot.sum.upper = backtransformer(x$upper.random, sm=effect, n=n.h.m.tot) } } else { if(meta.inf$sm == "IRFT"){ n.h.m = meta.inf$t.harmonic.mean[1:(length(meta.inf$t.harmonic.mean)-2)] n.h.m.es = n.h.m[order(sortdat.es$mean)] n.h.m.i2 = n.h.m[order(sortdat.i2$mean)] n.h.m.tot = meta.inf$t.harmonic.mean[length(meta.inf$t.harmonic.mean)] sortdat.es$mean = backtransformer(sortdat.es$mean, sm=effect, n=n.h.m.es) sortdat.es$lower = backtransformer(sortdat.es$lower, sm=effect, n=n.h.m.es) sortdat.es$upper = backtransformer(sortdat.es$upper, sm=effect, n=n.h.m.es) sortdat.i2$mean = backtransformer(sortdat.i2$mean, sm=effect, n=n.h.m.i2) sortdat.i2$lower = backtransformer(sortdat.i2$lower, sm=effect, n=n.h.m.i2) sortdat.i2$upper = backtransformer(sortdat.i2$upper, sm=effect, n=n.h.m.i2) if (method.meta == "fixed"){ plot.sum.effect = backtransformer(x$TE.fixed, sm=effect, n=n.h.m.tot) plot.sum.lower = backtransformer(x$lower.fixed, sm=effect, n=n.h.m.tot) plot.sum.upper = backtransformer(x$upper.fixed, sm=effect, n=n.h.m.tot) } else { plot.sum.effect = backtransformer(x$TE.random, sm=effect, n=n.h.m.tot) plot.sum.lower = backtransformer(x$lower.random, sm=effect, n=n.h.m.tot) plot.sum.upper = backtransformer(x$upper.random, sm=effect, n=n.h.m.tot) } } else { n.h.m.tot = meta.inf$n.harmonic.mean[length(meta.inf$n.harmonic.mean)] sortdat.es$mean = backtransformer(sortdat.es$mean, sm=effect, n=NULL) sortdat.es$lower = backtransformer(sortdat.es$lower, sm=effect, n=NULL) sortdat.es$upper = backtransformer(sortdat.es$upper, sm=effect, n=NULL) sortdat.i2$mean = backtransformer(sortdat.i2$mean, sm=effect, n=NULL) sortdat.i2$lower = backtransformer(sortdat.i2$lower, sm=effect, n=NULL) sortdat.i2$upper = backtransformer(sortdat.i2$upper, sm=effect, n=NULL) if (method.meta == "fixed"){ plot.sum.effect = backtransformer(x$TE.fixed, sm=effect, n=n.h.m.tot) plot.sum.lower = backtransformer(x$lower.fixed, sm=effect, n=n.h.m.tot) plot.sum.upper = backtransformer(x$upper.fixed, sm=effect, n=n.h.m.tot) } else { plot.sum.effect = backtransformer(x$TE.random, sm=effect, n=n.h.m.tot) plot.sum.lower = backtransformer(x$lower.random, sm=effect, n=n.h.m.tot) plot.sum.upper = backtransformer(x$upper.random, sm=effect, n=n.h.m.tot) } } } # Generate Forest Plots if (forest.lims[1] == "default") { if (class(x)[1] == "metacor"){ min = min(sortdat.es$mean)-0.2 } else { min = -0.2 } max = max(sortdat.es$mean) + 0.5 } else { min = forest.lims[1] max = forest.lims[2] } # Set ggtitles if (class(x)[1] == "metaprop"){ ggtitl = as.character("Proportion") } else if (class(x)[1] == "metacor"){ ggtitl = as.character("Correlation") } else { ggtitl = as.character("Rate") } ######################################## cat("===============") ######################################## title.es = with(sortdat.es, { paste0("hat(theta)['*']~'='~", paste0("'", format(round(mean, 2)), "'"), "~'['*", paste0("'", format(round(lower, 2), nsmall = 2), "'"), "*'-'*", paste0("'", format(round(upper, 2), nsmall = 2), "'"), "*']'*';'~italic(I)^2~'='~", paste0("'", format(round(i2, 2)*100,nsmall = 0), "'"), "*'%'")}) title.i2 = with(sortdat.i2,{ paste0("italic(I)^2~'='~", paste0("'", format(round(i2, 2)*100,nsmall = 0), "'"), "*'%'", "*';'~", "hat(theta)['*']~'='~", paste0("'", format(round(mean, 2)), "'"), "~'['*", paste0("'", format(round(lower, 2), nsmall = 2), "'"), "*'-'*", paste0("'", format(round(upper, 2), nsmall = 2), "'"), "*']'")}) forest.es = ggplot(sortdat.es, aes(x = studlab, y = mean, ymin = lower, ymax = upper)) + geom_text(aes(label = title.es, y = Inf), parse = T, hjust = "inward", size = 3 * text.scale) + geom_hline(yintercept = 0, color = "blue") + ylab(paste(ggtitl, " (", as.character(lastline$studlab), ")", sep = "")) + ggtitle(paste("Sorted by", ggtitl)) + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_text(color = "black", size = 12, face = "bold"), axis.text.y = element_text(color = "black", size = 9 * text.scale), plot.title = element_text(face = "bold", hjust = 0.5), axis.line.x = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.x = element_text(color = "black", size = 9 * text.scale)) + scale_y_continuous(limits = c(min, max)) + geom_rect(aes(ymin=plot.sum.lower, ymax=plot.sum.upper, xmin=0, xmax=Inf), alpha=0.08, fill="lightgreen", color=NA) + geom_hline(yintercept = plot.sum.effect, color = "darkgreen", linetype="dotted", size=0.5) + geom_point(shape = 15, size = 4.5, color = "grey") + geom_linerange(size = 0.9) + geom_pointrange(shape = 3, size = 0.3) forest.i2 = ggplot(sortdat.i2, aes(x = studlab, y = mean, ymin = lower, ymax = upper)) + geom_text(aes(label = title.i2, y = Inf), parse = T, hjust = "inward", size = 3 * text.scale) + geom_hline(yintercept = 0, color = "blue") + ylab(paste(ggtitl, " (", as.character(lastline$studlab), ")", sep = "")) + ggtitle(expression(Sorted~by~italic(I)^2)) + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_text(color = "black", size = 12, face = "bold"), axis.text.y = element_text(color = "black", size = 9 * text.scale), plot.title = element_text(face = "bold", hjust = 0.5), axis.line.x = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.x = element_text(color = "black", size = 9 * text.scale)) + scale_y_continuous(limits = c(min, max)) + geom_rect(aes(ymin=plot.sum.lower, ymax=plot.sum.upper, xmin=0, xmax=Inf), alpha=0.08, fill="lightgreen", color=NA) + geom_hline(yintercept = plot.sum.effect, color = "darkgreen", linetype="dotted", size=0.5) + geom_point(shape = 15, size = 4.5, color = "grey") + geom_linerange(size = 0.9) + geom_pointrange(shape = 3, size = 0.3) } else { # Create Sortdat data set for sorting sortdat = data.frame(studlab = meta.inf$studlab, mean = meta.inf$TE, lower = meta.inf$lower, upper = meta.inf$upper, i2 = meta.inf$I2) sortdat2 = sortdat[1:(nrow(sortdat) - 2), ] lastline = sortdat[nrow(sortdat), ] # Change summary label if (random == TRUE) { lastline[1] = "Random-Effects Model" } else { lastline[1] = "Fixed-Effect Model" } for (i in 2:4) { lastline[i] = format(round(lastline[i], 2), nsmall = 2) } # Sort sortdat.es = sortdat2[order(sortdat2$mean), ] sortdat.es$studlab = factor(sortdat.es$studlab, levels = sortdat.es$studlab[order(-sortdat.es$mean)]) sortdat.i2 = sortdat2[order(sortdat2$i2), ] sortdat.i2$studlab = factor(sortdat.i2$studlab, levels = sortdat.i2$studlab[order(-sortdat.i2$i2)]) # Generate Forest Plots if (forest.lims[1] == "default") { min = round(min(sortdat.es$lower) - 0.1, 2) max = round(max(sortdat.es$upper) + 0.3, 2) } else { min = forest.lims[1] max = forest.lims[2] } if (method.meta == "fixed"){ plot.sum.effect = x$TE.fixed plot.sum.lower = x$lower.fixed plot.sum.upper = x$upper.fixed } else { plot.sum.effect = x$TE.random plot.sum.lower = x$lower.random plot.sum.upper = x$upper.random } ######################################## cat("===============") ######################################## title.es = with(sortdat.es, { paste0("hat(theta)['*']~'='~", paste0("'", format(round(mean, 2)), "'"), "~'['*", paste0("'", format(round(lower, 2), nsmall = 2), "'"), "*'-'*", paste0("'", format(round(upper, 2), nsmall = 2), "'"), "*']'*';'~italic(I)^2~'='~", paste0("'", format(round(i2, 2)*100,nsmall = 0), "'"), "*'%'")}) title.i2 = with(sortdat.i2,{ paste0("italic(I)^2~'='~", paste0("'", format(round(i2, 2)*100,nsmall = 0), "'"), "*'%'", "*';'~", "hat(theta)['*']~'='~", paste0("'", format(round(mean, 2)), "'"), "~'['*", paste0("'", format(round(lower, 2), nsmall = 2), "'"), "*'-'*", paste0("'", format(round(upper, 2), nsmall = 2), "'"), "*']'")}) forest.es = ggplot(sortdat.es, aes(x = studlab, y = mean, ymin = lower, ymax = upper)) + geom_text(aes(label = title.es, y = Inf), parse = T, hjust = "inward", size = 3 * text.scale) + geom_hline(yintercept = 0, color = "blue") + ylim(min, max) + ylab(paste("Effect Size (", as.character(lastline$studlab), ")", sep = "")) + ggtitle("Sorted by Effect Size") + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_text(color = "black", size = 12, face = "bold"), axis.text.y = element_text(color = "black", size = 9 * text.scale), plot.title = element_text(face = "bold", hjust = 0.5), axis.line.x = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.x = element_text(color = "black", size = 9 * text.scale)) + geom_rect(aes(ymin=plot.sum.lower, ymax=plot.sum.upper, xmin=0, xmax=Inf), alpha=0.08, fill="lightgreen", color=NA) + geom_hline(yintercept = plot.sum.effect, color = "darkgreen", linetype="dotted", size=0.5) + geom_point(shape = 15, size = 4.5, color = "grey") + geom_linerange(size = 0.9) + geom_pointrange(shape = 3, size = 0.3) forest.i2 = ggplot(sortdat.i2, aes(x = studlab, y = mean, ymin = lower, ymax = upper)) + geom_text(aes(label = title.i2, y = Inf), parse = T, hjust = "inward", size = 3 * text.scale) + geom_hline(yintercept = 0, color = "blue") + ylim(min, max) + ylab(paste("Effect Size (", as.character(lastline$studlab), ")", sep = "")) + ggtitle(expression(Sorted~by~italic(I)^2)) + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_text(color = "black", size = 12, face = "bold"), axis.text.y = element_text(color = "black",size = 9 * text.scale), plot.title = element_text(face = "bold", hjust = 0.5), axis.line.x = element_line(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.x = element_text(color = "black", size = 9 * text.scale)) + geom_rect(aes(ymin=plot.sum.lower, ymax=plot.sum.upper, xmin=0, xmax=Inf), alpha=0.08, fill="lightgreen", color=NA) + geom_hline(yintercept = plot.sum.effect, color = "darkgreen", linetype="dotted", size=0.5) + geom_point(shape = 15, size = 4.5, color = "grey") + geom_linerange(size = 0.9) + geom_pointrange(shape = 3, size = 0.3) } # Generate baujat plot Define baujat.silent baujat.silent = function(x, yscale = 1, xlim, ylim, ...) { TE = x$TE seTE = x$seTE TE.fixed = metagen(TE, seTE, exclude = x$exclude)$TE.fixed k = x$k studlab = x$studlab SE = x$seTE m.inf = metainf(x, pooled = "fixed") TE.inf = m.inf$TE[1:length(TE)] seTE.inf = m.inf$seTE[1:length(TE)] ys = (TE.inf - TE.fixed)^2/seTE.inf^2 ys = ys * yscale xs = (TE - TE.fixed)^2/seTE^2 if (!is.null(x$exclude)) xs[x$exclude] = 0 res = data.frame(studlab = studlab, x = xs, y = ys, se = SE) return(res) } ######################################## cat("===============") ######################################## bjt = baujat.silent(x) BaujatPlot = ggplot(bjt, aes(x = x, y = y)) + geom_point(aes(size = (1/se)), color = "blue", alpha = 0.75) + geom_rug(color = "lightgray", alpha = 0.5) + theme(legend.position = "none") + xlab("Overall heterogeneity contribution") + ylab("Influence on pooled result") + geom_label_repel(label = bjt$studlab, color = "black", size = 1.5 * text.scale, alpha = 0.7) # Return ######################################## cat("===============] DONE \n") ######################################## # Prepare data for return return.data = cbind(sortdat2, cheungviechtdata[, 2:ncol(cheungviechtdata)], HetContrib = bjt$x, InfluenceEffectSize = bjt$y) if (x$sm %in% c("RR", "OR", "IRR")) {colnames(return.data)[1:2] = c("Author", effect)} else {colnames(return.data)[1:2] = c("Author", "effect")} returnlist = suppressWarnings(suppressMessages(list(BaujatPlot = BaujatPlot, InfluenceCharacteristics = rma.influence.plot, ForestEffectSize = forest.es, ForestI2 = forest.i2, Data = return.data))) if (return.separate.plots == T){class(returnlist) = c("InfluenceAnalysis", "rsp")} if (return.separate.plots == F){class(returnlist) = c("InfluenceAnalysis", "rsp.null")} returnlist }