[R] Factor loadings and principal component plots

[Jari Oksanen]

On Tue, 2004-05-04 at 09:56, Prof Brian Ripley wrote:
> On 4 May 2004, Jari Oksanen wrote:
> 
> > On Tue, 2004-05-04 at 09:34, Prof Brian Ripley wrote:
> > 
> > > 
> > > Yes, but princomp is the recommended way, not prcomp.
> > 
> > But the documentation seems to recommend prcomp:
> 
> For numerical accuracy, but not for flexibility.
> 
Wouldn't the best alternative be to combine flexibility and accuracy
into one alternative? I mean, I'd still use prcomp after reading the
help pages, and I'd put more weight on accuracy than on flexibility. A
quick exploitation of the princomp would yield the attached "flexible
prcomp" code.

prcomp is more flexible at least in one point: it can handle data with
less units than variables.

cheers, jari oksanen
-- 
Jari Oksanen <jarioksa at sun3.oulu.fi>
-------------- next part --------------
"prcomp.default" <-
    function (x, retx = TRUE, center = TRUE, scale. = FALSE, tol = NULL,
              subset = rep(TRUE, nrow(as.matrix(x))), ...) 
{
    x <- as.matrix(x)
    x <- x[subset, , drop = FALSE]
    x <- scale(x, center = center, scale = scale.)
    s <- svd(x, nu = 0)
    if (!is.null(tol)) {
        rank <- sum(s$d > (s$d[1] * tol))
        if (rank < ncol(x)) 
            s$v <- s$v[, 1:rank, drop = FALSE]
    }
    s$d <- s$d/sqrt(max(1, nrow(x) - 1))
    dimnames(s$v) <- list(colnames(x), paste("PC", seq(len = ncol(s$v)), 
                                             sep = ""))
    r <- list(sdev = s$d, rotation = s$v)
    if (retx) 
        r$x <- x %*% s$v
    class(r) <- "prcomp"
    r
}
"prcomp.formula" <-
    function (formula, data = NULL, subset, na.action, ...) 
{
    mt <- terms(formula, data = data)
    if (attr(mt, "response") > 0) 
        stop("response not allowed in formula")
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    mf$... <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval.parent(mf)
    if (any(sapply(mf, function(x) is.factor(x) || !is.numeric(x)))) 
        stop("PCA applies only to numerical variables")
    na.act <- attr(mf, "na.action")
    mt <- attr(mf, "terms")
    attr(mt, "intercept") <- 0
    x <- model.matrix(mt, mf)
    res <- prcomp.default(x, ...)
    cl[[1]] <- as.name("prcomp")
    res$call <- cl
    if (!is.null(na.act)) {
        res$na.action <- na.act
        if (!is.null(sc <- res$x)) 
            res$x <- napredict(na.act, sc)
    }
    res
}
"prcomp" <-
    function (x, ...) 
    UseMethod("prcomp")