R/metapsyInfluenceAnalysis.R
In MathiasHarrer/metapsyTools: Several R Helper Functions For The "Metapsy" Database
Defines functions
metapsyInfluenceAnalysis
Documented in
metapsyInfluenceAnalysis
#' 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 metapsyInfluenceAnalysis(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
#' @importFrom meta metainf
#' @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 rnorm punif qchisq qf qnorm qt reformulate reorder setNames uniroot model.matrix rstudent dffits model.matrix
#'
#' @export metapsyInfluenceAnalysis
#'
#' @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, 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")}
metapsyInfluenceAnalysis = function(x, random = FALSE, subplot.heights = c(30, 18),
subplot.widths = c(30, 30),
forest.lims = "default", return.separate.plots = FALSE,
text.scale = 1) {
# Validate
x = x
if (class(x)[1] %in% c("meta", "metabin", "metagen", "metacont", "metacor", "metainc", "metaprop", "metarate")) {
} 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(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.")
}
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)
}
}
# If study labels are only numeric: reset level indexing
if (sum(grepl("[A-Za-z]", levels(as.factor(cheungviechtdata$study)), perl = T)) == 0){
cheungviechtdata$study = factor(cheungviechtdata$study, levels = sort(as.numeric(levels(cheungviechtdata$study))))
}
scalefun = function(x) sprintf("%.1f", x)
cheungviechtdata = as.data.frame(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") +
scale_y_continuous(labels = scalefun)
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") +
scale_y_continuous(labels = scalefun)
# geom_hline(yintercept = dffits.thresh, linetype='dashed', color='black')
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") +
scale_y_continuous(labels = scalefun)
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") +
scale_y_continuous(labels = scalefun)
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)") +
scale_y_continuous(labels = scalefun)
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)") +
scale_y_continuous(labels = scalefun)
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") + scale_y_continuous(labels = scalefun)
# geom_hline(yintercept = hat.thresh, linetype='dashed', color='black')
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") +
scale_y_continuous(labels = scalefun)
rma.influence.plot = gridExtra::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)
}
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")
}
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
}
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)
}
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
# 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
}
MathiasHarrer/metapsyTools documentation built on May 1, 2022, 5:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions metapsyInfluenceAnalysis
Documented in metapsyInfluenceAnalysis
#' 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 metapsyInfluenceAnalysis(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
#' @importFrom meta metainf
#' @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 rnorm punif qchisq qf qnorm qt reformulate reorder setNames uniroot model.matrix rstudent dffits model.matrix
#'
#' @export metapsyInfluenceAnalysis
#'
#' @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, 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")}
metapsyInfluenceAnalysis = function(x, random = FALSE, subplot.heights = c(30, 18),
subplot.widths = c(30, 30),
forest.lims = "default", return.separate.plots = FALSE,
text.scale = 1) {
# Validate
x = x
if (class(x)[1] %in% c("meta", "metabin", "metagen", "metacont", "metacor", "metainc", "metaprop", "metarate")) {
} 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(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.")
}
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)
}
}
# If study labels are only numeric: reset level indexing
if (sum(grepl("[A-Za-z]", levels(as.factor(cheungviechtdata$study)), perl = T)) == 0){
cheungviechtdata$study = factor(cheungviechtdata$study, levels = sort(as.numeric(levels(cheungviechtdata$study))))
}
scalefun = function(x) sprintf("%.1f", x)
cheungviechtdata = as.data.frame(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") +
scale_y_continuous(labels = scalefun)
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") +
scale_y_continuous(labels = scalefun)
# geom_hline(yintercept = dffits.thresh, linetype='dashed', color='black')
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") +
scale_y_continuous(labels = scalefun)
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") +
scale_y_continuous(labels = scalefun)
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)") +
scale_y_continuous(labels = scalefun)
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)") +
scale_y_continuous(labels = scalefun)
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") + scale_y_continuous(labels = scalefun)
# geom_hline(yintercept = hat.thresh, linetype='dashed', color='black')
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") +
scale_y_continuous(labels = scalefun)
rma.influence.plot = gridExtra::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)
}
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")
}
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
}
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)
}
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
# 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
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.