|
1
|
|
#' @title cox2.display: table for coxph.object with model option: TRUE - allow "frailty" or "cluster" model
|
|
2
|
|
#' @description Table for coxph.object with model option: TRUE - allow "frailty" or "cluster" model
|
|
3
|
|
#' @param cox.obj.withmodel coxph.object with model option: TRUE
|
|
4
|
|
#' @param dec Decimal point, Default: 2
|
|
5
|
|
#' @return Table, cluster/frailty info, metrics, caption
|
|
6
|
|
#' @details GEE like - cluster, Mixed effect model like - frailty
|
|
7
|
|
#' @examples
|
|
8
|
|
#' library(survival);data(lung)
|
|
9
|
|
#' fit1 <- coxph(Surv(time, status) ~ ph.ecog + age + cluster(inst), data = lung, model = TRUE)
|
|
10
|
|
#' fit2 <- coxph(Surv(time, status) ~ ph.ecog + age + frailty(inst), data = lung, model = TRUE)
|
|
11
|
|
#' cox2.display(fit1)
|
|
12
|
|
#' cox2.display(fit2)
|
|
13
|
|
#' @rdname cox2.display
|
|
14
|
|
#' @export
|
|
15
|
|
#' @importFrom survival coxph cluster frailty
|
|
16
|
|
#' @importFrom stats formula update
|
|
17
|
|
|
|
18
|
|
cox2.display <- function (cox.obj.withmodel, dec = 2)
|
|
19
|
|
{
|
|
20
|
1
|
model <- cox.obj.withmodel
|
|
21
|
0
|
if(!any(class(model)=="coxph")){stop("Model not from Cox model")}
|
|
22
|
|
|
|
23
|
1
|
xf <- attr(model$terms, "term.labels") # Independent vars
|
|
24
|
1
|
xf.old <- xf
|
|
25
|
1
|
xc <- NULL
|
|
26
|
1
|
xc.vn <- NULL
|
|
27
|
1
|
mtype <- "normal"
|
|
28
|
|
|
|
29
|
1
|
if(length(grep("strata", xf)) > 0){
|
|
30
|
0
|
xf <- xf[-grep("strata",xf)]
|
|
31
|
1
|
} else if(length(grep("frailty\\(", xf)) > 0){
|
|
32
|
1
|
xf <- xf[-grep("frailty\\(",xf)]
|
|
33
|
1
|
mtype <- "frailty"
|
|
34
|
1
|
xc <- setdiff(xf.old, xf)
|
|
35
|
1
|
xc.vn <- strsplit(strsplit(xc, "frailty\\(")[[1]][2], "\\)")[[1]][1]
|
|
36
|
1
|
} else if(!(is.null(model$call$cluster))){
|
|
37
|
1
|
mtype <- "cluster"
|
|
38
|
|
#xfull <- strsplit(as.character(model$call[[2]][3]), " \\+ ")[[1]]
|
|
39
|
|
#xc <- setdiff(xfull, xf)
|
|
40
|
|
#xf <- ifelse(length(grep("cluster\\(",xf)) > 0, xf[-grep("cluster\\(",xf)], xf)
|
|
41
|
|
#xc <- setdiff(xf.old, xf)
|
|
42
|
1
|
xc <- as.character(model$call$cluster)
|
|
43
|
1
|
xc.vn <- xc
|
|
44
|
|
#xc.vn <- strsplit(strsplit(xc, "cluster\\(")[[1]][2], "\\)")[[1]][1]
|
|
45
|
|
|
|
46
|
|
}
|
|
47
|
|
|
|
48
|
1
|
formula.surv <- as.character(model$formula)[2]
|
|
49
|
1
|
formula.ranef <- paste(" + ", xc, sep = "")
|
|
50
|
1
|
mdata <- model$model
|
|
51
|
|
|
|
52
|
1
|
if (length(xc) == 0){
|
|
53
|
1
|
formula.ranef <- NULL
|
|
54
|
|
} else{
|
|
55
|
1
|
names(mdata)[names(mdata) == xc] <- xc.vn
|
|
56
|
|
}
|
|
57
|
|
|
|
58
|
|
#if (is.null(data)){
|
|
59
|
|
# mdata = data.frame(get(as.character(model$call)[3]))
|
|
60
|
|
#} else{
|
|
61
|
|
# mdata = data.frame(data)
|
|
62
|
|
#}
|
|
63
|
|
|
|
64
|
|
|
|
65
|
|
|
|
66
|
1
|
if(length(xf) == 1){
|
|
67
|
0
|
uni.res <- data.frame(summary(model)$coefficients)
|
|
68
|
|
#uni.res <- data.frame(summary(coxph(as.formula(paste("mdata[, 1]", "~", xf, formula.ranef, sep="")), data = mdata))$coefficients)
|
|
69
|
0
|
rn.uni <- lapply(list(uni.res), rownames)
|
|
70
|
0
|
names(uni.res)[ncol(uni.res)] <- "p"
|
|
71
|
0
|
uni.res2 <- NULL
|
|
72
|
0
|
if (mtype == "normal"){
|
|
73
|
0
|
uni.res2 <- uni.res[, c(1, 3, 4, 5)]
|
|
74
|
0
|
if (length(grep("robust.se", names(uni.res))) > 0){
|
|
75
|
0
|
uni.res2 <- uni.res[, c(1, 4, 5, 6)]
|
|
76
|
|
}
|
|
77
|
0
|
} else if (mtype == "cluster"){
|
|
78
|
0
|
uni.res2 <- uni.res[, c(1, 4, 5, 6)]
|
|
79
|
|
} else {
|
|
80
|
0
|
uni.res2 <- uni.res[-nrow(uni.res), c(1, 3, 4, 6)]
|
|
81
|
|
}
|
|
82
|
0
|
fix.all <- coxExp(uni.res2, dec = dec)
|
|
83
|
0
|
colnames(fix.all) <- c("HR(95%CI)", "P value")
|
|
84
|
|
#rownames(fix.all) = ifelse(mtype == "frailty", names(model$coefficients)[-length(model$coefficients)], names(model$coefficients))
|
|
85
|
0
|
rownames(fix.all) <- names(model$coefficients)
|
|
86
|
|
} else{
|
|
87
|
1
|
mdata2 <- cbind(matrix(sapply(mdata[, 1], `[[`, 1), ncol = 2), mdata[, -1])
|
|
88
|
1
|
names(mdata2)[1:2] <- as.character(model$formula[[2]][2:3])
|
|
89
|
1
|
if (!is.null(xc.vn)){
|
|
90
|
1
|
names(mdata2)[ncol(mdata2)] <- xc.vn
|
|
91
|
|
}
|
|
92
|
1
|
basemodel <- update(model, formula(paste(c(". ~ .", xf), collapse=' - ')), data = mdata2)
|
|
93
|
1
|
unis <- lapply(xf, function(x){
|
|
94
|
1
|
newfit <- update(basemodel, formula(paste0(". ~ . +", x)), data= mdata2)
|
|
95
|
1
|
uni.res <- data.frame(summary(newfit)$coefficients)
|
|
96
|
1
|
uni.res <- uni.res[grep(x, rownames(uni.res)), ]
|
|
97
|
|
#uni.res <- uni.res[c(2:nrow(uni.res), 1), ]
|
|
98
|
|
#uni.res <- data.frame(summary(coxph(as.formula(paste("mdata[, 1]", "~", x, formula.ranef, sep="")), data = mdata))$coefficients)
|
|
99
|
1
|
names(uni.res)[ncol(uni.res)] <- "p"
|
|
100
|
1
|
uni.res2 <- NULL
|
|
101
|
1
|
if (mtype == "normal"){
|
|
102
|
1
|
uni.res2 <- uni.res[, c(1, 3, 4, 5)]
|
|
103
|
1
|
} else if (mtype == "cluster"){
|
|
104
|
1
|
uni.res2 <- uni.res[, c(1, 4, 5, 6)]
|
|
105
|
|
} else {
|
|
106
|
1
|
uni.res2 <- uni.res[, c(1, 3, 4, 6)]
|
|
107
|
|
}
|
|
108
|
1
|
return(uni.res2)
|
|
109
|
|
})
|
|
110
|
1
|
rn.uni <- lapply(unis, rownames)
|
|
111
|
1
|
unis2 <- Reduce(rbind, unis)
|
|
112
|
1
|
uni.res <- unis2
|
|
113
|
1
|
mul.res <- data.frame(coefNA(model))
|
|
114
|
1
|
uni.res <- uni.res[rownames(uni.res) %in% rownames(mul.res), ]
|
|
115
|
1
|
colnames(mul.res)[ncol(mul.res)] <- "p"
|
|
116
|
1
|
fix.all <- cbind(coxExp(uni.res, dec = dec), coxExp(mul.res[rownames(uni.res), names(uni.res)], dec = dec))
|
|
117
|
1
|
colnames(fix.all) <- c("crude HR(95%CI)", "crude P value", "adj. HR(95%CI)", "adj. P value")
|
|
118
|
1
|
rownames(fix.all) <- rownames(uni.res)
|
|
119
|
|
}
|
|
120
|
|
|
|
121
|
|
## rownames
|
|
122
|
1
|
fix.all.list <- lapply(1:length(xf), function(x){fix.all[rownames(fix.all) %in% rn.uni[[x]], ]})
|
|
123
|
1
|
varnum.mfac <- which(lapply(fix.all.list, length) > ncol(fix.all))
|
|
124
|
0
|
lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])})
|
|
125
|
1
|
fix.all.unlist <- Reduce(rbind, fix.all.list)
|
|
126
|
|
|
|
127
|
1
|
rn.list <- lapply(1:length(xf), function(x){rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]})
|
|
128
|
1
|
varnum.2fac <- which(lapply(xf, function(x){length(sapply(mdata, levels)[[x]])}) == 2)
|
|
129
|
0
|
lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xf[x], ": ", levels(mdata[, xf[x]])[2], " vs ", levels(mdata[, xf[x]])[1], sep="")})
|
|
130
|
0
|
lapply(varnum.mfac, function(x){rn.list[[x]] <<- c(paste(xf[x],": ref.=", levels(mdata[, xf[x]])[1], sep=""), gsub(xf[x]," ", rn.list[[x]]))})
|
|
131
|
1
|
if (class(fix.all.unlist)[1] == "character"){
|
|
132
|
0
|
fix.all.unlist <- t(data.frame(fix.all.unlist))
|
|
133
|
|
}
|
|
134
|
1
|
rownames(fix.all.unlist) <- unlist(rn.list)
|
|
135
|
1
|
pv.colnum <- which(colnames(fix.all.unlist) %in% c("P value", "crude P value", "adj. P value"))
|
|
136
|
1
|
for (i in pv.colnum){
|
|
137
|
1
|
fix.all.unlist[, i] <- ifelse(as.numeric(fix.all.unlist[, i]) < 0.001, "< 0.001", round(as.numeric(fix.all.unlist[, i]), dec + 1))
|
|
138
|
|
}
|
|
139
|
|
|
|
140
|
|
|
|
141
|
|
## random effect
|
|
142
|
|
#ranef = unlist(model$vcoef)
|
|
143
|
|
#ranef.out = round(ranef, dec)
|
|
144
|
1
|
ranef.mat = NULL
|
|
145
|
1
|
if (mtype == "cluster"){
|
|
146
|
1
|
ranef.mat <- cbind(c(NA, NA), matrix(NA, length(xc) + 1, ncol(fix.all) - 1))
|
|
147
|
|
#clname = strsplit(xc, "\\(")[[1]]
|
|
148
|
1
|
clname <- c("cluster", xc)
|
|
149
|
1
|
cvname <- strsplit(paste(clname[-1], collapse = "("), "\\)")[[1]]
|
|
150
|
1
|
cvname <- paste(cvname[length(cvname)], collapse = ")")
|
|
151
|
1
|
rownames(ranef.mat) <- c(clname[1], cvname)
|
|
152
|
1
|
} else if (mtype == "frailty"){
|
|
153
|
1
|
ranef.mat <- cbind(c(NA, NA), matrix(NA, length(xc) + 1, ncol(fix.all) - 1))
|
|
154
|
1
|
clname <- strsplit(xc, "\\(")[[1]]
|
|
155
|
1
|
cvname <- strsplit(paste(clname[-1], collapse = "("), "\\)")[[1]]
|
|
156
|
1
|
cvname <- paste(cvname[length(cvname)], collapse = ")")
|
|
157
|
1
|
rownames(ranef.mat) <- c(clname[1], cvname)
|
|
158
|
|
}
|
|
159
|
|
|
|
160
|
|
|
|
161
|
|
## metric
|
|
162
|
|
#no.grp = unlist(lapply(model$frail, length))
|
|
163
|
1
|
no.obs = model$n
|
|
164
|
1
|
no.event = model$nevent
|
|
165
|
1
|
metric.mat = cbind(c(NA, no.obs, no.event), matrix(NA, 3, ncol(fix.all) - 1))
|
|
166
|
1
|
rownames(metric.mat) = c(NA, "No. of observations", "No. of events")
|
|
167
|
|
|
|
168
|
|
## Integrated ll
|
|
169
|
|
#ll = model$loglik[2]
|
|
170
|
|
#aic = -2 * ll -2*model$df[1]
|
|
171
|
|
|
|
172
|
|
## caption
|
|
173
|
1
|
surv.string <- as.character(attr(model$terms, "variables")[[2]])
|
|
174
|
1
|
time.var.name <- surv.string[2]
|
|
175
|
1
|
status.var.name <- surv.string[3]
|
|
176
|
1
|
intro <- paste("Cox model on time ('", time.var.name, "') to event ('", status.var.name, "')", sep="")
|
|
177
|
1
|
if (mtype == "cluster"){
|
|
178
|
1
|
intro <- paste("Marginal", intro, "- Group", cvname)
|
|
179
|
1
|
} else if (mtype == "frailty"){
|
|
180
|
1
|
intro <- paste("Frailty", intro, "- Group", cvname)
|
|
181
|
|
}
|
|
182
|
|
|
|
183
|
1
|
var.names0 <- attr(model$terms, "term.labels")
|
|
184
|
0
|
if(length(grep("strata", var.names0))>0) {intro <- paste(intro, " with '", var.names0[grep("strata", var.names0)], "'", sep="" )}
|
|
185
|
|
|
|
186
|
1
|
if (is.null(ranef.mat)){
|
|
187
|
1
|
return(list(table = fix.all.unlist, metric = metric.mat, caption = intro))
|
|
188
|
|
} else{
|
|
189
|
1
|
return(list(table = fix.all.unlist, ranef = ranef.mat, metric = metric.mat, caption = intro))
|
|
190
|
|
}
|
|
191
|
|
}
|
jinseob2kim