bandwidth.rq <- function(p, n, hs = TRUE, alpha = 0.05)
{
# Bandwidth selection for sparsity estimation two flavors:
# Hall and Sheather(1988, JRSS(B)) rate = O(n^{-1/3})
# Bofinger (1975, Aus. J. Stat) -- rate = O(n^{-1/5})
# Generally speaking, default method, hs=TRUE is preferred.
x0 <- qnorm(p)
f0 <- dnorm(x0)
if(hs)
n^(-1/3) * qnorm(1 - alpha/2)^(2/3) *
((1.5 * f0^2)/(2 * x0^2 + 1))^(1/3)
else n^-0.2 * ((4.5 * f0^4)/(2 * x0^2 + 1)^2)^ 0.2
}
plot.rq.process <- function(x, nrow = 3, ncol = 2, ...)
{
## Function to plot estimated quantile regression process
tdim <- dim(x$sol)
p <- tdim[1] - 3
m <- tdim[2]
oldpar <- par(no.readonly=TRUE)
par(mfrow = c(nrow, ncol))
ylab <- dimnames(x$sol)[[1]]
for(i in 1:p) {
plot(x$sol[1,], x$sol[3 + i, ], xlab = "tau",
ylab = ylab[3 + i], type = "l")
}
par(oldpar)
}
print.rqs <- function (x, ...)
{
if (!is.null(cl <- x$call)) {
cat("Call:\n")
dput(cl)
}
coef <- coef(x)
cat("\nCoefficients:\n")
print(coef, ...)
rank <- x$rank
nobs <- nrow(residuals(x))
p <- nrow(coef)
rdf <- nobs - p
cat("\nDegrees of freedom:", nobs, "total;", rdf, "residual\n")
if (!is.null(attr(x, "na.message")))
cat(attr(x, "na.message"), "\n")
invisible(x)
}
"print.rq" <-
function(x, ...)
{
if(!is.null(cl <- x$call)) {
cat("Call:\n")
dput(cl)
}
coef <- coef(x)
cat("\nCoefficients:\n")
print(coef, ...)
rank <- x$rank
nobs <- length(residuals(x))
if(is.matrix(coef))
p <- dim(coef)[1]
else p <- length(coef)
rdf <- nobs - p
cat("\nDegrees of freedom:", nobs, "total;", rdf, "residual\n")
if(!is.null(attr(x, "na.message")))
cat(attr(x, "na.message"), "\n")
invisible(x)
}
print.summary.rqs <- function(x, ...) {
lapply(x, print.summary.rq)
invisible(x)
}
"print.summary.rq" <-
function(x, digits = max(5, .Options$digits - 2), ...)
{
cat("\nCall: ")
dput(x$call)
coef <- x$coef
## df <- x$df
## rdf <- x$rdf
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)
}
"rq" <-
function (formula, tau = 0.5, data, subset, weights, na.action, method = "br",
model = TRUE, contrasts = NULL, ...)
{
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")
weights <- as.vector(model.weights(mf))
Y <- model.response(mf)
X <- model.matrix(mt, mf, contrasts)
eps <- .Machine$double.eps^(2/3)
Rho <- function(u,tau) u * (tau - (u < 0))
if(length(tau)>1){
if(any(tau < 0) || any(tau > 1))
stop("invalid tau: taus should be >= 0 and <= 1")
if(any(tau == 0)) tau[tau == 0] <- eps
if(any(tau == 1)) tau[tau == 1] <- 1 -eps
coef <- matrix(0,ncol(X),length(tau))
rho <- rep(0,length(tau))
fitted <- resid <- matrix(0,nrow(X),length(tau))
for(i in 1:length(tau)){
z <- {if (length(weights))
rq.wfit(X, Y, tau = tau[i], weights, method, ...)
else rq.fit(X, Y, tau = tau[i], method, ...)
}
coef[,i] <- z$coefficients
resid[,i] <- z$residuals
rho[i] <- sum(Rho(z$residuals,tau[i]))
fitted[,i] <- Y - z$residuals
}
taulabs <- paste("tau=",format(round(tau,3)))
dimnames(coef) <- list(dimnames(X)[[2]],taulabs)
dimnames(resid) <- list(dimnames(X)[[1]],taulabs)
fit <- z
fit$coefficients <- coef
fit$residuals <- resid
fit$fitted.values <- fitted
if(method == "lasso") class(fit) <- c("lassorqs","rqs")
else if(method == "scad") class(fit) <- c("scadrqs","rqs")
else class(fit) <- "rqs"
}
else{
process <- (tau < 0 || tau > 1)
if(tau == 0) tau <- eps
if(tau == 1) tau <- 1 -eps
fit <- {
if (length(weights))
rq.wfit(X, Y, tau = tau, weights, method, ...)
else rq.fit(X, Y, tau = tau, method, ...)
}
if(process)
rho <- list(x = fit$sol[1,],y = fit$sol[3,])
else {
dimnames(fit$residuals) <- list(dimnames(X)[[1]],NULL)
rho <- sum(Rho(fit$residuals,tau))
}
if(method == "lasso") class(fit) <- c("lassorq","rq")
else if(method == "scad") class(fit) <- c("scadrq","rq")
else class(fit) <- ifelse(process, "rq.process", "rq")
}
fit$na.action <- attr(mf, "na.action")
fit$formula <- formula
fit$terms <- mt
fit$xlevels <- .getXlevels(mt,mf)
fit$call <- call
fit$tau <- tau
fit$weights <- weights
fit$residuals <- drop(fit$residuals)
fit$rho <- rho
fit$method <- method
fit$fitted.values <- drop(fit$fitted.values)
attr(fit, "na.message") <- attr(m, "na.message")
if(model) fit$model <- mf
fit
}
"rq.fit" <-
function(x, y, tau = 0.5, method = "br", ...)
{
if(length(tau) > 1) {
tau <- tau[1]
warning("Multiple taus not allowed in rq.fit: solution restricted to first element")
}
fit <- switch(method,
fn = rq.fit.fnb(x, y, tau = tau, ...),
fnb = rq.fit.fnb(x, y, tau = tau, ...),
fnc = rq.fit.fnc(x, y, tau = tau, ...),
pfn = rq.fit.pfn(x, y, tau = tau, ...),
br = rq.fit.br(x, y, tau = tau, ...),
lasso = rq.fit.lasso(x, y, tau = tau, ...),
scad = rq.fit.scad(x, y, tau = tau, ...),
{
what <- paste("rq.fit.", method, sep = "")
if(exists(what, mode = "function"))
(get(what, mode = "function"))(x, y, ...)
else stop(paste("unimplemented method:", method))
}
)
fit$fitted.values <- y - fit$residuals
fit$contrasts <- attr(x, "contrasts")
fit
}
"rqs.fit"<-
function(x, y, tau = 0.5, tol = 0.0001)
{
# function to compute rq fits for multiple y's
x <- as.matrix(x)
p <- ncol(x)
n <- nrow(x)
m <- ncol(y)
z <- .Fortran("rqs",
as.integer(n),
as.integer(p),
as.integer(m),
as.integer(n + 5),
as.integer(p + 2),
as.double(x),
as.double(y),
as.double(tau),
as.double(tol),
flag = integer(m),
coef = double(p * m),
resid = double(n),
integer(n),
double((n + 5) * (p + 2)),
double(n),
PACKAGE = "quantreg")
if(sum(z$flag)>0){
if(any(z$flag)==2)
warning(paste(sum(z$flag==2),"out of",m,
"BS replications have near singular design"))
if(any(z$flag)==1)
warning(paste(sum(z$flag==1),"out of",m,"may be nonunique"))
}
return(t(matrix(z$coef, p, m)))
}
"formula.rq" <-
function (x, ...)
{
form <- x$formula
if (!is.null(form)) {
form <- formula(x$terms)
environment(form) <- environment(x$formula)
form
}
else formula(x$terms)
}
"predict.rq" <-
function (object, newdata, type = "none", interval = c("none",
"confidence"), level = 0.95, na.action = na.pass, ...)
{
if (missing(newdata))
return(object$fitted)
else {
tt <- terms(object)
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.action,
xlev = object$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses")))
.checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
}
pred <- drop(X %*% object$coefficients)
dots <- list(...)
if (length(dots$se))
boot <- (dots$se == "boot")
else boot <- FALSE
if (length(dots$mofn))
mofn <- dots$mofn
interval <- match.arg(interval)
if (!interval == "none") {
if (interval == "confidence") {
if (type == "percentile") {
if (boot) {
if(exists("mofn")) {# Rescale and recenter!!
n <- length(object$fitted)
factor <- ifelse(mofn < n, sqrt(mofn/n), 1)
XB <- X %*% t(summary(object, cov = TRUE, ...)$B)/factor
pl <- apply(XB, 1, function(x) quantile(x, (1 - level)/2))
pu <- apply(XB, 1, function(x) quantile(x, 1 - (1 - level)/2))
pl <- pred + factor * (pl - pred)
pu <- pred + factor * (pu - pred)
}
else {
XB <- X %*% t(summary(object, cov = TRUE, ...)$B)
pl <- apply(XB, 1, function(x) quantile(x, (1 - level)/2))
pu <- apply(XB, 1, function(x) quantile(x, 1 - (1 - level)/2))
}
pred <- cbind(pred, pl, pu)
colnames(pred) <- c("fit", "lower", "higher")
}
else stop("Percentile method requires se = \"boot\".")
}
else if (type == "direct") {
if (boot)
stop("Direct method incompatible with bootstrap covariance matrix estimation")
Z <- rq.fit(object$x, object$y, tau = -1)$sol
V <- summary(object, cov = TRUE, ...)
df <- V$rdf
tfrac <- qt(1 - (1 - level)/2, df)
Vun <- V$cov * V$scale^2
tau <- object$tau
bn <- tfrac * sqrt(diag(X %*% Vun %*% t(X)))
tauU <- pmin(tau + bn, 1 - 1/df)
tauL <- pmax(tau - bn, 1/df)
tauhat <- Z[1, ]
yhat <- X %*% Z[-(1:3), ]
n <- nrow(X)
pl <- yhat[cbind(1:n, cut(tauL, tauhat, label = FALSE))]
pu <- yhat[cbind(1:n, cut(tauU, tauhat, label = FALSE))]
pred <- cbind(pred, pl, pu)
colnames(pred) <- c("fit", "lower", "higher")
}
else {
V <- summary(object, cov = TRUE, ...)
df <- V$rdf
tfrac <- qt((1 - level)/2, df)
sdpred <- sqrt(diag(X %*% V$cov %*% t(X)))
pred <- cbind(pred, pred + tfrac * sdpred %o%
c(1, -1))
colnames(pred) <- c("fit", "lower", "higher")
}
}
else stop(paste("No interval method for", interval))
}
pred
}
"predict.rqs" <-
function (object, newdata, type = "Qhat", stepfun = FALSE, na.action = na.pass, ...)
{
## with all defaults
if(missing(newdata) & !stepfun) return(object$fitted)
## otherwise
tt <- delete.response(terms(object))
m <- if(missing(newdata)) model.frame(object) else model.frame(tt, newdata,
na.action = na.action, xlev = object$xlevels)
if(!is.null(cl <- attr(tt, "dataClasses"))) .checkMFClasses(cl, m)
X <- model.matrix(tt, m, contrasts = object$contrasts)
pred <- X %*% object$coefficients
## return stepfun or matrix
if(stepfun) {
pred <- t(pred)
taus <- object$tau
if(type == "Qhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) f[[i]] <- stepfun(taus, c(pred[1,i], pred[,i]))
} else
f <- stepfun(taus, c(pred[1,1], pred[,1]))
}
else if(type == "Fhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) {
o <- order(pred[,i])
f[[i]] <- stepfun(pred[o,i],c(taus[1],taus[o]))
}
} else
f <- stepfun(pred[,1],c(taus[1],taus))
}
else stop("type must be either 'Qhat' or 'Fhat'")
return(f)
} else {
return(pred)
}
}
"predict.rq.process" <-
function (object, newdata, type = "Qhat", stepfun = FALSE, na.action = na.pass, ...)
{
if (missing(newdata))
return(object$fitted)
else {
tt <- terms(object)
Terms <- delete.response(tt)
m <- model.frame(Terms, newdata, na.action = na.action,
xlev = object$xlevels)
if (!is.null(cl <- attr(Terms, "dataClasses")))
.checkMFClasses(cl, m)
X <- model.matrix(Terms, m, contrasts = object$contrasts)
}
pred <- t(X %*% object$sol[-(1:3),])
if(stepfun){
taus <- object$sol[1,]
if(type == "Qhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) f[[i]] <- stepfun(taus, c(pred[1,i], pred[,i]))
} else
f <- stepfun(taus, c(pred[1,1], pred[,1]))
}
else if(type == "Fhat"){
if(ncol(pred) > 1) {
f <- as.list(1:ncol(pred))
for(i in 1:ncol(pred)) {
o <- order(pred[,i])
f[[i]] <- stepfun(pred[,i],c(taus[1],taus))
}
} else
f <- stepfun(pred[,1],c(taus[1],taus))
}
else stop("type must be either 'Qhat' or 'Fhat'")
return(f)
}
return(pred)
}
"rearrange" <- function (f, xmin, xmax)
# Revised Version September 11 2007.
{
if (is.list(f))
lapply(f, rearrange)
else {
if (!is.stepfun(f))
stop("Only stepfuns can be rearranged.\n")
call <- attributes(f)$call;
right <- call[match("right",names(call))]=="TRUE()"
x <- knots(f)
n <- length(x)
if(missing(xmin)) xmin = x[1]
if(missing(xmax)) xmax = x[n]
x <- x[(x >= xmin) & (x <= xmax)]
x <- c(xmin, x, xmax)
n <- length(x)
y <- f(x)
o <- ifelse(rep(right,n-1), order(y[-1])+1, order(y[-n]))
x <- cumsum(c(x[1], diff(x)[o - right]))
y <- y[o]
y <- c(y[1], y, max(y))
stepfun(x, y, right = right)
}
}
# Function to compute regression quantiles using original simplex approach
# of Barrodale-Roberts/Koenker-d'Orey. There are several options.
# The options are somewhat different than those available for the Frisch-
# Newton version of the algorithm, reflecting the different natures of the
# problems typically solved. Succintly BR for "small" problems, FN for
# "large" ones. Obviously, these terms are conditioned by available hardware.
#
# Basically there are two modes of use:
# 1. For Single Quantiles:
#
# if tau is between 0 and 1 then only one quantile solution is computed.
#
# if ci = FALSE then just the point estimate and residuals are returned
# If the column dimension of x is 1 then ci is set to FALSE since
# since the rank inversion method has no proper null model.
# if ci = TRUE then there are two options for confidence intervals:
#
# 1. if iid = TRUE we get the original version of the rank
# inversion intervals as in Koenker (1994)
# 2. if iid = FALSE we get the new version of the rank inversion
# intervals which accounts for heterogeneity across
# observations in the conditional density of the response.
# The theory of this is described in Koenker-Machado(1999)
# Both approaches involve solving a parametric linear programming
# problem, the difference is only in the factor qn which
# determines how far the PP goes. In either case one can
# specify two other options:
# 1. interp = FALSE returns two intervals an upper and a
# lower corresponding to a level slightly
# above and slightly below the one specified
# by the parameter alpha and dictated by the
# essential discreteness in the test statistic.
# interp = TRUE returns a single interval based on
# linear interpolation of the two intervals
# returned: c.values and p.values which give
# the critical values and p.values of the
# upper and lower intervals. Default: interp = TRUE.
# 2. tcrit = TRUE uses Student t critical values while
# tcrit = FALSE uses normal theory ones.
# 2. For Multiple Quantiles:
#
# if tau < 0 or tau >1 then it is presumed that the user wants to find
# all of the rq solutions in tau, and the program computes the whole
# quantile regression solution as a process in tau, the resulting arrays
# containing the primal and dual solutions, betahat(tau), ahat(tau)
# are called sol and dsol. These arrays aren't printed by the default
# print function but they are available as attributes.
# It should be emphasized that this form of the solution can be
# both memory and cpu quite intensive. On typical machines it is
# not recommended for problems with n > 10,000.
# In large problems a grid of solutions is probably sufficient.
#
rq.fit.br <-
function (x, y, tau = 0.5, alpha = 0.1, ci = FALSE, iid = TRUE,
interp = TRUE, tcrit = TRUE)
{
tol <- .Machine$double.eps^(2/3)
eps <- tol
big <- .Machine$double.xmax
x <- as.matrix(x)
p <- ncol(x)
n <- nrow(x)
ny <- NCOL(y)
nsol <- 2
ndsol <- 2
# Check for Singularity of X since br fortran isn't very reliable about this
if (qr(x)$rank < p)
stop("Singular design matrix")
if (tau < 0 || tau > 1) {
nsol <- 3 * n
ndsol <- 3 * n
lci1 <- FALSE
qn <- rep(0, p)
cutoff <- 0
tau <- -1
}
else {
if (p == 1)
ci <- FALSE
if (ci) {
lci1 <- TRUE
if (tcrit)
cutoff <- qt(1 - alpha/2, n - p)
else cutoff <- qnorm(1 - alpha/2)
if (!iid) {
h <- bandwidth.rq(tau, n, hs = TRUE)
bhi <- rq.fit.br(x, y, tau + h, ci = FALSE)
bhi <- coefficients(bhi)
blo <- rq.fit.br(x, y, tau - h, ci = FALSE)
blo <- coefficients(blo)
dyhat <- x %*% (bhi - blo)
if (any(dyhat <= 0)) {
pfis <- (100 * sum(dyhat <= 0))/n
warning(paste(pfis, "percent fis <=0"))
}
f <- pmax(eps, (2 * h)/(dyhat - eps))
qn <- rep(0, p)
for (j in 1:p) {
qnj <- lm(x[, j] ~ x[, -j] - 1, weights = f)$resid
qn[j] <- sum(qnj * qnj)
}
}
else qn <- 1/diag(solve(crossprod(x)))
}
else {
lci1 <- FALSE
qn <- rep(0, p)
cutoff <- 0
}
}
z <- .Fortran("rqbr", as.integer(n), as.integer(p), as.integer(n +
5), as.integer(p + 3), as.integer(p + 4), as.double(x),
as.double(y), as.double(tau), as.double(tol), flag = as.integer(1),
coef = double(p), resid = double(n), integer(n), double((n +
5) * (p + 4)), double(n), as.integer(nsol), as.integer(ndsol),
sol = double((p + 3) * nsol), dsol = double(n * ndsol),
lsol = as.integer(0), h = integer(p * nsol), qn = as.double(qn),
cutoff = as.double(cutoff), ci = double(4 * p), tnmat = double(4 *
p), as.double(big), as.logical(lci1), PACKAGE = "quantreg")
if (z$flag != 0)
warning(switch(z$flag, "Solution may be nonunique", "Premature end - possible conditioning problem in x"))
if (tau < 0 || tau > 1) {
sol <- matrix(z$sol[1:((p + 3) * z$lsol)], p + 3)
dsol <- matrix(z$dsol[1:(n * z$lsol)], n)
vnames <- dimnames(x)[[2]]
dimnames(sol) <- list(c("tau", "Qbar", "Obj.Fun", vnames),
NULL)
return(list(sol = sol, dsol = dsol))
}
if (!ci) {
coef <- z$coef
dual <- z$dsol[1:n]
names(coef) <- dimnames(x)[[2]]
return(list(coefficients = coef, x = x, y = y, residuals = y - x %*% z$coef,
dual = dual))
}
if (interp) {
Tn <- matrix(z$tnmat, nrow = 4)
Tci <- matrix(z$ci, nrow = 4)
Tci[3, ] <- Tci[3, ] + (abs(Tci[4, ] - Tci[3, ]) * (cutoff -
abs(Tn[3, ])))/abs(Tn[4, ] - Tn[3, ])
Tci[2, ] <- Tci[2, ] - (abs(Tci[1, ] - Tci[2, ]) * (cutoff -
abs(Tn[2, ])))/abs(Tn[1, ] - Tn[2, ])
Tci[2, ][is.na(Tci[2, ])] <- -big
Tci[3, ][is.na(Tci[3, ])] <- big
coefficients <- cbind(z$coef, t(Tci[2:3, ]))
vnames <- dimnames(x)[[2]]
cnames <- c("coefficients", "lower bd", "upper bd")
dimnames(coefficients) <- list(vnames, cnames)
residuals <- y - drop(x %*% z$coef)
return(list(coefficients = coefficients, residuals = residuals))
}
else {
Tci <- matrix(z$ci, nrow = 4)
coefficients <- cbind(z$coef, t(Tci))
residuals <- y - drop(x %*% z$coef)
vnames <- dimnames(x)[[2]]
cnames <- c("coefficients", "lower bound", "Lower Bound",
"upper bd", "Upper Bound")
dimnames(coefficients) <- list(vnames, cnames)
c.values <- t(matrix(z$tnmat, nrow = 4))
c.values <- c.values[, 4:1]
dimnames(c.values) <- list(vnames, cnames[-1])
p.values <- if (tcrit)
matrix(pt(c.values, n - p), ncol = 4)
else matrix(pnorm(c.values), ncol = 4)
dimnames(p.values) <- list(vnames, cnames[-1])
list(coefficients = coefficients, residuals = residuals,
c.values = c.values, p.values = p.values)
}
}
"rq.fit.fnb" <-
function (x, y, tau = 0.5, beta = 0.99995, eps = 1e-06)
{
n <- length(y)
p <- ncol(x)
if (n != nrow(x))
stop("x and y don't match n")
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method. Set tau in (0,1)")
rhs <- (1 - tau) * apply(x, 2, sum)
d <- rep(1,n)
u <- rep(1,n)
wn <- rep(0,10*n)
wn[1:n] <- (1-tau) #initial value of dual solution
z <- .Fortran("rqfnb", as.integer(n), as.integer(p), a = as.double(t(as.matrix(x))),
c = as.double(-y), rhs = as.double(rhs), d = as.double(d),as.double(u),
beta = as.double(beta), eps = as.double(eps),
wn = as.double(wn), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1),PACKAGE= "quantreg")
if (z$info != 0)
stop(paste("Error info = ", z$info, "in stepy: singular design"))
coefficients <- -z$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
list(coefficients=coefficients, tau=tau, residuals=residuals)
}
"rq.fit.fnc" <-
function (x, y, R, r, tau = 0.5, beta = 0.9995, eps = 1e-06)
{
n1 <- length(y)
n2 <- length(r)
p <- ncol(x)
if (n1 != nrow(x))
stop("x and y don't match n1")
if (n2 != nrow(R))
stop("R and r don't match n2")
if (p != ncol(R))
stop("R and x don't match p")
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method. Set tau in (0,1)")
rhs <- (1 - tau) * apply(x, 2, sum)
u <- rep(1, max(n1,n2)) #upper bound vector and scratch vector
wn1 <- rep(0, 9 * n1)
wn1[1:n1] <- (1 - tau) #store the values of x1
wn2 <- rep(0, 6 * n2)
wn2[1:n2] <- 1 #store the values of x2
z <- .Fortran("rqfnc", as.integer(n1), as.integer(n2), as.integer(p),
a1 = as.double(t(as.matrix(x))), c1 = as.double(-y),
a2 = as.double(t(as.matrix(R))), c2 = as.double(-r),
rhs = as.double(rhs), d1 = double(n1), d2 = double(n2),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn1 = as.double(wn1), wn2 = as.double(wn2), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
if (z$info != 0)
stop(paste("Error info = ", z$info, "in stepy2: singular design"))
coefficients <- -z$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
it.count <- z$it.count
list(coefficients=coefficients, tau=tau, residuals=residuals, it = it.count)
}
"rq.fit.scad" <-
function (x, y, tau = 0.5, alpha = 3.2, lambda = 1, start = "rq", beta = 0.9995, eps = 1e-06)
{
n <- length(y)
p <- ncol(x)
if (n != nrow(x))
stop("x and y don't match n")
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method. Set tau in (0,1)")
if(length(lambda) == 1)
lambda <- c(0,rep(lambda,p-1))
if(length(lambda) != p)
stop(paste("lambda must be either of length ",p," or length one"))
if(any(lambda < 0))
stop("negative lambdas disallowed")
R <- diag(lambda,nrow = length(lambda))
R <- R[which(lambda != 0),, drop = FALSE]
r <- rep(0,nrow(R))
X <- rbind(x, R)
Y <- c(y, r)
N <- length(Y)
rhs <- (1 - tau) * apply(x, 2, sum) + apply(R,2,sum)
dscad <- function(x, a = 3.7, lambda = 2){
lambda * sign(x) * (abs(x) <= lambda) +
sign(x) * (a * lambda - abs(x)) / (a - 1) *
(abs(x) <= a * lambda) * (abs(x) > lambda)
}
binit <- switch(start,
rq = rq.fit.fnb(x, y, tau = tau)$coef[-1],
lasso = rq.fit.lasso(x, y, tau = tau, lambda = lambda)$coef[-1]
)
coef <- rep(.Machine$double.xmax,p)
vscad <- rhs - c(0,dscad(binit) * sign(binit))
it <- 0
while(sum(abs(binit - coef[-1])) > eps){
it <- it + 1
d <- rep(1, N)
u <- rep(1, N)
wn <- rep(0, 10 * N)
wn[1:N] <- c(rep((1 - tau),n),rep(.5,nrow(R)))
vrhs <- rhs - vscad
binit <- coef[-1]
z <- .Fortran("rqfnb", as.integer(N), as.integer(p), a = as.double(t(as.matrix(X))),
c = as.double(-Y), vrhs = as.double(vrhs), d = as.double(d),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn = as.double(wn), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
coef <- -z$wp[1:p]
vscad <- c(0,dscad(coef[2:p]) * sign(coef[2:p]))
}
if (z$info != 0)
stop(paste("Error info = ", z$info, "in stepy2: singular design"))
coefficients <- -z$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
it.count <- z$it.count
list(coefficients=coefficients, residuals=residuals, tau = tau,
lambda = lambda, it = it.count)
}
"rq.fit.lasso" <-
function (x, y, tau = 0.5, lambda = 1, beta = 0.9995, eps = 1e-06)
{
n <- length(y)
p <- ncol(x)
if (n != nrow(x))
stop("x and y don't match n")
if(length(lambda) == 1)
lambda <- c(0,rep(lambda,p-1))
if(length(lambda) != p)
stop(paste("lambda must be either of length ",p," or length one"))
if(any(lambda < 0))
stop("negative lambdas disallowed")
R <- diag(lambda,nrow = length(lambda))
R <- R[which(lambda != 0),, drop = FALSE]
r <- rep(0,nrow(R))
if (tau < eps || tau > 1 - eps)
stop("No parametric Frisch-Newton method. Set tau in (0,1)")
X <- rbind(x, R)
Y <- c(y, r)
N <- length(Y)
rhs <- (1 - tau) * apply(x, 2, sum) + 0.5 * apply(R,2,sum)
d <- rep(1, N)
u <- rep(1, N)
wn <- rep(0, 10 * N)
wn[1:N] <- rep(0.5, N)
z <- .Fortran("rqfnb", as.integer(N), as.integer(p), a = as.double(t(as.matrix(X))),
c = as.double(-Y), rhs = as.double(rhs), d = as.double(d),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn = as.double(wn), wp = double((p + 3) * p), aa = double(p *
p), it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
if (z$info != 0)
stop(paste("Error info = ", z$info, "in stepy2: singular design"))
coefficients <- -z$wp[1:p]
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% coefficients
it.count <- z$it.count
list(coefficients=coefficients, residuals=residuals, tau = tau,
lambda = lambda, it = it.count)
}
"rq.fit.pfn" <-
# This is an implementation (purely in R) of the preprocessing phase
# of the rq algorithm described in Portnoy and Koenker, Statistical
# Science, (1997) 279-300. In this implementation it can be used
# as an alternative method for rq() by specifying method="pfn"
# It should probably be used only on very large problems and then
# only with some caution. Very large in this context means roughly
# n > 100,000. The options are described in the paper, and more
# explicitly in the code. Again, it would be nice perhaps to have
# this recoded in a lower level language, but in fact this doesn't
# seem to make a huge difference in this case since most of the work
# is already done in the rq.fit.fnb calls.
#
function(x, y, tau = 0.5, Mm.factor = 0.8,
max.bad.fixup = 3, eps = 1e-6)
{
#rq function for n large --
n <- length(y)
if(nrow(x) != n)
stop("x and y don't match n")
if(tau < 0 | tau > 1)
stop("tau outside (0,1)")
p <- ncol(x)
m <- round(((p + 1) * n)^(2/3))
not.optimal <- TRUE
while(not.optimal) {
if(m < n)
s <- sample(n, m)
else {
b <- rq.fit.fnb(x, y, tau = tau, eps = eps)$coef
break
}
xx <- x[s, ]
yy <- y[s]
z <- rq.fit.fnb(xx, yy, tau = tau, eps = eps)
xxinv <- solve(chol(crossprod(xx)))
band <- sqrt(((x %*% xxinv)^2) %*% rep(1, p))
#sqrt(h<-ii)
r <- y - x %*% z$coef
M <- Mm.factor * m
lo.q <- max(1/n, tau - M/(2 * n))
hi.q <- min(tau + M/(2 * n), (n - 1)/n)
kappa <- quantile(r/pmax(eps, band), c(lo.q, hi.q))
sl <- r < band * kappa[1]
su <- r > band * kappa[2]
bad.fixup <- 0
while(not.optimal & (bad.fixup < max.bad.fixup)) {
xx <- x[!su & !sl, ]
yy <- y[!su & !sl]
if(any(sl)) {
glob.x <- c(t(x[sl, , drop = FALSE]) %*% rep(
1, sum(sl)))
glob.y <- sum(y[sl])
xx <- rbind(xx, glob.x)
yy <- c(yy, glob.y)
}
if(any(su)) {
ghib.x <- c(t(x[su, , drop = FALSE]) %*% rep(
1, sum(su)))
ghib.y <- sum(y[su])
xx <- rbind(xx, ghib.x)
yy <- c(yy, ghib.y)
}
z <- rq.fit.fnb(xx, yy, tau = tau, eps = eps)
b <- z$coef
r <- y - x %*% b
su.bad <- (r < 0) & su
sl.bad <- (r > 0) & sl
if(any(c(su.bad, sl.bad))) {
if(sum(su.bad | sl.bad) > 0.10000000000000001 *
M) {
warning("Too many fixups: doubling m")
m <- 2 * m
break
}
su <- su & !su.bad
sl <- sl & !sl.bad
bad.fixup <- bad.fixup + 1
}
else not.optimal <- FALSE
}
}
coefficients <- b
names(coefficients) <- dimnames(x)[[2]]
residuals <- y - x %*% b
return(list(coefficients=coefficients, tau=tau,
residuals=residuals))
}
"rq.wfit" <-
function(x, y, tau = 0.5, weights, method = "br", ...)
{
if(any(weights < 0))
stop("negative weights not allowed")
if(length(tau) > 1) {
tau <- tau[1]
warning("Multiple taus not allowed in rq.wfit: solution restricted to first element")
}
contr <- attr(x, "contrasts")
wx <- x * weights
wy <- y * weights
fit <- switch(method,
fn = rq.fit.fnb(wx, wy, tau = tau, ...),
fnb = rq.fit.fnb(wx, wy, tau = tau, ...),
br = rq.fit.br(wx, wy, tau = tau, ...),
fnc = rq.fit.fnc(wx, wy, tau = tau, ...),
pfn = rq.fit.pfn(wx, wy, tau = tau, ...), {
what <- paste("rq.fit.", method, sep = "")
if(exists(what, mode = "function"))
(get(what, mode = "function"))(x, y, ...)
else stop(paste("unimplemented method:", method))
}
)
if(length(fit$sol))
fit$fitted.values <- x %*% fit$sol[-(1:3),]
else
fit$fitted.values <- x %*% fit$coef
fit$residuals <- y - fit$fitted.values
fit$contrasts <- attr(x, "contrasts")
fit$weights <- weights
fit
}
"summary.rqs" <-
function (object, ...) {
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) <- "rq"
xsum[[i]] <- summary(xi, ...)
if(class(object)[1] == "dynrqs"){
class(xsum[[1]]) <- c("summary.dynrq", "summary.rq")
if(i == 1) xsum[[1]]$model <- object$model
}
}
class(xsum) <- "summary.rqs"
if(class(object)[1] == "dynrqs")
class(xsum) <- c("summary.dynrqs", "summary.rqs")
xsum
}
"logLik.rq" <- function(object, ...){
n <- length(object$residuals)
p <- length(object$coefficients)
pen <- (length(object$lambda) > 0)
tau <- object$tau
fid <- object$rho
val <- n * (log(tau * (1-tau)) - 1 - log(fid/n))
attr(val,"n") <- n
if(pen){
if(!hasArg(edfThresh)) edfThresh <- 0.0001
attr(val,"df") <- sum(abs(object$coefficients) > edfThresh)
}
else attr(val,"df") <- p
class(val) <- "logLik"
val
}
"logLik.rqs" <- function(object, ...){
n <- nrow(object$residuals)
p <- nrow(object$coefficients)
pen <- (length(object$lambda) > 0)
tau <- object$tau
fid <- object$rho
val <- n * (log(tau * (1-tau)) - 1 - log(fid/n))
attr(val,"n") <- n
if(pen){
if(!hasArg(edfThresh)) edfThresh <- 0.0001
attr(val,"df") <- sum(abs(object$coefficients) > edfThresh)
}
else attr(val,"df") <- p
class(val) <- "logLik"
val
}
"AIC.rq" <- function(object, ... , k = 2){
v <- logLik(object)
if(k <= 0)
k <- log(attr(v,"n"))
val <- AIC(v, k = k)
attr(val,"edf") <- attr(v,"df")
val
}
"extractAIC.rq" <- function(fit, scale, k=2, ...){
aic <- AIC(fit,k)
edf <- attr(aic, "edf")
c(edf, aic)
}
"AIC.rqs" <- function(object, ... , k = 2){
v <- logLik(object)
if(k <= 0)
k <- log(attr(v,"n"))
val <- AIC(v, k = k)
attr(val,"edf") <- attr(v,"df")
val
}
"summary.rq" <-
# This function provides methods for summarizing the output of the
# rq command. In this instance, "summarizing" means essentially provision
# of either standard errors, or confidence intervals for the rq coefficents.
# Since the preferred method for confidence intervals is currently the
# rank inversion method available directly from rq() by setting ci=TRUE, with br=TRUE.
# these summary methods are intended primarily for comparison purposes
# and for use on large problems where the parametric programming methods
# of rank inversion are prohibitively memory/time consuming. Eventually
# iterative versions of rank inversion should be developed that would
# employ the Frisch-Newton approach.
#
# Object is the result of a call to rq(), and the function returns a
# table of coefficients, standard errors, "t-statistics", and p-values, and, if
# covariance=TRUE a structure describing the covariance matrix of the coefficients,
# i.e. the components of the Huber sandwich.
#
# There are five options for "se":
#
# 1. "rank" strictly speaking this doesn't produce a "standard error"
# at all instead it produces a coefficient table with confidence
# intervals for the coefficients based on inversion of the
# rank test described in GJKP and Koenker (1994).
# 2. "iid" which presumes that the errors are iid and computes
# an estimate of the asymptotic covariance matrix as in KB(1978).
# 3. "nid" which presumes local (in tau) linearity (in x) of the
# the conditional quantile functions and computes a Huber
# sandwich estimate using a local estimate of the sparsity.
# 4. "ker" which uses a kernel estimate of the sandwich as proposed
# by Powell.
# 5. "boot" which uses a bootstrap method:
# "xy" uses xy-pair method
# "wxy" uses weighted (generalized) method
# "pwy" uses the parzen-wei-ying method
# "mcmb" uses the Markov chain marginal bootstrap method
#
#
function (object, se = NULL, covariance = FALSE, hs = TRUE, ...)
{
if(object$method == "lasso")
stop("no inference for lasso'd rq fitting: try rqss (if brave)")
mt <- terms(object)
m <- model.frame(object)
y <- model.response(m)
x <- model.matrix(mt, m, contrasts = object$contrasts)
wt <- as.vector(model.weights(object$model))
tau <- object$tau
eps <- .Machine$double.eps^(1/2)
coef <- coefficients(object)
if (is.matrix(coef))
coef <- coef[, 1]
vnames <- dimnames(x)[[2]]
resid <- object$residuals
n <- length(resid)
p <- length(coef)
rdf <- n - p
if (!is.null(wt)) {
resid <- resid * wt
x <- x * wt
y <- y * wt
}
if (is.null(se)) {
if (n < 1001 & covariance == FALSE)
se <- "rank"
else se <- "nid"
}
if (se == "rank") {
f <- rq.fit.br(x, y, tau = tau, ci = TRUE, ...)
}
if (se == "iid") {
xxinv <- diag(p)
xxinv <- backsolve(qr(x)$qr[1:p, 1:p,drop=FALSE], xxinv)
xxinv <- xxinv %*% t(xxinv)
pz <- sum(abs(resid) < eps)
h <- max(p + 1, ceiling(n * bandwidth.rq(tau, n, hs = hs)))
ir <- (pz + 1):(h + pz + 1)
ord.resid <- sort(resid[order(abs(resid))][ir])
xt <- ir/(n - p)
sparsity <- rq(ord.resid ~ xt)$coef[2]
cov <- sparsity^2 * xxinv * tau * (1 - tau)
scale <- 1/sparsity
serr <- sqrt(diag(cov))
}
else if (se == "nid") {
h <- bandwidth.rq(tau, n, hs = hs)
if (tau + h > 1)
stop("tau + h > 1: error in summary.rq")
if (tau - h < 0)
stop("tau - h < 0: error in summary.rq")
bhi <- rq.fit.fnb(x, y, tau = tau + h)$coef
blo <- rq.fit.fnb(x, y, tau = tau - h)$coef
dyhat <- x %*% (bhi - blo)
if (any(dyhat <= 0))
warning(paste(sum(dyhat <= 0), "non-positive fis"))
f <- pmax(0, (2 * h)/(dyhat - eps))
fxxinv <- diag(p)
fxxinv <- backsolve(qr(sqrt(f) * x)$qr[1:p, 1:p,drop=FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
fxxinv
scale <- mean(f)
serr <- sqrt(diag(cov))
}
else if (se == "ker") {
h <- bandwidth.rq(tau, n, hs = hs)
if (tau + h > 1)
stop("tau + h > 1: error in summary.rq")
if (tau - h < 0)
stop("tau - h < 0: error in summary.rq")
uhat <- c(y - x %*% coef)
h <- (qnorm(tau + h) - qnorm(tau - h))*
min(sqrt(var(uhat)), ( quantile(uhat,.75)- quantile(uhat, .25))/1.34 )
f <- dnorm(uhat/h)/h
fxxinv <- diag(p)
fxxinv <- backsolve(qr(sqrt(f) * x)$qr[1:p, 1:p,drop=FALSE], fxxinv)
fxxinv <- fxxinv %*% t(fxxinv)
cov <- tau * (1 - tau) * fxxinv %*% crossprod(x) %*%
fxxinv
scale <- mean(f)
serr <- sqrt(diag(cov))
}
else if (se == "boot") {
B <- boot.rq(x, y, tau, ...)
cov <- cov(B)
serr <- sqrt(diag(cov))
}
if( se == "rank"){
coef <- f$coef
}
else {
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
if(se == "iid") object$scale <- scale
if(se %in% c("nid", "ker")) {
object$Hinv <- fxxinv
object$J <- crossprod(x)
object$scale <- scale
}
else if (se == "boot") {
object$B <- B
}
}
object$coefficients <- coef
object$residuals <- resid
object$rdf <- rdf
object$tau <- tau
class(object) <- "summary.rq"
object
}
akj <- function(x,
z = seq(min(x), max(x), length = 2 * length(x)),
p = rep(1/length(x), length(x)),
h = -1, alpha = 0.5, kappa = 0.9, iker1 = 0)
{
nx <- length(x)
stopifnot(is.numeric(x),
length(p) == nx,
any((iker1 <- as.integer(iker1)) == 0:1))
nz <- length(z)
if(is.unsorted(x))
x <- sort(x)
.Fortran("akj",
as.double(x),
as.double(z),
as.double(p),
iker1,
dens = double(nz),
psi = double(nz),
score = double(nz),
as.integer(nx),
as.integer(nz),
h = as.double(h),
as.double(alpha),
as.double(kappa),
double(nx),
PACKAGE = "quantreg")[c("dens", "psi", "score", "h")]
}
"lm.fit.recursive" <-
function(X, y, int = TRUE)
{
if(int)
X <- cbind(1, X)
p <- ncol(X)
n <- nrow(X)
D <- qr(X[1:p, ])
A <- qr.coef(D, diag(p))
A[is.na(A)] <- 0
A <- crossprod(t(A))
Ax <- rep(0, p)
b <- matrix(0, p, n)
b[, p] <- qr.coef(D, y[1:p])
b[is.na(b)] <- 0
z <- .Fortran( "rls",
as.integer(n),
as.integer(p),
as.double(t(X)),
as.double(y),
b = as.double(b),
as.double(A),
as.double(Ax),
PACKAGE = "quantreg")
bhat <- matrix(z$b, p, n)
return(bhat)
}
"rq.fit.hogg" <-
function (x, y, taus = c(.1,.3,.5), weights = c(.7,.2,.1),
R= NULL, r = NULL, beta = 0.99995, eps = 1e-06)
{
n <- length(y)
n2 <- nrow(R)
m <- length(taus)
p <- ncol(x)+m
if (n != nrow(x))
stop("x and y don't match n")
if (m != length(weights))
stop("taus and weights differ in length")
if (any(taus < eps) || any(taus > 1 - eps))
stop("taus outside (0,1)")
W <- diag(weights)
if(m == 1) W <- weights
x <- as.matrix(x)
X <- cbind(kronecker(W,rep(1,n)),kronecker(weights,x))
y <- kronecker(weights,y)
rhs <- c(weights*(1 - taus)*n, sum(weights*(1-taus)) * apply(x, 2, sum))
if(n2!=length(r))
stop("R and r of incompatible dimension")
if(ncol(R)!=p)
stop("R and X of incompatible dimension")
d <- rep(1, m*n)
u <- rep(1, m*n)
if(length(r)){
wn1 <- rep(0, 10 * m*n)
wn1[1:(m*n)] <- .5
wn2 <- rep(0,6*n2)
wn2[1:n2] <- 1
z <- .Fortran("rqfnc", as.integer(m*n), as.integer(n2), as.integer(p),
a1 = as.double(t(as.matrix(X))), c1 = as.double(-y),
a2 = as.double(t(as.matrix(R))), c2 = as.double(-r),
rhs = as.double(rhs), d1 = double(m*n), d2 = double(n2),
as.double(u), beta = as.double(beta), eps = as.double(eps),
wn1 = as.double(wn1), wn2 = as.double(wn2), wp = double((p + 3) * p),
it.count = integer(3), info = integer(1), PACKAGE = "quantreg")
}
else{
wn <- rep(0, 10 * m*n)
wn[1:(m*n)] <- .5
z <- .Fortran("rqfnb", as.integer(m*n), as.integer(p), a = as.double(t(as.matrix(X))),
c = as.double(-y), rhs = as.double(rhs), d = as.double(d), as.double(u),
beta = as.double(beta), eps = as.double(eps), wn = as.double(wn),
wp = double((p + 3) * p), it.count = integer(2), info = integer(1),
PACKAGE = "quantreg")
}
if (z$info != 0)
warning(paste("Info = ", z$info, "in stepy: singular design: iterations ", z$it.count[1]))
coefficients <- -z$wp[1:p]
if(any(is.na(coefficients)))stop("NA coefs: infeasible problem?")
list(coefficients = coefficients, nit = z$it.count, flag = z$info)
}
