crq.fit.pow <- function(x, y, yc, tau=0.5, weights=NULL, start, left=TRUE,maxit = 500){
x <- as.matrix(x)
y <- as.vector(y)
n <- nrow(x)
p <- ncol(x)
if(missing(yc)) {
yc <- rep(0,n)
if(!left) warning("default right (!?) censoring at zero")
}
yc <- as.vector(yc)
if(length(weights)){
if (any(weights < 0))
stop("negative weights not allowed")
contr <- attr(x, "contrasts")
x <- x * weights
y <- y * weights
}
if(left) {
x <- -x
y <- -y
yc <- -yc
tau <- 1-tau
}
# Starting value of beta
if(missing(start)){
tol <- .Machine$double.eps^(2/3)
r <- rq.fit.br(x,y,tau)$residuals
h <- which(abs(r) < tol)[1:p]
}
else if(start[1] == "global"){
# Generate H matrix
H <- combos(n,p)
m <- ncol(H)
# Starting value of beta
U <- solve(x[H[,1],])
bh <- U %*% y[H[,1]]
f <- .Fortran("brutpow",
as.integer(n),
as.integer(p),
as.integer(m),
as.integer(H),
as.double(x),
as.double(y),
as.double(yc),
as.double(bh),
as.double(tau),
as.double(U),
double(p),
double(p),
kminz = integer(1),
nflag = as.integer(0),
PACKAGE = "quantreg")
if(f$nflag!=0)
warning(switch(f$nflag,
"Error in pivot: hout not in h",
"Error in pivot: hin in h",
"Error in pivot: hin out of bounds",
"Error in findk: k not found"
))
k <- f$kminz
bh <- solve(x[H[,k],],y[H[,k]])
residuals <- as.matrix(y - pmin(yc, x %*% bh))
return(list(coefficients= bh,residuals = residuals))
}
else if(is.numeric(start) && length(start) == p){
if(all(start %in% 1:n)) h <- start
else h <- order(abs(y - x %*% start))[1:p]
}
else
stop("Invalid starting value")
xhinv <- try(solve(x[h,]))
if(class(xhinv) == "try-error")
stop("Singular basic solution generated by 'start'")
f <- .Fortran("powell",
as.integer(n),
as.integer(p),
as.integer(2*p),
as.double(x),
as.double(y),
as.double(yc),
coef = double(p),
as.double(tau),
as.integer(h),
double(n),
as.double(xhinv),
double(n),
double(2*p),
double(p),
double(p),
as.integer(maxit),
nflag = as.integer(0),
PACKAGE = "quantreg")
if(f$nflag!=0)
warning(switch(f$nflag, "Max iterations reached",
"Solution may be nonunique",
"Error in pivot: hout not in h",
"Error in pivot: hin in h",
"Error in pivot: hin out of bounds"
))
coef <- f$coef
residuals <- as.matrix(y - pmin(yc, x %*% coef))
return (list(coefficients= coef, residuals = residuals))
}
Curv <- function (y, yc, ctype = c("left", "right"))
{
nn <- length(y)
if(length(yc) != nn)
stop("Event times and censoring times of different length")
ctype <- match.arg(ctype)
if(ctype == "right" && any(y > yc))
stop("Event times can not exceed ctimes for right censoring")
if(ctype == "left" && any(y < yc))
stop("Event times can not be less than ctimes for left censoring")
if(!(ctype == "left" || ctype == "right"))
stop("Invalid ctype for method Powell")
ss <- cbind(y, yc)
dimnames(ss) <- list(NULL, c("time", "ctime"))
attr(ss, "ctype") <- ctype
attr(ss, "type") <- ctype
class(ss) <- "Surv"
ss
}
QTECox <- function(x, smooth = TRUE){
# compute quantile treatment effect for a Cox PengHuang fit
g <- survival::survfit(x)
if(smooth)
g <- supsmu(1-g$surv,g$time)
taus <- (g$x[-1] + g$x[-length(g$x)])/2
Qhat <- (g$y[-1] + g$y[-length(g$y)])/2
dQ <- diff(Qhat)/diff(taus)
taus <- taus[-1]; Qhat <- Qhat[-1]
QTE <- outer(dQ * (1-taus) * log(1-taus) / Qhat, coef(x))
list(QTE = QTE , taus = taus)
}
boot.crq <- function(x, y, c, taus, method, ctype = "right", R=100,
mboot, bmethod = "jack", ...)
{
n <- length(y)
p <- ncol(x)
if(missing(mboot)) {
if(bmethod=="jack") mboot <- 2*ceiling(sqrt(n))
else mboot <- n }
A <- array(0,dim=c(p,length(taus),R))
for (i in 1:R){
if(bmethod == "jack") {
s <- sample(1:n,mboot)
yb <- y[-s]
xb <- x[-s, ]
cb <- c[-s]
w <- rep(1,n-mboot)
}
else if(bmethod == "xy-pair"){
w <- table(sample(1:n,mboot,replace=TRUE))
s <- as.numeric(names(w))
w <- as.numeric(w)
yb <- y[s]
xb <- x[s,]
cb <- c[s]
}
else if(bmethod == "Bose"){
w <- rexp(n)
yb <- y
xb <- x
cb <- c
}
else
stop("invalid bmethod for boot.crq")
if(method == "Portnoy")
a <- crq.fit.por(xb,yb,cb, weights = w, ctype = ctype, ... )
#else if(method == "Portnoy2")
#a <- crq.fit.por2(xb,yb,cb, weights = w, ctype = ctype, ... )
else if(method == "PengHuang")
a <- crq.fit.pen(xb,yb,cb, weights = w, ctype = ctype, ... )
else
stop("Invalid method for boot.crq")
if((i %% floor(R/10)) == 0 & n > 3000)
cat(paste("bootstrap roughly ",100*(i/R)," percent complete\n"))
A[,,i] <- coef(a,taus)
}
list(A = A, n = length(y), mboot = mboot, bmethod = bmethod)
}
coef.crq <- function(object, taus = 1:4/5, ...)
{
# Extract coefficients from the crq solution array
if(min(taus) < 0 || max(taus) > 1) stop("taus out of range [0,1]")
if(length(object$tau) == 1){
coef <- object$coefficients
return(coef)
}
taus <- sort(taus)
S <- object$sol
ctype <- object$ctype
r <- S[1, ]
r <- c(r[1],r)
if(is.unsorted(r)) r <- r[-length(r)] #kludge until why this happens is found
B <- S[-1,,drop = FALSE]
B <- t(cbind(B,B[,ncol(B),drop = FALSE]))
ts <- taus[taus > min(r) & taus < max(r)]
bin <- findInterval(ts,r)
wgt <- (ts - r[bin])/(r[bin + 1] - r[bin])
binlag <- bin - 1
binlag[binlag == 0] <- 1
coef <- t(wgt * B[bin, , drop = FALSE] + (1 - wgt) * B[binlag, , drop = FALSE])
nna <- length(taus) - length(ts)
if(nna > 0) {
if (ctype == "left")
coef <- cbind(matrix(NA, nrow(coef), nna), coef)
else
coef <- cbind(coef, matrix(NA, nrow(coef), nna))
}
taulabs <- paste("tau=", format(round(taus, 3)))
dimnames(coef)[[2]] <- taulabs
coef[-nrow(coef),] # Delete Qhat entries
}
crq <- function (formula, taus, data, subset, weights, na.action,
method = c("Powell","Portnoy","Portnoy2", "PengHuang"), contrasts = NULL, ...)
{
if(!requireNamespace("survival", quietly = TRUE))
stop("crq requires survival package to be installed")
#if(method == "Portnoy2") stop("Portnoy2 method not (yet) implemented")
Surv <- survival::Surv
call <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action"),
names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval.parent(mf)
if (method == "model.frame")
return(mf)
mt <- attr(mf, "terms")
X <- model.matrix(mt, mf, contrasts)
weights <- as.vector(model.weights(mf))
Y <- model.extract(mf, "response")
eps <- .Machine$double.eps^(2/3)
if(!inherits(Y,"Surv"))
stop("Response must be a survival object")
method <- match.arg(method)
if(method == "Powell") {
ctype <- attr(Y, "ctype")
if(!ctype %in% c("right","left"))
stop("Only right or left censoring Surv objects are allowed")
left <- (ctype == "left")
if (any(taus < -eps) || any(taus > 1 + eps))
stop("invalid taus: taus should be >= 0 and <= 1")
y <- Y[,1]
cen <- Y[,2]
if (length(taus) > 1) {
coef <- matrix(0, ncol(X), length(taus))
fitted <- resid <- matrix(0, nrow(X), length(taus))
for (i in 1:length(taus)) {
z <- crq.fit.pow(X, y, cen, tau = taus[i], weights, left = left, ...)
coef[, i] <- z$coefficients
resid[, i] <- z$residuals
fitted[, i] <- y - z$residuals
}
taulabs <- paste("tau=", format(round(taus, 3)))
dimnames(coef) <- list(dimnames(X)[[2]], taulabs)
dimnames(resid) <- list(dimnames(X)[[1]], taulabs)
fit <- list(coefficients = coef, residuals = resid, fitted.values = fitted)
fit$tau <- taus
class(fit) <- "crqs"
}
else {
fit <- crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...)
fit$tau <- taus
class(fit) <- "crq"
}
}
else if(method == "Portnoy"){
ctype <- "right"
if(attr(Y,"type") != "right"){
if(attr(Y,"type") == "left") ctype <- "left"
else stop("Only right censoring Surv objects are allowed for Portnoy method")
}
y <- Y[,1]
cen <- Y[,2]
fit <- crq.fit.por(X, y, cen, weights, ctype = ctype, ...)
class(fit) <- "crq"
}
else if(method == "Portnoy2"){
ctype <- "right"
if(attr(Y,"type") != "right"){
if(attr(Y,"type") == "left") ctype <- "left"
else stop("Only right censoring Surv objects are allowed for Portnoy method")
}
y <- Y[,1]
cen <- Y[,2]
fit <- crq.fit.por2(X, y, cen, weights, ctype = ctype, ...)
class(fit) <- "crq"
}
else if(method == "PengHuang"){
ctype <- "right"
if(attr(Y,"type") != "right"){
if(attr(Y,"type") == "left") ctype <- "left"
else stop("Only right censoring Surv objects are allowed for Portnoy method")
}
y <- Y[,1]
cen <- Y[,2]
fit <- crq.fit.pen(X, y, cen, weights, ctype = ctype, ...)
class(fit) <- "crq"
}
else
stop("Method not defined for crq")
fit$terms <- mt
fit$call <- call
fit$formula <- formula(mt)
fit$method <- method
fit$contrasts <- contrasts
fit$ctype <- ctype
attr(fit, "na.message") <- attr(m, "na.message")
fit
}
predict.crqs <-
function (object, newdata, type = NULL, ...) {
pred <- predict.rqs(object, newdata, ...)
}
predict.crq <-
function (object, newdata, ...) {
method <- object$method
if(method %in% c("Portnoy","Portnoy", "PengHuang")){ #Kludge to make crq sol matrix look like rq sol matrix
p2 <- nrow(object$sol)
object$sol <- object$sol[c(1,p2,p2,2:(p2-1)),]
}
pred <- switch(method,
"Powell" = predict.rq(object, newdata, ...),
"Portnoy" = predict.rq.process(object, newdata, ...),
"Portnoy2" = predict.rq.process(object, newdata, ...),
"PengHuang" = predict.rq.process(object, newdata, ...)
)
}
crq.fit.por <- function(x, y, cen, weights = NULL, grid, ctype = "right")
{
p <- ncol(x)
n <- length(y)
cen <- 1 - cen #NB: Fortran routine wants censoring indicator flipped (!!!)
mp <- max(n + 5 + max(1, sum(cen)),n+p)
eps <- 1e-04
kmax <- 10000 # perhaps this should be parameter in future
if(length(weights)){
if (any(weights < 0))
stop("negative weights not allowed")
contr <- attr(x, "contrasts")
x <- weights * x
y <- weights * y
}
if(ctype == "left") y <- -y
if(missing(grid)) grid <- seq(1/n,1-1/n,by=1/(5 + 3 * n^.4))
if(is.numeric(grid)){
ginit <- min(grid)
dgrid <- diff(grid)
gstep <- median(dgrid)
if(any(dgrid <0)) stop("grid is not monotonic")
if(gstep < eps) stop("grid stepsize too small")
nsol <- 3*n
mw = -1
}
else if(grid == "pivot"){
nsol <- 3*n
ginit <- 1/(2*n)
gstep <- 1/(2*n)
mw <- 20
}
else
stop("Invalid grid")
z <- .Fortran("crq",
as.integer(n),
as.integer(p),
as.integer(mp),
as.integer(p+2),
as.double(x),
as.double(y),
as.integer(cen),
as.double(ginit),
as.integer(mw),
as.double(gstep),
ift = integer(1),
h = integer(p),
xh = double(p*p),
wa = double(mp*p),
wb = double(mp),
wc = double(mp*(p+2)),
wd = double(mp),
we = double(mp),
wf = double(p),
iflag = integer(mp),
as.integer(nsol),
sol = double(nsol*(p+2)),
lsol = integer(1),
icen = integer(n),
tcen = double(n),
lcen = integer(1),
PACKAGE = "quantreg")
nw <- z$h[1]
flag <- z$ift
msg <- switch(flag,
paste("Error in input dimensions, n,p,mw "),
paste("Error in input dimensions, n,p,mw "),
paste("Error in input dimensions, n,p,mw "),
paste("Less than p=",p,"observations above tau = 0 solution"),
paste("Possible degeneracy at",nw,"tau values.",
"$tau.degen: first mp =", n + 5 + sum(cen)," such tau values"),
paste("Number of pivots to be saved in sol > nsol.",
"Redefine nsol: use nsol < n to save for tau = i/(nsol-1)"),
paste("Error with partial return: possible degeneracies",
"Max number of rq calls exceeded: dither x or increase mw"),
paste("Premature stop: defective conditional distribution"),
paste("Simplex iteration limit exceeded -- consider dithering y"))
#if(flag > 0 && flag != 5 && flag < 10)
if(flag %in% c(1:4,6,7))
ifelse(flag <= 3,stop(msg),warning(msg))
J <- z$lsol
B <- matrix(z$sol, nrow=p+2, ncol=nsol, byrow=FALSE)[,1:J, drop = FALSE]
if(B[1,J] < B[1,J-1]) B <- B[,-J] # SLP hack Oct 27 2014
ic <- z$icen
sp <- (1:n)[ic == 1]
tsp <- z$tcen[sp]
t1 <- z$wd[1:nw]
if(ctype == "left") {
B[1,] <- 1 - B[1,]
B[-1,] <- - B[-1,]
B <- B[,ncol(B):1]
}
dimnames(B) <- list(c("tau",dimnames(x)[[2]],"Qbar"),NULL)
a <- list(sol=B, Isplit=sp, tsp = tsp, status=ic, ctype = ctype)
fitted <- x %*% B[-c(1,dim(B)[1]),]
dimnames(fitted) <- list(NULL, paste("tau=",round(B[1,],4)))
a$fitted.values <- fitted
class(a) <- "crq"
return(a)
}
crq.fit.pen <- function(x, y, cen, weights=NULL,grid, ctype = "right" ){
p <- ncol(x)
n <- length(y)
if(missing(grid)) grid <- seq(1/n,1-1/n,by=min(0.01,1/(2*length(y)^.7)))
if(!is.numeric(grid)) stop("Invalid grid")
if(any(grid < 0) || any(grid > 1)) stop("Invalid grid")
m <- length(grid)
xbar <- apply(x,2,mean)
if(length(weights)){
if (any(weights < 0))
stop("negative weights not allowed")
contr <- attr(x, "contrasts")
x <- x * weights
y <- y * weights
}
if(ctype == "left") y <- -y
s <- rep(0,n)
u <- rep(1,n)
d <- rep(1,n)
r <- rep(1,p)
B <- matrix(0,p,m)
cc <- as.logical(cen)
y1 <- y[cc]
n1 <- length(y1)
x1 <- x[cc,]
z <- .Fortran("crqfnb", as.integer(n), as.integer(p),
a1 = as.double(t(as.matrix(x1))), c1 = as.double(-y1), n1=as.integer(n1),
as.double(x), as.double(y),as.double(cen),B =as.double(B),
g = as.double(grid),m = as.integer(m), as.double(r),
as.double(s), as.double(d), as.double(u),
wn = double(n1 * 9), wp = double((p + 3) * p),
info = integer(1), PACKAGE = "quantreg")
J <- z$m - 1
B <- matrix(-z$B, p, m)
B <- B[,1:J,drop = FALSE]
qhat <- t(xbar) %*% B
B <- rbind(grid[1:J],B,qhat)
dimnames(B) <- list(c("tau",dimnames(x)[[2]],"Qhat"),NULL)
if(ctype == "left") {
B[1,] <- 1 - B[1,]
B[-1,] <- - B[-1,]
B <- B[,ncol(B):1]
}
B <- list(sol=B, ctype = ctype)
class(B) <- "crq"
B
}
plot.summary.crqs <-
function (x, nrow = 3, ncol = 3, CoxPHit = NULL, ...) {
taus <- function(x) x$tau
xx <- unlist(lapply(x, taus))
coef <- lapply(x, coefficients)
p <- nrow(coef[[1]])
k <- ncol(coef[[1]])
if(k != 6) stop("summary.crqs object has wrong column dimension")
m <- length(xx)
blab <- dimnames(coef[[1]])[[1]]
a <- array(unlist(coef), c(p, k, m))
if(length(CoxPHit))
CoxQTE <- QTECox(CoxPHit)
oldpar <- par(no.readonly=TRUE)
par(mfrow = c(nrow, ncol))
for (i in 2:p) {
b <- a[i, 1, ]
bl <- a[i, 2, ]
bu <- a[i, 3, ]
plot(rep(xx, 2), c(bl, bu), xlab = "", ylab = "", type = "n")
title(paste(blab[i]), cex = 0.75)
polygon(c(xx, rev(xx)), c(bl, rev(bu)), col = "LightSkyBlue")
points(xx, b, cex = 0.5, pch = "o", col = "blue")
lines(xx, b, col = "blue")
abline(h = 0)
if(length(CoxPHit)) {
lines(CoxQTE$taus,CoxQTE$QTE[,i-1],col="red")
}
}
par(oldpar)
}
summary.crq <-
function (object, taus = 1:4/5, alpha = .05, se = "boot", covariance = TRUE, ...) {
mt <- terms(object)
m <- model.frame(object)
x <- model.matrix(mt, m, contrasts = object$contrasts)
Y <- model.response(m)
method <- object$method
ctype <- object$ctype
y <- Y[,1]
cen <- Y[,2]
wt <- as.vector(model.weights(object$model))
if (!is.null(wt)) {
resid <- resid * wt
x <- x * wt
y <- y * wt
cen <- cen * wt
}
if(method == "Powell"){
coef <- coefficients(object)
vnames <- dimnames(x)[[2]]
resid <- object$residuals
tau <- object$tau
n <- length(resid)
p <- length(coef)
rdf <- n - p
if(rdf < 1) stop("residual degrees of freedom nonpositive")
if (se == "boot") {
if(attr(Y,"type") == "left")
s <- cen < x %*% coef
else
s <- cen > x %*% coef
B <- boot.rq(x[s, ], y[s], tau, ...)
cov <- cov(B)
serr <- sqrt(diag(cov))
}
else stop("Only boot method is implemented for crq inference")
coef <- array(coef, c(p, 6))
fact <- qt(1 - alpha/2, rdf)
coef[, 2] <- coef[,2] - fact * serr
coef[, 3] <- coef[,3] + fact * serr
coef[, 4] <- serr
coef[, 5] <- coef[, 1]/coef[, 4]
coef[, 6] <- 2 * (1 - pt(abs(coef[, 5]),rdf))
cnames <- c("Value","Lower Bd","Upper Bd","Std Error","T Value","Pr(>|t|)")
dimnames(coef) <- list(vnames, cnames)
object <- object[c("call", "terms")]
if (covariance == TRUE)
object$cov <- cov
object$B <- B
object$coefficients <- coef
object$rdf <- rdf
object$tau <- tau
class(object) <- "summary.crq"
return(object)
}
else if(method == "Portnoy" || method == "Portnoy2" || method == "PengHuang") {
coef <- as.matrix(coef(object,taus))
coef <- coef[,apply(coef,2,function(x) any(!is.na(x))),drop = FALSE] # Delete NA columns if any
if(ctype == "right") taus <- taus[1:ncol(coef)]
else taus <- taus[(1 + length(taus)-ncol(coef)):length(taus)]
B <- boot.crq(x, y, cen, taus, method = method, ctype = ctype, ...)
bmethod <- B$bmethod
nas <- apply(is.na(B$A[1,,, drop = TRUE]),1,sum)
Brep <- dim(B$A)[3]
p <- dim(B$A)[1]
if(bmethod == "jack")
sqmn <- sqrt((B$n - B$mboot - p + 1)/B$mboot)
else
sqmn <- sqrt(B$mboot/B$n)
fact <- qnorm(1 - alpha/2)/qnorm(.75)
A <- apply(B$A, 1:2, quantile, probs = 1:3/4, na.rm = TRUE)
#D <- .5 * fact *(A[3,,]-A[1,,]) * sqmn
#L <- coef - D
#U <- coef + D
DU <- fact *(A[3,,]-A[2,,]) * sqmn
DL <- fact *(A[2,,]-A[1,,]) * sqmn
L <- coef - DL
U <- coef + DU
S <- (U - L)/(2 * qnorm(1 - alpha/2))
T <- coef/S
P <- 2 * (1 - pnorm(abs(T)))
G <- list()
cnames <- c("Value","Lower Bd","Upper Bd","Std Error","T Value","Pr(>|t|)")
for(i in 1:length(taus)){
tab <- cbind(coef[,i],L[,i],U[,i],S[,i],T[,i],P[,i])
dimnames(tab)[[2]] <- cnames
if(covariance)
cov <- var(t(B$A[,i,]), na.rm = TRUE) * sqmn^2
else
cov <- NULL
G[[i]] <- list(tau = taus[i], coefficients = tab, NAs = nas[i],
cov = cov, Brep = Brep, bmethod = bmethod)
}
class(G) <- "summary.crqs"
return(G)
}
else
stop("Invalid method for summary.crq")
}
"summary.crqs" <-
function (object, ...) {
if(!object$method == "Powell")
stop("Invalid method")
taus <- object$tau
xsum <- as.list(taus)
for(i in 1:length(taus)){
xi <- object
xi$coefficients <- xi$coefficients[,i]
xi$residuals <- xi$residuals[,i]
xi$tau <- xi$tau[i]
class(xi) <- "crq"
xsum[[i]] <- summary(xi, ...)
}
class(xsum) <- "summary.crqs"
xsum
}
print.crq <- function(x, ...){
if(!is.null(cl <- x$call)) {
cat("Call:\n")
dput(cl)
}
coef <- coef(x, ...)
cat("\nCoefficients:\n")
print.default(coef)
invisible(x)
}
print.summary.crqs <- function(x, ...)
lapply(x,print.summary.crq)
print.summary.crq <- function (x, digits = max(5, .Options$digits - 2), ...) {
coef <- x$coefficients
tau <- x$tau
if(length(x$NAs)) NAs <- x$NAs
else NAs <- 0
cat("\ntau: ")
print(format(round(tau, digits = digits)), quote = FALSE, ...)
if(NAs > 0)
cat(paste(" Number of NA Bootstrap Replications: ", NAs, "out of", x$Brep))
cat("\nCoefficients:\n")
print(format(round(coef, digits = digits)), quote = FALSE, ...)
invisible(x)
}
# SLP: Mar 31 2015 (Revised RWK: Aug 2015)
# An R version of crq.fit.por without any frills on a fixed grid.
crq.fit.por2 <- function(x, y, cen, weights = NULL, grid, ctype = "right") {
method <- "fn"
BIG <- 1e+6
ztol <- 1e-6
n <- length(y)
p <- ncol(x)
cen <- as.logical(cen)
uncen <- !cen
if(sum(uncen) < p)
stop("too few uncensored observations")
sol <- matrix(NA, p+2,1)
if (length(weights)) {
if (any(weights < 0))
stop("negative weights not allowed")
contr <- attr(x, "contrasts")
x <- weights * x
y <- weights * y
}
w <- rep(1,n)
if (ctype == "left") y <- -y
if(missing(grid)) { # Fixme: else if ...
g0 <- 1/(2*n)
dg <- 1/(5 + 3*n^.4)
grid <- seq(g0,1-g0,by=dg)
}
repeat { # Peel censored points below g0 fit
z <- rq.wfit(x, y, tau = g0, weights = w, method = method)
J <- which(abs(z$resid) < ztol & cen)
if(length(J)){
y <- y[-J]
x <- x[-J,]
w <- w[-J]
cen <- cen[-J]
}
else break
}
n <- length(y)
U <- which(cen)
if(length(cen) == 0) { # No censored points left!
z <- rq(y ~ x - 1, tau = grid, weights = w, method = method)
return(coef(z))
}
else { # Data Augmentation:
m <- sum(cen)
x <- rbind(x, x[cen, ])
y <- c(y, rep(BIG, m))
w[cen] <- 0
w <- c(w, rep(1, m))
}
taw <- rep(0, length(y))
K <- NULL
sol <- matrix(NA, p+2, length(grid))
b <- z$coef
xbar <- apply(x, 2, mean)
sol[,1] <- c(0, b, sum(xbar * b))
for(i in 2:length(grid)){
tau <- grid[i]
if(length(K)){ # Update w
w[K] <- (tau - taw[K])/(1 - taw[K])
idx <- which(which(cen) %in% K)[rank(K)]
w[n + idx] <- 1 - w[K]
}
z <- rq.wfit(x, y, tau = tau, weights = w, method = method)
b0 <- b
b <- z$coef
sol[,i] <- c(tau, b, sum(b * xbar))
if(any(abs(z$resid) < ztol & (y > BIG - 1))) break
J <- U[z$resid[U] < 0]
if(length(J)){
U <- setdiff(U, J)
K <- union(K, J)
y0 <- x[J,,drop = FALSE] %*% b0
y1 <- y[J] - z$resid[J]
gstep <- grid[i] - grid[i-1]
taw[J] <- tau - gstep * (y1 - y[J])/(y1 - y0)
}
}
if (ctype == "left") {
sol[1, ] <- 1 - sol[1, ]
sol[-1, ] <- -sol[-1, ]
}
dimnames(sol) <- list(c("tau", dimnames(x)[[2]], "Qbar"), NULL)
return(list(sol = sol, ctype = ctype))
}
