[R] Factor loadings and principal component plots
Jari Oksanen
jarioksa at sun3.oulu.fi
Tue May 4 10:02:42 CEST 2004
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")
More information about the R-help
mailing list