gx.robmva: Function to undertake a Robust Exploratory Multivariate Data Analysis

Description

The function carries out a robust Principal Components Analysis (PCA) and estimates the Mahalanobis distances for a non-compositional dataset and places them in an object to be saved and post-processed for display and further manipulation. Robust procedures are used, MCD, MVE or user supplied weights, for classical procedures see gx.mva. For results display see gx.rqpca.screeplot, gx.rqpca.plot, gx.rotate, gx.md.plot and gx.md.print. For closed compositional data sets, such as geochemical analyses, use gx.robmva.closed.

Usage

gx.robmva(xx, proc = "mcd", wts = NULL, main = deparse(substitute(xx)))

Arguments

a n by p data matrix to be processed.

proc

by default proc = "mcd" for the Minimimum Covariance Determinant (MCD) robust procedure. Setting proc = "mve" results in the Minimum Volume Ellipsoid (MVE) procedure being used. If p > 50 the MVE procedure is used.

wts

by default wts = NULL and the MCD or MVE estimation procedures will be used. If, however, a vector of n 0 or 1 weights are supplied these will be used for robust estimation and the value of proc

main

by default the name of the object xx, main = deparse(substitute(xx)), it may be replaced by the user, but this is not recommended, see Details below.

`Value`

The following are returned as an object to be saved for subsequent display, etc.:
mainby default (recommended) the input data matrix name.
inputthe data matrix name, input = deparse(substitute(xx)), retained to be used by post-processing display functions.
procthe robust procedure used, the value of proc will be "mcd", "mve" or "wts".
nthe total number of individuals (observations, cases or samples) in the input data matrix.
ncthe number of individuals remaining in the core data subset following the robust estimation, i.e. the sum of those individuals with wts = 1.
pthe number of variables on which the multivariate operations were based.
ifilrflag for gx.md.plot, set to FALSE.
matnamesthe row numbers and column headings of the input matrix.
wtsthe vector of weights for the n individuals arising from the robust estimation of the covariance matrix and means.
meanthe vector the weighted means for the p variables.
covthe p by p weighted covariance matrix for the n by p data matrix.
sdthe vector of weighted standard deviations for the p variables.
sndthe n by p matrix of weighted standard normal deviates.
rthe p by p matrix of weighted Pearson product moment correlation coefficients.
eigenvaluesthe vector of p eigenvalues of the scaled Pearson robust correlation matrix for RQ analysis, see Grunsky (2001).
econtribthe vector of p robustly estimated eigenvalues each expressed as a percentage of the sum of the eigenvalues.
eigenvectorsthe n by p matrix of robustly estimated eigenvectors.
rloadthe p by p matrix of robust Principal Component (PC) loadings.
rcrthe p by p matrix containing the percentages of the variability of each variable (columns) expressed in each robust PC (rows).
rqscorethe n by p matrix of the n individuals scores on the p robust PCs.
vcontriba vector of p variances of the columns of rqscore.
pvcontribthe vector of p variances of the columns of rqscore expressed as percentages.  This is a check on vector econtrib, the values should be identical for a classical PCA.  However, for robust PCAs this is not so as the trimmed individuals from the robust estimation have been re-introduced.  As a consequence pvcontrib can be very different from econtrib.  The plotting of PCs containing high proportions of the variance in robust PCAs can be useful for identifying outliers
cpvcontribthe vector of p cumulative sums of pvcontrib, see above.
mdthe vector of n robust Mahalanobis distances (MDs) for the n by p input matrix.
ppmthe vector of n robust predicted probabilities of population membership, see Garrett (1990).
epmthe vector of n robust empirical Chi-square probabilities for the MDs.
nrthe number of PCs that have been rotated.  At this stage of a data analysis nr = NULL in order to control PC plot axis labelling.

`Details`

If main is undefined the name of the matrix object passed to the function is used to identify the object.  This is the recommended procedure as it helps to track the progression of a data analysis.  Alternate plot titles are best defined when the saved object is passed to gx.rqpca.plot, gx.rqpca.screeplot or gx.md.plot for display.  If no plot title is required set main = "", or if a user defined plot title is required it may be defined, e.g., main = "Plot Title Text".

`References`

Garrett, R.G., 1990. A robust multivariate allocation procedure with applications to geochemical data. In Proc. Colloquium on Statistical Applications in the Earth Sciences (Eds F.P. Agterberg & G.F. Bonham-Carter). Geological Survey of Canada Paper 89-9, pp. 309-318.

Garrett, R.G., 1993. Another cry from the heart. Explore - Assoc. Exploration Geochemists Newsletter, 81:9-14.

Grunsky, E.C., 2001. A program for computinRQ-mode principal components analysis for S-Plus and R. Computers & Geosciences, 27(2):229-235.

Reimann, C., Filzmoser, P., Garrett, R. and Dutter, R., 2008. Statistical Data Analysis Explained: Applied Environmental Statistics with R. John Wiley & Sons, Ltd., 362 p.

`See Also`

ltdl.fix.df, remove.na, na.omit, gx.rqpca.screeplot, gx.rqpca.plot, 
gx.rotate, gx.md.plot, gx.md.print, gx.robmva

`Examples`

Run this code## Make test data available
data(sind)
attach(sind)
sind.mat <- as.matrix(sind[, -c(1:3)])
## Ensure all data are in the same units (mg/kg)
sind.mat2open <- sind.mat
sind.mat2open[, 2] <- sind.mat2open[, 2] * 10000

## Generate gx.robmva object after an clr transform using prior weights
sind.wts <- rep(1, 25)
## Set weights of all 6 high weighted sums, truly anomalous, individuals
## to zero
sind.wts[c(9, 12, 13, 15, 16, 21)] <- 0
sind.save.rob.clr <- gx.robmva(clr(sind.mat2open), wts = sind.wts)
## Display saved object with alternate main titles
gx.rqpca.screeplot(sind.save.rob.clr,
main = "Howarth & Sinding-Larsen
Stream Sediments, clr Transformed Data",
cex.main = 0.8)
gx.rqpca.plot(sind.save.rob.clr,
main = "Howarth & Sinding-Larsen
Stream Sediments, clr Transformed Data",
cex.main = 0.8)
gx.md.plot(sind.save.rob.clr,
main = "Howarth & Sinding-Larsen
Stream Sediments
clr Transformed Data",
cex.main = 0.8, cex.lab = 0.9)
gx.md.print(cbind(sind.save.rob.clr$md,sind.save.rob.clr$ppm,
ID, Zn, Cu, Cd, Fe, Mn), pcut = 0.5)

## Generate gx.robmva object after an ilr transform
## ** Function gx.robmva.closed should be used **
## ** This is for demonstration only **
sind.save.rob.ilr <- gx.robmva(ilr(sind.mat2open))
## Display saved object with alternate main titles
gx.rqpca.screeplot(sind.save.rob.ilr,
main = "Howarth & Sinding-Larsen
Stream Sediments, ilr Transformed Data",
cex.main = 0.8)
gx.rqpca.plot(sind.save.rob.ilr,
main = "Howarth & Sinding-Larsen
Stream Sediments, ilr Transformed Data",
cex.main = 0.8)
gx.md.plot(sind.save.rob.ilr,
main = "Howarth & Sinding-Larsen
Stream Sediments
ilr Transformed Data",
cex.main = 0.8, cex.lab = 0.9)
gx.md.print(cbind(sind.save.rob.ilr$md,sind.save.rob.ilr$ppm,
ID, Zn, Cu, Cd, Fe, Mn), pcut = 0.5)

## Clean-up and detach test data
rm(sind.mat)
rm(sind.mat2open)
rm(sind.wts)
rm(sind.save.rob.clr)
rm(sind.save.rob.ilr)
detach(sind)
Run the code above in your browser using DataLab