Vignette with example plots

Author: KatjaDanielzik

NAMESPACE

@@ -1,6 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
-export(extract_stan_fit)
+export(extract_fit)
 export(plot_volcano)
 export(prepare_volcano_df)
 import(ggplot2)

R/extract_stan_fit.R R/extract_fit.RR/extract_stan_fit.R renamed to R/extract_fit.R

@@ -1,4 +1,4 @@
-#' extract_stan_fit
+#' extract_fit
 #'
 #' Wrapper function to extract parameter draws from two common Stan interfaces.
 #' This function requires the respective stan interface (rstan, brms) 
@@ -21,7 +21,7 @@
 #'  # posterior <- extract_stan_fit(fit, "b_Intercept")
 #'  # End(Not run)
 
-extract_stan_fit <- function(fit, parameter_name) {
+extract_fit <- function(fit, parameter_name) {
   # rstan
   if (inherits(fit, "stanfit")) {
     if (!requireNamespace("rstan", quietly = TRUE)) {

R/plot_volcano.R

@@ -95,7 +95,7 @@ plot_volcano <- function(result,
   ## get threshold
   t <- result$meta$threshold
 
-  subtitle <- paste0("vertical black line: threshold for pi=",t)
+  subtitle <- paste0("vertical black line: threshold for pi = ",t)
 
   p <- ggplot(df,(aes(x=parameter.median,y=pi.value))) +
     geom_point()+
@@ -113,7 +113,7 @@ plot_volcano <- function(result,
 "errorbar: CrI ",result$meta$CrI.low,", ",result$meta$CrI.high)
 
     p <- p+
-      geom_errorbar(aes(xmin=parameter.low,xmax=parameter.high))+
+      geom_errorbar(aes(xmin=parameter.low,xmax=parameter.high),col="grey")+
       ggtitle(title, subtitle)
   }
 
@@ -149,7 +149,8 @@ plot_volcano <- function(result,
     }
 
     subtitle <- paste0(subtitle,'\n',
-                       "grey lines: label thresholds")
+                       paste0("grey lines: label thresholds, |parameter| > ",
+                              label.parameter.threshold,", pi > ",label.pi.threshold))
 
 
     p <- p+

R/prepare_volcano_df.R

@@ -5,7 +5,7 @@
 #' information (e.g., cell line, time point) based on parameter names and a user-provided
 #' annotation data frame, and returns a data frame that is ready for plotting.
 #'
-#' @param posterior A data frame of posterior draws (one row per draw) [extract_stan_fit()].
+#' @param posterior A data frame of posterior draws (one row per draw) [extract_fit()].
 #' @param annotation_df A data frame with at least one column:
 #'   \itemize{
 #'     \item \code{parameter}: the parameter name (e.g., `doubling.1`, `logOR.treatment`)
@@ -41,7 +41,7 @@
 #'   }
 #'}
 #'
-#' @seealso [extract_stan_fit()]
+#' @seealso [extract_fit()]
 #' 
 #' @import magrittr
 #' @importFrom tidyr pivot_longer

README.md

@@ -13,7 +13,7 @@ small sample sizes, and complex experimental designs.
 However, Bayesian posteriors of models with many parameters are often difficult
 to interpret at a glance.
 
-One way to quickly identify important biological changes of frequentist analysis 
+One way to quickly identify important biological changes based on frequentist analysis 
 are volcano plots (using fold-changes and p-values).
 
 Bayesian volcano plots bring together the uncertainty aware power of Bayesian 

data/annotation_df.rda

@@ -22,7 +22,8 @@ data("posterior")
 
 annotation_df <- as.data.frame(cbind(parameter=colnames(posterior),
 label=levels(as.factor(longitudinalMetabolomics_df$metabolite)),
-group=rep(c("A","B"))))
+group=rep(c("A","B")),value=rnorm(98,10,3)))
+annotation_df$value <- as.numeric(annotation_df$value)
 }
 
 \keyword{datasets}

man/annotation_df.Rd

@@ -1,6 +1,8 @@
 ---
 title: "BayesVolcano"
-author: Katja Danielzik (katja.danielzik@uni-due.de)
+package: BayesVolcano
+author: "Katja Danielzik (katja.danielzik@uni-due.de)"
+
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteIndexEntry{BayesVolcano}
@@ -15,6 +17,149 @@ knitr::opts_chunk$set(
 )
 ```
 
+# Data requirements
+
+BayesVolcano requires as input:
+- posteriors of a Bayesian model fit stored in a data frame, where columns are parameters and rows are iterations. A wrapper function to extract the posterior exists for fits generated by [rstan](https://mc-stan.org/rstan) or [brms](https://paulbuerkner.com/brms) with the helper function extract_fit.R
+
 ```{r setup}
 library(BayesVolcano)
 ```
+```{r}
+data("posterior")
+head(posterior[,1:4])
+```
+
+- a data frame mapping parameter names to biological entities with at least a column named "parameter" and a column named "label". Additional columns can be provided that might be helpful for later visualization (can be categorical (here: group) or numerical (here:value))
+
+```{r}
+data("annotation_df")
+head(annotation_df)
+```
+
+# Background
+Bayesian models are used to estimate effect sizes (e.g., gene expression changes,
+protein abundance differences, drug response effects) while accounting for uncertainty, 
+small sample sizes, and complex experimental designs.
+However, Bayesian posteriors of models with many parameters are often difficult to interpret at a glance.
+
+One way to quickly identify important biological changes based on frequentist analysis 
+are volcano plots (using fold-changes and p-values).
+
+Bayesian volcano plots bring together the uncertainty aware power of Bayesian 
+models and the familiar visualization of volcano plots.
+
+# Workflow
+## Preparation of plotting dataframe
+
+The prepare_volcano_df() function takes as input the posterior and the annotation
+data frame. For every parameter listed in the annotation data frame the function
+calculates one pi-value on basis of the following formula:
+
+With *i* being one entity that was modeled, *param* the estimated parameter
+and *t* the threshold for which to calculate the pi-value.
+
+![equation](https://latex.codecogs.com/svg.image?%5Cpi_%7Bi%7D=2*max%5Cleft(%5Cint_%7Bparam_%7Bi%7D=-%5Cinfty%7D%5E%7Bt%7Dp(param_%7Bi%7D)%5C,dparam_%7Bi%7D,%5Cint_%7Bparam_%7Bi%7D=t%7D%5E%7B%5Cinfty%7Dp(param_%7Bi%7D)%5C,dparam_%7Bi%7D%5Cright)-1)
+
+The user further can specify which credible interval (CrI) should be calculated
+from the posterior. This can later be visualized in the volcano plot.
+
+```{r}
+library(BayesVolcano)
+results <- prepare_volcano_df(posterior=posterior, 
+                   annotation_df = annotation_df, 
+                   threshold = 0.5, # for demonstration purposes, with centered parameters usually 0 is a valid choice
+                   CrI.low = 0.025, #2.5%
+                   CrI.high = 0.975) # 97.5%, in combination with CrI.low translates to 95% CrI
+
+```
+
+This returns a list with two elements:
+- a data frame with the parameter, pi.value, parameter.median, parameter.low (CrI.low),
+parameter.high (CrI.high), label (from annotation data frame) and additional
+columns (here "group") from the annotation data frame
+
+```{r}
+head(results$result)
+```
+- the meta information (set threshold, CrI.low, and CrI.high)
+```{r}
+results$meta
+```
+
+## Volcano Plot
+The plot_volcano() function is a wrapper function for convenient plotting
+of the results$results data frame with [ggplot2](https://ggplot2.tidyverse.org/)
+that automatically includes the
+meta information. It returns a ggplot object that can be further customized by
+the user.
+
+The simplest version is: 
+
+```{r}
+plot_volcano(results)
+```
+
+Where each point corresponds to one parameter of the Bayesian posterior. The
+set threshold for pi is displayed as vertical black line and described in the
+subtitle.
+
+Titles and x-axis labels can be set with:
+```{r}
+plot_volcano(results, 
+             title = "My amazing plot",
+             xlab = "My informative parameter")
+```
+
+
+To add the set CrI use:
+```{r}
+plot_volcano(results,
+             CrI = TRUE)
+```
+
+This adds the credible intervals as grey lines and adds a description of the errorbar
+to the subtitle.
+
+To color code the points the user can choose a character of numerical column
+from the annotation_df
+
+```{r}
+plot_volcano(results,
+             CrI = FALSE,
+             color = "group")
+
+plot_volcano(results,
+             CrI = FALSE,
+             color = "value")
+```
+
+To add labels the [ggrepel](https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html)
+is employed:
+
+```{r}
+plot_volcano(results, label = "label")
+```
+
+To only display labels above thresholds for the parameter value and pi set:
+
+```{r}
+plot_volcano(results, 
+             label = "label",
+             label.parameter.threshold = 1.5,
+             label.pi.threshold = 0.95)
+```
+
+Grey lines are added to visualize the threshold and a subtitle is added.
+
+Of course color and labels can be combined:
+```{r}
+plot_volcano(results, 
+             color = "group",
+             label = "label",
+             label.parameter.threshold = 1.5,
+             label.pi.threshold = 0.95,
+             title = "My plot",
+             xlab = "My informative parameter")
+```
+