Probably-Probably Plot

draw_pp.R

# # MIT License
#
# Copyright (c) 2024 Mickaël Canouil
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# library(data.table)
# library(ggplot2)
# library(ggtext)
# devtools::source_gist("d7cbc59a7fcbb6d231f432801ec7aa19") # pval_trans

draw_pp <- function(data, y, alpha = 0.05) {
  data <- data.table::as.data.table(data)#[, .SD, .SDcols = y]
  data.table::setnames(data, y, "pvalue", skip_absent = TRUE)
  ggplot2::ggplot(
    data = data[order(pvalue)][,
      c("exppval", "labels") := list(
        (1:.N - 0.5) / .N,
        paste0(
          "&lambda;<sub>gc</sub> = ", 
          format(median(qnorm(pvalue / 2)^2, na.rm = TRUE) / qchisq(0.5, df = 1), digits = 3, nsmall = 3)
        )
      )
    ][]
  ) +
    ggplot2::aes(x = .data[["exppval"]], y = .data[["pvalue"]], colour = .data[["labels"]], shape = .data[["labels"]]) +
    ggplot2::geom_abline(intercept = 0, slope = 1, colour = "black", linetype = 2) +
    ggplot2::geom_point(size = 0.60) +
    ggplot2::scale_x_continuous(trans = "pval", expand = ggplot2::expansion(c(0, 0.2)), limits = c(1, NA)) +
    ggplot2::scale_colour_viridis_d(begin = 0.5, end = 0.5) +
    ggplot2::scale_shape_discrete(solid = TRUE) +
    ggplot2::labs(x = "Expected P-value", y = "Observed P-value", colour = NULL, shape = NULL) +
    ggplot2::theme(
      legend.position = c(0.99, 0.01),
      legend.justification = c("right", "bottom"),
      legend.box.just = "right",
      legend.text = ggtext::element_markdown(),
      legend.margin = ggplot2::margin(1.5, 1.5, 1.5, 1.5),
      legend.spacing.x = ggplot2::unit(0, "pt"),
      legend.spacing.y = ggplot2::unit(0, "pt")
    ) +
    ggplot2::geom_hline(yintercept = alpha, linetype = 2, colour = "firebrick") +
    ggplot2::scale_y_continuous(
      trans = pval_trans(alpha = alpha, md = TRUE, colour = "firebrick"), 
      expand = ggplot2::expansion(mult = c(0, 0.2)), 
      limits = c(1, NA)
    ) +
    ggplot2::guides(shape = ggplot2::guide_legend(override.aes = list(size = 2)))
}