1

2

3
#' @title coxmeTable: Summary table of coxme.object(internal function)
4
#' @description Extract fixed effect table in coxme.object
5
#' @param mod coxme.object
6
#' @return beta, se, z, p of fixed effects
7
#' @details DETAILS
8
#' @examples 
9
#'  library(coxme)
10
#'  fit <- coxme(Surv(time, status) ~ ph.ecog + age + (1|inst), lung)
11
#'  jstable:::coxmeTable(fit)
12
#' @rdname coxmeTable
13
#' @importFrom coxme fixef
14
#' @importFrom stats pnorm
15

16

17
coxmeTable <- function (mod){
18 0
  if(!any(class(mod)=="coxme")){stop("Model not from mixed effects Cox model")}
19 1
  beta <- fixef(mod)
20 1
  nvar <- length(beta)
21 1
  nfrail <- nrow(mod$variance) - nvar
22 1
  se <- sqrt(diag(as.matrix(mod$variance))[nfrail + 1:nvar])
23 1
  z <- beta/se
24 1
  p <- 2*(1-pnorm(abs(z)))
25
  #p<- signif(1 - pchisq((beta/se)^2, 1), 2)
26 1
  table=data.frame(cbind(beta,se,z,p))
27 1
  return(table)
28
}
29

30

31

32

33
#' @title coxExp: transform the unit of coefficients in cox model(internal function)
34
#' @description Transform the unit of coefficients to "HR"
35
#' @param cox.coef cox model coefficients
36
#' @param dec Decimal point
37
#' @return The transforemed coefficients(95% CI), p-value
38
#' @details DETAILS
39
#' @examples 
40
#'  library(coxme)
41
#'  fit <- coxme(Surv(time, status) ~ ph.ecog + age + (1|inst), lung)
42
#'  jstable:::coxExp(jstable:::coxmeTable(fit))
43
#' @rdname coxExp
44
#' @importFrom stats pnorm
45

46
coxExp = function(cox.coef, dec){
47 1
  HR = paste(round(exp(cox.coef[,1]), dec), " (", round(exp(cox.coef[,1] - 1.96*cox.coef[,2]), dec), ",", round(exp(cox.coef[,1] + 1.96*cox.coef[,2]), dec), ")", sep="")
48 1
  pv = cox.coef[, "p"]
49
  #pv = 2*(1-pnorm(abs(cox.coef[, "z"])))
50 1
  return(cbind(HR, pv))
51
  
52
} 
53

54

55

56
#' @title extractAIC.coxme: Extract AIC from coxme.object
57
#' @description Extract AIC from coxme.object
58
#' @param fit coxme.object
59
#' @param scale NULL
60
#' @param k numeric specifying the 'weight' of the equivalent degrees of freedom (=: edf) part in the AIC formula.
61
#' @param ... further arguments (currently unused in base R).
62
#' @return AIC(Integreted, Penalized)
63
#' @details DETAILS
64
#' @examples 
65
#'  library(coxme)
66
#'  fit <- coxme(Surv(time, status) ~ ph.ecog + age + (1|inst), lung)
67
#'  extractAIC(fit)
68
#' @rdname extractAIC.coxme
69
#' @export 
70

71

72

73
extractAIC.coxme <- function(fit, scale = NULL, k = 2, ...){
74 0
  loglik <- fit$loglik + c(0, 0, fit$penalty)
75 0
  chi1 <- 2*diff(loglik[1:2]) 
76 0
  chi2 <- 2*diff(loglik[c(1,3)])
77 0
  c(chi1 - k*fit$df[1], chi2 - k*fit$df[2])
78
}
79

80

81

82

83
#' @title coxme.display: table for coxme.object (coxme package)
84
#' @description Make mixed effect model results from coxme.object (coxme package)
85
#' @param coxme.obj coxme.object
86
#' @param dec Decimal point, Default: 2
87
#' @return Fixed effect table, random effect, metrics, caption
88
#' @details DETAILS
89
#' @examples 
90
#'  library(coxme)
91
#'  fit <- coxme(Surv(time, status) ~ ph.ecog + age + (1|inst), lung)
92
#'  coxme.display(fit) 
93
#' @rdname coxme.display
94
#' @export
95
#' @importFrom coxme coxme 
96

97
coxme.display = function(coxme.obj, dec =2){
98 1
  model <- coxme.obj
99 0
  if(!any(class(model)=="coxme")){stop("Model not from mixed effects Cox model")}
100
  
101 1
  xf <- attr(model$terms, "term.labels") # Independent vars
102 1
  if(length(grep("strata", xf)) > 0){
103 0
    xf <- xf[-grep("strata",xf)]
104
  }
105
  
106 1
  formula.surv = as.character(model$formulaList$fixed)[2]
107 1
  formula.ranef = as.character(model$formulaList$random)
108 1
  mdata = data.frame(get(as.character(model$call)[3]))
109
  
110 1
  if(length(xf) == 1){
111 0
    uni.res = coxmeTable(coxme(as.formula(paste(formula.surv, "~", xf," + ", formula.ranef, sep="")), data = mdata))
112 0
    rn.uni <- lapply(list(uni.res), rownames)
113 0
    fix.all = coxExp(uni.res, dec = dec)
114 0
    colnames(fix.all) = c("HR(95%CI)", "P value")
115 0
    rownames(fix.all) = names(model$coefficients)
116
  } else{
117 1
    unis <- lapply(xf, function(x){coxmeTable(coxme(as.formula(paste(formula.surv, "~", x," + ", formula.ranef, sep="")), data = mdata))})
118 1
    rn.uni <- lapply(unis, rownames)
119 1
    unis2 <- Reduce(rbind, unis)
120 1
    uni.res <- unis2
121 1
    fix.all = cbind(coxExp(uni.res, dec = dec), coxExp(coxmeTable(model), dec = dec))
122 1
    colnames(fix.all) = c("crude HR(95%CI)", "crude P value", "adj. HR(95%CI)", "adj. P value")
123 1
    rownames(fix.all) = names(model$coefficients)
124
  }
125
  
126
  ## rownames
127 1
  fix.all.list = lapply(1:length(xf), function(x){fix.all[rownames(fix.all) %in% rn.uni[[x]], ]})        
128 1
  varnum.mfac = which(lapply(fix.all.list, length) > ncol(fix.all))
129 0
  lapply(varnum.mfac, function(x){fix.all.list[[x]] <<- rbind(rep(NA, ncol(fix.all)), fix.all.list[[x]])})
130 1
  fix.all.unlist = Reduce(rbind, fix.all.list)
131
  
132 1
  rn.list = lapply(1:length(xf), function(x){rownames(fix.all)[rownames(fix.all) %in% rn.uni[[x]]]})
133 1
  varnum.2fac = which(lapply(xf, function(x){length(sapply(mdata, levels)[[x]])}) == 2)
134 0
  lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xf[x], ": ", levels(mdata[, xf[x]])[2], " vs ", levels(mdata[, xf[x]])[1], sep="")})
135 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]]))})
136 1
  if (class(fix.all.unlist)[1] == "character"){
137 0
    fix.all.unlist = t(data.frame(fix.all.unlist))
138
  }
139 1
  rownames(fix.all.unlist) = unlist(rn.list)
140 1
  pv.colnum = which(colnames(fix.all.unlist) %in% c("P value", "crude P value", "adj. P value"))
141 1
  for (i in pv.colnum){
142 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))
143
  }
144

145

146
  ## random effect
147
  #ranef = unlist(model$vcoef)
148
  #ranef.out = round(ranef, dec)
149 1
  ranef.out = lapply(model$vcoef, function(x){round(x, dec)})
150 1
  ranef.mat = cbind(c(NA, unlist(ranef.out)), matrix(NA, length(unlist(ranef.out)) + 1, ncol(fix.all) - 1))
151 1
  rownames(ranef.mat)  = c("Random effect", paste(names(ranef.out), "(", unlist(lapply(ranef.out, names)), ")", sep=""))
152
  
153
  
154
  ## metric
155 1
  no.grp = unlist(lapply(model$frail, length))
156 1
  no.obs = model$n[2]
157 1
  no.event = model$n[1]
158 1
  metric.mat = cbind(c(NA, no.grp, no.obs, no.event), matrix(NA, length(no.grp) + 3, ncol(fix.all) - 1))
159 1
  rownames(metric.mat) = c(NA, paste("No. of groups(", names(no.grp), ")", sep=""), "No. of observations", "No. of events")
160
  
161
  ## Integrated ll
162
  #ll = model$loglik[2]
163
  #aic = -2 * ll -2*model$df[1] 
164
  
165
  ## caption
166 1
  surv.string <- as.character(attr(model$terms, "variables")[[2]])
167 1
  time.var.name <- surv.string[2]
168 1
  status.var.name <-  surv.string[3]
169 1
  intro <- paste("Mixed effects Cox model on time ('", time.var.name, "') to event ('", status.var.name, "')", " - Group ",paste(names(model$vcoef), collapse = ", ") ,sep="")
170 1
  var.names0 <- attr(model$terms, "term.labels")
171 0
  if(length(grep("strata", var.names0))>0) {intro <- paste(intro, " with '", var.names0[grep("strata", var.names0)], "'", sep="" )}
172
  
173 1
  return(list(table = fix.all.unlist, ranef = ranef.mat, metric = metric.mat, caption = intro))
174
}
175

176

177

178

179

180

Read our documentation on viewing source code .

Loading