R/plotPermResults.R

In yue-wang-biomath/RgnTX: Colocalization analysis of transcriptome elements in the presence of isoform heterogeneity and ambiguity

#' Plot permutation test results
#' @export plotPermResults
#' @importFrom graphics hist
#' @importFrom stats dnorm qnorm
#' @importFrom ggplot2 ggplot aes annotate arrow element_blank element_line element_text geom_histogram geom_line geom_polygon geom_segment labs stat_function theme unit xlab ylab ylim
#'
#' @description Show a graphical representation of permutation test.
#'
#' @usage plotPermResults(permTestTx_results = NULL, breaks = 15, alpha = 0.05,
#' test_type = "one-sided", binwidth = NULL)
#'
#' @param permTestTx_results A \code{permTestTx.results} list object, which can be generated by the \code{permTestTx} function.
#' @param breaks Histogram breaks. Default is 15.
#' @param alpha Significance level.
#' @param test_type The type of the test. This argument only receives either two options "one-sided" or "two-sided". Default is "one-sided".
#' @param binwidth Histogram binwidth.
#'
#' @return A plot object.
#'
#' @seealso \code{\link{permTestTx}}
#'
#' @examples
#' file <- system.file(package="RgnTX", "extdata", "permTestTx_results.rds")
#' permTestTx_results <- readRDS(file)
#' p_a <- plotPermResults(permTestTx_results, binwidth = 1)
#' p_a
plotPermResults <- function(permTestTx_results = NULL, breaks = 15, alpha = 0.05,
                            test_type = "one-sided", binwidth = NULL) {
    if(!is(permTestTx_results, 'permTestTx.results')){
        stop("Argument permTestTx_results must be a permTestTx.results object.")
    }
    rand.ev <- permTestTx_results$rand.ev
    orig.ev <- permTestTx_results$orig.ev
    zscore <- permTestTx_results$zscore
    pval <- permTestTx_results$pval
    ntimes <- length(rand.ev)
    type <- test_type

    if (orig.ev >= mean(rand.ev)) {
        if (test_type == "one-sided") {
            xcoords <- rand.ev
            rand.mean <- mean(xcoords, na.rm = TRUE)
            rand.sd <- sd(xcoords, na.rm = TRUE)
            xcoords <- xcoords[order(xcoords)]
            y <- dnorm(xcoords, rand.mean, sd = rand.sd)
            x.breaks <- hist(xcoords, breaks = breaks, plot = FALSE)$breaks # breaks, counts, density
            y.density <- hist(xcoords, breaks = breaks, plot = FALSE)$density # breaks, counts, density

            ymax <- max(max(y, na.rm = TRUE), max(y.density, na.rm = TRUE))

            # greater
            aux <- qnorm((1 - alpha), mean = rand.mean, sd = rand.sd)

            xmin <- min(orig.ev, min(xcoords, na.rm = TRUE), min(aux, na.rm = TRUE),
                na.rm = TRUE
            )
            xmax <- max(orig.ev, max(xcoords, na.rm = TRUE), max(aux, na.rm = TRUE),
                na.rm = TRUE
            )
            if (length(binwidth) == 0) {
                binwidth <- x.breaks[2] - x.breaks[1]
            }
            # plot histogram and normal
            data <- data.frame(x = xcoords)

            # plot polygon
            polygon1.x <- seq(aux, xmax, length = 50)
            polygon1.y <- dnorm(polygon1.x, rand.mean, rand.sd)

            d1 <- data.frame(x = c(aux, polygon1.x, xmax), y = c(
                0, polygon1.y, 0
            ))
            d2 <- data.frame(x = c(aux, xmax, xmax, aux), y = c(-0.05 *
                ymax, -0.05 * ymax, ymax, ymax))
            au2max <- xmax - aux
            # plot arrow
            darrow <- data.frame(x = c(aux + 0.4 * au2max, xmax - 0.4 *
                au2max), y = c(0.45 * max(y), 0.45 * max(y)), vx = c(-0.4 *
                au2max, 0.4 * au2max), vy = c(0, 0))

            dmean <- data.frame(x = rand.mean, y = 0, vx = 0, vy = ymax)
            dobv <- data.frame(x = orig.ev, y = 0, vx = 0, vy = ymax)

            # plot critical region
            text_size1 <- 11
            text_size2 <- 4

            p1 <- ggplot(data, aes(x)) +
                labs(subtitle = paste(
                    "p-value:",
                    pval, "\n", "ntimes:", ntimes, "\n", "zscore:", zscore,
                    "\n", "alpha:", alpha, "\n", "type:", paste(type, "test")
                )) +
                theme(
                    panel.background = element_blank(), plot.title = element_text(hjust = 0.5),
                    plot.subtitle = element_text(hjust = 1), panel.grid.major = element_line(
                        colour = "grey",
                        linetype = 9, size = 0.2
                    ), axis.ticks = element_blank(),
                    line = element_line(
                        colour = "white", size = 0.5, linetype = 9,
                        lineend = "butt"
                    ), legend.position = "bottom",
                    axis.text = element_text(size = text_size1),
                    legend.text = element_text(size = text_size1),
                    axis.title.y = element_text(size = text_size1, hjust = 0.6),
                    title = element_text(size = text_size1, face = "bold")
                ) +
                geom_histogram(aes(y = ..density..), binwidth = binwidth, col = "grey", fill = "grey") +
                xlab("") +
                ylab("Density
                    ") +
                stat_function(
                    fun = dnorm,
                    args = list(
                        mean = mean(data$x),
                        sd = sd(data$x)
                    ),
                    col = "black",
                    size = 1, n = 200
                ) +
                geom_polygon(data = d1, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .9) +
                geom_polygon(data = d2, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .1) +
                geom_segment(data = darrow, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), arrow = arrow(length = unit(0.20, "cm"), type = "closed"), size = 0.5, color = "black") +
                annotate("text", x = aux + 1 / 2 * au2max, y = 0.53 * max(y), label = "Critical Region", size = text_size2) +
                annotate("text", x = aux + 1 / 2 * au2max, y = 0.45 * max(y), label = paste0(alpha * 100, "%"), size = text_size2) +
                annotate("text", x = aux + 1 / 2 * au2max, y = 0.37 * max(y), label = paste("alpha =", alpha), size = text_size2) +
                annotate("text", x = rand.mean, y = 1.05 * max(y), label = "E.perm", size = text_size2, color = "blue", fontface = "bold") +
                annotate("text", x = orig.ev, y = 1.05 * max(y), label = "E.obs", size = text_size2, color = "red", fontface = "bold") +
                annotate("text", x = 0.99 * rand.mean, y = 0.3 * max(y), label = paste0(100 - alpha * 100, "%"), size = text_size2, color = "black") +
                geom_segment(data = dmean, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "blue") +
                geom_segment(data = dobv, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "red")
        }
        if (test_type == "two-sided") {
            xcoords <- rand.ev
            rand.mean <- mean(xcoords, na.rm = TRUE)
            rand.sd <- sd(xcoords, na.rm = TRUE)
            xcoords <- xcoords[order(xcoords)]
            y <- dnorm(xcoords, rand.mean, sd = rand.sd)
            x.breaks <- hist(xcoords, breaks = breaks, plot = FALSE)$breaks # breaks, counts, density
            y.density <- hist(xcoords, breaks = breaks, plot = FALSE)$density # breaks, counts, density
            ymax <- max(max(y, na.rm = TRUE), max(y.density, na.rm = TRUE))

            # greater
            aux <- qnorm(1 - alpha / 2, mean = rand.mean, sd = rand.sd)
            aux2 <- qnorm(alpha / 2, mean = rand.mean, sd = rand.sd)

            xmin <- min(orig.ev, min(xcoords, na.rm = TRUE), min(aux, na.rm = TRUE),
                na.rm = TRUE
            )
            xmax <- max(orig.ev, max(xcoords, na.rm = TRUE), max(aux, na.rm = TRUE),
                na.rm = TRUE
            )

            extendwidth <- max(rand.mean - xmin, xmax - rand.mean)
            xmax <- rand.mean + extendwidth
            xmin <- rand.mean - extendwidth

            if (length(binwidth) == 0) {
                binwidth <- x.breaks[2] - x.breaks[1]
            }
            # plot histogram and normal
            data <- data.frame(x = xcoords)

            # plot polygon
            polygon1.x <- seq(aux, xmax, length = 50)
            polygon1.y <- dnorm(polygon1.x, rand.mean, rand.sd)

            d1 <- data.frame(x = c(aux, polygon1.x, xmax), y = c(
                0, polygon1.y,
                0
            ))
            d2 <- data.frame(x = c(aux, xmax, xmax, aux), y = c(-0.05 *
                ymax, -0.05 * ymax, ymax, ymax))
            au2max <- xmax - aux

            # plot arrow
            darrow <- data.frame(x = c(aux + 0.45 * au2max, xmax - 0.45 *
                au2max), y = c(0.45 * max(y), 0.45 * max(y)), vx = c(-0.45 *
                au2max, 0.45 * au2max), vy = c(0, 0))

            # plot polygon2
            polygon2.x <- seq(xmin, aux2, length = 50)
            polygon2.y <- dnorm(polygon2.x, rand.mean, rand.sd)

            d3 <- data.frame(x = c(xmin, polygon2.x, aux2), y = c(
                0, polygon2.y,
                0
            ))
            d4 <- data.frame(x = c(xmin, aux2, aux2, xmin), y = c(-0.05 *
                ymax, -0.05 * ymax, ymax, ymax))
            au2max <- aux2 - xmin

            # plot arrow2
            darrow2 <- data.frame(x = c(xmin + 0.4 * au2max, aux2 - 0.4 *
                au2max), y = c(0.45 * max(y), 0.45 * max(y)), vx = c(-0.4 *
                au2max, 0.4 * au2max), vy = c(0, 0))


            dmean <- data.frame(x = rand.mean, y = 0, vx = 0, vy = ymax)
            dobv <- data.frame(x = orig.ev, y = 0, vx = 0, vy = ymax)
            text_size1 <- 11
            text_size2 <- 4
            p1 <- ggplot(data, aes(x)) +
                labs(subtitle = paste("p-value:", pval * 2, "\n", "ntimes:", ntimes, "\n", "zscore:", zscore, "\n", "alpha:", alpha, "\n", "type:", paste(type, "test"))) +
                theme(
                    panel.background = element_blank(),
                    plot.title = element_text(hjust = 0.5),
                    plot.subtitle = element_text(hjust = 1),
                    panel.grid.major = element_line(colour = "grey", linetype = 9, size = 0.2),
                    axis.ticks = element_blank(),
                    line = element_line(colour = "white", size = 0.5, linetype = 9, lineend = "butt"),
                    legend.position = "bottom",
                    axis.text = element_text(size = text_size1),
                    legend.text = element_text(size = text_size1),
                    axis.title.y = element_text(size = text_size1, hjust = 0.6),
                    title = element_text(size = text_size1, face = "bold")
                ) + # Draw histogram with densit
                geom_histogram(aes(y = ..density..), binwidth = binwidth, col = "grey", fill = "grey") +
                xlab("") +
                ylab("Density
                    ") +
                stat_function(
                    fun = dnorm,
                    args = list(
                        mean = mean(data$x),
                        sd = sd(data$x)
                    ),
                    col = "black",
                    size = 1, n = 200
                ) +
                geom_polygon(data = d1, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .9) +
                geom_polygon(data = d2, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .1) +
                geom_polygon(data = d3, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .9) +
                geom_polygon(data = d4, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .1) +
                geom_segment(data = darrow, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), arrow = arrow(length = unit(0.20, "cm"), type = "closed"), size = 0.5, color = "black") +
                geom_segment(data = darrow2, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), arrow = arrow(length = unit(0.20, "cm"), type = "closed"), size = 0.5, color = "black") +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.53 * max(y), label = "Critical Region", size = text_size2) +
                annotate("text", x = xmax - 1 / 2 * au2max, y = 0.53 * max(y), label = "Critical Region", size = text_size2) +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.45 * max(y), label = paste0(alpha / 2 * 100, "%"), size = text_size2) +
                annotate("text", x = xmax - 1 / 2 * au2max, y = 0.45 * max(y), label = paste0(alpha / 2 * 100, "%"), size = text_size2) +
                annotate("text", x = rand.mean, y = 1.05 * max(y), label = "E.perm", size = text_size2, color = "blue", fontface = "bold") +
                annotate("text", x = orig.ev, y = 1.05 * max(y), label = "E.obs", size = text_size2, color = "red", fontface = "bold") +
                annotate("text", x = 0.99 * rand.mean, y = 0.3 * max(y), label = paste0(100 - alpha * 100, "%"), size = text_size2, color = "black") +
                geom_segment(data = dmean, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "blue") +
                geom_segment(data = dobv, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "red")
        }
    }


    if (orig.ev < mean(rand.ev)) {
        if (test_type == "one-sided") {
            xcoords <- rand.ev
            rand.mean <- mean(xcoords, na.rm = TRUE)
            rand.sd <- sd(xcoords, na.rm = TRUE)
            xcoords <- xcoords[order(xcoords)]
            y <- dnorm(xcoords, rand.mean, sd = rand.sd)
            x.breaks <- hist(xcoords, breaks = breaks, plot = FALSE)$breaks # breaks, counts, density
            y.density <- hist(xcoords, breaks = breaks, plot = FALSE)$density # breaks, counts, density

            ymax <- max(max(y, na.rm = TRUE), max(y.density, na.rm = TRUE))

            # less
            aux <- qnorm(alpha, mean = rand.mean, sd = rand.sd)
            xmin <- min(orig.ev, min(xcoords, na.rm = TRUE), min(aux, na.rm = TRUE),
                na.rm = TRUE
            )
            xmax <- max(orig.ev, max(xcoords, na.rm = TRUE), max(aux, na.rm = TRUE),
                na.rm = TRUE
            )

            if (length(binwidth) == 0) {
                binwidth <- x.breaks[2] - x.breaks[1]
            }
            # plot histogram and normal
            data <- data.frame(x = xcoords)

            # plot polygon
            polygon1.x <- seq(xmin, aux, length = 50)
            polygon1.y <- dnorm(polygon1.x, rand.mean, rand.sd)

            d1 <- data.frame(x = c(xmin, polygon1.x, aux), y = c(
                0, polygon1.y,
                0
            ))
            d2 <- data.frame(x = c(aux, aux, xmin, xmin), y = c(-0.05 *
                ymax, ymax, ymax, -0.05 * ymax))
            au2max <- aux - xmin
            # plot arrow
            darrow <- data.frame(x = c(xmin + 0.4 * au2max, aux - 0.4 *
                au2max), y = c(0.45 * max(y), 0.45 * max(y)), vx = c(-0.4 *
                au2max, 0.4 * au2max), vy = c(0, 0))
            dmean <- data.frame(x = rand.mean, y = 0, vx = 0, vy = ymax)
            dobv <- data.frame(x = orig.ev, y = 0, vx = 0, vy = ymax)

            # plot critical region
            text_size1 <- 11
            text_size2 <- 4
            p1 <- ggplot(data, aes(x)) +
                labs(subtitle = paste("p-value:", pval, "\n", "ntimes:", ntimes, "\n", "zscore:", zscore, "\n", "alpha:", alpha, "\n", "type:", paste(type, "test"))) +
                theme(
                    panel.background = element_blank(),
                    plot.title = element_text(hjust = 0.5),
                    plot.subtitle = element_text(hjust = 1),
                    panel.grid.major = element_line(colour = "grey", linetype = 9, size = 0.2),
                    axis.ticks = element_blank(),
                    line = element_line(colour = "white", size = 0.5, linetype = 9, lineend = "butt"),
                    legend.position = "bottom",
                    axis.text = element_text(size = text_size1),
                    legend.text = element_text(size = text_size1),
                    axis.title.y = element_text(size = text_size1, hjust = 0.6),
                    title = element_text(size = text_size1, face = "bold")
                ) +
                geom_histogram(aes(y = ..density..), binwidth = binwidth,  col = "grey", fill = "grey") +
                xlab("") +
                ylab("Density
                    ") +
                stat_function(
                    fun = dnorm,
                    args = list(
                        mean = mean(data$x),
                        sd = sd(data$x)
                    ),
                    col = "black",
                    size = 1, n = 200
                ) +
                geom_polygon(data = d1, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .9) +
                geom_polygon(data = d2, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .1) +
                geom_segment(data = darrow, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), arrow = arrow(length = unit(0.20, "cm"), type = "closed"), size = 0.5, color = "black") +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.53 * max(y), label = "Critical Region", size = text_size2) +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.45 * max(y), label = paste0(alpha * 100, "%"), size = text_size2) +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.37 * max(y), label = paste("alpha =", alpha), size = text_size2) +
                annotate("text", x = rand.mean, y = 1.05 * max(y), label = "E.perm", size = text_size2, color = "blue", fontface = "bold") +
                annotate("text", x = orig.ev, y = 1.05 * max(y), label = "E.obs", size = text_size2, color = "red", fontface = "bold") +
                annotate("text", x = 0.99 * rand.mean, y = 0.3 * max(y), label = paste0(100 - alpha * 100, "%"), size = text_size2, color = "black") +
                geom_segment(data = dmean, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "blue") +
                geom_segment(data = dobv, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "red")
        }

        if (test_type == "two-sided") {
            xcoords <- rand.ev
            rand.mean <- mean(xcoords, na.rm = TRUE)
            rand.sd <- sd(xcoords, na.rm = TRUE)
            xcoords <- xcoords[order(xcoords)]
            y <- dnorm(xcoords, rand.mean, sd = rand.sd)
            x.breaks <- hist(xcoords, breaks = breaks, plot = FALSE)$breaks # breaks, counts, density
            y.density <- hist(xcoords, breaks = breaks, plot = FALSE)$density # breaks, counts, density

            ymax <- max(max(y, na.rm = TRUE), max(y.density, na.rm = TRUE))

            # less
            aux <- qnorm(alpha / 2, mean = rand.mean, sd = rand.sd)
            aux2 <- qnorm(1 - alpha / 2, mean = rand.mean, sd = rand.sd)

            xmin <- min(orig.ev, min(xcoords, na.rm = TRUE), min(aux, na.rm = TRUE),
                na.rm = TRUE
            )
            xmax <- max(orig.ev, max(xcoords, na.rm = TRUE), max(aux, na.rm = TRUE),
                na.rm = TRUE
            )

            extendwidth <- max(rand.mean - xmin, xmax - rand.mean)
            xmax <- rand.mean + extendwidth
            xmin <- rand.mean - extendwidth

            if (length(binwidth) == 0) {
                binwidth <- x.breaks[2] - x.breaks[1]
            }

            dmean <- data.frame(x = rand.mean, y = 0, vx = 0, vy = ymax)
            dobv <- data.frame(x = orig.ev, y = 0, vx = 0, vy = ymax)


            # plot histgram and normal
            data <- data.frame(x = xcoords)

            # plot polygon
            polygon1.x <- seq(xmin, aux, length = 50)
            polygon1.y <- dnorm(polygon1.x, rand.mean, rand.sd)

            polygon2.x <- seq(aux2, xmax, length = 50)
            polygon2.y <- dnorm(polygon2.x, rand.mean, rand.sd)

            d1 <- data.frame(x = c(xmin, polygon1.x, aux), y = c(
                0, polygon1.y,
                0
            ))
            d2 <- data.frame(x = c(aux, aux, xmin, xmin), y = c(-0.05 *
                ymax, ymax, ymax, -0.05 * ymax))
            au2max <- aux - xmin

            # plot arrow1
            darrow <- data.frame(x = c(xmin + 0.4 * au2max, aux - 0.4 *
                au2max), y = c(0.45 * max(y), 0.45 * max(y)), vx = c(-0.4 *
                au2max, 0.4 * au2max), vy = c(0, 0))


            d3 <- data.frame(x = c(aux2, polygon2.x, xmax), y = c(
                0, polygon2.y,
                0
            ))
            d4 <- data.frame(x = c(aux2, xmax, xmax, aux2), y = c(-0.05 *
                ymax, -0.05 * ymax, ymax, ymax))
            au2max2 <- xmax - aux2

            # plot arrow2
            darrow2 <- data.frame(x = c(aux2 + 0.4 * au2max2, xmax - 0.4 *
                au2max2), y = c(0.45 * max(y), 0.45 * max(y)), vx = c(-0.4 *
                au2max2, 0.4 * au2max2), vy = c(0, 0))
            text_size1 <- 11
            text_size2 <- 4

            p1 <- ggplot(data, aes(x)) +
                labs(subtitle = paste("p-value:", pval * 2, "\n", "ntimes:", ntimes, "\n", "zscore:", zscore, "\n", "alpha:", alpha, "\n", "type:", paste(type, "test"))) +
                theme(
                    panel.background = element_blank(),
                    plot.title = element_text(hjust = 0.5),
                    plot.subtitle = element_text(hjust = 1),
                    panel.grid.major = element_line(colour = "grey", linetype = 9, size = 0.2),
                    axis.ticks = element_blank(),
                    line = element_line(colour = "white", size = 0.5, linetype = 9, lineend = "butt"),
                    legend.position = "bottom",
                    axis.text = element_text(size = text_size1),
                    legend.text = element_text(size = text_size1),
                    axis.title.y = element_text(size = text_size1, hjust = 0.6),
                    title = element_text(size = text_size1, face = "bold")
                ) +
                geom_histogram(aes(y = ..density..), binwidth = binwidth, col = "grey", fill = "grey") +
                xlab("") +
                ylab("Density
                    ") +
                stat_function(
                    fun = dnorm,
                    args = list(
                        mean = mean(data$x),
                        sd = sd(data$x)
                    ),
                    col = "black",
                    size = 1, n = 200
                ) +
                geom_polygon(data = d1, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .9) +
                geom_polygon(data = d2, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .1) +
                geom_polygon(data = d3, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .9) +
                geom_polygon(data = d4, mapping = aes(x = x, y = y), fill = "#1b98e0", alpha = .1) +
                geom_segment(data = darrow, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), arrow = arrow(length = unit(0.20, "cm"), type = "closed"), size = 0.5, color = "black") +
                geom_segment(data = darrow2, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), arrow = arrow(length = unit(0.20, "cm"), type = "closed"), size = 0.5, color = "black") +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.53 * max(y), label = "Critical Region", size = text_size2) +
                annotate("text", x = xmax - 1 / 2 * au2max, y = 0.53 * max(y), label = "Critical Region", size = text_size2) +
                annotate("text", x = xmin + 1 / 2 * au2max, y = 0.45 * max(y), label = paste0(alpha / 2 * 100, "%"), size = text_size2) +
                annotate("text", x = xmax - 1 / 2 * au2max, y = 0.45 * max(y), label = paste0(alpha / 2 * 100, "%"), size = text_size2) +
                annotate("text", x = rand.mean, y = 1.05 * max(y), label = "E.perm", size = text_size2, color = "blue", fontface = "bold") +
                annotate("text", x = orig.ev, y = 1.05 * max(y), label = "E.obs", size = text_size2, color = "red", fontface = "bold") +
                annotate("text", x = 0.99 * rand.mean, y = 0.3 * max(y), label = paste0(100 - alpha * 100, "%"), size = text_size2, color = "black") +
                geom_segment(data = dmean, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "blue") +
                geom_segment(data = dobv, mapping = aes(x = x, y = y, xend = x + vx, yend = y + vy), size = 1.3, color = "red")
        }
    }
    return(p1)
}