R/pcaFit.R

In mvdalab: Multivariate Data Analysis Laboratory

pcaFit <- function (data, scale = TRUE, ncomp = NULL) {
X <- my.dummy.df(data)
if (is.null(ncomp)) {
  if (nrow(X) <= ncol(X)) {
    ncomp <- nrow(X)
  }
  else ncomp <- ncol(X)
} else ncomp <- ncomp
Xmeans <- colMeans(X)
Scale <- scale
nobj <- nrow(X)
npred <- ncol(X)
X <- X - rep(Xmeans, each = nobj)
if (Scale == TRUE) {
  scale. <- sqrt(colSums((X - rep(colMeans(X), each = nobj))^2)/(nobj -
                                                                   1))
  if (any(abs(scale.) < .Machine$double.eps^0.5))
    warning("Scaling with (near) zero standard deviation")
  X <- X/rep(scale., each = nobj)
} else X
SVD <- svd(X, nu = ncomp, nv = ncomp)
D.total <- diag((SVD$d^2))/(nrow(X) - 1)
if (ncomp == 1) {
  U <- as.matrix(SVD$u[, 1]) %*% diag(as.matrix(SVD$d[1:ncomp]))
  D <- diag(as.matrix(SVD$d[1:ncomp]^2))/(nrow(X) - 1)
} else {
  U <- SVD$u[, 1:ncomp] %*% diag((SVD$d[1:ncomp]))
  D <- diag((SVD$d[1:ncomp]^2))/(nrow(X) - 1)
}
loadings.a <- t(sqrt(D) %*% t(SVD$v[, 1:ncomp]))
loadings <- as.matrix(apply(loadings.a, 2, function(x) x/sqrt(crossprod(x)[1])))
scores <- as.matrix(U)
cumpress <- rep(0L, ncomp + 1)
  press <- matrix(0L, ncomp + 1, npred)
  s <- c(nobj, npred)
  for (i in c(0:ncomp)){
    if (i > 0){
      p2 <- as.matrix(loadings[, 1:i])
      srec <- scores[, 1:i] %*% t(p2)
      erec <- X - srec
      term3_p <- erec
      term1_p <- X * (rep(1, s[1]) %*% t(rowSums(p2^2)))
    }else{
      term1_p <- matrix(0L, nobj, npred)
      term3_p <- X
    }
    term1 <- term1_p^2
    term2 <- 2*term1_p*term3_p
    term3 <- term3_p^2
    press[i + 1,] <- colSums(rbind(colSums(term1), colSums(term2), colSums(term3)))
    cumpress[i + 1] <- sum(press[i + 1, ])
  }
  CV <- cumpress
  if(ncomp == 1) {
    Percents.Explained <- data.frame(Individual.Percent = (D / sum(diag(D.total))) *
                                       100, Cumulative.Percent = cumsum(D / sum(diag(D.total))) *
                                       100)
  } else {
    Percents.Explained <- data.frame(Individual.Percent = (diag(D)/sum(diag(D.total))) *
                                       100, Cumulative.Percent = cumsum(diag(D)/sum(diag(D.total))) *
                                       100)
  }
Percents.Explained$Comp <- 1:nrow(Percents.Explained)
row.names(loadings) <- names(X)
output <- list(loadings = loadings, scores = scores, D = D,
               Xdata = X, Percents.Explained = Percents.Explained, CV = CV,
               ncomp = ncomp, Xmeans = Xmeans, method = "PCA")
class(output$scores) <- "scores"
class(output$loadings) <- "loadings"
class(output) <- "mvdapca"
output
}