Nothing
# Copyright (C) 2014 Stephan Dlugosz, Limin Peng and Ruosha Li
# distributed under the terms of the GNU public license
crrQR <-
# function for quantile regression modeling of subdistribution functions
# arguments:
# ftime = vector of failure/censoring times (log-transformed)
# fstatus = vector with a unique code for each failure type and a
# separate code for censored observations
# X = (nobs x ncovs) matrix of fixed covariates
# failcode = code of fstatus that denotes the failure type of interest
# cencode = code of fstatus that denotes censored observations (default: 0)
# tau.range = subset of (0,1) denoting the range of qwuantiles to be computed
# tau.step = stepsize to produce a grid for covering the tau.range
# subset = logical vector length(ftime) indicating which cases to include
# na.action = function defining action to take for cases that have NA for
# any of ftime, fstatus, X, or subset.
# variance = should coefficient estimator variances be computed?
#
function(ftime, fstatus, X,
failcode = 1, cencode = 0, tau.range = c(0.01,0.99), tau.step = 0.01,
subset, na.action = na.omit, rq.method = "br", variance = TRUE) {
tau.range <- sort(tau.range)
if (length(tau.range) != 2 || tau.range[1] < 0 || tau.range[2] > 1) stop('invalid range for tau')
if (!missing(subset) && length(subset) != length(ftime)) stop('dimensions of subset and ftime do not match')
if (nrow(X) != length(ftime) | length(ftime) != length(fstatus)) stop('dimensions of ftime, fstatus and X do not match')
n <- length(ftime)
check <- data.frame(ftime = ftime, fstatus = fstatus)
if (!missing(subset)) check <- cbind(check, subset)
check <- cbind(check, X)
check <- na.action(check)
ftime <- check$ftime
fstatus <- check$fstatus
if (!missing(subset)) {
subset <- check$subset
X <- check[,-(1:3), drop = FALSE]
} else {
X <- check[,-(1:2), drop = FALSE]
}
rm(check)
if (!missing(subset)) {
ftime <- ftime[subset]
fstatus <- fstatus[subset]
X <- X[subset,, drop = FALSE]
}
X <- as.matrix(X)
ret <- crrQR.int(ftime, fstatus, X, tau.L = tau.range[1], tau.U = tau.range[2], tau.step = tau.step,
outcome = failcode, cencode = cencode, rq.method = rq.method, orig.num = n, variance = variance)
return(ret)
}
crrQR.int <- function(ftime, fstatus, X, tau.L, tau.U, tau.step,
outcome = 1, cencode = 0, noconst = FALSE, variance = TRUE, offset=0, max.lmt = 10^8,
rq.method = "br", orig.num = nrow(X), retry = TRUE, ...) {
eddcmp <- function(M) {
ev <- eigen(M)
return(list(evectors = ev$vectors, evalues = ev$values, im.evalues = NULL))
}
if (!noconst) {
X <- cbind(rep.int(1, dim(X)[1]), X)
dimnames(X)[[2]][1] <- "const"
}
o <- order(ftime, decreasing = TRUE)
ftime <- ftime[o]
fstatus <- fstatus[o]
X <- X[o,, drop = FALSE]
cvt.length <- ncol(X) ### number of covariates
tau.seq <- seq(tau.L, tau.U, tau.step)
L.tau.seq <- length(tau.seq)
num <- nrow(X)
FT0 <- as.numeric(fstatus == cencode)
cens <- (sum(FT0) > 0)
FT1 <- as.numeric(fstatus == outcome)
smallest <- .Machine$double.eps ^ 0.5
csurv.x <- 1.0
if (cens) {
#cat(' with censoring\n')
csurv.fit <- survfit(Surv(ftime, FT0) ~ 1)
#L.tpt <- length(csurv.fit$time)
csurv.indx <- unlist(lapply(ftime,
function(et) length(which(et > csurv.fit$time | (abs(et - csurv.fit$time) < smallest)))))
csurv.x <- csurv.fit$surv[csurv.indx] * FT1 + (1 - FT1)
rm(csurv.indx)
gc()
}
FT1.csurv.x <- FT1 / csurv.x
pseudo.resp <- c((ftime - offset) * FT1.csurv.x, max.lmt, max.lmt)
n.cvt.1 <- crossprod(-X, FT1.csurv.x)
est.beta.seq <- NULL
if (variance) {
est.var.seq <- NULL
inf.est.func <- list()
ind.conv.var <- rep(1, cvt.length)
est.var <- rep(0, cvt.length)
curr.inf <- array(rep(0, num*(cvt.length)), c(num, cvt.length))
compare.vec <- rank(-ftime, ties.method = "max")
}
L.bsq <- 1
ind.conv <- 1
X.FT1.csurv.x <- X * FT1.csurv.x
#cat('computing')
goon <- TRUE
while (goon) {
est.beta <- rep(0, cvt.length)
curr.tau <- tau.seq[L.bsq]
n.cvt.2 <- apply(2 * X * curr.tau, 2, sum)
pseudo.cvt <- rbind(X.FT1.csurv.x, as.vector(n.cvt.1), as.vector(n.cvt.2))
est.beta.obj <- try(switch(rq.method, fn = rq.fit.fnb(pseudo.cvt, pseudo.resp, ...),
fnb = rq.fit.fnb(pseudo.cvt, pseudo.resp, ...),
fnc = rq.fit.fnc(pseudo.cvt, pseudo.resp, ...),
pfn = rq.fit.pfn(pseudo.cvt, pseudo.resp, ...),
br = rq.fit.br(pseudo.cvt, pseudo.resp, ci = FALSE, ...),
lasso = rq.fit.lasso(pseudo.cvt, pseudo.resp, ...),
scad = rq.fit.scad(pseudo.cvt, pseudo.resp, ...),
{
what <- paste("rq.fit.", rq.method, sep = "")
if (exists(what, mode = "function")) (get(what, mode = "function"))(pseudo.cvt,
pseudo.resp, ...) else stop(paste("unimplemented method:",
rq.method))
}))
if (inherits(est.beta.obj, "try-error")) {
if (retry) {
cat('retry with method br ')
est.beta.obj <- rq.fit.br(pseudo.cvt, pseudo.resp, ci = FALSE, ...)
} else {
stop('estimation not successful, try different method')
}
}
est.beta <- est.beta.obj$coef
ind.conv <- as.numeric(abs(est.beta.obj$residuals[num + 1L]) > 1) * as.numeric(abs(est.beta.obj$residuals[num + 2L]) > 1)
if (ind.conv == 1) {
est.beta.seq <- rbind(est.beta.seq, est.beta)
if (variance) {
helper <- as.numeric((ftime) < X %*% as.vector(est.beta)) * FT1.csurv.x
tmp.1 <- ((helper - curr.tau) * X)
var.matrix.cmp.1 <- crossprod(tmp.1, tmp.1) / num
ZI.mat <- X * helper
tmp.2.1 <- .Call(comptmp21, num, ZI.mat, compare.vec)
tmp.2 <- FT0*(tmp.2.1/compare.vec)
var.matrix.cmp.2 <- crossprod(tmp.2, tmp.2) / num
var.matrix <- var.matrix.cmp.1 - var.matrix.cmp.2
var.dcmp <- eddcmp(var.matrix)
sigma.sqrt <- (var.dcmp$evector %*% diag(sqrt(var.dcmp$evalues),
length(var.dcmp$evalues),
length(var.dcmp$evalues)) %*% solve(var.dcmp$evector)) #eigen decomposition
est.beta.var <- NULL
ind.conv.var <- NULL
for (k in 1:cvt.length) {
pseudo.cvt.var <- rbind(X.FT1.csurv.x, as.vector(n.cvt.1), as.vector(n.cvt.2) + 2 * sigma.sqrt[,k] * sqrt(num))
est.beta.var.obj <- try(switch(rq.method, fn = rq.fit.fnb(pseudo.cvt.var, pseudo.resp, ...),
fnb = rq.fit.fnb(pseudo.cvt.var, pseudo.resp, ...),
fnc = rq.fit.fnc(pseudo.cvt.var, pseudo.resp, ...),
pfn = rq.fit.pfn(pseudo.cvt.var, pseudo.resp, ...),
br = rq.fit.br(pseudo.cvt.var, pseudo.resp, ci = FALSE, ...),
lasso = rq.fit.lasso(pseudo.cvt.var, pseudo.resp, ...),
scad = rq.fit.scad(pseudo.cvt.var, pseudo.resp, ...),
{
what <- paste("rq.fit.", rq.method, sep = "")
if (exists(what, mode = "function")) (get(what, mode = "function"))(pseudo.cvt,
pseudo.resp, ...) else stop(paste("unimplemented method:",
rq.method))
}))
if (inherits(est.beta.var.obj, "try-error")) {
if (retry) {
cat('retry with method br ')
est.beta.var.obj <- rq.fit.br(pseudo.cvt, pseudo.resp, ci = FALSE, ...)
} else {
stop('estimation not successful, try different method')
}
}
ind.conv.var <- c(ind.conv.var, as.numeric(abs(est.beta.var.obj$residuals[num + 1]) > 1) *
as.numeric(abs(est.beta.var.obj$residuals[num + 2]) > 1))
est.beta.var <- cbind(est.beta.var, est.beta.var.obj$coef - est.beta)
}
if (sum(ind.conv.var) == cvt.length) {
est.var <- diag(tcrossprod(est.beta.var, est.beta.var))
est.inv.deriv.matrx <- est.beta.var %*% solve(sigma.sqrt)*sqrt(num)
est.inv.deriv.matrx <- (est.inv.deriv.matrx + t(est.inv.deriv.matrx))/2
curr.inf <- tcrossprod(tmp.1 - tmp.2, est.inv.deriv.matrx)
} else {
goon <- FALSE
}
inf.est.func[[L.bsq]] <- curr.inf
est.var.seq <- rbind(est.var.seq, est.var)
}
L.bsq <- L.bsq + 1
} else {
goon <- FALSE
}
goon <- (L.bsq <= L.tau.seq & goon)
cat('.')
}
cat('done.\n')
if (L.bsq < L.tau.seq) {
cat('\n')
cat('stopping at tau=')
cat(tau.seq[L.bsq])
cat('\n')
L.tau.seq <- L.bsq - 1
if (L.bsq <= 1) return(list(tau.seq = NULL, L.bsq = 0))
tau.seq <- tau.seq[1:L.tau.seq]
tau.U <- tau.seq[L.tau.seq]
}
dimnames(est.beta.seq)[[1]] <- tau.seq
if (variance) dimnames(est.var.seq)[[1]] <- tau.seq
# return results
ret <- list(call = match.call(),
beta.seq = est.beta.seq,
tau.seq = tau.seq,
cvt.length = cvt.length,
n = num,
n.missing = orig.num - num)
if (variance) {
ret$var.seq = est.var.seq
ret$inf.func = inf.est.func
}
class(ret) <- c("crrQR")
return(ret)
}
summary.crrQR <-
# compute confidence intervals and tests
# as well as some summary statistics
# for a crrQR object
function(object, conf.int = 0.95, digits = max(options()$digits - 5, 2), ...) {
if (is.null(object$var.seq)) stop('please re-run crrQR with variance=TRUE')
L.tau.seq <- length(object$tau.seq)
tau.L <- object$tau.seq[1]
tau.U <- object$tau.seq[L.tau.seq]
tau.step <- {tau.U - tau.L}/L.tau.seq
tau.w1 <- (as.numeric(object$tau.seq < (tau.L + tau.U) / 2))[1:(L.tau.seq - 1)] ### weight for constancy test
inf.ave <- NULL
est.ave <- NULL
inf.cnst.test.w1 <- NULL
cnst.test.w1 <- NULL
for (k in 1:(object$cvt.length)) {
infs <- object$inf.func
infs[[length(infs)]] <- NULL
infs <- do.call(cbind, lapply(infs, function(x) x[,k,drop = FALSE]))
inf.ave <- cbind(inf.ave, apply(infs*tau.step, 1, sum) / (tau.U - tau.L))
est.ave <- c(est.ave, sum(object$beta.seq[1:(L.tau.seq - 1), k, drop = FALSE]*tau.step) / (tau.U - tau.L))
diff <- infs - matrix(rep.int(inf.ave[,k],ncol(infs)), ncol = (ncol(infs)), byrow = FALSE)
inf.cnst.test.w1 <- cbind(inf.cnst.test.w1, apply(t(t(diff)*tau.w1)*tau.step, 1, sum) / sum(tau.w1 * tau.step))
cnst.test.w1 <- c(cnst.test.w1,
sum((object$beta.seq[1:(L.tau.seq - 1), k, drop = FALSE] - est.ave[k]) * tau.w1*tau.step) /
sum(tau.w1*tau.step))
}
names(est.ave) <- dimnames(object$beta.seq)[[2]]
var.ave <- apply(inf.ave^2 / object$n, 2, mean)
names(var.ave) <- names(est.ave)
var.cnst.test.w1 <- apply(inf.cnst.test.w1^2 / object$n, 2, mean)
names(cnst.test.w1) <- names(est.ave)
names(var.cnst.test.w1) <- names(est.ave)
ret <- list(call = object$call,
ave.eff = est.ave,
var.ave.eff = var.ave,
p.signf.test = signif(1 - pchisq(est.ave^2/var.ave, 1), digits),
cnst.test = cnst.test.w1,
var.cnst.test = var.cnst.test.w1,
p.cnst.test = signif(1 - pchisq(cnst.test.w1^2/var.cnst.test.w1, 1), digits),
n = object$n,
n.missing = object$n.missing)
class(ret) <- "summary.crrQR"
return(ret)
}
print.summary.crrQR <-
# prints the summary object to a crrQR model
function(x, digits=max(options()$digits - 4, 3), ...) {
cat("Competing Risks Quantile Regression\n\n")
if (!is.null(x$call))
{ cat("Call:\n")
dput(x$call)
cat("\n")
}
savedig <- options(digits = digits)
on.exit(options(savedig))
out.bcr <- data.frame(x$ave.eff, sqrt(x$var.ave.eff), x$p.signf.test, x$cnst.test, x$var.cnst.test, x$p.cnst.test)
names(out.bcr) <- c("ave.eff", "se.ave.eff", "p.signf.test", "cnst.test", "var.cnst.test", "p.cnst.test")
print(out.bcr)
cat("\n")
cat("Num. cases =", x$n)
if (x$n.missing > 0) cat(" (", x$n.missing, " cases omitted due to missing values)", sep = "")
cat("\n")
invisible()
}
predict.crrQR <-
# for a crrQR object x, estimates subdistributions at covariate
# combinations given by rows of X. The terms in X
# must correspond exactly to the corresponding call to crrQR.
# add rearrangement=TRUE for using the rearrangement method to ensure increasing CiC predictions
# (cf. Chernozhukov and Fernandez-Val and Galichon, Quantile and probability curves without crossing, Econometrica 78, 1093--1125)
function(object, x, rearrangement = FALSE, ...) {
betas <- object$beta.seq
if (inherits(x, "data.frame")) x <- as.matrix(x)
if (inherits(x, "matrix")) {
if (ncol(x) != ncol(betas) - 1) stop('dimensions of X and estimated betas do not match')
x <- x[,dimnames(betas)[[2]][-1], drop = FALSE]
x <- cbind(1, x)
} else {
if (length(x) != ncol(betas) - 1) stop('dimensions of X and estimated betas do not match')
x <- x[dimnames(betas)[[2]][-1]]
x <- matrix(c(1, x), nrow = 1)
}
xb <- tcrossprod(betas, x)
exb <- exp(xb)
if (rearrangement) exb <- apply(exb, 2, sort)
lhat <- cbind(exb, object$tau.seq)
class(lhat) <- 'predict.crr'
return(lhat)
}
plot.predict.crrQR <-
# plots estimated subdistributions from predict.crrQR
function(x, lty=1:(ncol(x) - 1), color = 1, ylim = c(0,max(x[,ncol(x)])), xmin = 0, xmax = max(x[,-ncol(x)]), ...) {
if (length(lty) < ncol(x) - 1) lty <- rep(lty[1], ncol(x) - 1)
if (length(color) < ncol(x) - 1) color <- rep(color[1], ncol(x) - 1)
if (xmax < max(x[,-ncol(x)])) x <- x[apply(x, 1,max) < xmax,]
plot(c(xmin,xmax), ylim, type="n", ...)
for (j in 1:(ncol(x)-1))
lines(c(xmin,rep(x[,j], rep.int(2,nrow(x))), max(x[,j])),
c(0, 0, rep(x[,ncol(x)], rep.int(2,nrow(x)))),
lty = lty[j], col = color[j])
}
print.crrQR <-
# prints the essential parts of a crrQR Object
function(x, ...) {
cat('coefficients:\n')
print(x$beta.seq)
v <- sqrt(x$var.seq)
cat('standard errors:\n')
print(signif(v,4),...)
v <- 2 * (1 - pnorm(abs(x$beta.seq) / v))
cat('two-sided p-values:\n')
print(signif(v, 2),...)
invisible()
}
plot.crrQR <-
# plots the profile graph of the coefficients
# subset defines the subset of betas to plot
function(x, subset=NULL, main=NULL, ...) {
if (is.null(subset)) subset <- 1:(x$cvt.length)
if (is.null(main)) main <- ""
CI.lb.seq <- x$beta.seq - 1.96 * sqrt(x$var.seq)
CI.ub.seq <- x$beta.seq + 1.96 * sqrt(x$var.seq)
if (length(subset) > 1) op <- par(mfrow = c(ceiling((x$cvt.length + 1) / 2), 2))
#grp.frm.1<-apply(results[[risk]]$CI.lb.seq, 2, min)*0.99
#grp.frm.2<-apply(results[[risk]]$CI.ub.seq, 2, max)*1.01
grp.frm.1 <- apply(CI.lb.seq, 2, function(x) { quantile(x, probs = c(0.1)) } )*0.99
grp.frm.2 <- apply(CI.ub.seq, 2, function(x) { quantile(x, probs = c(0.9)) } )*1.01
char.cvt <- dimnames(x$beta.seq)[[2]]
for (k in subset) {
plot(x$tau.seq, x$beta.seq[, k], ylim = c(grp.frm.1[k],
grp.frm.2[k]), type = "s", xlab = "tau", ylab = "Coefficient",
main = paste(main, char.cvt[k], sep = ":"), ...)
lines(x$tau.seq, CI.lb.seq[,k], type = "s", lty = 3)
lines(x$tau.seq, CI.ub.seq[,k], type = "s", lty = 3)
lines(x$tau.seq, rep(0, length(x$tau.seq)), lty = 2)
lines(x$tau.seq, rep(x$est.ave[k], length(x$tau.seq)), col = 2)
}
if (length(subset) > 1) par(op)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.