From 7e5aeaa66b6a5a8c2cd34887c25eb050f24d2662 Mon Sep 17 00:00:00 2001 From: Jinseob Kim Date: Wed, 8 May 2024 01:22:01 +0900 Subject: [PATCH] Warning --- R/forestcox.R | 48 +++++++++++++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 13 deletions(-) diff --git a/R/forestcox.R b/R/forestcox.R index 768a916..1bc53ca 100644 --- a/R/forestcox.R +++ b/R/forestcox.R @@ -8,6 +8,7 @@ #' @param decimal.hr Decimal for hazard ratio, Default: 2 #' @param decimal.percent Decimal for percent, Default: 1 #' @param decimal.pvalue Decimal for pvalue, Default: 3 +#' @param cluster Cluster variable for coxph, Default: NULL #' @return Sub-group analysis table. #' @details This result is used to make forestplot. #' @examples @@ -50,7 +51,7 @@ #' @importFrom utils tail -TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, time_eventrate = 3 * 365, decimal.hr = 2, decimal.percent = 1, decimal.pvalue = 3) { +TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, time_eventrate = 3 * 365, decimal.hr = 2, decimal.percent = 1, decimal.pvalue = 3, cluster = NULL) { . <- NULL ### var_subgroup이 categorical variable이 아닌 경우 중단 ### @@ -113,7 +114,11 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, } } else { ### survey data가 아닌 경우 ### - model <- survival::coxph(formula, data = data, x = TRUE) + if (is.null(cluster)){ + model <- survival::coxph(formula, data = data, x = TRUE) + } else { + model <- survival::coxph(formula, data = data, x = TRUE, cluster = get(cluster)) + } # if (!is.null(model$xlevels) & length(model$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") # KM 구하기(categorical인 경우) @@ -236,18 +241,32 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, } } else { ### survey data가 아닌 경우 ### + if (is.null(cluster)){ + model <- data %>% + filter(!is.na(get(var_subgroup))) %>% + split(.[[var_subgroup]]) %>% + purrr::map(~ possible_coxph(formula, data = ., x = T)) + } else{ + model <- data %>% + filter(!is.na(get(var_subgroup))) %>% + split(.[[var_subgroup]]) %>% + purrr::map(~ tryCatch(coxph(formula, data = ., x = T, cluster = get(cluster)), error = function(e)NA)) + } - data %>% - filter(!is.na(get(var_subgroup))) %>% - split(.[[var_subgroup]]) %>% - purrr::map(~ possible_coxph(formula, data = ., x = T)) -> model + data %>% filter(!is.na(get(var_subgroup))) %>% select(dplyr::all_of(var_subgroup)) %>% table() %>% names() -> label_val xlev <- survival::coxph(formula, data = data)$xlevels - model.int <- possible_coxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), data = data) + + if (is.null(cluster)){ + model.int <- possible_coxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), data = data) + } else{ + model.int <- tryCatch(coxph(as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))), data = data, cluster = get(cluster)), error = function(e)NA) + } + # KM 구하기(categorical인 경우만) @@ -273,7 +292,8 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, for (i in rownames(prop)) { if (length(sub_xlev[[i]]) == 1) { - prop[i, paste0(xlabel, "=", sub_xlev[[i]])] <- res.kap.times[["0"]][["surv"]] + #prop[i, paste0(xlabel, "=", sub_xlev[[i]])] <- res.kap.times[["0"]][["surv"]] + prop[i, ] <- res.kap.times[["0"]][["surv"]] } else if (length(sub_xlev[[i]]) > 1) { surv.df <- data.frame(res.kap.times[[i]][c("strata", "surv")]) for (j in colnames(prop)) { @@ -304,8 +324,9 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, } else { model.int$call$formula <- as.formula(gsub(xlabel, paste0(xlabel, "*", var_subgroup), deparse(formula))) model.int$call$data <- as.name("data") - pv_anova <- anova(model.int) - pv_int <- round(pv_anova[nrow(pv_anova), 4], decimal.pvalue) + pv_anova <- tryCatch(anova(model.int), error = function(e)NA) + if (is.na(pv_anova) & !is.null(cluster)) {warning("Warning: Anova test is not available for cluster data. So Interaction P value is NA when 3 and more categorical subgroup variable.")} + pv_int <- tryCatch(round(pv_anova[nrow(pv_anova), 4], decimal.pvalue), error = function(e)NA) } } @@ -416,6 +437,7 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, #' @param decimal.percent Decimal for percent, Default: 1 #' @param decimal.pvalue Decimal for pvalue, Default: 3 #' @param line Include new-line between sub-group variables, Default: F +#' @param cluster Cluster variable for coxph, Default: NULL #' @return Multiple sub-group analysis table. #' @details This result is used to make forestplot. #' @examples @@ -449,17 +471,17 @@ TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, #' @importFrom magrittr %>% #' @importFrom dplyr bind_rows -TableSubgroupMultiCox <- function(formula, var_subgroups = NULL, var_cov = NULL, data, time_eventrate = 3 * 365, decimal.hr = 2, decimal.percent = 1, decimal.pvalue = 3, line = F) { +TableSubgroupMultiCox <- function(formula, var_subgroups = NULL, var_cov = NULL, data, time_eventrate = 3 * 365, decimal.hr = 2, decimal.percent = 1, decimal.pvalue = 3, line = F, cluster = NULL) { . <- NULL xlabel <- setdiff(as.character(formula)[[3]], "+")[1] - out.all <- TableSubgroupCox(formula, var_subgroup = NULL, var_cov = var_cov, data = data, time_eventrate = time_eventrate, decimal.hr = decimal.hr, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue) + out.all <- TableSubgroupCox(formula, var_subgroup = NULL, var_cov = var_cov, data = data, time_eventrate = time_eventrate, decimal.hr = decimal.hr, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue, cluster = cluster) out.all <- dplyr::mutate_all(out.all, as.character) if (is.null(var_subgroups)) { return(out.all) } else { - out.list <- purrr::map(var_subgroups, ~ TableSubgroupCox(formula, var_subgroup = ., var_cov = var_cov, data = data, time_eventrate = time_eventrate, decimal.hr = decimal.hr, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue)) + out.list <- purrr::map(var_subgroups, ~ TableSubgroupCox(formula, var_subgroup = ., var_cov = var_cov, data = data, time_eventrate = time_eventrate, decimal.hr = decimal.hr, decimal.percent = decimal.percent, decimal.pvalue = decimal.pvalue, cluster = cluster)) if (line) { out.newline <- out.list %>% purrr::map(~ rbind(NA, .)) result <- bind_rows(out.all, out.newline %>% dplyr::bind_rows() %>% dplyr::mutate_all(as.character))