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 .

Loading