princomp

[Paul Gilbert]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._