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

Usage Arguments Examples

robpca(x, pval = ncol(x), kmax = 10, alpha = 0.75, h, mcd = 1, plots = 1, labsd = 3, labod = 3, classic = 0, plotit = FALSE, pr = TRUE, SEED = TRUE, STAND = TRUE, est = tmean, varfun = winvar, scree = TRUE, xlab = "Principal Component", ylab = "Proportion of Variance")

`x`
`pval`
`kmax`
`alpha`
`h`
`mcd`
`plots`
`labsd`
`labod`
`classic`
`plotit`
`pr`
`SEED`
`STAND`
`est`
`varfun`
`scree`
`xlab`
`ylab`

##---- 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, pval = ncol(x), kmax = 10, alpha = 0.75, h, mcd = 1, 
    plots = 1, labsd = 3, labod = 3, classic = 0, plotit = FALSE, 
    pr = TRUE, SEED = TRUE, STAND = TRUE, est = tmean, varfun = winvar, 
    scree = TRUE, xlab = "Principal Component", ylab = "Proportion of Variance") 
{
    x <- elimna(x)
    if (pval != ncol(x)) 
        scree = FALSE
    if (STAND) 
        x = standm(x, est = est, scat = varfun)
    if (SEED) 
        set.seed(2)
    k <- pval
    if (pr) 
        print(paste("Number of principal components specified is", 
            pval))
    if (!plotit) 
        plots <- 0
    library(MASS)
    if (missing(x)) {
        stop("Error in robpca: You have to provide at least some data")
    }
    data <- as.matrix(x)
    n <- nrow(data)
    p <- ncol(data)
    if (n < p) {
        X.svd <- kernelEVD(data)
    }
    else {
        X.svd <- classSVD(data)
    }
    if (X.svd$rank == 0) {
        stop("All data points collapse!")
    }
    kmax <- max(min(floor(kmax), floor(n/2), X.svd$rank), 1)
    k <- floor(k)
    if (k < 0) {
        k <- 0
    }
    else if (k > kmax) {
        warning("Attention robpca: The number of principal components k = ", 
            k, " is larger then kmax = ", kmax, "; k is set to ", 
            kmax, ".")
        k <- kmax
    }
    if (!missing(h) & !missing(alpha)) {
        stop("Error in robpca: Both inputarguments alpha and h are provided. Only one is required.")
    }
    if (missing(h) & missing(alpha)) {
        h <- min(floor(2 * floor((n + kmax + 1)/2) - n + 2 * 
            (n - floor((n + kmax + 1)/2)) * alpha), n)
    }
    if (!missing(h) & missing(alpha)) {
        alpha <- h/n
        if (k == 0) {
            if (h < floor((n + kmax + 1)/2)) {
                h <- floor((n + kmax + 1)/2)
                alpha <- h/n
                warning("Attention robpca: h should be larger than (n+kmax+1)/2. It is set to its minimum value ", 
                  h, ".")
            }
        }
        else {
            if (h < floor((n + k + 1)/2)) {
                h <- floor((n + k + 1)/2)
                alpha <- h/n
                warning("Attention robpca: h should be larger than (n+k+1)/2. It is set to its minimum value ", 
                  h, ".")
            }
        }
        if (h > n) {
            alpha <- 0.75
            if (k == 0) {
                h <- floor(2 * floor((n + kmax + 1)/2) - n + 
                  2 * (n - floor((n + kmax + 1)/2)) * alpha)
            }
            else {
                h <- floor(2 * floor((n + k + 1)/2) - n + 2 * 
                  (n - floor((n + k + 1)/2)) * alpha)
            }
            warning("Attention robpca: h should be smaller than n = ", 
                n, ". It is set to its default value ", h, ".")
        }
    }
    if (missing(h) & !missing(alpha)) {
        if (alpha < 0.5) {
            alpha <- 0.5
            warning("Attention robpca: Alpha should be larger then 0.5. It is set to 0.5.")
        }
        if (alpha >= 1) {
            alpha <- 0.75
            warning("Attention robpca: Alpha should be smaller then 1. It is set to its default value 0.75.")
        }
        if (k == 0) {
            h <- floor(2 * floor((n + kmax + 1)/2) - n + 2 * 
                (n - floor((n + kmax + 1)/2)) * alpha)
        }
        else {
            h <- floor(2 * floor((n + k + 1)/2) - n + 2 * (n - 
                floor((n + k + 1)/2)) * alpha)
        }
    }
    labsd <- floor(max(0, min(labsd, n)))
    labod <- floor(max(0, min(labod, n)))
    out <- list()
    Xa <- X.svd$scores
    center <- X.svd$centerofX
    rot <- X.svd$loadings
    p1 <- ncol(Xa)
    if ((p1 <= min(floor(n/5), kmax)) & (mcd == 1)) {
        if (k != 0) {
            k <- min(k, p1)
        }
        else {
            k <- p1
        }
        if (h < floor((nrow(Xa) + ncol(Xa) + 1)/2)) {
            h <- floor((nrow(Xa) + ncol(Xa) + 1)/2)
            cat("Message from robpca: The number of non-outlying observations h is set to ", 
                h, " in order to make the mcd algorithm function.\n", 
                sep = "")
        }
        Xa.mcd <- cov.mcd(as.data.frame(Xa), quan = h)
        Xa.mcd.svd <- svd(Xa.mcd$cov)
        scores <- (Xa - matrix(data = rep(Xa.mcd$center, times = nrow(Xa)), 
            nrow = nrow(Xa), ncol = ncol(Xa), byrow = T)) %*% 
            Xa.mcd.svd$u
        out$M <- center + as.vector(Xa.mcd$center %*% t(rot))
        out$L <- Xa.mcd.svd$d[1:k]
        if (scree) {
            pv = out$L
            cs = pv/sum(pv)
            cm = cumsum(cs)
            plot(rep(c(1:ncol(x)), 2), c(cs, cm), type = "n", 
                xlab = xlab, ylab = ylab)
            points(c(1:ncol(x)), cs, pch = "*")
            lines(c(1:ncol(x)), cs, lty = 1)
            points(c(1:ncol(x)), cm, pch = ".")
            lines(c(1:ncol(x)), cm, lty = 2)
        }
        out$P <- X.svd$loadings %*% Xa.mcd.svd$u[, 1:k]
        out$T <- as.matrix(scores[, 1:k])
        if (is.list(dimnames(data))) {
            dimnames(out$T)[[1]] <- dimnames(data)[[1]]
        }
        out$h <- h
        out$k <- k
        out$alpha <- alpha
    }
    else {
        directions <- choose(n, 2)
        ndirect <- min(250, directions)
        all <- (ndirect == directions)
        seed <- 0
        B <- extradir(Xa, ndirect, seed, all)
        Bnorm <- vector(mode = "numeric", length = nrow(B))
        Bnorm <- apply(B, 1, vecnorm)
        Bnormr <- Bnorm[Bnorm > 1e-12]
        B <- B[Bnorm > 1e-12, ]
        A <- diag(1/Bnormr) %*% B
        Y <- Xa %*% t(A)
        Z <- matrix(data = 0, nrow = n, ncol = length(Bnormr))
        for (i in 1:ncol(Z)) {
            univ <- unimcd(Y[, i], quan = h)
            if (univ$smcd < 1e-12) {
                r2 <- qr(data[univ$weights == 1, ])$rank
                if (r2 == 1) {
                  stop("Error in robpca: At least ", sum(univ$weights), 
                    " observations are identical.")
                }
            }
            else {
                Z[, i] <- abs(Y[, i] - univ$tmcd)/univ$smcd
            }
        }
        H0 <- order(apply(Z, 1, max))
        Xh <- Xa[H0[1:h], ]
        Xh.svd <- classSVD(Xh)
        kmax <- min(Xh.svd$rank, kmax)
        if ((k == 0) & (plots == 0)) {
            test <- which((Xh.svd$eigenvalues/Xh.svd$eigenvalues[1]) <= 
                0.001)
            if (length(test) != 0) {
                k <- min(min(Xh.svd$rank, test[1]), kmax)
            }
            else {
                k <- min(Xh.svd$rank, kmax)
            }
            cumulative <- cumsum(Xh.svd$eigenvalues[1:k])/sum(Xh.svd$eigenvalues)
            if (cumulative[k] > 0.8) {
                k <- which(cumulative >= 0.8)[1]
            }
            cat("Message from robpca: The number of principal components is set by the algorithm. It is set to ", 
                k, ".\n", sep = "")
        }
        else {
            if ((k == 0) & (plots != 0)) {
                loc <- 1:kmax
                plot(loc, Xh.svd$eigenvalues[1:kmax], type = "b", 
                  axes = FALSE, xlab = "Component", ylab = "Eigenvalue")
                axis(2)
                axis(1, at = loc)
                cumv <- cumsum(Xh.svd$eigenvalues)/sum(Xh.svd$eigenvalues)
                text(loc, Xh.svd$eigenvalues[1:kmax] + par("cxy")[2], 
                  as.character(signif(cumv[1:kmax], 2)))
                box <- dialogbox(title = "ROBPCA", controls = list(), 
                  buttons = c("OK"))
                box <- dialogbox.add.control(box, where = 1, 
                  statictext.control(paste("How many principal components would you like to retain?\nMaximum = ", 
                    kmax, sep = ""), size = c(200, 20)))
                box <- dialogbox.add.control(box, where = 2, 
                  editfield.control(label = "Your choice:", size = c(30, 
                    10)))
                input <- as.integer(dialogbox.display(box)$values$"Your choice:")
                k <- max(min(min(Xh.svd$rank, input), kmax), 
                  1)
            }
            else {
                k <- min(min(Xh.svd$rank, k), kmax)
            }
        }
        if (k != X.svd$rank) {
            XRc <- Xa - matrix(data = rep(Xh.svd$centerofX, times = nrow(Xa)), 
                nrow = nrow(Xa), ncol = ncol(Xa), byrow = T)
            Xtilde <- XRc %*% Xh.svd$loadings[, 1:k] %*% t(Xh.svd$loadings[, 
                1:k])
            Rdiff <- XRc - Xtilde
            odh <- apply(Rdiff, 1, vecnorm)
            ms <- unimcd(odh^(2/3), h)
            cutoffodh <- sqrt(qnorm(0.975, ms$tmcd, ms$smcd)^3)
            indexset <- (odh <= cutoffodh)
            Xh.svd <- classSVD(Xa[indexset, ])
            kmax <- min(Xh.svd$rank, kmax)
        }
        center <- center + Xh.svd$centerofX %*% t(rot)
        rot <- rot %*% Xh.svd$loadings
        Xstar <- (Xa - matrix(data = rep(Xh.svd$centerofX, times = nrow(Xa)), 
            nrow = nrow(Xa), ncol = ncol(Xa), byrow = T)) %*% 
            Xh.svd$loadings
        Xstar <- as.matrix(Xstar[, 1:k])
        rot <- as.matrix(rot[, 1:k])
        mah <- mahalanobis(Xstar, center = rep(0, ncol(Xstar)), 
            cov = diag(Xh.svd$eigenvalues[1:k], nrow = k))
        oldobj <- prod(Xh.svd$eigenvalues[1:k])
        niter <- 100
        for (j in 1:niter) {
            mah.order <- order(mah)
            Xh <- as.matrix(Xstar[mah.order[1:h], ])
            Xh.svd <- classSVD(Xh)
            obj <- prod(Xh.svd$eigenvalues)
            Xstar <- (Xstar - matrix(data = rep(Xh.svd$centerofX, 
                times = nrow(Xstar)), nrow = nrow(Xstar), ncol = ncol(Xstar), 
                byrow = T)) %*% Xh.svd$loadings
            center <- center + Xh.svd$centerofX %*% t(rot)
            rot <- rot %*% Xh.svd$loadings
            mah <- mahalanobis(Xstar, center = rep(0, ncol(Xstar)), 
                cov = diag(x = Xh.svd$eigenvalues, nrow = length(Xh.svd$eigenvalues)))
            if ((Xh.svd$rank == k) & (abs(oldobj - obj) < 1e-12)) {
                break
            }
            else {
                oldobj <- obj
                if (Xh.svd$rank < k) {
                  j <- 1
                  k <- Xh.svd$rank
                }
            }
        }
        Xstar.mcd <- cov.mcd(as.data.frame(Xstar), quan = h)
        covf <- Xstar.mcd$cov
        centerf <- Xstar.mcd$center
        covf.eigen <- eigen(covf)
        covf.eigen.values.sort <- greatsort(covf.eigen$values)
        P6 <- covf.eigen$vectors
        P6 <- covf.eigen$vectors[, covf.eigen.values.sort$index]
        out$T <- (Xstar - matrix(data = rep(centerf, times = n), 
            nrow = n, ncol = ncol(Xstar), byrow = T)) %*% covf.eigen$vectors[, 
            covf.eigen.values.sort$index]
        if (is.list(dimnames(data))) {
            dimnames(out$T)[[1]] <- dimnames(data)[[1]]
        }
        out$P <- rot %*% covf.eigen$vectors[, covf.eigen.values.sort$index]
        out$M <- as.vector(center + centerf %*% t(rot))
        out$L <- as.vector(covf.eigen$values)
        out$k <- k
        out$h <- h
        out$alpha <- alpha
    }
    oldClass(out) <- "robpca"
    out <- CompRobustDist(data, X.svd$rank, out, classic)
    if (classic == 1) {
        out <- CompClassicDist(X.svd, out)
    }
    if (plots == 1) {
        plot(out, classic, labod = labod, labsd = labsd)
    }
    return(out)
  }