wlogreg:

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

Usage

wlogreg(x0, y, initwml = FALSE, const = 0.5, kmax = 1000, maxhalf = 10)

Arguments

x0
y
initwml
const
kmax
maxhalf

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 (x0, y, initwml = FALSE, const = 0.5, kmax = 1000, maxhalf = 10) 
{
    library(MASS)
    x0 = as.matrix(x0)
    p = ncol(x0) + 1
    p0 = p - 1
    sigmamin = 1e-04
    zz = elimna(cbind(x, y))
    x = as.matrix(zz[, 1:p0])
    y = zz[, p]
    n = nrow(x)
    print(n)
    x = as.matrix(cbind(rep(1, n), x))
    print(rep(1, n))
    y = as.numeric(y)
    if (initwml == T) {
        hp = floor(n * (1 - 0.25)) + 1
        mcdx = cov.mcd(x0, quantile.used = hp, method = "mcd")
        rdx = sqrt(mahalanobis(x0, center = mcdx$center, cov = mcdx$cov))
        vc = sqrt(qchisq(0.975, p - 1))
        wrd = (rdx <= vc)
        gstart = glm(y ~ x0, family = binomial, subset = wrd)$coef
    }
    else {
        gstart = glm(y ~ x0, family = binomial)$coef
    }
    sigmastart = 1/sqrt(sum(gstart^2))
    xistart = gstart * sigmastart
    stscores = x %*% xistart
    sigma1 = sigmastart
    oldobj = mean(phiBY3(stscores/sigmastart, y, const))
    kstep = jhalf = 1
    while ((kstep < kmax) & (jhalf < maxhalf)) {
        unisig <- function(sigma) {
            mean(phiBY3(stscores/sigma, y, const))
        }
        optimsig = nlminb(sigma1, unisig, lower = 0)
        sigma1 = optimsig$par
        if (sigma1 < sigmamin) {
            print("Explosion")
            kstep = kmax
        }
        else {
            gamma1 = xistart/sigma1
            scores = stscores/sigma1
            newobj = mean(phiBY3(scores, y, const))
            oldobj = newobj
            gradBY3 = colMeans((derphiBY3(scores, y, const) %*% 
                matrix(1, ncol = p)) * x)
            h = -gradBY3 + ((gradBY3 %*% xistart) * xistart)
            finalstep = h/sqrt(sum(h^2))
            xi1 = xistart + finalstep
            xi1 = xi1/(sum(xi1^2))
            scores1 = (x %*% xi1)/sigma1
            newobj = mean(phiBY3(scores1, y, const))
            hstep = jhalf = 1
            while ((jhalf <= maxhalf) & (newobj > oldobj)) {
                hstep = hstep/2
                xi1 = xistart + finalstep * hstep
                xi1 = xi1/sqrt(sum(xi1^2))
                scores1 = x %*% xi1/sigma1
                newobj = mean(phiBY3(scores1, y, const))
                jhalf = jhalf + 1
            }
            CONV = F
            if ((jhalf == maxhalf + 1) & (newobj > oldobj)) {
                CONV = T
            }
            else {
                jhalf = 1
                xistart = xi1
                oldobj = newobj
                stscores = x %*% xi1
                kstep = kstep + 1
            }
        }
    }
    if (kstep == kmax) {
        CONV = F
        result = list(convergence = FALSE, objective = 0, coef = t(rep(NA, 
            p)))
    }
    else {
        gammaest = xistart/sigma1
        stander = sterby3(x0, y, const, gammaest)
        result = list(convergence = CONV, coef = t(gammaest), 
            sterror = stander)
    }
    return(result)
  }

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