diff --git a/.gitignore b/.gitignore index 614100c..6ceaedc 100644 --- a/.gitignore +++ b/.gitignore @@ -20,8 +20,10 @@ bregr_*.tar.gz bregr_models .Rapp.history -*.png + demo_fix.R +*.png +*.pdf # Test artifacts tests/testthat/Rplots.pdf diff --git a/NAMESPACE b/NAMESPACE index cea357c..b9bc067 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,6 +28,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_survival_curves) diff --git a/R/04-show-nomogram-helpers.R b/R/04-show-nomogram-helpers.R new file mode 100644 index 0000000..622b53c --- /dev/null +++ b/R/04-show-nomogram-helpers.R @@ -0,0 +1,561 @@ +# Helper functions for nomogram creation +# +# Internal functions to support the br_show_nomogram function + +# Helper function to create Cox regression nomogram +.create_coxph_nomogram <- function(model, time_points, point_range, title, subtitle, model_name) { + # Extract model coefficients and terms + coefs <- stats::coef(model) + + # Check for intercept (Cox models don't have intercept coefficients but terms may indicate one) + model_terms <- stats::terms(model) + has_intercept_term <- attr(model_terms, "intercept") == 1 + has_intercept_coef <- "(Intercept)" %in% names(coefs) + + # Cox models are semi-parametric and don't include intercept coefficients + # even if the terms object indicates an intercept is present + if (has_intercept_term && !has_intercept_coef) { + cli::cli_inform("Cox model: intercept term present but no intercept coefficient (as expected for semi-parametric models)") + } else if (has_intercept_coef) { + # This would be unusual for a Cox model but handle it anyway + cli::cli_inform("removing intercept coefficient from Cox model") + coefs <- coefs[names(coefs) != "(Intercept)"] + } + + # Handle NA coefficients (for singular fits) while preserving coefficient-term correspondence + if (any(is.na(coefs))) { + na_coefs <- names(coefs)[is.na(coefs)] + cli::cli_inform("removing {length(na_coefs)} NA coefficient{?s} due to singular fit: {.val {na_coefs}}") + coefs <- coefs[!is.na(coefs)] + } + + if (length(coefs) == 0) { + cli::cli_abort("no valid coefficients found in the model") + } + + # Get model frame to understand variable ranges + model_frame <- broom.helpers::model_get_model_frame(model) + + # Get baseline survival for more accurate survival probability calculations + baseline_surv <- tryCatch( + { + survival::survfit(model) + }, + error = function(e) NULL + ) + + # Scale coefficients to point range + max_abs_coef <- max(abs(coefs)) + point_scale_factor <- diff(point_range) / (2 * max_abs_coef) + + # Create scales for each variable + nom_data <- list() + y_position <- length(coefs) + 1 # Start from top (removed redundant Points scale) + + # Variable scales - improved handling with proper line representations + for (i in seq_along(coefs)) { + var_name <- names(coefs)[i] + coef_val <- coefs[i] + + # Try to match coefficient name to original variable + # Handle factor variables with level suffixes + base_var_name <- gsub("^(.+?)[0-9]+$", "\\1", var_name) + + # Find the base variable in model frame + matching_vars <- names(model_frame)[grepl(paste0("^", base_var_name), names(model_frame))] + if (length(matching_vars) == 0) { + matching_vars <- names(model_frame)[grepl(var_name, names(model_frame), fixed = TRUE)] + } + + if (length(matching_vars) > 0) { + actual_var_name <- matching_vars[1] + var_data <- model_frame[[actual_var_name]] + + if (is.numeric(var_data)) { + # Continuous variable - create a proper scale with connecting line + var_range <- range(var_data, na.rm = TRUE) + # More points for a smooth connecting line + n_line_points <- 21 + var_values_line <- seq(var_range[1], var_range[2], length.out = n_line_points) + + # Create evenly spaced x-positions for the variable scale (NOT scaled by coefficient) + points_line <- seq(point_range[1] + diff(point_range) * 0.1, + point_range[2] - diff(point_range) * 0.1, + length.out = n_line_points + ) + + # Create labels and tick marks at meaningful intervals (every 5th point) + tick_indices <- seq(1, n_line_points, by = 5) + labels_line <- rep("", n_line_points) + labels_line[tick_indices] <- round(var_values_line[tick_indices], 1) + is_tick_line <- rep(FALSE, n_line_points) + is_tick_line[tick_indices] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = points_line, + label = labels_line, + var_name = actual_var_name, + type = "variable", + is_tick = is_tick_line, + stringsAsFactors = FALSE + ) + } else if (is.factor(var_data)) { + # Categorical variable - show line segment between reference and level + levels_found <- levels(var_data) + + # For factor variables, we need to map coefficient to the right level + if (grepl("[0-9]+$", var_name)) { + # Extract level number from coefficient name + level_num <- as.numeric(gsub(".*([0-9]+)$", "\\1", var_name)) + if (level_num <= length(levels_found)) { + # Reference level gets baseline points + ref_points <- point_range[1] + diff(point_range) * 0.2 + # Current level gets points based on coefficient + level_points <- ref_points + coef_val * point_scale_factor + + # Create line segment between the two points + n_line_points <- 11 + line_x <- seq(ref_points, level_points, length.out = n_line_points) + line_labels <- rep("", n_line_points) + line_labels[1] <- paste0(levels_found[1], " (ref)") + line_labels[n_line_points] <- levels_found[level_num] + line_ticks <- rep(FALSE, n_line_points) + line_ticks[c(1, n_line_points)] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = line_x, + label = line_labels, + var_name = actual_var_name, + type = "variable", + is_tick = line_ticks, + stringsAsFactors = FALSE + ) + } + } else { + # Simple two-level case + ref_points <- point_range[1] + diff(point_range) * 0.2 + level_points <- ref_points + coef_val * point_scale_factor + + # Create line segment between the two points + n_line_points <- 11 + line_x <- seq(ref_points, level_points, length.out = n_line_points) + line_labels <- rep("", n_line_points) + line_labels[1] <- paste0(levels_found[1], " (ref)") + line_labels[n_line_points] <- if (length(levels_found) > 1) levels_found[2] else "Other" + line_ticks <- rep(FALSE, n_line_points) + line_ticks[c(1, n_line_points)] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = line_x, + label = line_labels, + var_name = actual_var_name, + type = "variable", + is_tick = line_ticks, + stringsAsFactors = FALSE + ) + } + } + } else { + # If variable not found in model frame, create a generic scale + n_line_points <- 11 + points_vals <- seq(point_range[1] + diff(point_range) * 0.1, + point_range[2] - diff(point_range) * 0.1, + length.out = n_line_points + ) + labels_vals <- rep("", n_line_points) + labels_vals[c(1, 6, 11)] <- c("Low", "Medium", "High") + tick_vals <- rep(FALSE, n_line_points) + tick_vals[c(1, 6, 11)] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = points_vals, + label = labels_vals, + var_name = var_name, + type = "variable", + is_tick = tick_vals, + stringsAsFactors = FALSE + ) + } + + y_position <- y_position - 1 + } + + # Total points scale with proper line + y_position <- y_position - 0.5 + n_total_points <- 21 # More points for smoother line + total_points <- seq(point_range[1], point_range[2], length.out = n_total_points) + total_labels <- rep("", n_total_points) + # Show labels every 4th point + label_indices <- seq(1, n_total_points, by = 4) + total_labels[label_indices] <- total_points[label_indices] + total_ticks <- rep(FALSE, n_total_points) + total_ticks[label_indices] <- TRUE + + nom_data[[length(nom_data) + 1]] <- data.frame( + y = y_position, + x = total_points, + label = total_labels, + var_name = "Total Points", + type = "scale", + is_tick = total_ticks, + stringsAsFactors = FALSE + ) + + # Survival probability scales for each time point - improved calculation + if (length(time_points) > 0) { + for (j in seq_along(time_points)) { + y_position <- y_position - 1 + + # More accurate survival probability calculation + if (!is.null(baseline_surv)) { + # Convert months to days for proper time matching + time_in_days <- time_points[j] * 30.44 # Average days per month + time_idx <- which.min(abs(baseline_surv$time - time_in_days)) + + if (length(time_idx) > 0 && time_idx <= length(baseline_surv$surv)) { + baseline_surv_at_time <- baseline_surv$surv[time_idx] + + # Calculate survival probabilities based on linear predictor + # Linear predictor range corresponding to the point range + lp_range <- (total_points - mean(point_range)) / point_scale_factor + survival_probs <- baseline_surv_at_time^exp(lp_range) + survival_probs <- pmax(0.01, pmin(0.99, survival_probs)) + } else { + # Fallback: use average baseline hazard estimation + # Convert months to hazard time scale + hazard_time <- time_points[j] / 12 # Convert to years for hazard calculation + lp_range <- (total_points - mean(point_range)) / point_scale_factor + survival_probs <- exp(-0.5 * hazard_time * exp(lp_range)) # More realistic baseline + survival_probs <- pmax(0.01, pmin(0.99, survival_probs)) + } + } else { + # Fallback calculation with more realistic baseline hazard + hazard_time <- time_points[j] / 12 # Convert to years + lp_range <- (total_points - mean(point_range)) / point_scale_factor + survival_probs <- exp(-0.5 * hazard_time * exp(lp_range)) # More realistic baseline + survival_probs <- pmax(0.01, pmin(0.99, survival_probs)) + } + + # Create survival probability labels with fewer overlapping points + surv_labels <- rep("", length(total_points)) + surv_labels[label_indices] <- paste0(round(survival_probs[label_indices] * 100, 1), "%") + surv_ticks <- rep(FALSE, length(total_points)) + surv_ticks[label_indices] <- TRUE + + nom_data[[length(nom_data) + 1]] <- data.frame( + y = y_position, + x = total_points, + label = surv_labels, + var_name = paste0(time_points[j], "-month survival"), + type = "survival", + is_tick = surv_ticks, + stringsAsFactors = FALSE + ) + } + } + + # Combine all data + plot_data <- do.call(rbind, nom_data) + + # Create the plot with improved styling + if (is.null(title)) { + title <- paste("Nomogram for", model_name, "Model") + } + + p <- ggplot2::ggplot(plot_data, ggplot2::aes(x = .data$x, y = .data$y)) + + # Add subtle grid lines for easier reading + ggplot2::geom_vline( + xintercept = seq(point_range[1], point_range[2], by = 10), + color = "grey90", linewidth = 0.3 + ) + + # Add connecting lines for ALL scales + ggplot2::geom_line(ggplot2::aes(group = .data$y), + linewidth = 0.6, color = "black" + ) + + # Add tick marks only for labeled points + ggplot2::geom_point( + data = plot_data[plot_data$is_tick, ], + size = 1.8, color = "black" + ) + + # Add labels only for tick marks, with improved positioning to prevent overlap + ggplot2::geom_text( + data = plot_data[plot_data$is_tick & plot_data$label != "", ], + ggplot2::aes(label = .data$label), + vjust = -1.2, hjust = 0.5, size = 3, color = "black" + ) + + ggplot2::scale_y_continuous( + breaks = unique(plot_data$y), + labels = unique(plot_data$var_name)[order(unique(plot_data$y), decreasing = TRUE)], + limits = c(min(plot_data$y) - 0.5, max(plot_data$y) + 0.5) + ) + + ggplot2::scale_x_continuous( + limits = c(point_range[1] - 5, point_range[2] + 5) + ) + + ggplot2::labs( + title = title, + subtitle = subtitle, + x = "", + y = "" + ) + + ggplot2::theme_minimal() + + ggplot2::theme( + axis.text.x = ggplot2::element_blank(), + axis.ticks.x = ggplot2::element_blank(), + axis.text.y = ggplot2::element_text(size = 10, hjust = 1), + panel.grid = ggplot2::element_blank(), + plot.title = ggplot2::element_text(hjust = 0.5, size = 14, face = "bold"), + plot.subtitle = ggplot2::element_text(hjust = 0.5, size = 12), + panel.background = ggplot2::element_rect(fill = "white", color = "black") + ) + + return(p) +} + +# Helper function to create linear/GLM nomogram +.create_lm_nomogram <- function(model, fun_at, point_range, title, subtitle, model_name) { + # Extract model coefficients and terms + coefs <- stats::coef(model) + + # Check if model has intercept using proper method + model_terms <- stats::terms(model) + has_intercept <- attr(model_terms, "intercept") == 1 + + # Handle intercept removal if present + if (has_intercept) { + intercept <- coefs[1] + coefs <- coefs[-1] # Remove intercept + } else { + cli::cli_inform("model fitted without intercept") + } + + # Handle NA coefficients (for singular fits) while preserving coefficient-term correspondence + if (any(is.na(coefs))) { + na_coefs <- names(coefs)[is.na(coefs)] + cli::cli_inform("removing {length(na_coefs)} NA coefficient{?s} due to singular fit: {.val {na_coefs}}") + coefs <- coefs[!is.na(coefs)] + } + + if (length(coefs) == 0) { + cli::cli_abort("no valid coefficients found in the model") + } + + # Get model frame to understand variable ranges + model_frame <- broom.helpers::model_get_model_frame(model) + + # Set default prediction values if not provided + if (is.null(fun_at)) { + pred_range <- range(stats::predict(model), na.rm = TRUE) + fun_at <- seq(pred_range[1], pred_range[2], length.out = 5) + fun_at <- round(fun_at, 2) + } + + # Scale coefficients to point range + max_abs_coef <- max(abs(coefs)) + point_scale_factor <- diff(point_range) / (2 * max_abs_coef) + + # Create scales for each variable + nom_data <- list() + y_position <- length(coefs) + 1 # Start from top (removed redundant Points scale) + + # Variable scales - improved handling similar to Cox model + for (i in seq_along(coefs)) { + var_name <- names(coefs)[i] + coef_val <- coefs[i] + + # Get variable data from model frame + if (var_name %in% colnames(model_frame)) { + var_data <- model_frame[[var_name]] + + if (is.numeric(var_data)) { + # Continuous variable - create proper scale with connecting line + var_range <- range(var_data, na.rm = TRUE) + n_line_points <- 21 + var_values <- seq(var_range[1], var_range[2], length.out = n_line_points) + + # Create evenly spaced x-positions for the variable scale (NOT scaled by coefficient) + points <- seq(point_range[1] + diff(point_range) * 0.1, + point_range[2] - diff(point_range) * 0.1, + length.out = n_line_points + ) + + # Create labels at meaningful intervals + tick_indices <- seq(1, n_line_points, by = 5) + labels_line <- rep("", n_line_points) + labels_line[tick_indices] <- round(var_values[tick_indices], 1) + is_tick_line <- rep(FALSE, n_line_points) + is_tick_line[tick_indices] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = points, + label = labels_line, + var_name = var_name, + type = "variable", + is_tick = is_tick_line, + stringsAsFactors = FALSE + ) + } else { + # Categorical variable - improved handling + if (is.factor(var_data)) { + levels_found <- levels(var_data) + if (length(levels_found) > 1) { + # Reference level gets baseline points + ref_level <- levels_found[1] + ref_points <- point_range[1] + diff(point_range) * 0.2 + level_points <- ref_points + coef_val * point_scale_factor + + # Create line segment + n_line_points <- 11 + line_x <- seq(ref_points, level_points, length.out = n_line_points) + line_labels <- rep("", n_line_points) + line_labels[1] <- paste0(ref_level, " (ref)") + line_labels[n_line_points] <- if (length(levels_found) > 1) levels_found[2] else "Other" + line_ticks <- rep(FALSE, n_line_points) + line_ticks[c(1, n_line_points)] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = line_x, + label = line_labels, + var_name = var_name, + type = "variable", + is_tick = line_ticks, + stringsAsFactors = FALSE + ) + } + } + } + } else { + # Generic scale for unknown variables + n_line_points <- 11 + points_vals <- seq(point_range[1] + diff(point_range) * 0.1, + point_range[2] - diff(point_range) * 0.1, + length.out = n_line_points + ) + labels_vals <- rep("", n_line_points) + labels_vals[c(1, 6, 11)] <- c("Low", "Medium", "High") + tick_vals <- rep(FALSE, n_line_points) + tick_vals[c(1, 6, 11)] <- TRUE + + nom_data[[i]] <- data.frame( + y = y_position, + x = points_vals, + label = labels_vals, + var_name = var_name, + type = "variable", + is_tick = tick_vals, + stringsAsFactors = FALSE + ) + } + + y_position <- y_position - 1 + } + + # Total points scale with proper line + y_position <- y_position - 0.5 + n_total_points <- 21 + total_points <- seq(point_range[1], point_range[2], length.out = n_total_points) + total_labels <- rep("", n_total_points) + label_indices <- seq(1, n_total_points, by = 4) + total_labels[label_indices] <- total_points[label_indices] + total_ticks <- rep(FALSE, n_total_points) + total_ticks[label_indices] <- TRUE + + nom_data[[length(nom_data) + 1]] <- data.frame( + y = y_position, + x = total_points, + label = total_labels, + var_name = "Total Points", + type = "scale", + is_tick = total_ticks, + stringsAsFactors = FALSE + ) + + # Prediction scale - improved mapping + y_position <- y_position - 1 + + # Map prediction values to points more accurately + n_pred_points <- 21 + pred_points <- seq(point_range[1], point_range[2], length.out = n_pred_points) + + # Create prediction labels at the specified fun_at values + pred_labels <- rep("", n_pred_points) + pred_ticks <- rep(FALSE, n_pred_points) + + # Place fun_at values evenly across the prediction scale + if (length(fun_at) > 0) { + pred_indices <- round(seq(1, n_pred_points, length.out = length(fun_at))) + pred_labels[pred_indices] <- fun_at + pred_ticks[pred_indices] <- TRUE + } + + nom_data[[length(nom_data) + 1]] <- data.frame( + y = y_position, + x = pred_points, + label = pred_labels, + var_name = "Predicted Value", + type = "prediction", + is_tick = pred_ticks, + stringsAsFactors = FALSE + ) + + # Combine all data + plot_data <- do.call(rbind, nom_data) + + # Create the plot with improved styling + if (is.null(title)) { + title <- paste("Nomogram for", model_name, "Model") + } + + p <- ggplot2::ggplot(plot_data, ggplot2::aes(x = .data$x, y = .data$y)) + + # Add subtle grid lines for easier reading + ggplot2::geom_vline( + xintercept = seq(point_range[1], point_range[2], by = 10), + color = "grey90", linewidth = 0.3 + ) + + # Add connecting lines for ALL scales + ggplot2::geom_line(ggplot2::aes(group = .data$y), + linewidth = 0.6, color = "black" + ) + + # Add tick marks only for labeled points + ggplot2::geom_point( + data = plot_data[plot_data$is_tick, ], + size = 1.8, color = "black" + ) + + # Add labels only for tick marks, with improved positioning to prevent overlap + ggplot2::geom_text( + data = plot_data[plot_data$is_tick & plot_data$label != "", ], + ggplot2::aes(label = .data$label), + vjust = -1.2, hjust = 0.5, size = 3, color = "black" + ) + + ggplot2::scale_y_continuous( + breaks = unique(plot_data$y), + labels = unique(plot_data$var_name)[order(unique(plot_data$y), decreasing = TRUE)], + limits = c(min(plot_data$y) - 0.5, max(plot_data$y) + 0.5) + ) + + ggplot2::scale_x_continuous( + limits = c(point_range[1] - 5, point_range[2] + 5) + ) + + ggplot2::labs( + title = title, + subtitle = subtitle, + x = "", + y = "" + ) + + ggplot2::theme_minimal() + + ggplot2::theme( + axis.text.x = ggplot2::element_blank(), + axis.ticks.x = ggplot2::element_blank(), + axis.text.y = ggplot2::element_text(size = 10, hjust = 1), + panel.grid = ggplot2::element_blank(), + plot.title = ggplot2::element_text(hjust = 0.5, size = 14, face = "bold"), + plot.subtitle = ggplot2::element_text(hjust = 0.5, size = 12), + panel.background = ggplot2::element_rect(fill = "white", color = "black") + ) + + return(p) +} diff --git a/R/04-show.R b/R/04-show.R index db6e9de..8bc7212 100644 --- a/R/04-show.R +++ b/R/04-show.R @@ -867,3 +867,93 @@ br_show_residuals <- function(breg, idx = NULL, plot_type = "fitted") { return(p) } + + +#' Show nomogram for regression models +#' +#' @description +#' `r lifecycle::badge('experimental')` +#' +#' Creates a nomogram (graphical calculator) for regression models, particularly +#' useful for Cox proportional hazards models. A nomogram allows visual calculation +#' of predicted outcomes by assigning points to variable values and summing them +#' to get total points that correspond to predicted probabilities. +#' +#' @param breg A `breg` object with fitted regression models. +#' @param idx Index or name of the model to use for the nomogram. +#' If NULL, uses the first model. +#' @param time_points For Cox models, time points at which to show survival probabilities. +#' Default is c(12, 24, 36) representing months. +#' @param fun_at For non-survival models, the function values at which to show predictions. +#' @param point_range Range of points to use in the nomogram scale. Default is c(0, 100). +#' @param title Plot title. If NULL, generates automatic title. +#' @param subtitle Plot subtitle. +#' @returns A ggplot2 object showing the nomogram. +#' @export +#' @family br_show +#' @examples +#' \donttest{ +#' # Cox regression nomogram +#' +#' 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" +#' ) +#' p <- br_show_nomogram(mds) +#' p +#' +#' +#' # Linear regression nomogram +#' mds_lm <- br_pipeline( +#' mtcars, +#' y = "mpg", +#' x = c("hp", "wt"), +#' x2 = "vs", +#' method = "gaussian" +#' ) +#' p2 <- br_show_nomogram(mds_lm, fun_at = c(15, 20, 25, 30)) +#' p2 +#' } +#' @testexamples +#' expect_s3_class(p, "ggplot") +#' expect_s3_class(p2, "ggplot") +br_show_nomogram <- function(breg, + idx = NULL, + time_points = c(12, 24, 36), + fun_at = NULL, + point_range = c(0, 100), + title = NULL, + subtitle = NULL) { + assert_breg_obj_with_results(breg) + + # Get the model to use + if (is.null(idx)) { + cli::cli_inform("{.arg idx} not set, use the first model") + idx <- 1 + } else { + if (length(idx) != 1) { + cli::cli_abort("please specify one model") + } + } + + model <- br_get_models(breg, idx) + model_name <- if (!rlang::is_string(idx)) { + br_get_model_names(breg)[idx] + } else { + idx + } + + # Check model type and dispatch to appropriate function + if (inherits(model, "coxph")) { + .create_coxph_nomogram(model, time_points, point_range, title, subtitle, model_name) + } else if (inherits(model, c("lm", "glm"))) { + .create_lm_nomogram(model, fun_at, point_range, title, subtitle, model_name) + } else { + cli::cli_abort("Nomograms are currently supported for Cox regression (coxph) and linear/generalized linear models (lm/glm)") + } +} diff --git a/man/br_show_fitted_line.Rd b/man/br_show_fitted_line.Rd index f7a9abc..b88b253 100644 --- a/man/br_show_fitted_line.Rd +++ b/man/br_show_fitted_line.Rd @@ -45,6 +45,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/man/br_show_fitted_line_2d.Rd b/man/br_show_fitted_line_2d.Rd index d9bbbc5..9ebe13e 100644 --- a/man/br_show_fitted_line_2d.Rd +++ b/man/br_show_fitted_line_2d.Rd @@ -42,6 +42,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/man/br_show_forest.Rd b/man/br_show_forest.Rd index 15aa8fa..0971750 100644 --- a/man/br_show_forest.Rd +++ b/man/br_show_forest.Rd @@ -71,6 +71,7 @@ Other br_show: \code{\link{br_show_fitted_line_2d}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/man/br_show_forest_ggstats.Rd b/man/br_show_forest_ggstats.Rd index c8014d5..5d1ded3 100644 --- a/man/br_show_forest_ggstats.Rd +++ b/man/br_show_forest_ggstats.Rd @@ -39,6 +39,7 @@ Other br_show: \code{\link{br_show_fitted_line_2d}()}, \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/man/br_show_forest_ggstatsplot.Rd b/man/br_show_forest_ggstatsplot.Rd index 66d139b..a6ef319 100644 --- a/man/br_show_forest_ggstatsplot.Rd +++ b/man/br_show_forest_ggstatsplot.Rd @@ -41,6 +41,7 @@ Other br_show: \code{\link{br_show_fitted_line_2d}()}, \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/man/br_show_nomogram.Rd b/man/br_show_nomogram.Rd new file mode 100644 index 0000000..ee5fecf --- /dev/null +++ b/man/br_show_nomogram.Rd @@ -0,0 +1,87 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/04-show.R +\name{br_show_nomogram} +\alias{br_show_nomogram} +\title{Show nomogram for regression models} +\usage{ +br_show_nomogram( + breg, + idx = NULL, + time_points = c(12, 24, 36), + fun_at = NULL, + point_range = c(0, 100), + title = NULL, + subtitle = NULL +) +} +\arguments{ +\item{breg}{A \code{breg} object with fitted regression models.} + +\item{idx}{Index or name of the model to use for the nomogram. +If NULL, uses the first model.} + +\item{time_points}{For Cox models, time points at which to show survival probabilities. +Default is c(12, 24, 36) representing months.} + +\item{fun_at}{For non-survival models, the function values at which to show predictions.} + +\item{point_range}{Range of points to use in the nomogram scale. Default is c(0, 100).} + +\item{title}{Plot title. If NULL, generates automatic title.} + +\item{subtitle}{Plot subtitle.} +} +\value{ +A ggplot2 object showing the nomogram. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +Creates a nomogram (graphical calculator) for regression models, particularly +useful for Cox proportional hazards models. A nomogram allows visual calculation +of predicted outcomes by assigning points to variable values and summing them +to get total points that correspond to predicted probabilities. +} +\examples{ +\donttest{ +# Cox regression nomogram + +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" +) +p <- br_show_nomogram(mds) +p + + +# Linear regression nomogram +mds_lm <- br_pipeline( + mtcars, + y = "mpg", + x = c("hp", "wt"), + x2 = "vs", + method = "gaussian" +) +p2 <- br_show_nomogram(mds_lm, fun_at = c(15, 20, 25, 30)) +p2 +} +} +\seealso{ +Other br_show: +\code{\link{br_show_fitted_line}()}, +\code{\link{br_show_fitted_line_2d}()}, +\code{\link{br_show_forest}()}, +\code{\link{br_show_forest_ggstats}()}, +\code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_residuals}()}, +\code{\link{br_show_risk_network}()}, +\code{\link{br_show_survival_curves}()}, +\code{\link{br_show_table}()}, +\code{\link{br_show_table_gt}()} +} +\concept{br_show} diff --git a/man/br_show_residuals.Rd b/man/br_show_residuals.Rd index fe54d0e..4a313de 100644 --- a/man/br_show_residuals.Rd +++ b/man/br_show_residuals.Rd @@ -55,6 +55,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, \code{\link{br_show_table}()}, diff --git a/man/br_show_risk_network.Rd b/man/br_show_risk_network.Rd index f1af6df..b8c7845 100644 --- a/man/br_show_risk_network.Rd +++ b/man/br_show_risk_network.Rd @@ -37,6 +37,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_survival_curves}()}, \code{\link{br_show_table}()}, diff --git a/man/br_show_survival_curves.Rd b/man/br_show_survival_curves.Rd index 4eaed7c..2f4705a 100644 --- a/man/br_show_survival_curves.Rd +++ b/man/br_show_survival_curves.Rd @@ -62,6 +62,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_table}()}, diff --git a/man/br_show_table.Rd b/man/br_show_table.Rd index c989d3d..fec93f1 100644 --- a/man/br_show_table.Rd +++ b/man/br_show_table.Rd @@ -50,6 +50,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/man/br_show_table_gt.Rd b/man/br_show_table_gt.Rd index d9fe4de..bdbfc31 100644 --- a/man/br_show_table_gt.Rd +++ b/man/br_show_table_gt.Rd @@ -48,6 +48,7 @@ Other br_show: \code{\link{br_show_forest}()}, \code{\link{br_show_forest_ggstats}()}, \code{\link{br_show_forest_ggstatsplot}()}, +\code{\link{br_show_nomogram}()}, \code{\link{br_show_residuals}()}, \code{\link{br_show_risk_network}()}, \code{\link{br_show_survival_curves}()}, diff --git a/tests/testthat/test-nomogram.R b/tests/testthat/test-nomogram.R new file mode 100644 index 0000000..35e2b87 --- /dev/null +++ b/tests/testthat/test-nomogram.R @@ -0,0 +1,152 @@ +test_that("br_show_nomogram works for Cox models", { + skip_if_not_installed("survival") + + # Create Cox 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" + ) + + # Test basic nomogram + p <- br_show_nomogram(mds) + expect_s3_class(p, "ggplot") + + # Test with custom time points + p2 <- br_show_nomogram(mds, time_points = c(6, 12)) + expect_s3_class(p2, "ggplot") + + # Test with specific model index + p3 <- br_show_nomogram(mds, idx = 1) + expect_s3_class(p3, "ggplot") +}) + +test_that("br_show_nomogram works for linear models", { + # Create linear model + mds_lm <- br_pipeline( + mtcars, + y = "mpg", + x = c("hp", "wt"), + x2 = "vs", + method = "gaussian" + ) + + # Test basic nomogram + p <- br_show_nomogram(mds_lm) + expect_s3_class(p, "ggplot") + + # Test with custom prediction values + p2 <- br_show_nomogram(mds_lm, fun_at = c(15, 20, 25, 30)) + expect_s3_class(p2, "ggplot") +}) + +test_that("br_show_nomogram handles unsupported models", { + # This would test an unsupported model type if we had one + # For now, we can test error handling with invalid idx + mds_lm <- br_pipeline( + mtcars, + y = "mpg", + x = c("hp", "wt"), + x2 = "vs", + method = "gaussian" + ) + + # Test error for multiple indices + expect_error(br_show_nomogram(mds_lm, idx = c(1, 2))) +}) + +test_that("br_show_nomogram produces correct plot structure", { + # Create simple model for structure testing + mds_lm <- br_pipeline( + mtcars[1:10, ], # Small dataset for faster testing + y = "mpg", + x = "hp", + x2 = "vs", + method = "gaussian" + ) + + p <- br_show_nomogram(mds_lm) + + # Test that plot has the expected structure + expect_s3_class(p, "ggplot") + expect_true("data" %in% names(p)) + expect_true("layers" %in% names(p)) + expect_true("theme" %in% names(p)) +}) + +test_that("br_show_nomogram handles models without intercepts", { + # Create model data that can handle no-intercept fitting + test_data <- mtcars[1:15, ] + test_data$vs <- factor(test_data$vs) + + # Fit model without intercept manually to test coefficient handling + no_int_model <- lm(mpg ~ hp + wt - 1, data = test_data) + + # Test that our coefficient handling logic works + coefs <- stats::coef(no_int_model) + model_terms <- stats::terms(no_int_model) + has_intercept <- attr(model_terms, "intercept") == 1 + + expect_false(has_intercept) + expect_false("(Intercept)" %in% names(coefs)) + expect_true(length(coefs) >= 2) + expect_false(any(is.na(coefs))) +}) + +test_that("br_show_nomogram handles singular coefficient matrices", { + # Create data with collinear variables to test NA coefficient handling + singular_data <- data.frame( + y = 1:10, + x1 = 1:10, + x2 = 2 * (1:10), # x2 = 2 * x1, creating collinearity + x3 = rnorm(10) + ) + + # Fit model that will have singular coefficients + singular_model <- lm(y ~ x1 + x2 + x3, data = singular_data) + coefs <- stats::coef(singular_model) + + # Check that we can handle NA coefficients + if (any(is.na(coefs))) { + # Test that our NA handling preserves coefficient-term correspondence + non_na_coefs <- coefs[!is.na(coefs)] + expect_true(length(non_na_coefs) > 0) + expect_true(all(!is.na(non_na_coefs))) + } +}) + +test_that("br_show_nomogram handles Cox model intercept behavior correctly", { + skip_if_not_installed("survival") + + # Create Cox model to test intercept handling + 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" + ) + + model <- br_get_models(mds, 1) + + # Test that Cox model behavior is as expected + coefs <- stats::coef(model) + model_terms <- stats::terms(model) + has_intercept_term <- attr(model_terms, "intercept") == 1 + has_intercept_coef <- "(Intercept)" %in% names(coefs) + + # Cox models have intercept in terms but not in coefficients + expect_true(has_intercept_term) + expect_false(has_intercept_coef) + + # Test that nomogram creation works correctly (may include informative messages) + # The function should work without errors regardless of messages + suppressMessages(p <- br_show_nomogram(mds)) + expect_s3_class(p, "ggplot") +}) diff --git a/tests/testthat/test-roxytest-testexamples-04-show.R b/tests/testthat/test-roxytest-testexamples-04-show.R index cfaef7a..a87e5b6 100644 --- a/tests/testthat/test-roxytest-testexamples-04-show.R +++ b/tests/testthat/test-roxytest-testexamples-04-show.R @@ -141,3 +141,37 @@ test_that("Function br_show_residuals() @ L741", { expect_s3_class(br_show_residuals(m, idx = 1), "ggplot") }) + +test_that("Function br_show_nomogram() @ L925", { + + + # Cox regression nomogram + + 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" + ) + p <- br_show_nomogram(mds) + p + + + # Linear regression nomogram + mds_lm <- br_pipeline( + mtcars, + y = "mpg", + x = c("hp", "wt"), + x2 = "vs", + method = "gaussian" + ) + p2 <- br_show_nomogram(mds_lm, fun_at = c(15, 20, 25, 30)) + p2 + + expect_s3_class(p, "ggplot") + expect_s3_class(p2, "ggplot") +}) +