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
}

Read our documentation on viewing source code .

Loading