mbmreg:

In musto101/wilcox_R: What the package does (short line)

Usage

mbmreg(x, y, iter = 20, bend = 2 * sqrt(ncol(x) + 1)/nrow(x), xout = FALSE, outfun = outpro, ...)

Arguments

x
y
iter
bend
xout
outfun
...

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (x, y, iter = 20, bend = 2 * sqrt(ncol(x) + 1)/nrow(x), 
    xout = FALSE, outfun = outpro, ...) 
{
    x <- as.matrix(x)
    p1 <- ncol(x) + 1
    p <- ncol(x)
    xy <- cbind(x, y)
    xy <- elimna(xy)
    x <- xy[, 1:p]
    y <- xy[, p1]
    if (xout) {
        m <- cbind(x, y)
        flag <- outfun(x, plotit = FALSE, ...)$keep
        m <- m[flag, ]
        x <- m[, 1:p]
        y <- m[, p1]
    }
    x <- as.matrix(x)
    if (is.matrix(y)) {
        if (ncol(y) == 1) 
            y = as.vector(y)
    }
    x1 <- cbind(1, x)
    library(MASS)
    reslms <- lmsreg(x, y)$resid
    sighat <- sqrt(median(reslms^2))
    sighat <- 1.4826 * (1 + (5/(length(y) - ncol(x) - 1))) * 
        sighat
    if (sighat == 0) 
        warning("The estimated measure of scale, based on the residuals using lms regression, is zero")
    temp <- ifelse(sighat * reslms > 0, abs(reslms)/sighat, 0 * 
        reslms)
    wt <- ifelse(temp <= 2.5, 1, 0)
    init <- lsfit(x, y, wt)
    resid <- init$residuals
    nu <- sqrt(1 - hat(x1))
    low <- ncol(x) + 1
    for (it in 1:iter) {
        ev <- sort(abs(resid))
        scale <- median(ev[c(low:length(y))])/qnorm(0.75)
        rov <- (resid/scale)/nu
        psi <- ifelse(abs(rov) <= bend, rov, bend * sign(rov))
        wt <- nu * psi/(resid/scale)
        wt <- ifelse(temp <= 2.5, wt, 0)
        new <- lsfit(x, y, wt)
        if (abs(max(new$coef - init$coef) < 1e-04)) 
            break
        init$coef <- new$coef
        resid <- new$residuals
    }
    resid <- y - x1 %*% new$coef
    if (abs(max(new$coef - init$coef) >= 1e-04)) 
        paste("failed to converge in", iter, "steps")
    list(coef = new$coef, residuals = resid, w = wt)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.