1

2
#' @title svycoxph.display: table for svycoxph.object in survey package.
3
#' @description Table for complex design cox model.
4
#' @param svycoxph.obj svycoxph.object
5
#' @param decimal digit, Default: 2
6
#' @return List including table, metric, caption
7
#' @details DETAILS
8
#' @examples 
9
#'  library(survival);data(pbc)
10
#'  pbc$sex = factor(pbc$sex)
11
#'  pbc$stage = factor(pbc$stage)
12
#'  pbc$randomized<-with(pbc, !is.na(trt) & trt>0)
13
#'  biasmodel<-glm(randomized~age*edema,data=pbc,family=binomial)
14
#'  pbc$randprob<-fitted(biasmodel)
15
#'  
16
#'  if (is.null(pbc$albumin)) pbc$albumin<-pbc$alb ##pre2.9.0
17
#'  
18
#'  dpbc <- survey::svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
19
#'  
20
#'  model <- survey::svycoxph(Surv(time,status>0)~ sex + protime + albumin + stage,design=dpbc)
21
#'  svycox.display(model)
22
#' @seealso 
23
#'  \code{\link[survey]{svycoxph}}
24
#'  \code{\link[stats]{AIC}}
25
#' @rdname svycox.display
26
#' @export 
27
#' @importFrom survey svycoxph
28
#' @importFrom stats AIC
29
 
30

31
svycox.display <- function(svycoxph.obj, decimal = 2){
32
  
33 1
  model <- svycoxph.obj
34 0
  if(!any(class(model)=="svycoxph")){stop("Model not from Survey cox model")}
35 1
  xf <- attr(model$terms, "term.labels")
36 1
  formula.surv <- as.character(model$formula)[2]
37 1
  design.model <- model$survey.design
38
  
39 1
  if(length(xf) == 1){
40 0
    uni.res <- data.frame(summary(model)$coefficients)
41
    #uni.res <- data.frame(summary(survey::svycoxph(as.formula(paste(formula.surv, "~", xf, sep="")), design = design.model))$coefficients)
42 0
    names(uni.res)[ncol(uni.res)] <- "p"
43 0
    uni.res2 <- uni.res[, c("coef", grep("se", colnames(uni.res), value = T)[length(grep("se", colnames(uni.res)))], "z", "p")]
44
    
45 0
    fix.all = coxExp(uni.res2, dec = decimal)
46 0
    colnames(fix.all) <- c("HR(95%CI)", "P value")
47
    #rownames(fix.all) = ifelse(mtype == "frailty", names(model$coefficients)[-length(model$coefficients)], names(model$coefficients))
48 0
    rownames(fix.all) <-  names(model$coefficients)
49
    } else {
50 1
    unis <- lapply(xf, function(x){
51 1
    uni.res <- data.frame(summary(survey::svycoxph(as.formula(paste(formula.surv, "~", x, sep="")), design = design.model))$coefficients)
52 1
    names(uni.res)[ncol(uni.res)] <- "p"
53 1
    uni.res2 <- uni.res[, c("coef", grep("se", colnames(uni.res), value = T)[length(grep("se", colnames(uni.res)))], "z", "p")]
54 1
    return(uni.res2)
55
    })
56
    
57 1
   unis2 <- Reduce(rbind, unis)
58 1
   uni.res <- unis2
59 1
    mul.res <- data.frame(summary(model)$coefficients)
60 1
    uni.res <- uni.res[rownames(uni.res) %in% rownames(mul.res), ]  ## set
61 1
    colnames(mul.res)[ncol(mul.res)] <- "p"
62 1
    fix.all <- cbind(coxExp(uni.res, dec = decimal), coxExp(mul.res[rownames(uni.res), names(uni.res)], dec = decimal))
63 1
    colnames(fix.all) <- c("crude HR(95%CI)", "crude P value", "adj. HR(95%CI)", "adj. P value")
64 1
    rownames(fix.all) <- rownames(uni.res)
65
    }
66
  
67
  ## rownames
68 1
  fix.all.list <- lapply(xf, function(x){fix.all[grepl(x, rownames(fix.all)),]})      
69 1
  varnum.mfac <- which(lapply(fix.all.list, length) > ncol(fix.all))
70 1
  lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])})
71 1
  fix.all.unlist <- Reduce(rbind, fix.all.list)
72
  
73 1
  rn.list <- lapply(xf, function(x){rownames(fix.all)[grepl(x, rownames(fix.all))]})
74 1
  varnum.2fac <- which(xf %in% names(model$xlevels)[lapply(model$xlevels, length) == 2])
75 1
  lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xf[x], ": ", model$xlevels[[xf[x]]][2], " vs ", model$xlevels[[xf[x]]][1], sep="")})
76 1
  lapply(varnum.mfac, function(x){rn.list[[x]] <<- c(paste(xf[x],": ref.=", model$xlevels[[xf[x]]][1], sep=""), gsub(xf[x],"   ", rn.list[[x]]))})
77 1
  if (class(fix.all.unlist)[1] == "character"){
78 0
    fix.all.unlist <- t(data.frame(fix.all.unlist))
79
  }
80 1
  rownames(fix.all.unlist) <- unlist(rn.list)
81 1
  pv.colnum <- which(colnames(fix.all.unlist) %in% c("P value", "crude P value", "adj. P value"))
82 1
  for (i in pv.colnum){
83 1
    fix.all.unlist[, i] <- ifelse(as.numeric(fix.all.unlist[, i]) < 0.001, "< 0.001", round(as.numeric(fix.all.unlist[, i]), decimal + 1))
84
  }
85
  
86

87
  ## metric
88 1
  no.obs <- model$n
89 1
  no.event <- model$nevent
90 1
  aic <- round(stats::AIC(model)[2], decimal)
91 1
  metric.mat <- cbind(c(NA, no.obs, no.event, aic), matrix(NA, 4, ncol(fix.all) - 1))
92 1
  rownames(metric.mat) <- c(NA, "No. of observations", "No. of events", "AIC")
93
  
94
  ## caption
95 1
  surv.string <- as.character(attr(model$terms, "variables")[[2]])
96 1
  time.var.name <- surv.string[2]
97 1
  status.var.name <-  surv.string[3]
98 1
  intro <- paste("Survey cox model on time ('", time.var.name, "') to event ('", status.var.name, "')", sep="")
99
  
100 1
  return(list(table = fix.all.unlist, metric = metric.mat, caption = intro))
101
  
102
  
103
}

Read our documentation on viewing source code .

Loading