diff --git a/NAMESPACE b/NAMESPACE index a6cb523..ae4698f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -27,6 +27,7 @@ export(br_show_fitted_line_2d) export(br_show_forest) export(br_show_forest_ggstats) export(br_show_forest_ggstatsplot) +export(br_show_nomogram) export(br_show_residuals) export(br_show_risk_network) export(br_show_table) diff --git a/R/04-show.R b/R/04-show.R index f943b7c..32a366b 100644 --- a/R/04-show.R +++ b/R/04-show.R @@ -527,6 +527,562 @@ br_show_table_gt <- function( t } +#' Show a nomogram for survival models +#' +#' @description +#' `r lifecycle::badge('stable')` +#' +#' This function creates a nomogram visualization for survival models, particularly +#' Cox proportional hazards models. A nomogram is a graphical calculating device +#' that provides a visual representation of a regression model to calculate +#' individualized predictions and risk scores. For survival models, it includes +#' survival probability predictions at specified time points. +#' +#' The nomogram displays: +#' - A points scale (0-100) for scoring variable contributions +#' - Individual variable scales with their ranges and point contributions +#' - A total points scale for summing individual variable points +#' - Survival probability scales for clinical predictions at specified time points +#' +#' @param breg A regression object with results (must pass `assert_breg_obj_with_results()`). +#' @param idx Length-1 vector. Index or name (focal variable) of the model. +#' Only one model is supported for nomogram visualization. +#' @param funlabel Character string for the function label on the nomogram. +#' Default is "Linear Predictor" for Cox models. +#' @param fun Function to be applied to the linear predictor. For Cox models, +#' this could be the survival function at a specific time point. +#' @param fun.at Numeric vector of points where the function should be evaluated. +#' @param lp Logical indicating whether to include linear predictor axis. +#' Default is `TRUE`. +#' @param points Logical indicating whether to include points axis for calculating +#' total score. Default is `TRUE`. +#' @param total.points Logical indicating whether to include total points axis. +#' Default is `TRUE`. +#' @param surv.at Numeric vector of time points (in years) for survival probability +#' prediction. Default is c(3, 5, 10) for 3, 5, and 10-year survival. +#' @param time.inc Time increment for baseline survival calculation. +#' Default is 365.25 (days per year). +#' @param ... Additional arguments passed to the plotting function. +#' @returns A plot +#' @export +#' @family br_show +#' @examples +#' # Cox proportional hazards model +#' lung <- survival::lung |> +#' dplyr::filter(ph.ecog != 3) +#' lung$ph.ecog <- factor(lung$ph.ecog) +#' +#' mds <- br_pipeline( +#' lung, +#' y = c("time", "status"), +#' x = c("age", "ph.ecog"), +#' x2 = "sex", +#' method = "coxph" +#' ) +#' +#' # Create nomogram for the first model (age) with default survival times +#' br_show_nomogram(mds, idx = 1) +#' +#' # Create nomogram for the second model (ph.ecog) +#' br_show_nomogram(mds, idx = 2) +#' +#' # Use with custom survival time points and labels +#' br_show_nomogram(mds, idx = "age", surv.at = c(1, 2, 5), +#' funlabel = "Log Hazard") +#' +#' @testexamples +#' expect_true(inherits(br_show_nomogram(mds, idx = 1), "ggplot")) +br_show_nomogram <- function( + breg, + idx = 1, + funlabel = "Linear Predictor", + fun = NULL, + fun.at = NULL, + lp = TRUE, + points = TRUE, + total.points = TRUE, + surv.at = c(3, 5, 10), + time.inc = 365.25, + ...) { + + assert_breg_obj_with_results(breg) + if (length(idx) != 1) { + cli_abort("length-1 {.arg idx} (integer index or a focal variable name) is required") + } + + # Get the specific model + mod <- br_get_models(breg, idx) + # br_get_models returns the model directly when idx has length 1 + + # Check if it's a supported model type + model_name <- insight::model_name(mod) + if (!model_name %in% c("coxph", "survreg")) { + cli_abort("nomogram visualization is currently only supported for survival models (coxph, survreg)") + } + + # Get model results and data + results <- br_get_results(breg, tidy = FALSE) + model_data <- br_get_data(breg) + + # For nomogram, we want to show all variables in the specific model + # The idx parameter selects which model to visualize + # We need to get the focal variable name for this model + model_names <- br_get_model_names(breg) + if (is.character(idx)) { + focal_var <- idx + } else { + focal_var <- model_names[idx] + } + + # Filter results for the specific model's focal variable + model_results <- results |> + dplyr::filter(.data$Focal_variable == focal_var) + + # Create a nomogram using ggplot2 following medical nomogram conventions + .create_nomogram_plot(mod, model_results, model_data, funlabel, fun, fun.at, lp, points, total.points, surv.at, time.inc, ...) +} + +# Helper function to create nomogram plot following medical nomogram conventions +.create_nomogram_plot <- function(model, results, data, funlabel, fun, fun.at, lp, points, total.points, surv.at, time.inc, ...) { + + # Extract model coefficients (use log scale for Cox models) + coefs <- stats::coef(model) + + # Get variable information from results - handle all variables in the model + all_vars <- results |> + dplyr::select(.data$variable, .data$estimate, .data$var_type, .data$var_class, .data$label, .data$reference_row) |> + dplyr::distinct() |> + dplyr::arrange(.data$variable) + + if (nrow(all_vars) == 0) { + cli_abort("No variables found for nomogram creation") + } + + # Calculate nomogram structure like rms package + nomogram_scales <- .build_nomogram_scales(model, all_vars, data, coefs) + + if (length(nomogram_scales$variables) == 0) { + cli_abort("Unable to create nomogram scales from model results") + } + + # Calculate survival probabilities for Cox models + survival_scales <- NULL + if (insight::model_name(model) == "coxph" && !is.null(surv.at) && length(surv.at) > 0) { + baseline_surv <- tryCatch({ + survival::survfit(model) + }, error = function(e) { + cli_warn("Could not calculate baseline survival: {e$message}") + NULL + }) + + if (!is.null(baseline_surv)) { + survival_scales <- .build_survival_scales(model, baseline_surv, surv.at, time.inc, nomogram_scales) + } + } + + # Create the nomogram plot + .plot_nomogram(nomogram_scales, survival_scales, points, total.points, lp, funlabel, surv.at) +} + +# Build nomogram scales following rms conventions +.build_nomogram_scales <- function(model, all_vars, data, coefs) { + + # Initialize scales + variable_scales <- list() + max_points <- 100 + + # Get unique variables (excluding reference rows for factors) + unique_vars <- unique(all_vars$variable) + + # Process each variable + for (var_name in unique_vars) { + var_data <- all_vars[all_vars$variable == var_name, ] + + if (!var_name %in% names(data)) { + next + } + + var_type <- var_data$var_type[1] + var_class <- var_data$var_class[1] + + if (var_type == "continuous") { + # For continuous variables, create a range and calculate points + var_range <- range(data[[var_name]], na.rm = TRUE) + + # Get coefficient for this variable + if (var_name %in% names(coefs)) { + coef_val <- coefs[var_name] + + # Create sequence of values across the range + var_seq <- seq(var_range[1], var_range[2], length.out = 11) + + # Calculate points for each value (relative to minimum) + points_contrib <- coef_val * (var_seq - var_range[1]) + + variable_scales[[var_name]] <- list( + variable = var_name, + type = "continuous", + values = var_seq, + points = points_contrib, + range = var_range, + coefficient = coef_val + ) + } + + } else if (var_class == "factor") { + # For factor variables, use levels and their estimates + + # Get all levels including reference + factor_data <- var_data[!is.na(var_data$estimate) | var_data$reference_row, ] + + if (nrow(factor_data) > 0) { + # Include reference level (0 points) + ref_row <- factor_data[factor_data$reference_row %in% TRUE, ] + non_ref_rows <- factor_data[!factor_data$reference_row %in% TRUE & !is.na(factor_data$estimate), ] + + if (nrow(ref_row) > 0) { + level_names <- c(ref_row$label[1], non_ref_rows$label) + point_values <- c(0, non_ref_rows$estimate) + } else { + level_names <- non_ref_rows$label + point_values <- non_ref_rows$estimate + } + + variable_scales[[var_name]] <- list( + variable = var_name, + type = "factor", + values = level_names, + points = point_values, + levels = level_names + ) + } + } + } + + # Normalize all points to 0-100 scale + all_points <- unlist(lapply(variable_scales, function(x) x$points)) + if (length(all_points) > 0) { + max_abs_point <- max(abs(all_points), na.rm = TRUE) + if (max_abs_point > 0) { + for (i in seq_along(variable_scales)) { + variable_scales[[i]]$points_norm <- (variable_scales[[i]]$points / max_abs_point) * max_points + } + } + } + + # Calculate total points range + total_range <- range(unlist(lapply(variable_scales, function(x) x$points_norm)), na.rm = TRUE) + + return(list( + variables = variable_scales, + total_range = total_range, + max_points = max_points + )) +} + +# Build survival probability scales +.build_survival_scales <- function(model, baseline_surv, surv.at, time.inc, nomogram_scales) { + + # Convert survival times from years to model time units + surv_times <- surv.at * time.inc + + survival_scales <- list() + + for (i in seq_along(surv.at)) { + surv_time <- surv_times[i] + + # Find closest time point in baseline survival + time_idx <- which.min(abs(baseline_surv$time - surv_time)) + if (length(time_idx) > 0 && time_idx <= length(baseline_surv$surv)) { + baseline_surv_t <- baseline_surv$surv[time_idx] + + # Create a proper range of linear predictors for this model + # Use the actual range of point contributions, scaled appropriately + total_points_range <- nomogram_scales$total_range + + # Convert points back to linear predictor scale + # Assuming the points are normalized to 0-100, we need to scale back + max_abs_coef <- max(abs(unlist(lapply(nomogram_scales$variables, function(x) { + if (!is.null(x$coefficient)) x$coefficient else max(abs(x$points), na.rm = TRUE) + })), na.rm = TRUE)) + + if (max_abs_coef > 0) { + # Scale points back to approximate linear predictor values + lp_scale_factor <- max_abs_coef / nomogram_scales$max_points + lp_min <- total_points_range[1] * lp_scale_factor + lp_max <- total_points_range[2] * lp_scale_factor + } else { + # Fallback if we can't determine scale + lp_min <- -2 + lp_max <- 2 + } + + # Create sequence of linear predictor values + lp_seq <- seq(lp_min, lp_max, length.out = 20) + + # Calculate survival probabilities: S(t|x) = S0(t)^exp(β'x) + # Ensure baseline survival is valid (between 0 and 1) + if (baseline_surv_t > 0 && baseline_surv_t <= 1) { + # Calculate survival probabilities + surv_probs <- pmax(0.01, pmin(0.99, baseline_surv_t^exp(lp_seq))) + + # Map back to points scale for positioning + points_seq <- lp_seq / lp_scale_factor + + survival_scales[[paste0("surv_", surv.at[i])]] <- list( + time_years = surv.at[i], + probabilities = surv_probs, + points_scale = points_seq, + baseline_surv = baseline_surv_t + ) + } + } + } + + return(survival_scales) +} + +# Plot the nomogram using ggplot2 +.plot_nomogram <- function(nomogram_scales, survival_scales, points, total.points, lp, funlabel, surv.at) { + + # Calculate layout parameters + n_vars <- length(nomogram_scales$variables) + current_y <- 1 + + # Determine the x-axis range for consistent alignment + x_min <- min(nomogram_scales$total_range) - 25 + x_max <- max(nomogram_scales$total_range) + 10 + axis_start <- nomogram_scales$total_range[1] + axis_end <- nomogram_scales$total_range[2] + + # Create base plot with no default data + p <- ggplot2::ggplot() + + ggplot2::theme_void() + + ggplot2::theme( + plot.title = ggplot2::element_text(hjust = 0.5, size = 14, face = "bold"), + plot.margin = ggplot2::margin(20, 20, 20, 20) + ) + + ggplot2::labs( + title = "Nomogram for Survival Prediction", + x = NULL, + y = NULL + ) + + # Track all y positions and labels for proper axis setup + y_positions <- c() + y_labels_list <- c() + + # Add variable scales + for (i in seq_along(nomogram_scales$variables)) { + var_scale <- nomogram_scales$variables[[i]] + y_pos <- current_y + y_positions <- c(y_positions, y_pos) + y_labels_list <- c(y_labels_list, var_scale$variable) + + # Add axis line for this variable (use full axis length for consistency) + p <- p + ggplot2::annotate("segment", + x = axis_start, xend = axis_end, + y = y_pos, yend = y_pos, + linewidth = 0.8, color = "black") + + # Add tick marks and labels for this variable + if (var_scale$type == "continuous") { + # For continuous variables, show selected tick marks + tick_indices <- seq(1, length(var_scale$points_norm), by = 2) + for (j in tick_indices) { + x_pos <- var_scale$points_norm[j] + label_val <- round(var_scale$values[j], 1) + + # Only add ticks within the axis range + if (x_pos >= axis_start && x_pos <= axis_end) { + # Tick mark + p <- p + ggplot2::annotate("segment", + x = x_pos, xend = x_pos, + y = y_pos - 0.1, yend = y_pos + 0.1, + linewidth = 0.5, color = "black") + + # Label above the tick + p <- p + ggplot2::annotate("text", + x = x_pos, y = y_pos + 0.25, + label = as.character(label_val), + hjust = 0.5, vjust = 0, + size = 2.8) + } + } + } else { + # For factors, show all levels + for (j in seq_along(var_scale$points_norm)) { + x_pos <- var_scale$points_norm[j] + label_val <- var_scale$values[j] + + # Only add ticks within the axis range + if (x_pos >= axis_start && x_pos <= axis_end) { + # Tick mark + p <- p + ggplot2::annotate("segment", + x = x_pos, xend = x_pos, + y = y_pos - 0.1, yend = y_pos + 0.1, + linewidth = 0.5, color = "black") + + # Label above the tick + p <- p + ggplot2::annotate("text", + x = x_pos, y = y_pos + 0.25, + label = as.character(label_val), + hjust = 0.5, vjust = 0, + size = 2.8) + } + } + } + + current_y <- current_y + 1 + } + + # Add points scale if requested + if (points) { + y_pos <- current_y + y_positions <- c(y_positions, y_pos) + y_labels_list <- c(y_labels_list, "Points") + + # Points axis line + p <- p + ggplot2::annotate("segment", + x = axis_start, xend = axis_end, + y = y_pos, yend = y_pos, + linewidth = 1, color = "black") + + # Points tick marks and labels + point_ticks <- seq(0, nomogram_scales$max_points, by = 20) + for (pt in point_ticks) { + if (pt >= axis_start && pt <= axis_end) { + # Tick mark + p <- p + ggplot2::annotate("segment", + x = pt, xend = pt, + y = y_pos - 0.1, yend = y_pos + 0.1, + linewidth = 0.5, color = "black") + + # Label + p <- p + ggplot2::annotate("text", + x = pt, y = y_pos + 0.25, + label = as.character(pt), + hjust = 0.5, vjust = 0, + size = 2.8) + } + } + + current_y <- current_y + 1 + } + + # Add total points scale if requested + if (total.points) { + y_pos <- current_y + y_positions <- c(y_positions, y_pos) + y_labels_list <- c(y_labels_list, "Total Points") + + # Total points axis line + p <- p + ggplot2::annotate("segment", + x = axis_start, xend = axis_end, + y = y_pos, yend = y_pos, + linewidth = 1, color = "black") + + # Total points ticks + point_ticks <- seq(0, nomogram_scales$max_points, by = 20) + for (pt in point_ticks) { + if (pt >= axis_start && pt <= axis_end) { + # Tick mark + p <- p + ggplot2::annotate("segment", + x = pt, xend = pt, + y = y_pos - 0.1, yend = y_pos + 0.1, + linewidth = 0.5, color = "black") + + # Label + p <- p + ggplot2::annotate("text", + x = pt, y = y_pos + 0.25, + label = as.character(pt), + hjust = 0.5, vjust = 0, + size = 2.8) + } + } + + current_y <- current_y + 1 + } + + # Add survival probability scales if available + if (!is.null(survival_scales) && length(survival_scales) > 0) { + for (i in seq_along(survival_scales)) { + surv_scale <- survival_scales[[i]] + y_pos <- current_y + y_positions <- c(y_positions, y_pos) + y_labels_list <- c(y_labels_list, paste0(surv_scale$time_years, "-Year Survival")) + + # Survival axis line + p <- p + ggplot2::annotate("segment", + x = axis_start, xend = axis_end, + y = y_pos, yend = y_pos, + linewidth = 1, color = "red") + + # Add probability tick marks - spread them out to avoid overlapping + prob_values <- c(0.9, 0.7, 0.5, 0.3, 0.1) + + # Calculate well-spaced positions for probability labels + prob_positions <- seq(axis_start, axis_end, length.out = length(prob_values)) + + for (j in seq_along(prob_values)) { + prob <- prob_values[j] + x_pos <- prob_positions[j] + + # Find the closest probability value from the scale for this position + # This ensures the labels correspond roughly to the actual survival function + if (length(surv_scale$probabilities) > 0) { + actual_prob_idx <- round(j * length(surv_scale$probabilities) / length(prob_values)) + actual_prob_idx <- max(1, min(actual_prob_idx, length(surv_scale$probabilities))) + display_prob <- surv_scale$probabilities[actual_prob_idx] + + # Only show if it's a reasonable probability (between 0.05 and 0.95) + if (!is.na(display_prob) && display_prob >= 0.05 && display_prob <= 0.95) { + # Tick mark + p <- p + ggplot2::annotate("segment", + x = x_pos, xend = x_pos, + y = y_pos - 0.1, yend = y_pos + 0.1, + linewidth = 0.5, color = "red") + + # Label with proper spacing + p <- p + ggplot2::annotate("text", + x = x_pos, y = y_pos + 0.25, + label = paste0(round(display_prob * 100), "%"), + hjust = 0.5, vjust = 0, + size = 2.8, color = "red") + } + } + } + + current_y <- current_y + 1 + } + } + + # Add y-axis labels on the left side (avoiding duplication) + for (i in seq_along(y_positions)) { + p <- p + ggplot2::annotate("text", + x = x_min, + y = y_positions[i], + label = y_labels_list[i], + hjust = 1, vjust = 0.5, + size = 3.5, fontface = "plain") + } + + # Set axis limits properly (no y-axis labels to avoid duplication) + p <- p + + ggplot2::scale_y_continuous( + limits = c(0.5, current_y - 0.5), + breaks = NULL, + labels = NULL + ) + + ggplot2::scale_x_continuous( + limits = c(x_min, x_max), + breaks = NULL, + labels = NULL + ) + + return(p) +} + #' Show residuals vs fitted plot for regression models #' #' @description @@ -695,3 +1251,4 @@ br_show_residuals <- function(breg, idx = NULL, plot_type = "fitted") { return(p) } +>>>>>>> origin/main diff --git a/fixed_nomogram_age.png b/fixed_nomogram_age.png new file mode 100644 index 0000000..a21f838 Binary files /dev/null and b/fixed_nomogram_age.png differ diff --git a/fixed_nomogram_custom.png b/fixed_nomogram_custom.png new file mode 100644 index 0000000..319aad2 Binary files /dev/null and b/fixed_nomogram_custom.png differ diff --git a/fixed_nomogram_ecog.png b/fixed_nomogram_ecog.png new file mode 100644 index 0000000..0a72ec5 Binary files /dev/null and b/fixed_nomogram_ecog.png differ diff --git a/nomogram_example.png b/nomogram_example.png new file mode 100644 index 0000000..d2c8787 Binary files /dev/null and b/nomogram_example.png differ diff --git a/nomogram_test.png b/nomogram_test.png new file mode 100644 index 0000000..a21f838 Binary files /dev/null and b/nomogram_test.png differ diff --git a/nomogram_test2.png b/nomogram_test2.png new file mode 100644 index 0000000..0a72ec5 Binary files /dev/null and b/nomogram_test2.png differ diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf new file mode 100644 index 0000000..20a25d7 Binary files /dev/null and b/tests/testthat/Rplots.pdf differ diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R index 8849056..454c61f 100644 --- a/tests/testthat/test-basic.R +++ b/tests/testthat/test-basic.R @@ -1,3 +1,55 @@ test_that("multiplication works", { expect_equal(2 * 2, 4) }) + +test_that("br_show_nomogram works with Cox regression models", { + # Create a simple Cox regression model using survival data + lung <- survival::lung + if (!is.null(lung)) { + lung <- lung[!is.na(lung$ph.ecog), ] + lung <- lung[lung$ph.ecog != 3, ] + lung$ph.ecog <- factor(lung$ph.ecog) + + # Create breg object with Cox regression + mds <- br_pipeline( + lung, + y = c("time", "status"), + x = c("age", "ph.ecog"), + x2 = "sex", + method = "coxph" + ) + + # Test that nomogram function works + expect_no_error( + plot <- br_show_nomogram(mds, idx = 1) + ) + + # Test that it returns a ggplot object + plot <- br_show_nomogram(mds, idx = 1) + expect_true(inherits(plot, "ggplot")) + } +}) + +test_that("br_show_nomogram validates input parameters", { + # Create a simple model for testing + lung <- survival::lung + if (!is.null(lung)) { + lung <- lung[!is.na(lung$ph.ecog), ] + lung <- lung[lung$ph.ecog != 3, ] + lung$ph.ecog <- factor(lung$ph.ecog) + + mds <- br_pipeline( + lung, + y = c("time", "status"), + x = c("age"), + x2 = "sex", + method = "coxph" + ) + + # Test that multiple indices are rejected + expect_error( + br_show_nomogram(mds, idx = c(1, 2)), + "length-1.*idx.*is required" + ) + } +})