Source

mulligan rubinstein bounds / myrq.R

Full commit
######################################################################
## rq that drops collinear x's automatically, largely copied from rq in
## the quantreg package
myrq <- 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 <- model.weights(mf)
  Y <- model.response(mf)
  X <- model.matrix(mt, mf, contrasts)
    ols <- lm(Y ~ X-1,x=TRUE,weights=weights)
  if (any(is.na(ols$coef))) {
    usex <- sapply(colnames(X),
                   function(xi) {!is.na(ols$coef[paste("X",xi,sep="")]) })
    addcoefs <- names(ols$coef)[is.na(ols$coef)]
  } else {
    usex <- (1:ncol(X)>0)
    addcoefs <- NULL
  }    
  eps <- .Machine$double.eps^(2/3)
  Rho <- function(u, tau) u * (tau - (u < 0))
  if (length(tau) > 1) {
    if (any(tau < -eps) || any(tau > 1 + eps)) 
      stop("invalid tau:  taus should be >= 0 and <= 1")
        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)==length(Y)) 
          rq.wfit(X[,usex], Y, tau = tau[i], weights, method, 
                          ...)
        else rq.fit(X[,usex], Y, tau = tau[i], method, ...)
      }
      coef[usex, i] <- z$coefficients
            if (length(addcoefs)>0) coef[!usex, i] <- NA
      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)
    fit <- {
      if (length(weights)) 
        rq.wfit(X[,usex], Y, tau = tau, weights, method, ...)
      else rq.fit(X[,usex], Y, tau = tau, method, ...)
    }
    if (length(addcoefs)>0) {
      coef <- fit$coef
      fit$coefficients <- (1:ncol(X))*0
      fit$coefficients[usex] <- coef
      fit$coefficients[!usex] <- NA
      names(fit$coefficients) <- dimnames(X)[[2]]
    }
    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$fitted.values <- drop(fit$fitted.values)
  attr(fit, "na.message") <- attr(m, "na.message")
  if (model) 
    fit$model <- mf
  return(fit)
}

################################################################################
# myrq fix for dropping colinear x's for rq breaks predict afterward,
# this fixes predict
mypredict <- function(object,newdata) {
  object$coefficients[is.na(object$coefficients)] <- 0
  predict(object,newdata=newdata)
}