diff --git a/DESCRIPTION b/DESCRIPTION index 7956cdb..62b06a0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -38,24 +38,34 @@ Imports: DBI, dplyr, duckdb, + circlize, + forcats, ggplot2, ggrepel, + ggridges, glmnet, glue, + graphics, + grDevices, jsonlite, hardhat, methods, parsnip, purrr, + patchwork, + RColorBrewer, readr, recipes, rlang, rsample, + scales, sgof, + stats, stringr, tibble, tidyr, tune, + utils, vip, withr, workflows, diff --git a/NAMESPACE b/NAMESPACE index 0b9d60e..789117a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,9 +19,17 @@ export(getConfusionMatrix) export(getNumFeat) export(loadMLInputTibble) export(plotBaselineComparison) -export(plotDefaultEval) +export(plotCM) +export(plotCrossDrug) +export(plotDrugDist) +export(plotDrugPerf) export(plotFishers) +export(plotMDR) export(plotPRC) +export(plotROC) +export(plotShuffleVsReal) +export(plotStratifiedPerf) +export(plotTopClusters) export(plotTopFeatsVI) export(predictML) export(removeTopFeats) @@ -87,6 +95,7 @@ importFrom(ggplot2,xlab) importFrom(ggplot2,ylim) importFrom(glmnet,glmnet) importFrom(glue,glue) +importFrom(grDevices,colorRampPalette) importFrom(graphics,barplot) importFrom(hardhat,tune) importFrom(jsonlite,fromJSON) @@ -124,7 +133,9 @@ importFrom(rsample,validation_set) importFrom(rsample,vfold_cv) importFrom(stats,coef) importFrom(stats,fisher.test) +importFrom(stats,median) importFrom(stats,reformulate) +importFrom(stats,reorder) importFrom(stringr,str_remove) importFrom(stringr,str_split) importFrom(tibble,add_column) diff --git a/R/arg_check_ml.R b/R/arg_check_ml.R index 2056323..a0b8af8 100644 --- a/R/arg_check_ml.R +++ b/R/arg_check_ml.R @@ -41,11 +41,11 @@ NULL "genome_drug.resistant_phenotype" %in% colnames(tibble), "resistant_classes" %in% colnames(tibble) )) { - stop(paste( - "The input tibble must have a target variable column named", - "either `genome_drug.resistant_phenotype` or `resistant_classes`, but", + stop( + "The input tibble must have a target variable column named ", + "either `genome_drug.resistant_phenotype` or `resistant_classes`, but ", "not both." - )) + ) } } } @@ -107,11 +107,11 @@ NULL is_train_val_test <- (train_prop > 0 && val_prop > 0 && test_prop > 0) if (!(is_pure_cv || is_train_val_test)) { - stop(paste0( + stop( "Unsupported split proportions. Use either:\n", " - pure CV: split = c(1, 0)\n", " - classical splits (train/val/test): c(train_prop, val_prop) with train > 0, val > 0, and test = 1 - train - val > 0." - )) + ) } # Sum must be <= 1 for all valid cases @@ -192,10 +192,10 @@ NULL } if (smallest_n_obs_rs < 1 || smallest_n_obs_rs %% 1 != 0) { - stop(paste( - "The `smallest_n_obs_rs` argument must be a whole, positive", + stop( + "The `smallest_n_obs_rs` argument must be a whole, positive ", "number greater than 0." - )) + ) } } @@ -207,10 +207,10 @@ NULL #' .checkArgUsePCA <- function(use_pca) { if (!is.logical(use_pca)) { - stop(paste( - "The `use_pca` argument can only take logical values:", + stop( + "The `use_pca` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -240,10 +240,10 @@ NULL #' .checkArgMultiClass <- function(multi_class) { if (!is.logical(multi_class)) { - stop(paste( - "The `multi_class` argument can only take logical values:", + stop( + "The `multi_class` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -281,10 +281,10 @@ NULL #' .checkArgModel <- function(model) { if (!(model %in% c("LR", "RF", "BT"))) { - stop(paste( - "The `model` argument must be one of:", + stop( + "The `model` argument must be one of: ", "'LR' (logistic regression), 'RF' (random forest), 'BT' (boosted tree)." - )) + ) } } @@ -434,10 +434,10 @@ NULL #' .checkArgSelectBestMetric <- function(select_best_metric) { if (!(select_best_metric %in% c("f_meas", "pr_auc", "mcc", "bal_accuracy"))) { - stop(paste( - "The `select_best_metric` argument must be one of: 'f_meas',", + stop( + "The `select_best_metric` argument must be one of: 'f_meas', ", "'pr_auc', 'mcc', 'bal_accuracy'." - )) + ) } } @@ -453,10 +453,10 @@ NULL .checkArgTibble(test_data_plus_predictions, ml = TRUE) if (!(".pred_class" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` is missing a `.pred_class`", + stop( + "`test_data_plus_predictions` is missing a `.pred_class` ", "column. Try using the output of `predictML()`." - )) + ) } } @@ -489,10 +489,10 @@ NULL #' .checkArgPropVITopFeats <- function(prop_vi_top_feats) { if (!is.numeric(prop_vi_top_feats) || length(prop_vi_top_feats) != 2) { - stop(paste( - "The `prop_vi_top_feats` argument must be a vector of length", + stop( + "The `prop_vi_top_feats` argument must be a vector of length ", "2 with numeric values." - )) + ) } if (any(prop_vi_top_feats > 1) || any(prop_vi_top_feats < 0)) { @@ -515,10 +515,10 @@ NULL #' .checkArgShuffleLabels <- function(shuffle_labels) { if (!is.logical(shuffle_labels)) { - stop(paste( - "The `shuffle_labels` argument can only take logical values:", + stop( + "The `shuffle_labels` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -530,10 +530,10 @@ NULL #' .checkArgReturnTuneRes <- function(return_tune_res) { if (!is.logical(return_tune_res)) { - stop(paste( - "The `return_tune_res` argument can only take logical values:", + stop( + "The `return_tune_res` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -545,10 +545,10 @@ NULL #' .checkArgReturnFit <- function(return_fit) { if (!is.logical(return_fit)) { - stop(paste( - "The `return_fit` argument can only take logical values:", + stop( + "The `return_fit` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -561,10 +561,10 @@ NULL #' .checkArgReturnPred <- function(return_pred) { if (!is.logical(return_pred)) { - stop(paste( - "The `return_pred` argument can only take logical values:", + stop( + "The `return_pred` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -576,10 +576,10 @@ NULL #' .checkArgVerbose <- function(verbose) { if (!is.logical(verbose)) { - stop(paste( - "The `verbose` argument can only take logical values:", + stop( + "The `verbose` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -613,10 +613,10 @@ NULL #' .checkArgByVI <- function(by_vi) { if (!is.logical(by_vi)) { - stop(paste( - "The `by_vi` argument can only take logical values:", + stop( + "The `by_vi` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -629,10 +629,10 @@ NULL #' .checkArgByNum <- function(by_num) { if (!is.logical(by_num)) { - stop(paste( - "The `by_num` argument can only take logical values:", + stop( + "The `by_num` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -657,10 +657,10 @@ NULL } if (any(diff(percent_removal_vec) <= 0)) { - stop(paste( - "The elements of `percent_removal_vec` must be in ascending", + stop( + "The elements of `percent_removal_vec` must be in ascending ", "order with no repeats." - )) + ) } } @@ -673,10 +673,10 @@ NULL #' .checkArgReturnFeats <- function(return_feats) { if (!is.logical(return_feats)) { - stop(paste( - "The `return_feats` argument can only take logical values:", + stop( + "The `return_feats` argument can only take logical values: ", "`TRUE` or `FALSE`." - )) + ) } } @@ -713,10 +713,10 @@ NULL if (!(y_default_eval %in% c("avg_f1_score", "avg_log2_apop", "avg_bal_acc", "avg_mcc", "avg_nmcc")) ) { - stop(paste( - "`y_default_eval` must be one of:", + stop( + "`y_default_eval` must be one of: ", "'avg_f1_score', 'avg_log2_apop', 'avg_bal_acc', 'avg_mcc', 'avg_nmcc'." - )) + ) } } diff --git a/R/core_ml.R b/R/core_ml.R index ec29d41..7e398f2 100644 --- a/R/core_ml.R +++ b/R/core_ml.R @@ -531,10 +531,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } f1 <- test_data_plus_predictions |> @@ -561,10 +561,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } auprc <- test_data_plus_predictions |> @@ -590,10 +590,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } auprc <- .calculateAUPRC(test_data_plus_predictions) @@ -603,15 +603,15 @@ getConfusionMatrix <- function(test_data_plus_predictions) { nrow(test_data_plus_predictions) if (prior > 0.3 && prior < 0.7) { - warning(paste( - "Classes are roughly balanced.", + warning( + "Classes are roughly balanced. ", "Calculation of log2(AUPRC/prior) may be inappropriate." - )) + ) } else if (prior >= 0.7) { - warning(paste( - "Classes are imbalanced for this model.", + warning( + "Classes are imbalanced for this model. ", "The use of the log2(AUPRC/prior) metric may be more informative in this imbalanced model." - )) + ) } log2_apop <- log2(auprc / prior) |> round(2) @@ -631,10 +631,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } bal_acc <- test_data_plus_predictions |> @@ -660,10 +660,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } sens <- test_data_plus_predictions |> @@ -690,10 +690,10 @@ getConfusionMatrix <- function(test_data_plus_predictions) { if (!("genome_drug.resistant_phenotype" %in% colnames(test_data_plus_predictions))) { - stop(paste( - "`test_data_plus_predictions` does not have a column for", + stop( + "`test_data_plus_predictions` does not have a column for ", "`genome_drug.resistant_phenotype`." - )) + ) } spec <- test_data_plus_predictions |> @@ -799,7 +799,7 @@ extractTopFeats <- function( dplyr::arrange(dplyr::desc(Importance)) if (!is.na(n_top_feats)) { - top_feats_and_VIs <- feats_arranged |> dplyr::slice(1:n_top_feats) + top_feats_and_VIs <- feats_arranged |> dplyr::slice(seq_len(n_top_feats)) } else if (any(!is.na(prop_vi_top_feats))) { cum_vi_lower <- prop_vi_top_feats[1] * sum(feats_arranged$Importance) cum_vi_upper <- prop_vi_top_feats[2] * sum(feats_arranged$Importance) @@ -818,10 +818,10 @@ extractTopFeats <- function( # Take a different approach if using multi-class (the previous code would give # a less meaningful result). if (is(fit$fit$actions$model$spec, "multinom_reg")) { - warning(paste( - "Extracting top features from a multi-class model.", + warning( + "Extracting top features from a multi-class model. ", "The `prop_vi_top_feats` and `n_top_feats` arguments do not apply." - )) + ) fit_penalty <- .getFitHps(fit)["penalty"] |> as.numeric() glmnet_fit <- parsnip::extract_fit_engine(fit) diff --git a/R/data.R b/R/data.R index 0f06f74..07cbed9 100644 --- a/R/data.R +++ b/R/data.R @@ -7,6 +7,9 @@ #' @format A tibble with 60 rows and 82 columns: `genome_id`, #' `genome_drug.resistant_phenotype`, and 80 binary feature columns. #' @source `inst/scripts/make_demo_data.R`. +#' @examples +#' data(demo_ml_tibble) +#' dim(demo_ml_tibble) "demo_ml_tibble" #' Demo LR fit @@ -15,4 +18,7 @@ #' #' @format A fitted `workflow` object (output of [fitBestModel()]). #' @source `inst/scripts/make_demo_data.R`. +#' @examples +#' data(demo_fit) +#' class(demo_fit) "demo_fit" diff --git a/R/generate_matrices_ml.R b/R/generate_matrices_ml.R index 55e04c6..1f4de04 100644 --- a/R/generate_matrices_ml.R +++ b/R/generate_matrices_ml.R @@ -944,10 +944,10 @@ skipImbalancedMatrix <- function(genome_ids, } }), prefix = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { - if (is.na(i) || i < 3) { + if (is.na(i) || i < 5) { NA_character_ } else { - paste(x[1:(i - 4)], collapse = "_") + paste(x[seq_len(i - 4)], collapse = "_") } }), drug = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { @@ -986,7 +986,7 @@ skipImbalancedMatrix <- function(genome_ids, ## -------------------------------------------------------------- ## Leave-one-drug-out with genome blocking ## -------------------------------------------------------------- - purrr::walk(names(grouped), function(leave_one_out) { + new_files <- purrr::map(names(grouped), function(leave_one_out) { ## Read held-out drug to get its genome IDs heldout_tbl <- arrow::read_parquet(grouped[[leave_one_out]]$file) @@ -1036,7 +1036,6 @@ skipImbalancedMatrix <- function(genome_ids, )) arrow::write_parquet(combined, out_file) - created <<- c(created, out_file) log( "debug", @@ -1045,7 +1044,13 @@ skipImbalancedMatrix <- function(genome_ids, " | removed ", length(heldout_genomes), " genomes" ) ) - }) + + out_file + }) |> + purrr::compact() |> + unlist() + + created <- c(created, new_files) } log("info", "All LOO-drug matrices generated with genome-level blocking") @@ -1105,10 +1110,10 @@ skipImbalancedMatrix <- function(genome_ids, } }), prefix = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { - if (is.na(i) || i < 3) { + if (is.na(i) || i < 5) { NA_character_ } else { - paste(x[1:(i - 4)], collapse = "_") + paste(x[seq_len(i - 4)], collapse = "_") } }), drug = purrr::pmap_chr(list(parts, idx_sparse), \(x, i) { diff --git a/R/globals.R b/R/globals.R index ef6b85e..3c440fe 100644 --- a/R/globals.R +++ b/R/globals.R @@ -4,6 +4,9 @@ # Package-level imports for packages used with :: notation #' @importFrom jsonlite fromJSON write_json #' @importFrom glmnet glmnet +#' @importFrom grDevices colorRampPalette +#' @importFrom stats median +#' @importFrom stats reorder #' @keywords internal "_PACKAGE" @@ -22,10 +25,19 @@ utils::globalVariables(c( "cum_imp", # Data columns + "abs_imp", "adj_p_value", "antibiotic", "bal_acc", + "category", + "class_abbr", + "cluster", + "contribution", + "diff_top2", "drug", + "drug_abbr", + "drug_class", + "drug_label", "drug_or_class", "feature", "feature_id", @@ -36,6 +48,7 @@ utils::globalVariables(c( "fit_mtry", "fit_trees", "gene", + "genome_drug.antibiotic", "genome_drug.genome_id", "genome_drug.resistant_phenotype", "genome_id", @@ -43,9 +56,13 @@ utils::globalVariables(c( "i_strat", "idx_sparse", "idx_strat", + "label", "model", "neg_log10_adj_p", "mcc", + "mean_margin", + "median_mcc", + "n_feat_types", "nmcc", "num_obs", "seed", @@ -56,6 +73,7 @@ utils::globalVariables(c( "pair_id", "parts", "phenotype", + "proteinName", "precision", "prefix", "prefix_key", @@ -66,14 +84,22 @@ utils::globalVariables(c( "res_prop", "resistant_classes", "run_time_sec", + "shuffled", + "shuffled_label", "sig_after_bh", "significance", "strat_value", "strat_value_test", "stratification", + "test_country", "test_drug", "test_file", + "test_year", + "tested_on", + "total", + "train_country", "train_prop", + "train_year", "value", # Base R functions that need explicit import diff --git a/R/ife_ml.R b/R/ife_ml.R index e6e38c0..f94a2a3 100644 --- a/R/ife_ml.R +++ b/R/ife_ml.R @@ -34,10 +34,10 @@ removeTopFeats <- function(ml_input_tibble, top_feat_tibble) { feats_to_remove <- top_feat_tibble |> dplyr::pull(Variable) if (length(feats_to_remove) == 0) { - stop(paste( - "No more feats to remove; try adjusting `percent_removal_vec`", + stop( + "No more feats to remove; try adjusting `percent_removal_vec` ", "in `runIFE()` to include a wider range between percentages." - )) + ) } ml_input_tibble_top_feats_removed <- ml_input_tibble |> @@ -82,7 +82,7 @@ removeTopFeats <- function(ml_input_tibble, top_feat_tibble) { #' @export runIFE <- function( ml_input_tibble, by_num = TRUE, by_vi = FALSE, - percent_removal_vec = 10 * 1:9, mix_vec = 0, return_feats = FALSE, + percent_removal_vec = 10 * seq_len(9), mix_vec = 0, return_feats = FALSE, verbose = TRUE ) { .checkArgByNum(by_num) @@ -92,10 +92,10 @@ runIFE <- function( .checkArgVerbose(verbose) if (sum(c(by_num, by_vi)) != 1) { - stop(paste( - "Set either `by_num` or `by_vi` to `TRUE` and the other to", + stop( + "Set either `by_num` or `by_vi` to `TRUE` and the other to ", "`FALSE` (not both `TRUE` or both `FALSE`)." - )) + ) } # Create empty objects to store results. @@ -119,7 +119,7 @@ runIFE <- function( message("Training with all features") } - for (i in 1:(length(percent_removal_vec) + 1)) { + for (i in seq_len(length(percent_removal_vec) + 1)) { if (i == 1) { n_feats_removed <- 0 } @@ -227,19 +227,19 @@ runIFE <- function( # Removal of top features if (by_vi) { if (verbose) { - message(paste("Cumulatively removed top ", - percent_removal_vec[i], "% of total variable importance", - sep = "" - )) + message( + "Cumulatively removed top ", + percent_removal_vec[i], "% of total variable importance" + ) } ml_input_tibble <- ml_input_tibble |> removeTopFeats(feats_to_remove_tibble) } else if (by_num) { if (verbose) { - message(paste("Cumulatively removed top ", - percent_removal_vec[i], "% of all features", - sep = "" - )) + message( + "Cumulatively removed top ", + percent_removal_vec[i], "% of all features" + ) } feats_to_remove_tibble <- ml_res$top_feat_tibble @@ -264,7 +264,7 @@ runIFE <- function( # Adjust `percent_removal_vec` to appropriate length (if IFE stopped running # early due to encountering an error). - percent_removal_vec <- percent_removal_vec[1:(length(num_obs_vec) - 1)] + percent_removal_vec <- percent_removal_vec[seq_len(length(num_obs_vec) - 1)] ife_performance_tibble <- tibble::tibble( percent_removed = c(0, percent_removal_vec), diff --git a/R/plot_ml.R b/R/plot_ml.R index a986394..3fde725 100644 --- a/R/plot_ml.R +++ b/R/plot_ml.R @@ -16,41 +16,44 @@ #' @importFrom ggplot2 theme #' @importFrom ggplot2 xlab #' @importFrom ggplot2 ylim +#' @importFrom graphics barplot #' @importFrom tune extract_fit_parsnip #' @importFrom vip vip #' @importFrom yardstick pr_curve -#' @importFrom graphics barplot NULL -#' plotPRC() +#' Plot a Precision-Recall Curve +#' +#' Generates a precision-recall curve (PRC) for AMR phenotype prediction results. +#' @param test_data_plus_predictions A tibble of test data with added prediction +#' columns (e.g. the output of `predictML()` or `runMLmodels(return_pred=TRUE)`), +#' or a path to a TSV file containing the same. #' -#' Plots the precision-recall curve given a set of test data plus predicted AMR -#' phenotypes. +#' @return A `ggplot2` object showing the precision-recall curve. +#' +#' @details +#' The function uses `yardstick::pr_curve()` to compute the PR curve and then +#' visualizes it using `ggplot2`. #' -#' @param test_data_plus_predictions Test data (tibble) with an added column for -#' predicted phenotype labels, such as the output of `predict()`. -#' @return A precision-recall curve as a `ggplot2` object #' @examples -#' preds <- tibble::tibble( -#' genome_id = paste0("g", 1:10), -#' genome_drug.resistant_phenotype = factor( -#' rep(c("Resistant", "Susceptible"), each = 5), -#' levels = c("Resistant", "Susceptible") -#' ), -#' .pred_class = factor( -#' c( -#' "Resistant", "Resistant", "Susceptible", "Resistant", "Susceptible", -#' "Susceptible", "Resistant", "Susceptible", "Susceptible", "Resistant" -#' ), -#' levels = c("Resistant", "Susceptible") -#' ), -#' .pred_Resistant = c(0.9, 0.8, 0.4, 0.7, 0.3, 0.2, 0.6, 0.1, 0.2, 0.55), -#' .pred_Susceptible = c(0.1, 0.2, 0.6, 0.3, 0.7, 0.8, 0.4, 0.9, 0.8, 0.45) -#' ) +#' data(demo_fit) +#' data(demo_ml_tibble) +#' preds <- predictML(demo_fit, demo_ml_tibble) #' plotPRC(preds) +#' #' @export plotPRC <- function(test_data_plus_predictions) { + if (is.character(test_data_plus_predictions)) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions) + } .checkArgTestDataPlusPredictions(test_data_plus_predictions) + test_data_plus_predictions <- test_data_plus_predictions |> + dplyr::mutate( + genome_drug.resistant_phenotype = factor( + genome_drug.resistant_phenotype, + levels = c("Resistant", "Susceptible") + ) + ) prc <- yardstick::pr_curve( test_data_plus_predictions, @@ -64,120 +67,160 @@ plotPRC <- function(test_data_plus_predictions) { return(prc) } -#' plotTopFeatsVI() +#' Plot a Receiver Operating Characteristic (ROC) Curve +#' +#' Generates a ROC curve for AMR phenotype prediction results. +#' +#' @param test_data_plus_predictions A tibble of test data with prediction +#' columns (e.g. the output of `predictML()` or `runMLmodels(return_pred=TRUE)`), +#' or a path to a TSV file containing the same. #' -#' Generates a plot showing the top features and their variable importance -#' scores. +#' @return A ROC curve plotted using `ggplot2::autoplot()`. #' -#' @param fit Best model fit, such as the output of `fitBestModel()` -#' @param n_top_feats [num] Number of top features to plot -#' @return Variable importance plot (a `ggplot2` object) #' @examples #' data(demo_fit) -#' plotTopFeatsVI(demo_fit, n_top_feats = 10) +#' data(demo_ml_tibble) +#' preds <- predictML(demo_fit, demo_ml_tibble) +#' plotROC(preds) +#' #' @export -plotTopFeatsVI <- function(fit, n_top_feats = 10) { - .checkArgWflow(fit) - .checkArgNTopFeats(n_top_feats) +plotROC <- function(test_data_plus_predictions) { + if (is.character(test_data_plus_predictions)) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions) + } + .checkArgTestDataPlusPredictions(test_data_plus_predictions) + test_data_plus_predictions <- test_data_plus_predictions |> + dplyr::mutate( + genome_drug.resistant_phenotype = factor( + genome_drug.resistant_phenotype, + levels = c("Resistant", "Susceptible") + ) + ) - vip <- fit |> - tune::extract_fit_parsnip() |> - vip::vip(num_features = n_top_feats) + - ggplot2::xlab("Top features") + + roc <- yardstick::roc_curve( + test_data_plus_predictions, + genome_drug.resistant_phenotype, .pred_Resistant + ) |> + ggplot2::autoplot(type = "se") + ggplot2::theme(panel.grid = ggplot2::element_blank()) - return(vip) + return(roc) } -#' plotDefaultEval() -#' -#' Plots performance metric or runtime vs. training data proportion or number -#' of cross-validation folds, colored by model. -#' -#' @param default_eval_tibble Output of `findOptimalMLDefaults()` -#' @param x_default_eval [chr] x value of default evaluation plot: "train_prop" -#' or "n_fold" -#' @param y_default_eval [chr] y value of default evaluation plot. It can be -#' "avg_runtime_sec" or one of the following performance metrics: -#' "avg_f1_score", "avg_log2_apop", "avg_bal_acc", "avg_mcc", or "avg_nmcc" -#' @param xlab [chr] Label for x axis -#' @param ylab [chr] Label for y axis -#' @return A `ggplot2` scatterplot (performance metric or runtime vs. -#' `train_prop` or `n_fold`), colored by model +#' Plot a Confusion Matrix Heatmap +#' +#' Produces a heatmap visualization of the confusion matrix for AMR predictions. +#' +#' @param test_data_plus_predictions A tibble containing true and predicted +#' phenotype labels (e.g. the output of `predictML()`), or a path to a TSV file +#' containing the same. +#' +#' @return A heatmap (`ggplot2` object) showing the confusion matrix. +#' #' @examples -#' default_eval <- tibble::tibble( -#' train_prop = c(0.5, 0.6, 0.7, 0.5, 0.6, 0.7), -#' n_fold = rep(5, 6), -#' model = rep(c("LR", "RF"), each = 3), -#' avg_f1_score = c(0.72, 0.78, 0.83, 0.70, 0.75, 0.80) -#' ) -#' plotDefaultEval(default_eval, -#' x_default_eval = "train_prop", -#' y_default_eval = "avg_f1_score" -#' ) +#' data(demo_fit) +#' data(demo_ml_tibble) +#' preds <- predictML(demo_fit, demo_ml_tibble) +#' plotCM(preds) +#' #' @export -plotDefaultEval <- function( - default_eval_tibble, x_default_eval = "train_prop", - y_default_eval = "avg_f1_score", xlab = "Train data proportion", - ylab = "Average F1 score" -) { - .checkArgTibble(default_eval_tibble) - .checkArgXDefaultEval(x_default_eval) - .checkArgYDefaultEval(y_default_eval) - .checkArgXYLabs(xlab = xlab, ylab = ylab) - - if (x_default_eval == "n_fold") { - default_eval_tibble <- default_eval_tibble |> - dplyr::filter(train_prop == 0.8) - } else { - default_eval_tibble <- default_eval_tibble |> - dplyr::filter(train_prop != 0.8) +plotCM <- function(test_data_plus_predictions) { + if (is.character(test_data_plus_predictions)) { + test_data_plus_predictions <- readr::read_tsv(test_data_plus_predictions) } - - default_eval_plot <- ggplot2::ggplot( - default_eval_tibble, - ggplot2::aes( - x = unlist(default_eval_tibble[x_default_eval]), - y = unlist(default_eval_tibble[y_default_eval]), color = model + .checkArgTestDataPlusPredictions(test_data_plus_predictions) + test_data_plus_predictions <- test_data_plus_predictions |> + dplyr::mutate( + genome_drug.resistant_phenotype = factor( + genome_drug.resistant_phenotype, + levels = c("Resistant", "Susceptible") + ), + .pred_class = factor( + .pred_class, + levels = c("Resistant", "Susceptible") + ) ) - ) + - ggplot2::geom_line(size = 1.5) + - ggplot2::geom_point(size = 3) + - ggplot2::theme( - axis.line = ggplot2::element_line(linewidth = 1.5), - axis.ticks = ggplot2::element_line(linewidth = 1.5, colour = "black"), - axis.text = ggplot2::element_text(size = 16, colour = "black"), - axis.title = ggplot2::element_text(size = 16, face = "bold"), - panel.grid = ggplot2::element_blank(), - panel.background = ggplot2::element_blank(), - legend.text = ggplot2::element_text(size = 16), - legend.title = ggplot2::element_text(size = 16, face = "bold") + test_data_plus_predictions |> + yardstick::conf_mat( + truth = genome_drug.resistant_phenotype, + estimate = .pred_class + ) |> + ggplot2::autoplot(type = "heatmap") +} + +#' Plot Top Feature Importances +#' +#' Creates a bar plot showing the most important features affecting +#' AMR phenotype predictions. +#' +#' @param topfeat A tibble of feature importance scores (e.g. the output of +#' `extractTopFeats()` or `runMLmodels()`), with `Variable`, `Importance`, and +#' `Sign` columns, or a path to a TSV file containing the same. +#' @param n_top_feats Number of top features to display (default: 10). +#' +#' @return A bar plot of variable importance (`ggplot2` object). +#' +#' @examples +#' data(demo_fit) +#' top_feats <- extractTopFeats(demo_fit, n_top_feats = 10) +#' plotTopFeatsVI(top_feats, n_top_feats = 10) +#' +#' @export +plotTopFeatsVI <- function(topfeat, n_top_feats = 10) { + if (is.character(topfeat)) { + topfeat <- readr::read_tsv(topfeat) + } + .checkArgNTopFeats(n_top_feats) + + vip <- topfeat |> + dplyr::slice_max(order_by = Importance, n = n_top_feats) |> + dplyr::mutate( + Variable = factor(Variable, levels = rev(Variable)), # preserve order as shown in table + Sign = factor(Sign, levels = c("POS", "NEG")) + ) |> + ggplot2::ggplot(ggplot2::aes(x = Importance, y = Variable, fill = Sign)) + + ggplot2::geom_col() + + ggplot2::scale_fill_manual( + values = c( + "POS" = "#c6d8d3", + "NEG" = "#f6c9a1" + ) ) + - ggplot2::labs(x = xlab, y = ylab, color = "Model") + ggplot2::labs( + x = "Importance", + y = "Features" + ) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::theme( + panel.grid.minor = ggplot2::element_blank(), + axis.text.y = ggplot2::element_text(size = 10) + ) - return(default_eval_plot) + return(vip) } -#' plotBaselineComparison() +#' Compare Baseline Performance With and Without Shuffled Labels #' -#' Generates a bar plot that compares model performance with and without -#' randomly shuffled AMR phenotype labels. +#' Produces a bar plot comparing balanced accuracy for each antibiotic using +#' true AMR labels vs. randomly shuffled labels. #' -#' @importFrom graphics barplot +#' @param non_shuffled_label_results Output of `runMLPipeline(shuffle_labels = FALSE)` +#' @param shuffled_label_results Output of `runMLPipeline(shuffle_labels = TRUE)` +#' +#' @return A barplot comparing balanced accuracy across models. #' -#' @param non_shuffled_label_results Output of `runMLPipeline()` -#' (`shuffle_labels = FALSE`) -#' @param shuffled_label_results Output of `runMLPipeline()` -#' (`shuffle_labels = TRUE`) -#' @return A bar plot with balanced accuracy comparisons per antibiotic #' @examples -#' non_shuffled <- tibble::tibble( -#' antibiotic = c("AMP", "CIP", "CRO"), -#' bal_acc = c(0.88, 0.81, 0.92) +#' non_shuffled <- list( +#' performance_tibble = tibble::tibble( +#' antibiotic = c("AMP", "CIP", "CRO"), +#' bal_acc = c(0.88, 0.81, 0.92) +#' ) #' ) -#' shuffled <- tibble::tibble( -#' antibiotic = c("AMP", "CIP", "CRO"), -#' bal_acc = c(0.52, 0.49, 0.55) +#' shuffled <- list( +#' performance_tibble = tibble::tibble( +#' antibiotic = c("AMP", "CIP", "CRO"), +#' bal_acc = c(0.52, 0.49, 0.55) +#' ) #' ) #' plotBaselineComparison(non_shuffled, shuffled) #' @export @@ -185,18 +228,14 @@ plotBaselineComparison <- function( non_shuffled_label_results, shuffled_label_results ) { - .checkArgTibble(non_shuffled_label_results) - .checkArgTibble(shuffled_label_results) + .checkArgTibble(non_shuffled_label_results$performance_tibble) + .checkArgTibble(shuffled_label_results$performance_tibble) - drugs <- non_shuffled_label_results |> - dplyr::select(antibiotic) |> - dplyr::pull() - - non_shuffled_bal_acc <- non_shuffled_label_results |> + non_shuffled_bal_acc <- non_shuffled_label_results$performance_tibble |> dplyr::select(bal_acc) |> dplyr::pull() - shuffled_bal_acc <- shuffled_label_results |> + shuffled_bal_acc <- shuffled_label_results$performance_tibble |> dplyr::select(bal_acc) |> dplyr::pull() @@ -204,13 +243,12 @@ plotBaselineComparison <- function( nrow = 2, byrow = TRUE ) - colnames(bal_acc_matrix) <- drugs - rownames(bal_acc_matrix) <- c("Non-shuffled labels", "Shuffled labels") + rownames(bal_acc_matrix) <- c("Non-Shuffled Labels", "Shuffled Labels") baseline_comparison_barplot <- barplot(bal_acc_matrix, beside = TRUE, legend.text = TRUE, col = c("skyblue", "lightpink"), - ylab = "Balanced accuracy", xlab = "Antibiotic" + ylab = "Balanced accuracy" ) return(baseline_comparison_barplot) @@ -275,8 +313,8 @@ plotFishers <- function( stop("Missing required columns: ", paste(missing_cols, collapse = ", ")) } - plot_df <- fisher_df %>% - dplyr::arrange(adj_p_value) %>% + plot_df <- fisher_df |> + dplyr::arrange(adj_p_value) |> dplyr::mutate( rank = dplyr::row_number(), neg_log10_adj_p = -log10(adj_p_value), @@ -308,7 +346,7 @@ plotFishers <- function( ) if (label_top_n > 0) { - label_df <- plot_df %>% + label_df <- plot_df |> dplyr::slice_head(n = label_top_n) p <- p + @@ -325,3 +363,936 @@ plotFishers <- function( return(p) } + +#' Plot drug phenotype distribution +#' +#' Reads metadata and generates a stacked bar plot showing counts of resistant +#' and susceptible phenotypes per antibiotic. +#' +#' @param metadata_path Character. Path to directory containing `metadata.parquet`. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' plotDrugDist(metadata_path = system.file("extdata", package = "amRml")) +plotDrugDist <- function(metadata_path = ".") { + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) + + ##################### phenotype distribution (drugs) ######################### + drug_dist <- metadata |> + dplyr::distinct( + genome.genome_id, + genome_drug.antibiotic, + drug_abbr, + genome_drug.resistant_phenotype + ) |> + dplyr::count( + genome_drug.antibiotic, + drug_abbr, + genome_drug.resistant_phenotype + ) |> + dplyr::group_by(genome_drug.antibiotic, drug_abbr) |> + dplyr::mutate(total = sum(n)) |> + dplyr::ungroup() |> + dplyr::mutate( + label = paste0(genome_drug.antibiotic, " (", drug_abbr, ")"), + label = forcats::fct_reorder(label, total) + ) + + p <- ggplot2::ggplot( + drug_dist, + ggplot2::aes( + x = label, + y = n, + fill = genome_drug.resistant_phenotype + ) + ) + + ggplot2::geom_col(color = "black", width = 0.8) + + ggplot2::coord_flip() + + ggplot2::scale_fill_manual( + values = c( + "Resistant" = "#d4872a", + "Susceptible" = "#5b8db8" + ), + name = "Phenotype" + ) + + ggplot2::labs( + x = "Antibiotic", + y = "Number of unique genomes" + ) + + ggplot2::theme_classic(base_size = 14) + + p +} + +#' Plot drug-level model performance +#' +#' Generates heatmaps and ridge plots summarizing model performance (MCC) +#' across drugs and feature types. +#' +#' @param metadata_path Character. Path to a directory containing `metadata.parquet`. +#' @param performance_path A performance tibble (with `drug_label`, `shuffled`, +#' `drug_or_class`, `feature_type`, `feature_subtype`, and `mcc` columns), or a +#' directory path containing `all_perf.parquet`. +#' +#' @return A patchwork ggplot object combining multiple panels. +#' @export +#' +#' @examples +#' # Several models per drug x feature type so the ridge densities can be drawn. +#' performance <- tidyr::expand_grid( +#' drug_or_class = c("AMP", "CIP", "NAL"), +#' feature_type = c("genes", "proteins"), +#' feature_subtype = c("binary", "counts"), +#' replicate = 1:8 +#' ) +#' performance$drug_label <- "drug" +#' performance$shuffled <- FALSE +#' performance$mcc <- 0.6 + 0.3 * sin(seq_len(nrow(performance))) +#' plotDrugPerf( +#' metadata_path = system.file("extdata", package = "amRml"), +#' performance_path = performance +#' ) +plotDrugPerf <- function(metadata_path = ".", performance_path = ".") { + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) + + if (!is.data.frame(performance_path)) { + performance_path <- arrow::read_parquet( + file.path(performance_path, "all_perf.parquet") + ) + } + performance <- performance_path + + plot_df <- metadata |> + dplyr::distinct(genome.genome_id, genome_drug.antibiotic, drug_abbr) |> + dplyr::count(genome_drug.antibiotic, drug_abbr, name = "total") + + ######################## drug performances ################################# + median_drug <- performance |> + dplyr::filter( + drug_label == "drug", + !shuffled # keep real models; remove if you want both + ) |> + dplyr::group_by(drug_or_class, feature_type, feature_subtype) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> + dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) + + drug_p1 <- ggplot2::ggplot( + median_drug, + ggplot2::aes( + x = feature_type, + y = drug_or_class, + fill = median_mcc + ) + ) + + ggplot2::geom_tile(color = "grey90", width = 0.9) + + ggplot2::scale_fill_gradientn( + colors = c( + "#EAF2FF", # very light blue + "#5F84C9", # medium blue + "#0F2A5A" # dark blue + ), + values = scales::rescale(c(-1, 0, 1)), + breaks = c(-1, -0.5, 0, 0.5, 1), + labels = scales::label_number(accuracy = 0.01), + name = "Median MCC" +) + + # ggplot2::scale_fill_gradientn( + # colors = c( + # "white", + # "#BDD7E7", + # "#6BAED6", + # "#2171B5", + # "#08306B" + # ), + # values = scales::rescale(c(0, 0.4, 0.6, 0.8, 1)), + # labels = scales::label_number(accuracy = 0.1), + # name = "Median MCC" + # ) + + ggplot2::labs(x = "Feature type") + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme( + axis.text = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(size = 12), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "bottom" + ) + + ggplot2::coord_fixed() + + drug_p1 + + median_feature <- performance |> + dplyr::filter( + drug_label == "drug", + !shuffled + ) |> + dplyr::group_by(drug_or_class, feature_type) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::left_join(plot_df, by = c("drug_or_class" = "drug_abbr")) |> + dplyr::mutate(drug_or_class = reorder(drug_or_class, total)) + + + feat_pal <- c( + "args" = "#56B4E9", # sky blue + "cogs" = "#E69F00", # orange + "genes" = "#009E73", # bluish green + "domains" = "#F0E442", # yellow + "proteins" = "#CC79A7", # reddish purple + "struct" = "#D55E00" # vermillion + ) + + rc_perf <- ggplot2::ggplot( + median_feature |> + dplyr::distinct( + drug_or_class, + feature_type, median_mcc + ), + ggplot2::aes(x = median_mcc, y = drug_or_class) + ) + + ggridges::geom_density_ridges( + scale = 0.75, + rel_min_height = 0.01, + alpha = 0.4, + fill = "grey90", + colour = "grey70" + ) + + ggplot2::geom_point( + position = ggplot2::position_jitter(height = 0.1), + size = 1.5, + alpha = 0.8, + ggplot2::aes(color = feature_type) + ) + + ggplot2::scale_color_manual(values = feat_pal, name = "Feature type") + + ggplot2::stat_summary( + fun = median, + geom = "point", + size = 2, + color = "black" + ) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::theme( + axis.text = ggplot2::element_text(size = 10, colour = "black"), + axis.title = ggplot2::element_text(size = 12), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1), + legend.position = "right", + panel.grid.minor = ggplot2::element_blank(), + axis.line = ggplot2::element_line(color = "black") + ) + + rc_perf + + final_plot <- patchwork::wrap_plots( + drug_p1, rc_perf, + widths = c(2, 2), # adjust proportions + guides = "collect" + ) & + ggplot2::theme( + legend.position = "bottom" + ) + + final_plot +} + +#' Plot cross-drug prediction as a heatmap +#' +#' Creates a heatmap showing cross-drug model performance (MCC), where models +#' trained on one drug are evaluated on another. +#' +#' @param cross_test_performance_path A cross-drug performance tibble (with +#' `drug`, `test_drug`, and `mcc` columns), or a directory path +#' containing `cross_drug_perf.parquet`. +#' @param drug_performance_path A performance tibble (with `drug_label`, +#' `drug_or_class`, and `mcc` columns), or a directory path containing +#' `all_perf.parquet`. +#' @param metadata_path Character. Path to a directory containing `metadata.parquet`. +#' +#' @return A ComplexHeatmap object. +#' @export +#' +#' @examples +#' cross_drug <- tibble::tibble( +#' drug = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), +#' test_drug = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), +#' mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) +#' ) +#' performance <- tibble::tibble( +#' drug_label = "drug", +#' drug_or_class = c("AMP", "CIP", "NAL"), +#' mcc = c(0.8, 0.7, 0.6) +#' ) +#' plotCrossDrug( +#' cross_test_performance_path = cross_drug, +#' drug_performance_path = performance, +#' metadata_path = system.file("extdata", package = "amRml") +#' ) +plotCrossDrug <- function( + cross_test_performance_path = ".", + drug_performance_path = ".", + metadata_path = "." +) { + if (!is.data.frame(cross_test_performance_path)) { + cross_test_performance_path <- arrow::read_parquet( + file.path(cross_test_performance_path, "cross_drug_perf.parquet") + ) + } + cross_drug <- cross_test_performance_path + if (!is.data.frame(drug_performance_path)) { + drug_performance_path <- arrow::read_parquet( + file.path(drug_performance_path, "all_perf.parquet") + ) + } + performance <- drug_performance_path + metadata <- arrow::read_parquet(file.path(metadata_path, "metadata.parquet")) + + + heatmap_df <- cross_drug |> + dplyr::filter(!is.na(drug), !is.na(test_drug)) |> + # dplyr::filter(test_drug %in% (cross_drug |> dplyr::pull(drug))) |> + dplyr::group_by(drug, test_drug) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") + + same_drugs <- performance |> + dplyr::filter( + drug_label == "drug", + drug_or_class %in% (cross_drug |> + dplyr::distinct(drug) |> + dplyr::pull()) + ) |> + dplyr::group_by(drug_or_class) |> + dplyr::summarise(median_mcc = median(mcc, na.rm = TRUE), .groups = "drop") |> + dplyr::mutate(test_drug = drug_or_class) |> + dplyr::distinct(drug_or_class, test_drug, median_mcc) |> + dplyr::rename(drug = drug_or_class) + + heatmap_df <- heatmap_df |> + dplyr::add_row(same_drugs) |> + dplyr::left_join( + metadata |> + dplyr::distinct(drug_abbr, class_abbr), + by = c("drug" = "drug_abbr") + ) |> + dplyr::rename("drug_class" = "class_abbr") |> + dplyr::left_join( + metadata |> + dplyr::distinct(drug_abbr, class_abbr), + by = c("test_drug" = "drug_abbr") + ) + + # Row annotation (already similar to what you did) + annotation_row <- heatmap_df |> + dplyr::distinct(drug, drug_class) |> + tibble::column_to_rownames("drug") + + # Column annotation + annotation_col <- heatmap_df |> + dplyr::distinct(test_drug, class_abbr) |> + tibble::column_to_rownames("test_drug") + + mat <- heatmap_df |> + dplyr::select(drug, test_drug, median_mcc) |> + tidyr::pivot_wider(names_from = test_drug, values_from = median_mcc) |> + tibble::column_to_rownames("drug") |> + as.matrix() + + row_order <- heatmap_df |> + dplyr::distinct(drug, drug_class) |> + dplyr::arrange(drug_class, drug) |> + dplyr::pull(drug) + + col_order <- heatmap_df |> + dplyr::distinct(test_drug, class_abbr) |> + dplyr::arrange(class_abbr, test_drug) |> + dplyr::pull(test_drug) + + # mat[is.na(mat)] <- 0 + mat <- mat[row_order, col_order] + + # Align annotations + annotation_row <- annotation_row[row_order, , drop = FALSE] + annotation_col <- annotation_col[col_order, , drop = FALSE] + + + # Collect all classes from both row and column + classes <- base::union( + annotation_row$drug_class, + annotation_col$class_abbr + ) + + # Create ONE named color vector + # class_colors <- stats::setNames( + # scales::hue_pal()(length(classes)), + # classes + # ) + + class_colors <- stats::setNames( + colorRampPalette(RColorBrewer::brewer.pal(9, "Pastel1"))(length(classes)), + classes + ) + + + heat_colors <- colorRampPalette(RColorBrewer::brewer.pal(11, "RdBu"))(100) + max_val <- max(abs(mat), na.rm = TRUE) + + # ---- Convert annotations ---- + ha_row <- ComplexHeatmap::rowAnnotation( + drug_class = annotation_row$drug_class, + col = list(drug_class = class_colors), + show_annotation_name = FALSE, + show_legend = FALSE + ) + + ha_col <- ComplexHeatmap::HeatmapAnnotation( + class_abbr = annotation_col$class_abbr, + col = list(class_abbr = class_colors), + show_annotation_name = FALSE, na_col = "grey3", + annotation_legend_param = list(labels_gp = grid::gpar(fontsize = 10)) + ) + + # ---- Color function (instead of breaks + palette) ---- + col_fun <- circlize::colorRamp2( + seq(-max_val, max_val, length.out = length(heat_colors)), + heat_colors + ) + # ---- Heatmap ---- + cross_drug_hm <- ComplexHeatmap::Heatmap( + mat, + name = "median_mcc", + col = col_fun, + cluster_rows = FALSE, + cluster_columns = FALSE, + row_order = row_order, + column_order = col_order, + left_annotation = ha_row, + top_annotation = ha_col, + width = grid::unit(ncol(mat) * 1.5, "cm"), + height = grid::unit(nrow(mat) * 1.5, "cm"), + show_row_names = TRUE, + show_column_names = TRUE, + column_title = "tested on", + row_title = "trained on", + column_title_side = "bottom", + row_title_side = "right", + row_names_gp = grid::gpar(fontsize = 14), + column_names_gp = grid::gpar(fontsize = 14), + column_names_rot = 45, + + # remove borders like pheatmap + rect_gp = grid::gpar(col = NA), + + # legends + show_heatmap_legend = TRUE, + use_raster = FALSE + ) + + cross_drug_hm +} + +#' Plot stratified model performance +#' +#' Visualizes model performance (MCC) stratified by year or country, +#' comparing within-group vs cross-group evaluation. +#' +#' @param year_or_country Character. Either "year" or "country". +#' @param stratified_performance_path Character. Path to stratified performance files. +#' @param stratified_cross_performance_path Character. Path to cross-stratified performance files. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' \dontrun{ +#' plotStratifiedPerf("year", +#' stratified_performance_path = "data/Campylobacter/ML_year_performance", +#' stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance" +#' ) +#' } +plotStratifiedPerf <- function(year_or_country = "year", + stratified_performance_path = ".", + stratified_cross_performance_path = ".") { + perf <- arrow::read_parquet(file.path( + stratified_performance_path, + paste0(year_or_country, "_perf.parquet") + )) + + cross_test <- arrow::read_parquet(file.path( + stratified_cross_performance_path, + paste0( + "cross_", + year_or_country, + "_perf.parquet" + ) + )) + # if (year_or_country == "year") { + all <- perf |> + # dplyr::rename("train_year" = "strat_value") |> + dplyr::mutate(strat_value_test = strat_value) |> + dplyr::select( + drug_label, drug_or_class, + strat_value, strat_value_test, feature_type, feature_subtype, mcc + ) |> + dplyr::bind_rows(cross_test |> + dplyr::select( + drug_label, drug_or_class, + strat_value, strat_value_test, feature_type, + feature_subtype, mcc + )) |> + dplyr::mutate(category = dplyr::if_else( + strat_value == strat_value_test, "same", "different" + )) + # } else { + # all <- perf |> + # dplyr::rename("train_country" = "strat_value") |> + # dplyr::mutate(test_country = train_country) |> + # dplyr::select( + # drug_label, drug_or_class, + # train_country, test_country, + # feature_type, feature_subtype, mcc + # ) |> + # dplyr::bind_rows(cross_test |> + # dplyr::select( + # drug_label, drug_or_class, + # train_country, test_country, + # feature_type, feature_subtype, mcc + # )) |> + # dplyr::mutate(category = dplyr::if_else( + # train_country == test_country, "same country", "different country" + # )) + # } + + # fill_vals <- if (year_or_country == "year") { + # c( + # "same year bin" = "#b3cde3", + # "different year bin" = "#fbb4ae" + # ) + # } else { + # c( + # "same country" = "#b3cde3", + # "different country" = "#fbb4ae" + # ) + # } + + fill_vals <- c( + "same" = "#b3cde3", + "different" = "#fbb4ae" + ) + + plot <- ggplot2::ggplot( + all |> + dplyr::filter(drug_label == "drug", !is.na(mcc)), + ggplot2::aes(x = mcc, y = drug_or_class, fill = category) + ) + + ggridges::geom_density_ridges( + alpha = 0.5, + scale = 1, + rel_min_height = 0.01, + position = "identity" + ) + + ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "black") + + ggplot2::scale_fill_manual(values = fill_vals) + + ggplot2::theme_minimal(base_size = 14) + + ggplot2::labs( + title = if (year_or_country == "year") { + "Performance of AMR drug models with temporal holdouts" + } else { + "Performance of AMR drug models with geographical holdouts" + }, + x = "MCC", + y = "Drug", + fill = "Tested on" + ) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 14), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 14, colour = "black"), + axis.text.y = ggplot2::element_text(size = 14, colour = "black"), + axis.title.y = ggplot2::element_blank(), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "bottom", + plot.title = ggplot2::element_text(face = "bold"), + panel.grid.minor = ggplot2::element_blank(), + plot.margin = ggplot2::margin(0, 0, 0, 0) + ) + plot +} + +#' Plot multi-drug resistance (MDR) model performance +#' +#' Generates violin plots of performance, feature importance summaries, +#' and prediction confusion-style visualizations for MDR models. +#' +#' @param MDR_performance_path Character. Path to `MDR_perf.parquet`. +#' @param MDR_top_feature_path Character. Path to `MDR_top_features.parquet`. +#' @param MDR_pred_path Character. Path to `MDR_pred.parquet`. +#' +#' @return A list of ggplot objects. +#' @export +#' +#' @examples +#' \dontrun{ +#' plotMDR( +#' MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", +#' MDR_pred_path = "data/Campylobacter/MDR_ML_pred" +#' ) +#' } +plotMDR <- function(MDR_performance_path = ".", MDR_top_feature_path = ".", + MDR_pred_path = ".") { + MDR_perf <- arrow::read_parquet(file.path(MDR_performance_path, "MDR_perf.parquet")) + + # ---- Violin plot ---- + perf_plot <- ggplot2::ggplot( + MDR_perf, + ggplot2::aes(x = feature_type, y = mcc) + ) + + + # violins (overall distribution per feature type) + ggplot2::geom_violin(fill = "grey85", color = NA, alpha = 0.8) + + + # points (colored by binary vs counts) + ggplot2::geom_jitter( + ggplot2::aes(color = feature_subtype), + width = 0.12, size = 2, alpha = 0.8 + ) + + ggplot2::scale_color_manual(values = c( + "binary" = "#7B9CB5", + "counts" = "#CC8644" + )) + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::labs( + # title = "MDR model performances", + # subtitle = "Violin = distribution per feature type; points = binary vs counts", + x = "Feature type", + y = "MCC", + color = "Feature\nsubtype" + ) + + ggplot2::theme( + legend.position = "right", + plot.title = ggplot2::element_text(face = "bold") + ) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 14, colour = "black"), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "none", + title = ggplot2::element_text(face = "bold"), + panel.background = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major.x = ggplot2::element_blank(), # remove vertical lines + panel.grid.major.y = ggplot2::element_line(color = "grey80"), # keep horizontal lines + + axis.line = ggplot2::element_line(color = "black") + ) + + ggplot2::scale_y_continuous(limits = c(0, 1)) + + perf_plot + + MDR_pred <- arrow::read_parquet(file.path(MDR_pred_path, "MDR_pred.parquet")) |> + dplyr::mutate( + diff_top2 = purrr::pmap_dbl(dplyr::across(dplyr::contains(".pred") & dplyr::where(is.numeric)), function(...) { + x <- c(...) + sx <- sort(x, decreasing = TRUE) + sx[1] - sx[2] + }) # Difference between prediction probabilities of top two classes + ) |> + dplyr::select( + genome_id, resistant_classes, .pred_class, diff_top2, + feature_type, feature_subtype, seed + ) |> + dplyr::group_by(resistant_classes, .pred_class, feature_type) |> + dplyr::summarise(mean_margin = mean(diff_top2), n = n(), .groups = "drop") |> + dplyr::group_by(resistant_classes, feature_type) |> # normalize within true class + dplyr::mutate(sum = sum(n), prop = n / sum(n)) |> + dplyr::ungroup() + + MDR_pred_plot <- ggplot2::ggplot( + MDR_pred, + ggplot2::aes( + x = resistant_classes, + y = .pred_class + ) + ) + + ggplot2::geom_tile(ggplot2::aes(fill = prop)) + + ggplot2::geom_point(ggplot2::aes(size = mean_margin), color = "black") + + ggplot2::facet_wrap(~feature_type) + + ggplot2::scale_fill_distiller( + palette = "RdBu", + direction = 1, # flip with -1 if needed + name = "Prediction proportion" + ) + + ggplot2::labs(x = "true class", y = "predicted class") + + ggplot2::scale_size(range = c(1, 6), name = "Mean margin") + + ggplot2::coord_equal() + + ggplot2::theme_minimal() + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + + ggplot2::theme( + axis.title = ggplot2::element_text(colour = "black", size = 10), + axis.text.x = ggplot2::element_text(angle = 45, hjust = 1, size = 10, colour = "black"), + axis.text.y = ggplot2::element_text(size = 10, colour = "black"), + legend.title = ggplot2::element_text(size = 12), + legend.text = ggplot2::element_text(size = 10), + legend.position = "right", + title = ggplot2::element_text(face = "bold") + ) + + MDR_pred_plot + + # MDR_feat <- arrow::read_parquet(file.path( + # MDR_top_feature_path,"MDR_top_features.parquet")) |> + # pivot_longer(-c(Variable, feature_type, feature_subtype, seed), + # values_to = "Importance", + # names_to = "Resistant_classes") |> + # filter(!Importance == 0) + # + # MDR_feat_clean <- MDR_feat |> + # dplyr::filter(feature_type != "struct") |> + # dplyr::group_by(Resistant_classes, feature_type, feature_subtype, seed) |> + # dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> + # dplyr::ungroup() |> + # dplyr::mutate(Variable = gsub( ".NCBIFAM", "", Variable)) |> + # dplyr::mutate(Variable = gsub("^X", "", Variable)) |> + # dplyr::mutate(Variable = dplyr::if_else( + # feature_type == "domains", gsub("_.*", "", Variable), Variable)) |> + # dplyr::mutate(Variable = dplyr::if_else( + # feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable)) |> + # dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> + # dplyr::mutate( + # cluster = dplyr::coalesce(cluster, Variable) + # ) + # + # cluster_df <- MDR_feat_clean |> + # dplyr::group_by(Resistant_classes, cluster) |> + # dplyr::summarise( + # Importance = median(Importance, na.rm = TRUE), + # .groups = "drop" + # ) + # + # top_clusters <- cluster_df |> + # group_by(Resistant_classes) |> + # group_modify(~{ + # + # df <- .x + # + # top_pos <- df |> + # arrange(desc(Importance)) |> + # slice_head(n = 10) + # + # top_neg <- df |> + # arrange(Importance) |> + # slice_head(n = 10) + # + # bind_rows(top_pos, top_neg) + # }) |> + # ungroup() + # + # top_clusters <- top_clusters |> + # dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> + # dplyr::mutate( + # proteinName = dplyr::coalesce(proteinName, cluster), # fallback + # proteinName = stringr::str_trunc(proteinName, 50) + # ) |> + # dplyr::distinct(Resistant_classes, proteinName, Importance) |> + # # ✅ reorder AFTER naming + # dplyr::group_by(Resistant_classes) |> + # dplyr::mutate( + # proteinName = forcats::fct_reorder(proteinName, Importance) + # ) |> + # dplyr::ungroup() + # + # ggplot(top_clusters, + # aes(x = Importance, y = proteinName)) + + # + # # line (lollipop stem) + # geom_segment( + # aes(x = 0, xend = Importance, + # y = proteinName, yend = proteinName), + # color = "grey60" + # ) + + # + # # dot + # geom_point( + # aes(color = Importance > 0), + # size = 3 + # ) + + # + # facet_wrap(~ Resistant_classes, scales = "free_y") + + # + # scale_color_manual( + # values = c("TRUE" = "#5b8db8", # positive + # "FALSE" = "#d4872a"), # negative + # guide = "none" + # ) + + # + # theme_minimal(base_size = 13) + + # labs( + # x = "Median importance", + # y = "Cluster" + # ) + + # theme( + # panel.grid.minor = element_blank(), + # strip.text = element_text(face = "bold") + # ) +} + +#' Compare shuffled vs real model performance +#' +#' Creates boxplots comparing performance (MCC) between real and shuffled labels +#' across feature types. +#' +#' @param metadata_path Character. Unused; retained for backward compatibility. +#' @param performance_path A performance tibble (with `feature_type`, +#' `feature_subtype`, `mcc`, and `shuffled` columns), or a directory path +#' containing `all_perf.parquet`. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' performance <- tibble::tibble( +#' feature_type = rep(c("genes", "proteins"), each = 4), +#' feature_subtype = rep(c("binary", "counts"), times = 4), +#' mcc = c(0.7, 0.6, 0.65, 0.55, 0.1, 0.05, 0.08, 0.02), +#' shuffled = rep(c(FALSE, TRUE), each = 4) +#' ) +#' plotShuffleVsReal(performance_path = performance) +plotShuffleVsReal <- function(metadata_path = ".", performance_path = ".") { + if (!is.data.frame(performance_path)) { + performance_path <- arrow::read_parquet( + file.path(performance_path, "all_perf.parquet") + ) + } + performance <- performance_path + + performance |> + dplyr::mutate( + shuffled_label = dplyr::if_else(shuffled, "shuffled", "real") + ) |> + ggplot2::ggplot(ggplot2::aes(x = feature_subtype, y = mcc, fill = shuffled_label)) + + ggplot2::geom_boxplot( + width = 0.55, outlier.size = 0.8, outlier.alpha = 0.4, + outlier.color = "grey50", linewidth = 0.4 + ) + + ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "grey60", linewidth = 0.4) + + ggplot2::scale_fill_manual( + values = c("real" = "#7B9CB5", "shuffled" = "#C4B8A8"), + name = NULL + ) + + ggplot2::scale_y_continuous(limits = c(-0.2, 1), breaks = seq(-0.2, 1, 0.2)) + + ggplot2::facet_wrap(~feature_type, nrow = 1) + + ggplot2::theme_minimal(base_size = 12) + + theme( + panel.grid.major.x = element_blank(), + panel.grid.minor = element_blank(), + panel.grid.major.y = element_line(color = "#E5E2D9", linewidth = 0.4), + strip.text = element_text(color = "grey30", face = "bold", size = 10), + # strip.background = element_rect(fill = "#EEEAE0", color = NA), + axis.text = element_text(color = "grey45"), + legend.position = "top", + legend.text = element_text(color = "grey40", size = 10) + ) + + labs( + x = NULL, y = "MCC" + ) +} + +#' Plot top contributing feature clusters +#' +#' Identifies top contributing clusters across feature types and drugs, +#' and visualizes their relative contributions. +#' +#' @param top_feat_path Character. Path to `all_top_features.parquet`. +#' @param cluster_feature_path Character. Path to `cluster_feature.parquet`. +#' @param protein_names_path Character. Path to `protein_names.parquet`. +#' @param top_n Integer. Number of top features to retain per model. +#' +#' @return A ggplot object. +#' @export +#' +#' @examples +#' \dontrun{ +#' plotTopClusters( +#' top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", +#' protein_names_path = "data/Campylobacter/", top_n = 10 +#' ) +#' } +plotTopClusters <- function(top_feat_path = ".", cluster_feature_path = ".", + protein_names_path = ".", top_n = 10) { + top_feat <- arrow::read_parquet(file.path(top_feat_path, "all_top_features.parquet")) + cluster_feature <- arrow::read_parquet(file.path(cluster_feature_path, "cluster_feature.parquet")) + protein_names <- arrow::read_parquet(file.path(protein_names_path, "protein_names.parquet")) + + # which clusters appear in top n across feature types per drug + # join top features with cluster mapping, filter out struct and shuffled + top_feat_clean <- top_feat |> + dplyr::filter(!shuffled, feature_type != "struct", drug_label == "drug") |> + dplyr::group_by(drug_or_class, feature_type, feature_subtype, seed) |> + dplyr::slice_max(Importance, n = top_n, with_ties = FALSE) |> + dplyr::ungroup() |> + dplyr::mutate(Variable = gsub(".NCBIFAM", "", Variable)) |> + dplyr::mutate(Variable = gsub("^X", "", Variable)) |> + dplyr::mutate(Variable = dplyr::if_else( + feature_type == "domains", gsub("_.*", "", Variable), Variable + )) |> + dplyr::mutate(Variable = dplyr::if_else( + feature_type == "proteins", gsub("fig.", "fig|", Variable), Variable + )) |> + dplyr::left_join(cluster_feature, by = c("Variable" = "feature")) |> + dplyr::mutate( + cluster = dplyr::coalesce(cluster, Variable), # fallback to Variable if no match + Importance_signed = dplyr::if_else(Sign == "NEG", -Importance, Importance) + ) + + shared_mat <- top_feat_clean |> + dplyr::group_by(drug_or_class, feature_type, cluster) |> + dplyr::summarise(abs_imp = median(Importance, na.rm = TRUE), .groups = "drop") |> + # convert to contribution within each feature_type + dplyr::group_by(drug_or_class, feature_type) |> + dplyr::mutate(contribution = abs_imp / sum(abs_imp, na.rm = TRUE)) |> + # pick top n contributors + dplyr::slice_max(contribution, n = top_n, with_ties = FALSE) |> + dplyr::ungroup() |> + dplyr::add_count(drug_or_class, cluster, name = "n_feat_types") |> + dplyr::left_join(protein_names, by = c("cluster" = "proteinID")) |> + dplyr::mutate( + proteinName = stringr::str_trunc(proteinName, 50), + proteinName = forcats::fct_reorder(proteinName, n_feat_types) + ) + + feat_plot <- ggplot2::ggplot( + shared_mat, + ggplot2::aes( + x = feature_type, + y = proteinName, + fill = contribution + ) + ) + + ggplot2::geom_tile(color = "#FAFAF7", linewidth = 0.5, width = 0.9, height = 0.9) + + # coord_fixed() + + ggplot2::scale_fill_distiller( + palette = "RdPu", + direction = 1, + name = "contribution", + na.value = "#EEEAE0" + ) + + ggplot2::facet_wrap(~drug_or_class, scales = "free_y") + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme( + panel.grid.major.x = ggplot2::element_blank(), + panel.grid.minor = ggplot2::element_blank(), + panel.grid.major.y = ggplot2::element_line(color = "#E5E2D9", linewidth = 0.4), + strip.text = ggplot2::element_text(color = "grey30", face = "bold", size = 10), + strip.background = ggplot2::element_rect(fill = "#EEEAE0", color = NA), + axis.title.y = ggplot2::element_blank(), + axis.text.x = ggplot2::element_text(color = "black", angle = 30, hjust = 1, size = 6), + axis.text.y = ggplot2::element_text(color = "black", size = 6), + legend.position = "bottom", + legend.text = ggplot2::element_text(color = "grey40", size = 10), + legend.title = ggplot2::element_text(color = "grey40", size = 10) + ) + + feat_plot +} diff --git a/R/prep_ml.R b/R/prep_ml.R index 9ce2ef2..efd4a97 100644 --- a/R/prep_ml.R +++ b/R/prep_ml.R @@ -96,18 +96,18 @@ loadMLInputTibble <- function(parquet_path) { "genome_drug.resistant_phenotype" %in% colnames(long_tibble), "resistant_classes" %in% colnames(long_tibble) )) { - stop(paste( - "The tibble loaded from `parquet_path` must have a target", - "variable column named either `genome_drug.resistant_phenotype` or", + stop( + "The tibble loaded from `parquet_path` must have a target ", + "variable column named either `genome_drug.resistant_phenotype` or ", "`resistant_classes`, but not both." - )) + ) } if (!all(c("genome_id", "feature_id", "value") %in% colnames(long_tibble))) { - stop(paste( - "The data is missing one or more of the following required", + stop( + "The data is missing one or more of the following required ", "columns: 'genome_id', 'feature_id', 'value'." - )) + ) } target_var <- .getTargetVarName(long_tibble) diff --git a/R/run_ML.R b/R/run_ML.R index e85d42e..bd0fe0a 100644 --- a/R/run_ML.R +++ b/R/run_ML.R @@ -326,7 +326,7 @@ createMLinputList <- function(path, # Case A: stratified -> prefix before the stratify label if (!is.na(i)) { if (i - 1 >= 1) { - return(paste(x[1:(i - 1)], collapse = "_")) + return(paste(x[seq_len(i - 1)], collapse = "_")) } return("") } @@ -334,11 +334,11 @@ createMLinputList <- function(path, # Case B: unstratified -> prefix is first two tokens if (x[2] == "drug" && x[3] != "class") { # Case A: Cje_drug_X - return(paste(x[1:2], collapse = "_")) + return(paste(x[seq_len(2)], collapse = "_")) } if (x[2] == "drug" && x[3] == "class") { # Case A: Cje_drug_X - return(paste(x[1:3], collapse = "_")) + return(paste(x[seq_len(3)], collapse = "_")) } }) ) @@ -659,7 +659,7 @@ runMDRmodels <- function(path, if (inherits(res, "try-error")) { warning( "Model failed for: ", output_prefix, - "\n Error: ", attr(res, "condition")$message + "\n ", attr(res, "condition")$message ) return(NULL) } @@ -964,7 +964,7 @@ runMLmodels <- function(path, if (inherits(res, "try-error")) { warning( "Model failed for: ", output_prefix, - "\n Error: ", attr(res, "condition")$message + "\n ", attr(res, "condition")$message ) return(NULL) } diff --git a/R/run_ml_pipeline.R b/R/run_ml_pipeline.R index 5c47271..db0fc64 100644 --- a/R/run_ml_pipeline.R +++ b/R/run_ml_pipeline.R @@ -136,18 +136,16 @@ runMLPipeline <- function( } if (model != "LR" & multi_class) { - stop(paste( - "Multi-class classification can be performed only with", + stop( + "Multi-class classification can be performed only with ", "logistic regression in our current implementation." - )) + ) } if (verbose) { message( - paste( - "`ml_input_tibble` has", num_obs_ml_input_tibble, "observations of", - getNumFeat(ml_input_tibble), "features." - ) + "`ml_input_tibble` has ", num_obs_ml_input_tibble, " observations of ", + getNumFeat(ml_input_tibble), " features." ) } @@ -158,7 +156,11 @@ runMLPipeline <- function( dplyr::mutate(prop = n / sum(n)) if (verbose) { - print(rs_props_original_ml_input_tibble) + rs_props_summary <- paste( + utils::capture.output(rs_props_original_ml_input_tibble), + collapse = "\n" + ) + message(rs_props_summary) } res_prop_original_ml_input_tibble <- rs_props_original_ml_input_tibble |> @@ -183,19 +185,19 @@ runMLPipeline <- function( external_test_data <- TRUE if (multi_class) { - stop(paste( - "Multi-class classification cannot be performed with", + stop( + "Multi-class classification cannot be performed with ", "external test data with the current implementation." - )) + ) } num_obs_test_data <- nrow(test_data) if (verbose) { - message(paste( - "`test_data` has", num_obs_test_data, - "observations of", getNumFeat(test_data), "features." - )) + message( + "`test_data` has ", num_obs_test_data, + " observations of ", getNumFeat(test_data), " features." + ) } rs_props_test_data <- test_data |> @@ -203,7 +205,11 @@ runMLPipeline <- function( dplyr::mutate(prop = n / sum(n)) if (verbose) { - print(rs_props_test_data) + test_data_props_summary <- paste( + utils::capture.output(rs_props_test_data), + collapse = "\n" + ) + message(test_data_props_summary) } # Get the proportion of resistant genomes in the external `test_data`. @@ -248,12 +254,12 @@ runMLPipeline <- function( # Stratify by R/S phenotype res_train_indices <- sample(res_indices, n_train_res) - sus_train_indices <- (1:n_train)[!(1:n_train %in% res_indices)] + sus_train_indices <- seq_len(n_train)[!(seq_len(n_train) %in% res_indices)] train_indices <- c(res_train_indices, sus_train_indices) - val_indices <- (1:nrow(original_ml_input_tibble))[!(1:nrow( + val_indices <- seq_len(nrow(original_ml_input_tibble))[!(seq_len(nrow( original_ml_input_tibble - ) %in% train_indices)] + )) %in% train_indices)] data_split[[2]] <- as.integer(train_indices) data_split[[3]] <- as.integer(val_indices) @@ -261,7 +267,7 @@ runMLPipeline <- function( n_train <- nrow(original_ml_input_tibble) n_val <- 0 - data_split[[2]] <- 1:nrow(original_ml_input_tibble) + data_split[[2]] <- seq_len(nrow(original_ml_input_tibble)) data_split[[3]] <- NA data_split[[4]] <- tibble::tibble(id = "Resample1") } @@ -314,7 +320,7 @@ runMLPipeline <- function( nmcc <- .calculatenMCC(test_data_plus_predictions) if (verbose) { - message(paste("Matthews correlation coefficient:", mcc, "| nMCC:", nmcc)) + message("Matthews correlation coefficient: ", mcc, " | nMCC: ", nmcc) } top_feat_tibble <- extractTopFeats(fit, @@ -342,7 +348,7 @@ runMLPipeline <- function( # `top_feats`. n_random_feats <- n_top_feats - length(top_feats) - random_indices <- sample(1:length(zero_vi_feats), n_random_feats) + random_indices <- sample(seq_along(zero_vi_feats), n_random_feats) random_feats <- zero_vi_feats[random_indices] top_feats <- c(top_feats, random_feats) diff --git a/man/calculateEvalMets.Rd b/man/calculateEvalMets.Rd index 0506645..5b61252 100644 --- a/man/calculateEvalMets.Rd +++ b/man/calculateEvalMets.Rd @@ -15,8 +15,8 @@ F1 score, AUPRC, balanced accuracy, nMCC, and log2(AUPRC/prior) } \description{ Returns the F1 score, area under the precision-recall curve (AUPRC), balanced -accuracy, normalized (to a 0 to 1 scale instead of -1 to 1) Matthews -correlation coefficient (nMCC), and log2(AUPRC/prior) based on the AMR +accuracy, Matthews correlation coefficient (MCC), normalized MCC (nMCC), +and log2(AUPRC/prior) based on the AMR phenotype predictions by an ML model compared against the actual values. } \examples{ diff --git a/man/demo_fit.Rd b/man/demo_fit.Rd index 0810296..9f29307 100644 --- a/man/demo_fit.Rd +++ b/man/demo_fit.Rd @@ -16,4 +16,8 @@ data(demo_fit) \description{ A tuned logistic-regression workflow fitted on \link{demo_ml_tibble}. } +\examples{ +data(demo_fit) +class(demo_fit) +} \keyword{datasets} diff --git a/man/demo_ml_tibble.Rd b/man/demo_ml_tibble.Rd index 021834d..79f685a 100644 --- a/man/demo_ml_tibble.Rd +++ b/man/demo_ml_tibble.Rd @@ -19,4 +19,8 @@ Stratified subset (30 Resistant + 30 Susceptible) of the AMP-genes-binary matrix from the bundled \code{Sfl_parquet.duckdb}, restricted to 80 feature columns. } +\examples{ +data(demo_ml_tibble) +dim(demo_ml_tibble) +} \keyword{datasets} diff --git a/man/dot-calculateMCC.Rd b/man/dot-calculateMCC.Rd new file mode 100644 index 0000000..edcc550 --- /dev/null +++ b/man/dot-calculateMCC.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/core_ml.R +\name{.calculateMCC} +\alias{.calculateMCC} +\title{.calculateMCC()} +\usage{ +.calculateMCC(test_data_plus_predictions) +} +\arguments{ +\item{test_data_plus_predictions}{Test data (tibble) with an added column for +predicted phenotype labels, such as the output of \code{predictML()}} +} +\value{ +Matthews correlation coefficient (MCC), range -1 to 1 +} +\description{ +Returns the Matthews correlation coefficient (MCC) +based on the AMR phenotype predictions by an +ML model compared against the actual values. +} diff --git a/man/dot-calculatenMCC.Rd b/man/dot-calculatenMCC.Rd index de13978..8763081 100644 --- a/man/dot-calculatenMCC.Rd +++ b/man/dot-calculatenMCC.Rd @@ -11,11 +11,10 @@ predicted phenotype labels, such as the output of \code{predictML()}} } \value{ -Normalized (to a 0 to 1 scale instead of -1 to 1) Matthews -correlation coefficient (nMCC) +Normalized Matthews correlation coefficient (nMCC), range 0 to 1 } \description{ -Returns the normalized (to a 0 to 1 scale instead of -1 to 1) Matthews -correlation coefficient (nMCC) based on the AMR phenotype predictions by an -ML model compared against the actual values. +Returns the normalized (0 to 1) Matthews correlation coefficient (nMCC) +based on the AMR phenotype predictions by an ML model compared against +the actual values. } diff --git a/man/dot-parquet2LOODrugMatrix.Rd b/man/dot-parquet2LOODrugMatrix.Rd new file mode 100644 index 0000000..9ad18ce --- /dev/null +++ b/man/dot-parquet2LOODrugMatrix.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generate_matrices_ml.R +\name{.parquet2LOODrugMatrix} +\alias{.parquet2LOODrugMatrix} +\title{Build leave-one-out (LOO) merged parquet matrices from drug parquet files.} +\usage{ +.parquet2LOODrugMatrix(path, verbosity = c("minimal", "debug")) +} +\arguments{ +\item{path}{Character. Base directory containing stratified parquet matrices. +Expected subdirs: matrix/..} + +\item{verbosity}{Character. "minimal" or "debug"; when "debug", prints detailed steps.} +} +\value{ +Invisibly returns a tibble with paths of created LOO parquet files. +} +\description{ +Build leave-one-out (LOO) merged parquet matrices from drug parquet files. +} diff --git a/man/plotBaselineComparison.Rd b/man/plotBaselineComparison.Rd index b212ab2..8d6403d 100644 --- a/man/plotBaselineComparison.Rd +++ b/man/plotBaselineComparison.Rd @@ -2,32 +2,34 @@ % Please edit documentation in R/plot_ml.R \name{plotBaselineComparison} \alias{plotBaselineComparison} -\title{plotBaselineComparison()} +\title{Compare Baseline Performance With and Without Shuffled Labels} \usage{ plotBaselineComparison(non_shuffled_label_results, shuffled_label_results) } \arguments{ -\item{non_shuffled_label_results}{Output of \code{runMLPipeline()} -(\code{shuffle_labels = FALSE})} +\item{non_shuffled_label_results}{Output of \code{runMLPipeline(shuffle_labels = FALSE)}} -\item{shuffled_label_results}{Output of \code{runMLPipeline()} -(\code{shuffle_labels = TRUE})} +\item{shuffled_label_results}{Output of \code{runMLPipeline(shuffle_labels = TRUE)}} } \value{ -A bar plot with balanced accuracy comparisons per antibiotic +A barplot comparing balanced accuracy across models. } \description{ -Generates a bar plot that compares model performance with and without -randomly shuffled AMR phenotype labels. +Produces a bar plot comparing balanced accuracy for each antibiotic using +true AMR labels vs. randomly shuffled labels. } \examples{ -non_shuffled <- tibble::tibble( - antibiotic = c("AMP", "CIP", "CRO"), - bal_acc = c(0.88, 0.81, 0.92) +non_shuffled <- list( + performance_tibble = tibble::tibble( + antibiotic = c("AMP", "CIP", "CRO"), + bal_acc = c(0.88, 0.81, 0.92) + ) ) -shuffled <- tibble::tibble( - antibiotic = c("AMP", "CIP", "CRO"), - bal_acc = c(0.52, 0.49, 0.55) +shuffled <- list( + performance_tibble = tibble::tibble( + antibiotic = c("AMP", "CIP", "CRO"), + bal_acc = c(0.52, 0.49, 0.55) + ) ) plotBaselineComparison(non_shuffled, shuffled) } diff --git a/man/plotCM.Rd b/man/plotCM.Rd new file mode 100644 index 0000000..673bca0 --- /dev/null +++ b/man/plotCM.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotCM} +\alias{plotCM} +\title{Plot a Confusion Matrix Heatmap} +\usage{ +plotCM(test_data_plus_predictions) +} +\arguments{ +\item{test_data_plus_predictions}{A tibble containing true and predicted +phenotype labels (e.g. the output of \code{predictML()}), or a path to a TSV file +containing the same.} +} +\value{ +A heatmap (\code{ggplot2} object) showing the confusion matrix. +} +\description{ +Produces a heatmap visualization of the confusion matrix for AMR predictions. +} +\examples{ +data(demo_fit) +data(demo_ml_tibble) +preds <- predictML(demo_fit, demo_ml_tibble) +plotCM(preds) + +} diff --git a/man/plotCrossDrug.Rd b/man/plotCrossDrug.Rd new file mode 100644 index 0000000..4c9163d --- /dev/null +++ b/man/plotCrossDrug.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotCrossDrug} +\alias{plotCrossDrug} +\title{Plot cross-drug prediction as a heatmap} +\usage{ +plotCrossDrug( + cross_test_performance_path = ".", + drug_performance_path = ".", + metadata_path = "." +) +} +\arguments{ +\item{cross_test_performance_path}{A cross-drug performance tibble (with +\code{drug}, \code{test_drug}, and \code{mcc} columns), or a directory path +containing \code{cross_drug_perf.parquet}.} + +\item{drug_performance_path}{A performance tibble (with \code{drug_label}, +\code{drug_or_class}, and \code{mcc} columns), or a directory path containing +\code{all_perf.parquet}.} + +\item{metadata_path}{Character. Path to a directory containing \code{metadata.parquet}.} +} +\value{ +A ComplexHeatmap object. +} +\description{ +Creates a heatmap showing cross-drug model performance (MCC), where models +trained on one drug are evaluated on another. +} +\examples{ +cross_drug <- tibble::tibble( + drug = c("AMP", "AMP", "CIP", "CIP", "NAL", "NAL"), + test_drug = c("CIP", "NAL", "AMP", "NAL", "AMP", "CIP"), + mcc = c(0.3, 0.2, 0.4, 0.25, 0.15, 0.35) +) +performance <- tibble::tibble( + drug_label = "drug", + drug_or_class = c("AMP", "CIP", "NAL"), + mcc = c(0.8, 0.7, 0.6) +) +plotCrossDrug( + cross_test_performance_path = cross_drug, + drug_performance_path = performance, + metadata_path = system.file("extdata", package = "amRml") +) +} diff --git a/man/plotDefaultEval.Rd b/man/plotDefaultEval.Rd deleted file mode 100644 index ca2dcfe..0000000 --- a/man/plotDefaultEval.Rd +++ /dev/null @@ -1,48 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot_ml.R -\name{plotDefaultEval} -\alias{plotDefaultEval} -\title{plotDefaultEval()} -\usage{ -plotDefaultEval( - default_eval_tibble, - x_default_eval = "train_prop", - y_default_eval = "avg_f1_score", - xlab = "Train data proportion", - ylab = "Average F1 score" -) -} -\arguments{ -\item{default_eval_tibble}{Output of \code{findOptimalMLDefaults()}} - -\item{x_default_eval}{\link[rlang:chr]{chr} x value of default evaluation plot: "train_prop" -or "n_fold"} - -\item{y_default_eval}{\link[rlang:chr]{chr} y value of default evaluation plot. It can be -"avg_runtime_sec" or one of the following performance metrics: -"avg_f1_score", "avg_log2_apop", "avg_bal_acc", or "avg_nmcc"} - -\item{xlab}{\link[rlang:chr]{chr} Label for x axis} - -\item{ylab}{\link[rlang:chr]{chr} Label for y axis} -} -\value{ -A \code{ggplot2} scatterplot (performance metric or runtime vs. -\code{train_prop} or \code{n_fold}), colored by model -} -\description{ -Plots performance metric or runtime vs. training data proportion or number -of cross-validation folds, colored by model. -} -\examples{ -default_eval <- tibble::tibble( - train_prop = c(0.5, 0.6, 0.7, 0.5, 0.6, 0.7), - n_fold = rep(5, 6), - model = rep(c("LR", "RF"), each = 3), - avg_f1_score = c(0.72, 0.78, 0.83, 0.70, 0.75, 0.80) -) -plotDefaultEval(default_eval, - x_default_eval = "train_prop", - y_default_eval = "avg_f1_score" -) -} diff --git a/man/plotDrugDist.Rd b/man/plotDrugDist.Rd new file mode 100644 index 0000000..9188dfe --- /dev/null +++ b/man/plotDrugDist.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotDrugDist} +\alias{plotDrugDist} +\title{Plot drug phenotype distribution} +\usage{ +plotDrugDist(metadata_path = ".") +} +\arguments{ +\item{metadata_path}{Character. Path to directory containing \code{metadata.parquet}.} +} +\value{ +A ggplot object. +} +\description{ +Reads metadata and generates a stacked bar plot showing counts of resistant +and susceptible phenotypes per antibiotic. +} +\examples{ +plotDrugDist(metadata_path = system.file("extdata", package = "amRml")) +} diff --git a/man/plotDrugPerf.Rd b/man/plotDrugPerf.Rd new file mode 100644 index 0000000..8fd65df --- /dev/null +++ b/man/plotDrugPerf.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotDrugPerf} +\alias{plotDrugPerf} +\title{Plot drug-level model performance} +\usage{ +plotDrugPerf(metadata_path = ".", performance_path = ".") +} +\arguments{ +\item{metadata_path}{Character. Path to a directory containing \code{metadata.parquet}.} + +\item{performance_path}{A performance tibble (with \code{drug_label}, \code{shuffled}, +\code{drug_or_class}, \code{feature_type}, \code{feature_subtype}, and \code{mcc} columns), or a +directory path containing \code{all_perf.parquet}.} +} +\value{ +A patchwork ggplot object combining multiple panels. +} +\description{ +Generates heatmaps and ridge plots summarizing model performance (MCC) +across drugs and feature types. +} +\examples{ +# Several models per drug x feature type so the ridge densities can be drawn. +performance <- tidyr::expand_grid( + drug_or_class = c("AMP", "CIP", "NAL"), + feature_type = c("genes", "proteins"), + feature_subtype = c("binary", "counts"), + replicate = 1:8 +) +performance$drug_label <- "drug" +performance$shuffled <- FALSE +performance$mcc <- 0.6 + 0.3 * sin(seq_len(nrow(performance))) +plotDrugPerf( + metadata_path = system.file("extdata", package = "amRml"), + performance_path = performance +) +} diff --git a/man/plotMDR.Rd b/man/plotMDR.Rd new file mode 100644 index 0000000..6573e4d --- /dev/null +++ b/man/plotMDR.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotMDR} +\alias{plotMDR} +\title{Plot multi-drug resistance (MDR) model performance} +\usage{ +plotMDR( + MDR_performance_path = ".", + MDR_top_feature_path = ".", + MDR_pred_path = "." +) +} +\arguments{ +\item{MDR_performance_path}{Character. Path to \code{MDR_perf.parquet}.} + +\item{MDR_top_feature_path}{Character. Path to \code{MDR_top_features.parquet}.} + +\item{MDR_pred_path}{Character. Path to \code{MDR_pred.parquet}.} +} +\value{ +A list of ggplot objects. +} +\description{ +Generates violin plots of performance, feature importance summaries, +and prediction confusion-style visualizations for MDR models. +} +\examples{ +\dontrun{ +plotMDR( + MDR_performance_path = "data/Campylobacter/MDR_ML_performance", MDR_top_feature_path = "data/Campylobacter/MDR_ML_top_features", + MDR_pred_path = "data/Campylobacter/MDR_ML_pred" +) +} +} diff --git a/man/plotPRC.Rd b/man/plotPRC.Rd index 3bad222..7581ca3 100644 --- a/man/plotPRC.Rd +++ b/man/plotPRC.Rd @@ -2,37 +2,29 @@ % Please edit documentation in R/plot_ml.R \name{plotPRC} \alias{plotPRC} -\title{plotPRC()} +\title{Plot a Precision-Recall Curve} \usage{ plotPRC(test_data_plus_predictions) } \arguments{ -\item{test_data_plus_predictions}{Test data (tibble) with an added column for -predicted phenotype labels, such as the output of \code{predict()}.} +\item{test_data_plus_predictions}{A tibble of test data with added prediction +columns (e.g. the output of \code{predictML()} or \code{runMLmodels(return_pred=TRUE)}), +or a path to a TSV file containing the same.} } \value{ -A precision-recall curve as a \code{ggplot2} object +A \code{ggplot2} object showing the precision-recall curve. } \description{ -Plots the precision-recall curve given a set of test data plus predicted AMR -phenotypes. +Generates a precision-recall curve (PRC) for AMR phenotype prediction results. +} +\details{ +The function uses \code{yardstick::pr_curve()} to compute the PR curve and then +visualizes it using \code{ggplot2}. } \examples{ -preds <- tibble::tibble( - genome_id = paste0("g", 1:10), - genome_drug.resistant_phenotype = factor( - rep(c("Resistant", "Susceptible"), each = 5), - levels = c("Resistant", "Susceptible") - ), - .pred_class = factor( - c( - "Resistant", "Resistant", "Susceptible", "Resistant", "Susceptible", - "Susceptible", "Resistant", "Susceptible", "Susceptible", "Resistant" - ), - levels = c("Resistant", "Susceptible") - ), - .pred_Resistant = c(0.9, 0.8, 0.4, 0.7, 0.3, 0.2, 0.6, 0.1, 0.2, 0.55), - .pred_Susceptible = c(0.1, 0.2, 0.6, 0.3, 0.7, 0.8, 0.4, 0.9, 0.8, 0.45) -) +data(demo_fit) +data(demo_ml_tibble) +preds <- predictML(demo_fit, demo_ml_tibble) plotPRC(preds) + } diff --git a/man/plotROC.Rd b/man/plotROC.Rd new file mode 100644 index 0000000..45e0f54 --- /dev/null +++ b/man/plotROC.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotROC} +\alias{plotROC} +\title{Plot a Receiver Operating Characteristic (ROC) Curve} +\usage{ +plotROC(test_data_plus_predictions) +} +\arguments{ +\item{test_data_plus_predictions}{A tibble of test data with prediction +columns (e.g. the output of \code{predictML()} or \code{runMLmodels(return_pred=TRUE)}), +or a path to a TSV file containing the same.} +} +\value{ +A ROC curve plotted using \code{ggplot2::autoplot()}. +} +\description{ +Generates a ROC curve for AMR phenotype prediction results. +} +\examples{ +data(demo_fit) +data(demo_ml_tibble) +preds <- predictML(demo_fit, demo_ml_tibble) +plotROC(preds) + +} diff --git a/man/plotShuffleVsReal.Rd b/man/plotShuffleVsReal.Rd new file mode 100644 index 0000000..3ed40f0 --- /dev/null +++ b/man/plotShuffleVsReal.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotShuffleVsReal} +\alias{plotShuffleVsReal} +\title{Compare shuffled vs real model performance} +\usage{ +plotShuffleVsReal(metadata_path = ".", performance_path = ".") +} +\arguments{ +\item{metadata_path}{Character. Unused; retained for backward compatibility.} + +\item{performance_path}{A performance tibble (with \code{feature_type}, +\code{feature_subtype}, \code{mcc}, and \code{shuffled} columns), or a directory path +containing \code{all_perf.parquet}.} +} +\value{ +A ggplot object. +} +\description{ +Creates boxplots comparing performance (MCC) between real and shuffled labels +across feature types. +} +\examples{ +performance <- tibble::tibble( + feature_type = rep(c("genes", "proteins"), each = 4), + feature_subtype = rep(c("binary", "counts"), times = 4), + mcc = c(0.7, 0.6, 0.65, 0.55, 0.1, 0.05, 0.08, 0.02), + shuffled = rep(c(FALSE, TRUE), each = 4) +) +plotShuffleVsReal(performance_path = performance) +} diff --git a/man/plotStratifiedPerf.Rd b/man/plotStratifiedPerf.Rd new file mode 100644 index 0000000..1fe5854 --- /dev/null +++ b/man/plotStratifiedPerf.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotStratifiedPerf} +\alias{plotStratifiedPerf} +\title{Plot stratified model performance} +\usage{ +plotStratifiedPerf( + year_or_country = "year", + stratified_performance_path = ".", + stratified_cross_performance_path = "." +) +} +\arguments{ +\item{year_or_country}{Character. Either "year" or "country".} + +\item{stratified_performance_path}{Character. Path to stratified performance files.} + +\item{stratified_cross_performance_path}{Character. Path to cross-stratified performance files.} +} +\value{ +A ggplot object. +} +\description{ +Visualizes model performance (MCC) stratified by year or country, +comparing within-group vs cross-group evaluation. +} +\examples{ +\dontrun{ +plotStratifiedPerf("year", + stratified_performance_path = "data/Campylobacter/ML_year_performance", + stratified_cross_performance_path = "data/Campylobacter/cross_test_ML_year_performance" +) +} +} diff --git a/man/plotTopClusters.Rd b/man/plotTopClusters.Rd new file mode 100644 index 0000000..0420d67 --- /dev/null +++ b/man/plotTopClusters.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_ml.R +\name{plotTopClusters} +\alias{plotTopClusters} +\title{Plot top contributing feature clusters} +\usage{ +plotTopClusters( + top_feat_path = ".", + cluster_feature_path = ".", + protein_names_path = ".", + top_n = 10 +) +} +\arguments{ +\item{top_feat_path}{Character. Path to \code{all_top_features.parquet}.} + +\item{cluster_feature_path}{Character. Path to \code{cluster_feature.parquet}.} + +\item{protein_names_path}{Character. Path to \code{protein_names.parquet}.} + +\item{top_n}{Integer. Number of top features to retain per model.} +} +\value{ +A ggplot object. +} +\description{ +Identifies top contributing clusters across feature types and drugs, +and visualizes their relative contributions. +} +\examples{ +\dontrun{ +plotTopClusters( + top_feat_path = "data/Campylobacter/ML_top_features", cluster_feature_path = "data/Campylobacter/", + protein_names_path = "data/Campylobacter/", top_n = 10 +) +} +} diff --git a/man/plotTopFeatsVI.Rd b/man/plotTopFeatsVI.Rd index 9bd920b..21fbb48 100644 --- a/man/plotTopFeatsVI.Rd +++ b/man/plotTopFeatsVI.Rd @@ -2,23 +2,27 @@ % Please edit documentation in R/plot_ml.R \name{plotTopFeatsVI} \alias{plotTopFeatsVI} -\title{plotTopFeatsVI()} +\title{Plot Top Feature Importances} \usage{ -plotTopFeatsVI(fit, n_top_feats = 10) +plotTopFeatsVI(topfeat, n_top_feats = 10) } \arguments{ -\item{fit}{Best model fit, such as the output of \code{fitBestModel()}} +\item{topfeat}{A tibble of feature importance scores (e.g. the output of +\code{extractTopFeats()} or \code{runMLmodels()}), with \code{Variable}, \code{Importance}, and +\code{Sign} columns, or a path to a TSV file containing the same.} -\item{n_top_feats}{\link[pillar:num]{num} Number of top features to plot} +\item{n_top_feats}{Number of top features to display (default: 10).} } \value{ -Variable importance plot (a \code{ggplot2} object) +A bar plot of variable importance (\code{ggplot2} object). } \description{ -Generates a plot showing the top features and their variable importance -scores. +Creates a bar plot showing the most important features affecting +AMR phenotype predictions. } \examples{ data(demo_fit) -plotTopFeatsVI(demo_fit, n_top_feats = 10) +top_feats <- extractTopFeats(demo_fit, n_top_feats = 10) +plotTopFeatsVI(top_feats, n_top_feats = 10) + } diff --git a/man/runIFE.Rd b/man/runIFE.Rd index 2e39c9e..31afb6b 100644 --- a/man/runIFE.Rd +++ b/man/runIFE.Rd @@ -4,13 +4,13 @@ \alias{runIFE} \title{runIFE Removes top features identified by ML models and retrains iteratively; -returns nMCC at each iteration.} +returns MCC at each iteration.} \usage{ runIFE( ml_input_tibble, by_num = TRUE, by_vi = FALSE, - percent_removal_vec = 10 * 1:9, + percent_removal_vec = 10 * seq_len(9), mix_vec = 0, return_feats = FALSE, verbose = TRUE @@ -49,7 +49,7 @@ list along with top features removed per iteration if \code{return_feats = TRUE} \description{ runIFE Removes top features identified by ML models and retrains iteratively; -returns nMCC at each iteration. +returns MCC at each iteration. } \examples{ data(demo_ml_tibble) diff --git a/vignettes/intro.Rmd b/vignettes/intro.Rmd index cbe4f2b..b0c5182 100644 --- a/vignettes/intro.Rmd +++ b/vignettes/intro.Rmd @@ -196,7 +196,7 @@ plotPRC(results$pred) ``` ```{r plot-vi, fig.width = 5, fig.height = 4} -plotTopFeatsVI(results$fit, n_top_feats = 10) +plotTopFeatsVI(results$top_feat_tibble, n_top_feats = 10) ``` For a baseline comparison against random labels, fit a shuffled-label pipeline and compare: