1 ```#' @title TableSubgroupCox: Sub-group analysis table for Cox/svycox model. ``` 2 ```#' @description Sub-group analysis table for Cox/svycox model. ``` 3 ```#' @param formula formula with survival analysis. ``` 4 ```#' @param var_subgroup 1 sub-group variable for analysis, Default: NULL ``` 5 ```#' @param var_cov Variables for additional adjust, Default: NULL ``` 6 ```#' @param data Data or svydesign in survey package. ``` 7 ```#' @param time_eventrate Time for kaplan-meier based event rate calculation, Default = 365 * 3 ``` 8 ```#' @param decimal.hr Decimal for hazard ratio, Default: 2 ``` 9 ```#' @param decimal.percent Decimal for percent, Default: 1 ``` 10 ```#' @param decimal.pvalue Decimal for pvalue, Default: 3 ``` 11 ```#' @return Sub-group analysis table. ``` 12 ```#' @details This result is used to make forestplot. ``` 13 ```#' @examples ``` 14 ```#' library(survival);library(dplyr) ``` 15 ```#' lung %>% ``` 16 ```#' mutate(status = as.integer(status == 1), ``` 17 ```#' sex = factor(sex), ``` 18 ```#' kk = factor(as.integer(pat.karno >= 70))) -> lung ``` 19 ```#' TableSubgroupCox(Surv(time, status) ~ sex, data = lung, time_eventrate = 100) ``` 20 ```#' TableSubgroupCox(Surv(time, status) ~ sex, var_subgroup = "kk", data = lung, ``` 21 ```#' time_eventrate = 100) ``` 22 ```#' ``` 23 ```#' ## survey design ``` 24 ```#' library(survey) ``` 25 ```#' data.design <- svydesign(id = ~1, data = lung) ``` 26 ```#' TableSubgroupCox(Surv(time, status) ~ sex, data = data.design, time_eventrate = 100) ``` 27 ```#' TableSubgroupCox(Surv(time, status) ~ sex, var_subgroup = "kk", data = data.design, ``` 28 ```#' time_eventrate = 100) ``` 29 ```#' @seealso ``` 30 ```#' \code{\link[purrr]{safely}},\code{\link[purrr]{map}},\code{\link[purrr]{map2}} ``` 31 ```#' \code{\link[survival]{coxph}} ``` 32 ```#' \code{\link[survey]{svycoxph}} ``` 33 ```#' \code{\link[stats]{confint}} ``` 34 ```#' @rdname TableSubgroupCox ``` 35 ```#' @export ``` 36 ```#' @importFrom purrr possibly map_dbl map map2 ``` 37 ```#' @importFrom dplyr group_split select filter mutate bind_cols ``` 38 ```#' @importFrom magrittr %>% ``` 39 ```#' @importFrom tibble tibble ``` 40 ```#' @importFrom survival coxph ``` 41 ```#' @importFrom survey svycoxph ``` 42 ```#' @importFrom stats confint coefficients ``` 43 ```#' @importFrom car Anova ``` 44 ```#' @importFrom utils tail ``` 45 46 ```TableSubgroupCox <- function(formula, var_subgroup = NULL, var_cov = NULL, data, time_eventrate = 3 * 365, decimal.hr = 2, decimal.percent = 1, decimal.pvalue = 3){ ``` 47 ``` ``` 48 1 ``` . <- NULL ``` 49 ``` ``` 50 1 ``` if (any(class(data) == "survey.design" & !is.null(var_subgroup))){ ``` 51 0 ``` if (is.numeric(data\$variables[[var_subgroup]])) stop("var_subgroup must categorical.") ``` 52 0 ``` if (length(levels(data\$variables[[as.character(formula[[3]])]])) != 2) stop("Independent variable must have 2 levels.") ``` 53 1 ``` } else if(any(class(data) == "data.frame" & !is.null(var_subgroup))){ ``` 54 0 ``` if (is.numeric(data[[var_subgroup]])) stop("var_subgroup must categorical.") ``` 55 0 ``` if (length(levels(data[[as.character(formula[[3]])]])) != 2) stop("Independent variable must have 2 levels.") ``` 56 ``` } ``` 57 ``` ``` 58 ``` ``` 59 ``` ## functions with error ``` 60 1 ``` possible_table <- purrr::possibly(table, NA) ``` 61 0 ``` possible_prop.table <- purrr::possibly(function(x){prop.table(x, 1)[2, ] * 100}, NA) ``` 62 1 ``` possible_pv <- purrr::possibly(function(x){summary(x)[["coefficients"]][1, ] %>% tail(1)}, NA) ``` 63 1 ``` possible_coxph <- purrr::possibly(survival::coxph, NA) ``` 64 1 ``` possible_svycoxph <- purrr::possibly(survey::svycoxph, NA) ``` 65 1 ``` possible_confint <- purrr::possibly(stats::confint, NA) ``` 66 0 ``` possible_modely <- purrr::possibly(function(x){purrr::map_dbl(x, .[["y"]], 1)}, NA) ``` 67 1 ``` possible_rowone <- purrr::possibly(function(x){x[1, ]}, NA) ``` 68 ``` ``` 69 1 ``` formula.km <- formula ``` 70 1 ``` var_cov <- setdiff(var_cov, c(as.character(formula[[3]]), var_subgroup)) ``` 71 1 ``` if (is.null(var_subgroup)){ ``` 72 1 ``` if (!is.null(var_cov)){ ``` 73 0 ``` formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+"))) ``` 74 ``` } ``` 75 1 ``` if (any(class(data) == "survey.design")){ ``` 76 1 ``` model <- survey::svycoxph(formula, design = data, x= T) ``` 77 0 ``` if (!is.null(model\$xlevels) & length(model\$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") ``` 78 1 ``` res.kap <- survey::svykm(formula.km, design = data) ``` 79 1 ``` prop <- round(100 * sapply(res.kap, function(x){1 - x[["surv"]][which.min(abs(x[["time"]] - time_eventrate))]}), decimal.percent) ``` 80 1 ``` names(prop) <- model\$xlevels[[1]] ``` 81 ``` } else{ ``` 82 1 ``` model <- survival::coxph(formula, data = data, x= TRUE) ``` 83 0 ``` if (!is.null(model\$xlevels) & length(model\$xlevels[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") ``` 84 1 ``` res.kap <- survival::survfit(formula.km, data = data) ``` 85 1 ``` res.kap.times <- summary(res.kap, times = time_eventrate, extend = T) ``` 86 1 ``` prop <- round(100 * (1 - res.kap.times[["surv"]]), decimal.percent) ``` 87 1 ``` names(prop) <- model\$xlevels[[1]] ``` 88 ``` #out.kap <- paste(res.kap.times[["n.event"]], " (", round(100 * (1 - res.kap.times[["surv"]]), decimal.percent), ")", sep = "") ``` 89 ``` } ``` 90 ``` ``` 91 1 ``` Point.Estimate <- round(exp(coef(model)), decimal.hr)[1] ``` 92 ``` ``` 93 ``` ``` 94 ``` #if (length(Point.Estimate) > 1){ ``` 95 ``` # stop("Formula must contain 1 independent variable only.") ``` 96 ``` #} ``` 97 ``` ``` 98 ``` ``` 99 1 ``` CI <- round(exp(confint(model)[1, ]), decimal.hr) ``` 100 1 ``` event <- purrr::map_dbl(model\$y, 1) %>% tail(model\$n) ``` 101 ``` #prop <- round(prop.table(table(event, model\$x[, 1]), 2)[2, ] * 100, decimal.percent) ``` 102 1 ``` pv <- round(summary(model)\$coefficients[1, "Pr(>|z|)"], decimal.pvalue) ``` 103 ``` ``` 104 1 ``` tibble::tibble(Variable = "Overall", Count = model\$n, Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2]) %>% ``` 105 1 ``` cbind(t(prop)) %>% ``` 106 1 ``` mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out ``` 107 ``` ``` 108 1 ``` return(out) ``` 109 1 ``` } else if (length(var_subgroup) > 1 | any(grepl(var_subgroup, formula))){ ``` 110 0 ``` stop("Please input correct subgroup variable.") ``` 111 ``` } else{ ``` 112 1 ``` if (!is.null(var_cov)){ ``` 113 0 ``` formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+"))) ``` 114 ``` } ``` 115 1 ``` if (any(class(data) == "survey.design")){ ``` 116 1 ``` data\$variables[[var_subgroup]] %>% table %>% names -> label_val ``` 117 1 ``` label_val %>% purrr::map(~possible_svycoxph(formula, design = subset(data, get(var_subgroup) == .), x = TRUE)) -> model ``` 118 1 ``` xlev <- survey::svycoxph(formula, design = data)\$xlevels ``` 119 1 ``` xlabel <- names(attr(model[[which(!is.na(model))[1]]]\$x, "contrast"))[1] ``` 120 1 ``` pvs_int <- possible_svycoxph(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data) %>% summary %>% coefficients ``` 121 1 ``` pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) ``` 122 0 ``` if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") ``` 123 ``` ``` 124 ``` ``` 125 1 ``` if (length(label_val) > 2){ ``` 126 0 ``` data.int <- data\$variables ``` 127 0 ``` model.int <- survival::coxph(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), data = data.int, weights = get(names(data\$allprob)), robust =T) ``` 128 0 ``` model.int\$call\$formula <- as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))) ``` 129 0 ``` model.int\$call\$data <- as.name("data.int") ``` 130 0 ``` pv_anova <- car::Anova(model.int, test.statistics = "Wald") ``` 131 0 ``` pv_int <- pv_anova[nrow(pv_anova), 3] ``` 132 ``` } ``` 133 1 ``` res.kap <- purrr::map(label_val, ~survey::svykm(formula.km, design = subset(data, get(var_subgroup) == . ))) ``` 134 1 ``` mkz <- function(reskap){ ``` 135 1 ``` round(100 * sapply(reskap, function(x){1 - x[["surv"]][which.min(abs(x[["time"]] - time_eventrate))]}), decimal.percent) ``` 136 ``` } ``` 137 1 ``` prop <- purrr::map(res.kap, mkz) %>% dplyr::bind_cols() %>% t ``` 138 ``` #prop <- purrr::map(res.kap, ~round(100 * sapply(., function(x){1 - x[["surv"]][which.min(abs(x[["time"]] - time_eventrate))]}), decimal.percent)) ``` 139 1 ``` colnames(prop) <- xlev[[1]] ``` 140 ``` ``` 141 ``` } else{ ``` 142 1 ``` data %>% filter(!is.na(get(var_subgroup))) %>% group_split(get(var_subgroup)) %>% purrr::map(~possible_coxph(formula, data = ., x= T)) -> model ``` 143 1 ``` data %>% filter(!is.na(get(var_subgroup))) %>% select(var_subgroup) %>% table %>% names -> label_val ``` 144 1 ``` xlev <- survival::coxph(formula, data = data)\$xlevels ``` 145 1 ``` xlabel <- names(attr(model[[which(!is.na(model))[1]]]\$x, "contrast"))[1] ``` 146 1 ``` model.int <- possible_coxph(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), data = data) ``` 147 1 ``` pvs_int <- model.int %>% summary %>% coefficients ``` 148 1 ``` pv_int <- round(pvs_int[nrow(pvs_int), ncol(pvs_int)], decimal.pvalue) ``` 149 0 ``` if (!is.null(xlev) & length(xlev[[1]]) != 2) stop("Categorical independent variable must have 2 levels.") ``` 150 ``` ``` 151 1 ``` if (length(label_val) > 2){ ``` 152 0 ``` model.int\$call\$formula <- as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))) ``` 153 0 ``` model.int\$call\$data <- as.name("data") ``` 154 0 ``` pv_anova <- car::Anova(model.int, test.statistics = "Wald") ``` 155 0 ``` pv_int <- round(pv_anova[nrow(pv_anova), 3], decimal.pvalue) ``` 156 ``` } ``` 157 1 ``` res.kap.times <- data %>% filter(!is.na(get(var_subgroup))) %>% group_split(get(var_subgroup)) %>% purrr::map(~survival::survfit(formula.km, data = .)) %>% purrr::map(~summary(., times = time_eventrate, extend = T)) ``` 158 1 ``` prop <- res.kap.times %>% purrr::map(~round(100 * (1 - .[["surv"]]), decimal.percent)) %>% dplyr::bind_cols() %>% t ``` 159 1 ``` colnames(prop) <- xlev[[1]] ``` 160 ``` } ``` 161 ``` ``` 162 1 ``` model %>% purrr::map_dbl("n", .default = NA) -> Count ``` 163 1 ``` model %>% purrr::map("coefficients", .default = NA) %>% purrr::map_dbl(1) %>% exp %>% round(decimal.hr) -> Point.Estimate ``` 164 1 ``` model %>% purrr::map(possible_confint) %>% purrr::map(possible_rowone) %>% Reduce(rbind, .) %>% exp %>% round(decimal.hr) -> CI ``` 165 ``` #model %>% purrr::map("y") %>% purrr::map(~purrr::map_dbl(., 1)) %>% purrr::map(~tail(., length(.)/2)) -> event ``` 166 ``` #purrr::map2(event, model, ~possible_table(.x, .y[["x"]][, 1])) %>% purrr::map(possible_prop.table) %>% purrr::map(~round(., decimal.percent)) %>% Reduce(rbind, .) -> prop ``` 167 1 ``` model %>% purrr::map(possible_pv) %>% purrr::map_dbl(~round(., decimal.pvalue)) -> pv ``` 168 ``` ``` 169 ``` ``` 170 1 ``` tibble::tibble(Variable = paste(" ", label_val) , Count = Count, Percent = round(Count/sum(Count) * 100, decimal.percent), `Point Estimate` = Point.Estimate, Lower = CI[, 1], Upper = CI[, 2]) %>% ``` 171 1 ``` cbind(prop) %>% ``` 172 1 ``` mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out ``` 173 ``` ``` 174 1 ``` return(rbind(c(var_subgroup, rep(NA, ncol(out) - 2), ifelse(pv_int >= 0.001, pv_int, "<0.001")), out)) ``` 175 ``` } ``` 176 ```} ``` 177 178 179 180 181 ```#' @title TableSubgroupMultiCox: Multiple sub-group analysis table for Cox/svycox model. ``` 182 ```#' @description Multiple sub-group analysis table for Cox/svycox model. ``` 183 ```#' @param formula formula with survival analysis. ``` 184 ```#' @param var_subgroups Multiple sub-group variables for analysis, Default: NULL ``` 185 ```#' @param var_cov Variables for additional adjust, Default: NULL ``` 186 ```#' @param data Data or svydesign in survey package. ``` 187 ```#' @param time_eventrate Time for kaplan-meier based event rate calculation, Default = 365 * 3 ``` 188 ```#' @param decimal.hr Decimal for hazard ratio, Default: 2 ``` 189 ```#' @param decimal.percent Decimal for percent, Default: 1 ``` 190 ```#' @param decimal.pvalue Decimal for pvalue, Default: 3 ``` 191 ```#' @param line Include new-line between sub-group variables, Default: F ``` 192 ```#' @return Multiple sub-group analysis table. ``` 193 ```#' @details This result is used to make forestplot. ``` 194 ```#' @examples ``` 195 ```#' library(survival);library(dplyr) ``` 196 ```#' lung %>% ``` 197 ```#' mutate(status = as.integer(status == 1), ``` 198 ```#' sex = factor(sex), ``` 199 ```#' kk = factor(as.integer(pat.karno >= 70)), ``` 200 ```#' kk1 = factor(as.integer(pat.karno >= 60))) -> lung ``` 201 ```#' TableSubgroupMultiCox(Surv(time, status) ~ sex, var_subgroups = c("kk", "kk1"), ``` 202 ```#' data=lung, time_eventrate = 100, line = TRUE) ``` 203 ```#' ``` 204 ```#' ## survey design ``` 205 ```#' library(survey) ``` 206 ```#' data.design <- svydesign(id = ~1, data = lung) ``` 207 ```#' TableSubgroupMultiCox(Surv(time, status) ~ sex, var_subgroups = c("kk", "kk1"), ``` 208 ```#' data = data.design, time_eventrate = 100) ``` 209 ```#' @seealso ``` 210 ```#' \code{\link[purrr]{map}} ``` 211 ```#' \code{\link[dplyr]{bind}} ``` 212 ```#' @rdname TableSubgroupMultiCox ``` 213 ```#' @export ``` 214 ```#' @importFrom purrr map ``` 215 ```#' @importFrom magrittr %>% ``` 216 ```#' @importFrom dplyr bind_rows ``` 217 218 219 ```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){ ``` 220 ``` ``` 221 1 ``` . <- NULL ``` 222 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) ``` 223 ``` ``` 224 1 ``` if (is.null(var_subgroups)){ ``` 225 1 ``` return(out.all) ``` 226 ``` } else { ``` 227 1 ``` 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)) ``` 228 1 ``` if (line){ ``` 229 1 ``` out.newline <- out.list %>% purrr::map(~rbind(NA, .)) ``` 230 1 ``` return(rbind(out.all, out.newline %>% dplyr::bind_rows())) ``` 231 ``` } else{ ``` 232 1 ``` return(rbind(out.all, out.list %>% dplyr::bind_rows())) ``` 233 ``` } ``` 234 ``` } ``` 235 ``` } ``` 236 237 238

Read our documentation on viewing source code .