#' 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{metamean}, \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
#' library(meta)
#' meta = metagen(TE, seTE, studlab = 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", "metamean", "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', 'metamean', '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)
        if (x$sm == "MLN") {
          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") {
            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

}