From 98a0ff82c031559ecb2a645bdbaf671d6c53270f Mon Sep 17 00:00:00 2001 From: Nick Kennedy Date: Thu, 31 Aug 2023 17:08:15 +0100 Subject: [PATCH] Fix labels for `coxph` models --- R/default_forest_panels.R | 23 ++++++++++++++--- R/forest_model.R | 49 +++++++++++++++++------------------- R/model_frame_coxph_simple.R | 21 ++++++++++++++++ 3 files changed, 63 insertions(+), 30 deletions(-) create mode 100644 R/model_frame_coxph_simple.R diff --git a/R/default_forest_panels.R b/R/default_forest_panels.R index ea85170..fd0bc09 100644 --- a/R/default_forest_panels.R +++ b/R/default_forest_panels.R @@ -58,6 +58,24 @@ default_forest_panels <- function(model = NULL, factor_separate_line = FALSE, me measure <- "Estimate" } } + + max_interaction_terms <- 0 + if (!is.null(model)) { + forest_terms_basic <- make_forest_terms_basic(model) + if ( + any( + factor_separate_line * + forest_terms_basic$is_interaction * vapply(forest_terms_basic$interaction_terms_are_factors, sum, integer(1)) > 0 + ) + ) { + warning("`factor_separate_line = TRUE` is not supported when there are interacting terms that include one or more factors") + factor_separate_line <- FALSE + } + max_interaction_terms <- max( + vapply(forest_terms_basic$interaction_terms_are_factors, sum, integer(1)) + ) + } + is_na <- function(x) sapply(x, is.na) panels <- list( forest_panel(width = 0.03), @@ -83,10 +101,6 @@ default_forest_panels <- function(model = NULL, factor_separate_line = FALSE, me ), forest_panel(width = 0.03) ) - forest_terms_basic <- make_forest_terms_basic(model) - max_interaction_terms <- max( - vapply(forest_terms_basic$interaction_terms_are_factors, sum, integer(1)) - ) if (max_interaction_terms > 1) { interacting_level_panels <- lapply( seq_len(max_interaction_terms), @@ -98,6 +112,7 @@ default_forest_panels <- function(model = NULL, factor_separate_line = FALSE, me ) panels <- c(panels[1:2], interacting_level_panels, panels[-(1:3)]) } + if (factor_separate_line) { panels[[2]][c("width", "width_group", "width_fixed")] <- list(0.02, 1, TRUE) panels[[3]][c("width", "width_group", "width_fixed")] <- list(0.02, 1, TRUE) diff --git a/R/forest_model.R b/R/forest_model.R index b2bc3d3..9ccc6fb 100644 --- a/R/forest_model.R +++ b/R/forest_model.R @@ -160,9 +160,6 @@ forest_model <- function(model, exponentiate <- inherits(model_list[[1]], "coxph") || (inherits(model_list[[1]], "glm") && model_list[[1]]$family$link == "logit") } - if (missing(panels)) { - panels <- default_forest_panels(model_list[[1]], factor_separate_line = factor_separate_line) - } } else { if (is.null(exponentiate)) { exponentiate <- inherits(model, "coxph") || @@ -172,33 +169,28 @@ forest_model <- function(model, if (exponentiate) trans <- exp else trans <- I - forest_terms_basic <- make_forest_terms_basic(model) - - if ( - any( - factor_separate_line * - forest_terms_basic$is_interaction * vapply(forest_terms_basic$interaction_terms_are_factors, sum, integer(1)) > 0 - ) - ) { - warn("`factor_separate_line = TRUE` is not supported when there are interacting terms that include one or more factors") - factor_separate_line <- FALSE - } - - if (is.null(panels)) { - panels <- default_forest_panels(model, factor_separate_line = factor_separate_line) + if (!is.null(panels) && !is.list(panels)) { + stop("panels should be a list") } - stopifnot(is.list(panels)) + make_forest_terms <- function(model) { + forest_terms_basic <- make_forest_terms_basic(model) - make_forest_terms <- function(model, forest_terms_basic) { tidy_model <- broom::tidy(model, conf.int = TRUE) data <- stats::model.frame(model) mmtrx <- stats::model.matrix(model) + if (inherits(model, "coxph")) { + # Cope with the way survival:::model.frame.coxph drops attributes from factors + data_for_labels <- model_frame_coxph_simple(model) + } else { + data_for_labels <- data + } + forest_labels <- tibble::tibble( - variable = names(data), + variable = names(data_for_labels), label = vapply( - data, + data_for_labels, function(x) attr(x, "label", exact = TRUE) %||% NA_character_, character(1) ) %>% @@ -249,7 +241,7 @@ forest_model <- function(model, ) if (inherits(model, "coxph")) { event_detail_tab <- lapply( - seq_len(ncol(model$y)), + setNames(seq_len(ncol(model$y)), colnames(model$y)), function(y_col) { apply( cols, @@ -371,7 +363,7 @@ forest_model <- function(model, variable = coalesce(label, variable) ) if (!is.null(covariates)) { - forest_terms <- filter(forest_terms, term_label %in% covariates) + forest_terms <- filter(forest_terms, term %in% covariates) } if (show_global_p != "none") { @@ -403,9 +395,14 @@ forest_model <- function(model, forest_terms$model_name <- NULL } } else { - forest_terms <- make_forest_terms(model, forest_terms_basic) + forest_terms <- make_forest_terms(model) } + if (is.null(panels)) { + panels <- default_forest_panels(model, factor_separate_line = factor_separate_line) + } + + # #use_exp <- grepl("exp", deparse(trans)) if (!is.null(limits)) { forest_terms <- forest_terms %>% @@ -447,13 +444,13 @@ make_forest_terms_basic <- function(model) { mdl_terms <- stats::terms(model) term_labels <- attr(mdl_terms, "term.labels") mdl_factors <- attr(mdl_terms, "factors") - mdl_data_classes <- attr(mdl_terms, "dataClasses") + mdl_data_classes <- attr(mdl_terms, "dataClasses")[remove_backticks(rownames(mdl_factors))] mdl_data_classes_factors <- mdl_data_classes %in% c("logical", "factor", "character") names(mdl_data_classes_factors) <- names(mdl_data_classes) mdl_terms_inc_factors <- colSums((mdl_factors == 1) & (mdl_data_classes_factors[remove_backticks(rownames(mdl_factors))])) > 0 mdl_terms_are_logical <- colMeans((mdl_factors == 0) | (mdl_data_classes == "logical")) == 1 - tibble::tibble( + forest_terms_basic <- tibble::tibble( term_number = 1:length(term_labels), term_label = term_labels, variable = remove_backticks(term_label), diff --git a/R/model_frame_coxph_simple.R b/R/model_frame_coxph_simple.R new file mode 100644 index 0000000..566cc0d --- /dev/null +++ b/R/model_frame_coxph_simple.R @@ -0,0 +1,21 @@ +model_frame_coxph_simple <- function(model) { + cph_terms <- terms(model) + fcall <- model$call + indx <- match(c("formula", "data", "weights", "subset", + "na.action", "cluster", "id", "istate"), names(fcall), + nomatch = 0) + if (indx[1] == 0) + stop("The coxph call is missing a formula!") + temp <- fcall[c(1, indx)] + temp[[1]] <- quote(stats::model.frame) + temp$formula <- cph_terms + if (!is.null(attr(temp$formula, "specials")$tt)) { + coxenv <- new.env(parent = environment(temp$formula)) + assign("tt", function(x) x, envir = coxenv) + environment(temp$formula) <- coxenv + } + if (is.null(environment(model$terms))) + mf <- eval(temp, parent.frame()) + else mf <- eval(temp, environment(model$terms), parent.frame()) + mf +}