From f3009d68a458dffd652e944ace0d6b1f534a68de Mon Sep 17 00:00:00 2001 From: aduvermy <arnaud.duvermy@ens-lyon.fr> Date: Tue, 24 Oct 2023 16:20:26 +0200 Subject: [PATCH] update flat_full --- dev/flat_full.Rmd | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/dev/flat_full.Rmd b/dev/flat_full.Rmd index 2130617..cb68a25 100644 --- a/dev/flat_full.Rmd +++ b/dev/flat_full.Rmd @@ -3672,10 +3672,10 @@ subset_glance <- function(glance_df, focus){ #' values include "AIC", "BIC", "logLik", "deviance", "df.resid", and #' "dispersion". If \code{NULL}, all available metrics will be plotted. #' -#' @return A ggplot object displaying density plots for the specified metrics. +#' @return A ggplot object displaying histogram plots for the specified metrics. #' #' @importFrom reshape2 melt -#' @importFrom ggplot2 aes facet_wrap geom_density theme_bw theme ggtitle +#' @importFrom ggplot2 aes facet_wrap geom_histogram theme_bw theme ggtitle #' #' @export #' @@ -3691,7 +3691,7 @@ metrics_plot <- function(l_tmb, focus = NULL) { } long_glance_df <- reshape2::melt(glance_df, variable.name = "metric") p <- ggplot2::ggplot(long_glance_df) + - ggplot2::geom_density(ggplot2::aes(x = value, col = metric, fill = metric), alpha = 0.4) + + ggplot2::geom_histogram(ggplot2::aes(x = value, col = metric, fill = metric)) + ggplot2::facet_wrap(~metric, scales = "free") + ggplot2::theme_bw() + ggplot2::theme(legend.position = 'null') + @@ -3767,8 +3767,7 @@ test_that("metrics_plot returns a ggplot object", { #' } evaluateDispersion <- function(TMB_dispersion_df, DESEQ_dispersion_df, color2use) { disp_comparison_dtf <- rbind(TMB_dispersion_df, DESEQ_dispersion_df) - disp_plot <- dispersion_plot(disp_comparison_dtf, col = "from") + ggplot2::scale_color_manual(values = color2use) - + disp_plot <- dispersion_plot(disp_comparison_dtf, col = "from", pch = "from") + ggplot2::scale_color_manual(values = color2use) return(list(disp_plot = disp_plot, data = disp_comparison_dtf)) } @@ -5611,7 +5610,7 @@ getLabelExpected <- function(comparison_df, coeff_threshold, alt_hypothesis) { #' @param ... additional params to pass ggplot2::aes #' @return A ggplot object representing the ROC curve plot. #' @importFrom plotROC geom_roc -#' @importFrom ggplot2 ggtitle theme_bw aes sym +#' @importFrom ggplot2 ggtitle theme_bw aes sym xlab ylab #' #' @examples #' comparison_data <- data.frame( @@ -5639,7 +5638,9 @@ roc_plot <- function(comparison_df, ...) { p <- ggplot2::ggplot(comparison_df, ggplot2::aes(d = !isDE , m = p.adj, !!!args )) + plotROC::geom_roc(n.cuts = 0, labels = FALSE) + ggplot2::theme_bw() + - ggplot2::ggtitle("ROC curve") + ggplot2::ggtitle("ROC curve") + + ggplot2::xlab("False positive rate") + + ggplot2::ylab("True positive rate") ## -- annotation AUC df_AUC <- subset(plotROC::calc_auc(p) , select = -c(PANEL, group)) @@ -5788,7 +5789,7 @@ identity_plot <- function(comparison_df, ...){ ggplot2::ggplot(comparison_df) + - ggplot2::geom_point(ggplot2::aes(x = actual, y = estimate, !!!args), alpha = 0.6) + + ggplot2::geom_point(ggplot2::aes(x = actual, y = estimate, !!!args), alpha = 0.6, size = 2) + ggplot2::geom_abline(intercept = 0, slope = 1, lty = 3, col = 'red', linewidth = 1) + ggplot2::facet_wrap(~description, scales = "free") + ggplot2::theme_bw() + @@ -5892,15 +5893,15 @@ getGrobTable <- function(df){ #' #' @param mock_obj A list containing simulation data and ground truth. #' @param list_tmb A list of model results. -#' @param dds_obj a DESeq2 object #' @param coeff_threshold A threshold for comparing estimates. #' @param alt_hypothesis The alternative hypothesis for comparisons ("greater", "less", "greaterAbs"). #' @param report_file File name to save the generated report. If NULL, the report will not be exported. +#' @param dds_obj a DESeq2 object #' @importFrom ggplot2 aes geom_point geom_abline facet_wrap theme_bw ggtitle #' @return A list containing settings, plots, and evaluation results. #' @export -simulationReport <- function(mock_obj, list_tmb = NULL, dds_obj = NULL , - coeff_threshold = 0, alt_hypothesis = "greaterAbs", report_file = NULL){ +simulationReport <- function(mock_obj, list_tmb = NULL, + coeff_threshold = 0, alt_hypothesis = "greaterAbs", report_file = NULL, dds_obj = NULL ){ #-- init TMB_comparison_df <- data.frame() @@ -5935,12 +5936,13 @@ simulationReport <- function(mock_obj, list_tmb = NULL, dds_obj = NULL , comparison_df <- rbind( DESEQ_comparison_df, TMB_comparison_df ) - color2use <- c("#D2B4DE", "#A2D9CE") - color2use <- color2use[c(!is.null(list_tmb), !is.null(dds_obj))] + #color2use <- c("#D2B4DE", "#A2D9CE") + color2use <- c("#500472", "#79cbb8") + color2use <- color2use[c(!is.null(dds_obj), !is.null(list_tmb))] # -- plotting roc_curve <- roc_plot(comparison_df, col = "from" ) + ggplot2::scale_color_manual(values = color2use) - id_plot <- identity_plot(comparison_df, col = "from") + ggplot2::scale_color_manual(values = color2use) + id_plot <- identity_plot(comparison_df, col = "from", pch = "from") + ggplot2::scale_color_manual(values = color2use) #metrics_plot <- metrics_plot(list_tmb) evalDisp <- evaluateDispersion(TMB_dispersion_df, DESEQ_dispersion_df, color2use ) dispersion_plot <- evalDisp$disp_plot @@ -7383,6 +7385,8 @@ l_tmb <- fitModelParallel(formula = kij ~ varB + (varB | varA), family = glmmTMB::nbinom2(link = "log"), log_file = "log.txt", n.cores = 1) +## -- output +l_tmb$gene1 ## -- evaluation resSimu <- simulationReport(mock_data, list_tmb = l_tmb, @@ -7400,6 +7404,11 @@ resSimu$roc_plot ``` +# About mixed model evaluation + +**HTRfit** offers a versatile simulation framework capable of fitting various types of models involving mixed effects, thanks to its implementation of **glmmTMB**. By combining the functionalities of `init_variable()` and `add_interaction()`, **HTRfit** enables the simulation of even the most complex experimental designs. However, it's important to note that as of now, HTRfit supports the evaluation of only *Type I* mixed models. In this context, *Type I* models are defined as those that follow the structure: `~ varA + (1 | varB)` or `~ varA + (varA | varB)`. Models not conforming to this specific form cannot be evaluated using **HTRfit's** current implementation. Nonetheless, you are welcome to extend its capabilities by implementing support for additional model types. + + ```{r development-inflate, eval=FALSE} setwd("/Users/ex_dya/Documents/LBMC/HTRfit/") #usethis::create_package(path = "/Users/ex_dya/Documents/LBMC/HTRfit/") @@ -7409,5 +7418,5 @@ usethis::use_pipe(export = TRUE) devtools::document() # Keep eval=FALSE to avoid infinite loop in case you hit the knit button # Execute in the console directly -fusen::inflate(flat_file = "dev/flat_full.Rmd", vignette_name = "HTRfit") +fusen::inflate(flat_file = "dev/flat_full.Rmd", vignette_name = "HTRfit", open_vignette = T, overwrite = T) ``` -- GitLab