Nothing
#' Estimate a quantile regression estimator of residual lifetime from survival data
#'
#' Using three estimation methods
#' (1) L1-minimization(non-smooth estimating equation
#' (2) Induced smoothing approach (smooth estimating equation
#' (3) Iterative procedure with induced smoothing approach (smooth estimating equation
#'
#' @param formula is a formula expression, of the form \code{response ~ predictors}.
#' The \code{response} is a \code{Surv} object with right censoring.
#' @param data is an optional \code{data.frame} in which to interpret the variables occurring in the \code{formula}.
#' @param t0 is the follow-up time (or basetime of analysis). The default followup time is set to 0.
#' @param Q is the quantile. The default quantile is set to 0.5.
#' @param nB is number of multiplier bootstrapping for V matrix estimation. The default number of bootstrapping is set to 100.
#' @param method is an option for specifying the methods of parameters estimation;
#' \code{smooth} is the default in which parameters estimates and their standard errors are obtained via induced smoothed estimating equations.
#' \code{nonsmooth} uses a L1-minimization method for non-smooth object functions in coefficient estimation.
#' \code{iterative} simultaneously estimates parameters and their standard errors based on the iterative updates for the parameter estimates.
#' @param se is an option for specifying the methods of standard errors estimation;
#' The available options are \code{pmb} and \code{fmb}.
#' The \code{pmb} is the default option.
#' It estimates the standard errors via partial multiplier bootstrapping and is only available
#' when \code{method = "smooth"} or \code{method = "iterative"}.
#' The \code{fmb} uses a full multiplier bootstrapping in standard errors estimation.
#' @param init is an option for specifying the initial values of the parameters estimates.
#' Available options are \code{rq} and \code{noeffect}.
#' These options correspond to the estimates from the \code{quantreg::rq()} and a zero vector, respectively.
#' Alternatively, user defined numerical vector is also allowed.
#' @param verbose Shows computation status.
#' @param control controls maximum number of iteration, tolerance of convergence and whether to display output for each iteration when
#' \code{method = "iterative"}.
#'
#' @return An object of class "\code{qris}" contains model fitting results.
#' The "\code{qris}" object is a list containing at least the following components:
#' \describe{
#' \item{coefficient}{a vector of point estimates}
#' \item{stderr}{a vector of standard error of point estimates}
#' \item{vcov}{a matrix of the estimated variance-covariance matrix}
#' \item{maxiter}{a number of iteration until convergence (only for iterative procedure)}
#' }
#'
#' @export
#' @importFrom survival Surv survfit is.Surv
#' @importFrom quantreg rq.wfit
#' @importFrom nleqslv nleqslv
#' @importFrom stats pnorm rnorm complete.cases quantile
#' @importFrom stringr str_replace
#' @import Rcpp
#'
#' @example inst/examples/ex_qris.R
#' @useDynLib qris
qris <- function(formula, data, t0 = 0, Q = 0.5, nB = 100,
method = c("smooth", "iterative", "nonsmooth"),
se = c("fmb","pmb"),
init = c("rq", "noeffect"),
verbose = FALSE, control = qris.control()) {
scall <- match.call()
mnames <- c("", "formula", "data")
cnames <- names(scall)
cnames <- cnames[match(mnames, cnames, 0)]
mcall <- scall[cnames]
mcall[[1]] <- as.name("model.frame")
m <- eval(mcall, parent.frame())
mterms <- attr(m, "terms")
allnames <- all.vars(formula)
obj <- unclass(m[,1])
method <- match.arg(method)
se <- match.arg(se)
formula0 <- formula
if (!is.Surv(m[[1]]) || attr(obj, "type") != "right")
stop("qris only supports Surv object with right censoring.", call. = FALSE)
formula[[2]] <- NULL
## Create data; the first 2 columns are from Surv(), e.g., time, status, x1, x2, ...
intModel <- FALSE
n <- nrow(obj)
X <- covariate <- model.matrix(mterms, m)
nc <- ncol(covariate)
if (nc == 1) intModel <- TRUE
## Some checks
if(t0 < 0) stop("basetime must be 0 or positive number")
if(length(Q) > 1) stop("Multiple taus not allowed in qris")
if(Q <= 0 | Q >= 1) stop("Tau must be scalar number between 0 and 1")
data <- data.frame(Z = obj[,1])
data$"log(Z-t0)" <- log(pmax(data$Z - t0, 1e-4))
data$"I[Z>t0]" <- as.numeric(data$Z >= t0)
data$delta <- c(obj[1:(n - 1),2], 1)
data <- na.omit(data)
n <- nrow(data)
## Rcpp IPCW with jump weight
# Use survfit
sv <- survfit(Surv(Z, 1 - delta) ~ 1, data)
if (t0 <= sv$time[1]) {
ghatt0 <- 1
} else {
ghatt0 <- sv$surv[min(which(sv$time>t0))-1]
}
W <- data$delta / sv$surv[findInterval(data$Z, sv$time)] * ghatt0
W[is.na(W)] <- max(W, na.rm = TRUE)
data$weight <- W
logZ <- data[,2]
I <- data[,3]
H <- diag(1 / n, nc, nc)
## Intercept model print message and output a qris object
if(nc < 2 | intModel) {
out <- NULL
s <- survfit(Surv(time, status) ~ 1, data = as.data.frame(obj), subset = time > t0)
tmp <- quantile(s, .5)
out$coefficient <- as.numeric(tmp$quantile)
cat(paste0("Intercept-only model detected; the ", round(Q, 2),
" quantile residual time is ", round(out$coefficient, 2), "\n"))
out$varNames <- colnames(covariate)
out$ci <- c(tmp$lower, tmp$upper)
out$call <- scall
out$formula <- formula0
out$data <- as.data.frame(do.call(cbind, m))
out$para <- list(method = method, Q = Q, t0 = t0, nB = nB)
class(out) <- "qris"
return(out)
}
## Define initial value
if (is.character(init)) {
init <- tryCatch(match.arg(init),
error = function(e) {
tmp <- conditionMessage(e)
tmp <- str_replace(tmp, "'arg'", "'init'")
stop(tmp, " or a numerical vector \n")
## on.exit(options(show.error.messages = F))
})
if (init == "rq") betastart <- as.vector(rq.wfit(X, data[,2], tau = Q, weights = I * W)$coef)
if (init == "noeffect") betastart <- rep(0, nc)
} else {
if (!is.numeric(init)) stop("User specified initial value must be a numerical vector")
if (length(init) != nc) stop("User specified initial value must match the number of covariates")
betastart <- as.vector(init)
}
## collect all useful information
info <- list(X = X, I = I, W = W, Q = Q, nB = nB, nc = nc, n = n, H = H, t0 = t0,
logZ = logZ, data = data, betastart = betastart, se = se,
verbose = verbose, control = control)
## pass to fit
out <- qris.fit(info, method)
out$call <- scall
out$formula <- formula0
out$data <- as.data.frame(do.call(cbind, m))
names(out$data) <- allnames
out$varNames <- colnames(covariate)
out$para <- list(method = method, Q = Q, t0 = t0, nB = nB)
if (info$control$trace) {
rownames(out$trace.coefficient) <- rownames(out$trace.stderr) <- NULL
colnames(out$trace.coefficient) <- colnames(out$trace.stderr) <- out$varNames
}
out <- out[order(names(out))]
class(out) <- "qris"
return(out)
}
#' Calculate the weighted Kaplan-Meier estimate
#'
#' @param Time is a vector of observed time, which is minimum of failure time and censored time
#' @param event is a vector of event indicator (not censored = 1, censored = 0)
#' @param wgt is a vector of weight
#'
#' @return
#' A data frame containing the following components:
#' \describe{
#' \item{deathtime}{the observed time}
#' \item{ndeath}{a vector of number of subject who experienced event at deathtime}
#' \item{nrisk}{a vector of number of subject who are possible to experience event at deathtime}
#' \item{survp}{a vector of survival probability at deathtime}
#' }
#' @noRd
ghat <- function(Time, event, wgt = 1) {
deathtime <- c(0, sort(unique(Time[event > 0])))
ndeath <- colSums(outer(Time, deathtime, "==") * event * wgt)
nrisk <- colSums(outer(Time, deathtime, ">=") * wgt)
survp <- cumprod(1 - ndeath / nrisk)
data.frame(deathtime, ndeath, nrisk, survp)
}
#' Auxiliary for Controlling qris
#'
#' Auxiliary function as user interface for \code{qris} fitting.
#'
#' When \code{trace} is TRUE, output for each iteration is printed to the screen.
#'
#' @param maxiter max number of iteration.
#' @param tol tolerance of convergence
#' @param trace a binary variable, determine whether to save output for each iteration.
#'
#' @export
#' @return A list with the arguments as components.
#' @seealso \code{\link{qris}}
qris.control <- function(maxiter = 10, tol = 1e-3, trace = FALSE) {
list(maxiter = maxiter, tol = tol, trace = trace)
}
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.