add plotting functions

DESCRIPTION

@@ -10,15 +10,20 @@ Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
 Imports: 
+    bayesplot,
     cli,
     cmdstanr,
     dplyr,
     ggplot2,
+    gridExtra,
     loo,
     magrittr,
+    patchwork,
     ppclust,
     rlang,
+    scales,
     stats,
+    stringr,
     TailRank,
     tibble,
     tidyr,

NAMESPACE

@@ -9,6 +9,10 @@ export(get_clonal_peaks)
 export(get_initialization)
 export(karyo_to_int)
 export(model_selection_h)
+export(plot_elbo_h)
+export(plot_model_selection_h)
+export(plot_posterior_clocks_h)
+export(plot_posterior_weights_h)
 export(plot_timing)
 export(plot_timing_h)
 export(prepare_input_data)

R/model_selection_plot_h.R

@@ -0,0 +1,45 @@
+#' Plot the behavior of the model selection scores vs the number of clusters
+#'
+#' @description Plot the BIC, Log Likelihood, ICL scores obtatined after the fit for each number of cluster.
+#'
+#' @param model_selection_tibble a tibble with 3 scores and k_max values, one for each inference 
+#' @param best_K the integer corresponding to the best number of components
+#' 
+#' @return model_selection_plot 
+#' @export
+
+plot_model_selection_h <- function(model_selection_tibble, best_K) {
+  create_plot <- function(data, x, y, best_K, y_label, score_name) {
+    ggplot2::ggplot(data, ggplot2::aes(x = !!ggplot2::sym(x), y = !!ggplot2::sym(y))) +
+      ggplot2::geom_line(color = "steelblue", linewidth = 1) +
+      ggplot2::geom_point(color = "steelblue", size = 3) +
+      ggplot2::geom_point(data = data[data[[x]] == best_K, ],
+                          ggplot2::aes(x = !!ggplot2::sym(x), y = !!ggplot2::sym(y)), color = "firebrick", size = 4) +
+      ggplot2::labs(y = y_label, x = "Number of Clusters (K)", title = score_name) +
+      ggplot2::theme_minimal(base_size = 14) +
+      ggplot2::theme(
+        legend.position = "none",
+        plot.title = ggplot2::element_text(face = "bold", size = 14, hjust = 0.5),
+        axis.title.x = ggplot2::element_text(margin = ggplot2::margin(t = 10)),
+        axis.title.y = ggplot2::element_text(margin = ggplot2::margin(r = 10)),
+        panel.grid.minor = ggplot2::element_blank()
+      )
+  }
+
+  bic_plot <- create_plot(model_selection_tibble, "K", "BIC", best_K, "BIC", "BIC vs K")
+  log_lik_plot <- create_plot(model_selection_tibble, "K", "Log_lik", best_K, "Log-Likelihood", "Log-Likelihood vs K")
+  icl_plot <- create_plot(model_selection_tibble, "K", "ICL", best_K, "ICL", "ICL vs K")
+
+  model_selection_plot <- (bic_plot) / icl_plot / (log_lik_plot)  +
+    patchwork::plot_annotation(
+      title = "Model Selection Graphs: Scores vs Number of Clusters",
+      caption = "Source: Your Data",
+      theme = ggplot2::theme_minimal(base_size = 12) +
+        ggplot2::theme(
+          plot.title = ggplot2::element_text(size = 16, hjust = 0.5),
+          plot.caption = ggplot2::element_text(size = 10, hjust = 0.5, margin = ggplot2::margin(t = 10))
+        )
+    )
+
+  return(model_selection_plot)
+}

R/plot_elbo_h.R

@@ -0,0 +1,34 @@
+#' Plot elbo Function
+#'
+#' @description Plot the ELBO behaviour for each fit fixing the number of cluster.
+#' @param elbo_iteration data.frame':	#iterations until convergence of elbo obs. of  2 variables: iteration: int, elbo     : num  
+#' one of the element list elbo_iteration obtained from fit_h "results"
+#' 
+#' @return p 
+#' @export
+plot_elbo_h <- function(elbo_iteration){
+
+  data <- elbo_iteration
+
+  p <- ggplot2::ggplot(data, ggplot2::aes(x = .data$iteration, y = .data$elbo)) +
+    ggplot2::geom_line(color = "blue", size = 0.2) +  
+    # geom_point(color = "black", size = 1.5) +  
+    ggplot2::labs(
+      title = paste0("ELBO Convergence for K = ", nrow(data) - 1),  
+      x = "Iteration",
+      y = "ELBO"
+    ) +
+    # theme_classic(base_size = 12) +  
+    ggplot2::theme(
+      axis.title = ggplot2::element_text(size = 14), 
+      axis.text = ggplot2::element_text(size = 12), 
+      panel.grid = ggplot2::element_line(color = "gray100"),
+    ) +
+    ggplot2::scale_x_continuous(
+      breaks = seq(min(data$iteration), max(data$iteration), by = 1)  # Ensure only integers appear
+    ) +
+    ggplot2::scale_y_continuous(labels = scales::comma)  
+
+
+  return(p)
+}

R/plot_posterior_clocks_h.R

@@ -0,0 +1,30 @@
+#' Plot posterior clocks distributions obtained from the hierarchical model fit
+#'
+#' @param results list of 4: $data, $draws_and_summary, $log_lik_matrix_list and $elbo_iterations 
+#' @param K index of inference whose results want to be plotted   
+#' 
+#' @return areas_tau 
+#' @export
+#' 
+plot_posterior_clocks_h <- function(results, K){
+  draws <- results$draws_and_summary[[K]]$draws
+
+  names <- paste("tau[", 1:K, "]", sep = "")
+
+  areas_tau <- bayesplot::mcmc_areas(
+    draws,
+    pars = names,
+    prob = 0.8, # 80% intervals
+    prob_outer = 0.95, # 99%
+    point_est = "median"
+  )+
+    ggplot2::labs(
+      title = "Approximate Posterior distributions",
+      subtitle = "With median and 80% and 95% intervals"
+    )+
+    ggplot2::xlim(0, 1) 
+
+   return(areas_tau)
+}
+
+

R/plot_posterior_weights_h.R

@@ -0,0 +1,29 @@
+#' Plot posterior weights distributions obtained from the hierarchical model fit
+#'
+#' @param results list of 4: $data, $draws_and_summary, $log_lik_matrix_list and $elbo_iterations 
+#' @param K index of inference whose results want to be plotted  
+#' 
+#' @return areas_tau 
+#' @export
+
+plot_posterior_weights_h <- function(results, K){
+
+  draws <- results$draws_and_summary[[K]]$draws
+  S = nrow(results$data$accepted_cna)
+
+  intervals_weigths_per_tau <- list()
+  for (k in 1:K){
+    names_weights <- paste("w[",1:S,",", k, "]", sep = "")
+    p <- bayesplot::mcmc_intervals(draws, pars = names_weights, point_est = "median", prob = 0.8, prob_outer = 0.95)+
+      ggplot2::labs(
+        title =  stringr::str_wrap( paste0("Posterior distributions of the weigths for tau ",k), width = 30 + K + sqrt(S)),
+        subtitle = "With median and 80% and 95% intervals"
+      )
+    intervals_weigths_per_tau[[k]] <- ggplot2::ggplotGrob(p)
+  }
+  p <- gridExtra::grid.arrange(grobs = intervals_weigths_per_tau, ncol = K) #add global title
+
+  return(p)
+}
+
+