Nothing
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))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.