|
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
|
2
|
. <- NULL
|
|
49
|
|
|
|
50
|
2
|
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
|
2
|
} 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
|
2
|
possible_table <- purrr::possibly(table, NA)
|
|
61
|
0
|
possible_prop.table <- purrr::possibly(function(x){prop.table(x, 1)[2, ] * 100}, NA)
|
|
62
|
2
|
possible_pv <- purrr::possibly(function(x){summary(x)[["coefficients"]][1, ] %>% tail(1)}, NA)
|
|
63
|
2
|
possible_coxph <- purrr::possibly(survival::coxph, NA)
|
|
64
|
2
|
possible_svycoxph <- purrr::possibly(survey::svycoxph, NA)
|
|
65
|
2
|
possible_confint <- purrr::possibly(stats::confint, NA)
|
|
66
|
0
|
possible_modely <- purrr::possibly(function(x){purrr::map_dbl(x, .[["y"]], 1)}, NA)
|
|
67
|
2
|
possible_rowone <- purrr::possibly(function(x){x[1, ]}, NA)
|
|
68
|
|
|
|
69
|
2
|
formula.km <- formula
|
|
70
|
2
|
var_cov <- setdiff(var_cov, c(as.character(formula[[3]]), var_subgroup))
|
|
71
|
2
|
if (is.null(var_subgroup)){
|
|
72
|
2
|
if (!is.null(var_cov)){
|
|
73
|
0
|
formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+")))
|
|
74
|
|
}
|
|
75
|
2
|
if (any(class(data) == "survey.design")){
|
|
76
|
2
|
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
|
2
|
res.kap <- survey::svykm(formula.km, design = data)
|
|
79
|
2
|
prop <- round(100 * sapply(res.kap, function(x){1 - x[["surv"]][which.min(abs(x[["time"]] - time_eventrate))]}), decimal.percent)
|
|
80
|
2
|
names(prop) <- model$xlevels[[1]]
|
|
81
|
|
} else{
|
|
82
|
2
|
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
|
2
|
res.kap <- survival::survfit(formula.km, data = data)
|
|
85
|
2
|
res.kap.times <- summary(res.kap, times = time_eventrate, extend = T)
|
|
86
|
2
|
prop <- round(100 * (1 - res.kap.times[["surv"]]), decimal.percent)
|
|
87
|
2
|
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
|
2
|
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
|
2
|
CI <- round(exp(confint(model)[1, ]), decimal.hr)
|
|
100
|
2
|
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
|
2
|
pv <- round(summary(model)$coefficients[1, "Pr(>|z|)"], decimal.pvalue)
|
|
103
|
|
|
|
104
|
2
|
tibble::tibble(Variable = "Overall", Count = model$n, Percent = 100, `Point Estimate` = Point.Estimate, Lower = CI[1], Upper = CI[2]) %>%
|
|
105
|
2
|
cbind(t(prop)) %>%
|
|
106
|
2
|
mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out
|
|
107
|
|
|
|
108
|
2
|
return(out)
|
|
109
|
2
|
} else if (length(var_subgroup) > 1 | any(grepl(var_subgroup, formula))){
|
|
110
|
0
|
stop("Please input correct subgroup variable.")
|
|
111
|
|
} else{
|
|
112
|
2
|
if (!is.null(var_cov)){
|
|
113
|
0
|
formula <- as.formula(paste0(deparse(formula), " + ", paste(var_cov, collapse = "+")))
|
|
114
|
|
}
|
|
115
|
2
|
if (any(class(data) == "survey.design")){
|
|
116
|
2
|
data$variables[[var_subgroup]] %>% table %>% names -> label_val
|
|
117
|
2
|
label_val %>% purrr::map(~possible_svycoxph(formula, design = subset(data, get(var_subgroup) == .), x = TRUE)) -> model
|
|
118
|
2
|
xlev <- survey::svycoxph(formula, design = data)$xlevels
|
|
119
|
2
|
xlabel <- names(attr(model[[which(!is.na(model))[1]]]$x, "contrast"))[1]
|
|
120
|
2
|
pvs_int <- possible_svycoxph(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), design = data) %>% summary %>% coefficients
|
|
121
|
2
|
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
|
2
|
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
|
2
|
res.kap <- purrr::map(label_val, ~survey::svykm(formula.km, design = subset(data, get(var_subgroup) == . )))
|
|
134
|
2
|
mkz <- function(reskap){
|
|
135
|
2
|
round(100 * sapply(reskap, function(x){1 - x[["surv"]][which.min(abs(x[["time"]] - time_eventrate))]}), decimal.percent)
|
|
136
|
|
}
|
|
137
|
2
|
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
|
2
|
colnames(prop) <- xlev[[1]]
|
|
140
|
|
|
|
141
|
|
} else{
|
|
142
|
2
|
data %>% filter(!is.na(get(var_subgroup))) %>% group_split(get(var_subgroup)) %>% purrr::map(~possible_coxph(formula, data = ., x= T)) -> model
|
|
143
|
2
|
data %>% filter(!is.na(get(var_subgroup))) %>% select(var_subgroup) %>% table %>% names -> label_val
|
|
144
|
2
|
xlev <- survival::coxph(formula, data = data)$xlevels
|
|
145
|
2
|
xlabel <- names(attr(model[[which(!is.na(model))[1]]]$x, "contrast"))[1]
|
|
146
|
2
|
model.int <- possible_coxph(as.formula(gsub(xlabel, paste(xlabel, "*", var_subgroup, sep=""), deparse(formula))), data = data)
|
|
147
|
2
|
pvs_int <- model.int %>% summary %>% coefficients
|
|
148
|
2
|
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
|
2
|
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
|
2
|
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
|
2
|
prop <- res.kap.times %>% purrr::map(~round(100 * (1 - .[["surv"]]), decimal.percent)) %>% dplyr::bind_cols() %>% t
|
|
159
|
2
|
colnames(prop) <- xlev[[1]]
|
|
160
|
|
}
|
|
161
|
|
|
|
162
|
2
|
model %>% purrr::map_dbl("n", .default = NA) -> Count
|
|
163
|
2
|
model %>% purrr::map("coefficients", .default = NA) %>% purrr::map_dbl(1) %>% exp %>% round(decimal.hr) -> Point.Estimate
|
|
164
|
2
|
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
|
2
|
model %>% purrr::map(possible_pv) %>% purrr::map_dbl(~round(., decimal.pvalue)) -> pv
|
|
168
|
|
|
|
169
|
|
|
|
170
|
2
|
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
|
2
|
cbind(prop) %>%
|
|
172
|
2
|
mutate(`P value` = ifelse(pv >= 0.001, pv, "<0.001"), `P for interaction` = NA) -> out
|
|
173
|
|
|
|
174
|
2
|
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
|
2
|
. <- NULL
|
|
222
|
2
|
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
|
2
|
if (is.null(var_subgroups)){
|
|
225
|
2
|
return(out.all)
|
|
226
|
|
} else {
|
|
227
|
2
|
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
|
2
|
if (line){
|
|
229
|
2
|
out.newline <- out.list %>% purrr::map(~rbind(NA, .))
|
|
230
|
2
|
return(rbind(out.all, out.newline %>% dplyr::bind_rows()))
|
|
231
|
|
} else{
|
|
232
|
2
|
return(rbind(out.all, out.list %>% dplyr::bind_rows()))
|
|
233
|
|
}
|
|
234
|
|
}
|
|
235
|
|
}
|
|
236
|
|
|
|
237
|
|
|
|
238
|
|
|
jinseob2kim