R/gx.mva.R

In rgr: Applied Geochemistry EDA

gx.mva <-
function(xx, main = deparse(substitute(xx)))
{
     # Procedure to undertake non-robust multivariate data analyses; the
     # object generated is identical to that of gx.robmva so that the
     # saved list may be passed to other rotation and display functions.
     # Thus weights are set to 1, and other variables are set to 
     # appropriate defaults.  The estimation of Mahalanobis distances is
     # only undertaken if x is non-singular, i.e., the lowest eigenvalue
     # is > 10e-4.
     #
     # Note this procedure uses svd() rather than the classic solve().
     #
     # PCA output may be plotted with gx.rqpca.plot and gx.rqpca.screeplot, 
     # and Mahalanobis distances may be plotted with gx.md.plot.
     #
     if(!is.matrix(xx)) stop(deparse(substitute(xx)), " is not a Matrix")
     # Remove any rows containing NAs
     temp.x <- remove.na(xx)
     x <- temp.x$x; n <- temp.x$n; p <- temp.x$m
     matnames <- dimnames(xx)
     if(is.null(matnames[[1]])) matnames[[1]] <- c(1:n)
     wts <- numeric(n); wts[1:n] <- 1
     nc <- n
     cat("  n =", n, "\tnc =", n, "\tp =", p, "\t\tnc/p =", round(nc/p, 2), "\n")
     if(nc < 5 * p)
         cat("  *** Proceed with Care, Core Size is < 5p ***\n")
     if(nc < 3 * p)
         cat("  *** Proceed With Great Care, Core Size = ", n, ", which is < 3p ***\n")
     # Compute means & SDs, and standardize the data set.  Note cov.wt()
     # is used in order to have cov and r entries for the saved object
     save <- cov.wt(x, wt = wts, cor = TRUE)
     xmean <- save$center
     xsd <- sqrt(diag(save$cov))
     # Compute SNDs
     temp <- sweep(x, 2, xmean, "-")
     snd <- sweep(temp, 2, xsd, "/")
     # The following procedure duplicates that in rqpca() for iwght = "r"
     # Standardize centred x for R-mode PCA and compute
     w <- sweep(temp, 2, sqrt(n) * xsd, "/")
     wt <- t(as.matrix(w))
     a <- wt %*% as.matrix(w)
     b <- svd(a)
     cat("  Eigenvalues:", signif(b$d, 4), "\n")
     sumc <- sum(b$d)
     econtrib <- 100 * (b$d/sumc)
     cat("     as %ages:", round(econtrib, 1), "\n")
     rqscore <- w %*% b$v
     vcontrib <- numeric(p)
     for (j in 1:p) vcontrib[j] <- var(rqscore[, j])
     cat("  Score S^2s :", signif(vcontrib, 4), "\n")
     b1 <- b$v * 0
     diag(b1) <- sqrt(b$d)
     rload <- b$v %*% b1
     rcr <- rload[,  ] * 0
     rcr1 <- apply(rload^2, 1, sum)
     rcr <- 100 * sweep(rload^2, 1, rcr1, "/")
     dimnames(rload)[[1]] <- dimnames(rcr)[[1]] <- matnames[[2]]
     #
     # Test for non-singularity and compute Mahalanobis distances
     if(b$d[p] > 0.001) {
         md <- mahalanobis(x, save$center, save$cov)
         temp <- (nc - p)/(p * (nc + 1))
         ppm <- 1 - pf(temp * md, p, nc - p)
         epm <- 1 - pchisq(md, p)
     }
     else {
         cat("  Lowest eigenvalue < 10^-4, Mahalanobis distances not computed\n")
         md <- NULL
         ppm <- NULL
         epm <- NULL
     }
     invisible(list(main = main, input = deparse(substitute(xx)), proc = "cov", 
         n = n, nc = nc, p = p, ifilr = FALSE, matnames = matnames, wts = wts,
         mean = xmean, cov = save$cov, sd = xsd, snd = snd, r = save$cor,
         eigenvalues = b$d, econtrib = econtrib, eigenvectors = b$v,
         rload = rload, rcr = rcr, rqscore = rqscore, md = md, ppm = ppm,
         epm = epm, nr = NULL))
}