Revision e71c994d653adbe4d49f0fed4cb9d80ae22d8ce1 authored by Roger Koenker on 24 June 2012, 00:00:00 UTC, committed by Gabor Csardi on 24 June 2012, 00:00:00 UTC
1 parent ff94aa3
crq.R
``````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",
))
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
class(ss) <- "Surv"
ss
}

QTECox <- function(x, smooth = TRUE){
# compute quantile treatment effect for a Cox PengHuang fit
g <- 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 = "Bose", ...)
{
n <- length(y)
p <- ncol(x)
if(missing(mboot)) mboot <- n
A <- array(0,dim=c(p,length(taus),R))

for (i in 1:R){
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 == "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 )
cat(paste("bootstrap roughly ",100*(i/R)," percent complete\n"))
A[,,i] <- coef(a,taus)
}
list(A = A, n = length(y), mboot = mboot)
}

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
r <- S[1, ]
r <- c(r[1],r)
r <- (r[-1]+r[-length(r)])/2
r <- c(0,r)
if(is.unsorted(r)) r <- r[-length(r)] #kludge until why this happens is found
B <- S[-1,]
J <- length(r)
B <- t(cbind(B[,1],B))
ts <- taus[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)
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","PengHuang"), contrasts = NULL, ...)
{
require(survival)
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 == "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\$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","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,  ...),
"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 #!!! Fortran routine wants censoring indicator flipped.
mp <- n + 5 + max(1, sum(cen))
eps <- 1e-04
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"))
if(flag > 0 && flag != 5 && flag < 8)
ifelse(flag <= 3,stop(msg),warning(msg))
J <- z\$lsol
B <- matrix(z\$sol, nrow=p+2, ncol=nsol, byrow=FALSE)[,1:J]
dimnames(B) <- list(c("tau",dimnames(x)[[2]],"Qbar"),NULL)
ic <- z\$icen
sp <- (1:n)[ic == 1]
tsp <- z\$tcen[sp]
t1 <- z\$wd[1:nw]
if(ctype == "left") {
sol[1,] <- 1 - sol[1,]
sol[-1,] <- - sol[-1,]
}
a<-list(sol=B, Isplit=sp, tsp = tsp,status=ic)
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]
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  <- list(sol=B)
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)
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 (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, 4))
dimnames(coef) <- list(vnames, c("Value", "Std. Error", "t value",
"Pr(>|t|)"))
coef[, 2] <- serr
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- if(rdf > 0)
2 * (1 - pt(abs(coef[, 3]),rdf))
else NA
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 == "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
taus <- taus[1:ncol(coef)]
B <- boot.crq(x, y, cen, taus, method = method, ctype = ctype, ...)
sqmn <- sqrt(B\$mboot/B\$n)
fact <-   qnorm(1 - alpha/2)/qnorm(.75)
B <- apply(B\$A, 1:2, quantile, probs = 1:3/4, na.rm = TRUE)
D <- .5 * fact *(B[3,,]-B[1,,]) * sqmn
L <- coef - D
U <- coef + D
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
G[[i]] <- list(tau = taus[i], coefficients = tab)
}
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
cat("\ntau: ")
print(format(round(tau, digits = digits)), quote = FALSE, ...)
cat("\nCoefficients:\n")
print(format(round(coef, digits = digits)), quote = FALSE, ...)
invisible(x)

}
``````

Computing file changes ...