|
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
|
|
|
jinseob2kim