Below is a draft of a function "princomp" which reproduces Splus 3.3
results. In
writing this I have found there is an certain inconsistency in the scaling (N vs
N-1) used when cor is T or F. The net result is that for data which has already
been scaled (so cor=cov) the function gives slightly different results depending
on whether the argument cor is T or F. I have some reluctance in reproducing
this behavior, and perhaps someone could check if it has been changed in a more
recent version of Splus:
> z <- scale(iris[,,2], scale=T, center=T)
> z1 <- princomp(z, cor=F)
> z2 <- princomp(z, cor=T)
> z1$sdev
Comp.1 Comp.2 Comp.3 Comp.4
1.6934621 0.7316756 0.6221717 0.3601935> z2$sdev
Comp.1 Comp.2 Comp.3 Comp.4
1.7106550 0.7391040 0.6284883 0.3638504
I have at leasted tried to make the inconsistency obvious in the code, and
document it. Any other thoughts would be appreciated.
Paul Gilbert
______
princomp <- function(x, cor=FALSE, scores=TRUE,
subset=rep(TRUE, nrow(as.matrix(x)))) {
z<- as.matrix(x)[subset,, drop=F]
N <- nrow(z)
if(cor) cv <- get("cor",envir=.GlobalEnv)(z)
else cv <- cov(z)
# (svd can be used but gives different signs for some vectors)
edc <- eigen(cv)
cn <- paste("Comp.", 1:ncol(cv), sep="")
names(edc$values) <- cn
dimnames(edc$vectors) <- list(dimnames(x)[[2]], cn)
scr<- NULL
if (cor)
{sdev <- sqrt(edc$values)
sc <- (apply(z,2,var)*(N-1)/N)^0.5
if (scores)
scr<-(scale(z,center=T,scale=T) %*% edc$vectors)*sqrt(N/(N-1))
}
else
{sdev <- sqrt(edc$values*(N-1)/N)
sc <- rep(1, ncol(z))
if (scores)
scr<- (scale(z,center=T,scale=F) %*% edc$vectors)
}
names(sc) <- dimnames(x)[[2]]
edc <-list(sdev=sdev, loadings=edc$vectors,
center=apply(z,2,mean), scale=sc, n.obs=N, scores=scr)
# The splus function also return list elements factor.sdev, correlations
# and coef, but these are not documented in the help. coef seems to equal
# load. The splus function also return list elements call and terms which
# are not supported here.
class(edc) <- "princomp"
edc
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._