EMMAgeo (version 0.9.4)

test.robustness: Function to test model robustness.

Description

All possible combinations of number of end-members and weight transformation limits are used to perform EMMA. The resulting loadings are written to an output matrix and mode positions (i.e. class with maximum loading) of all loadings are evaluated and returned.

Usage

test.robustness(X, q, l, P, c, classunits, ID, rotation = "Varimax", ol.rej,
  mRt.rej, plot = FALSE, ..., pm = FALSE)

Arguments

X

Numeric matrix with m samples (rows) and n variables (columns).

q

Numeric vector with number of end-members to be modelled.

l

Numeric vector specifying the weight tranformation limits, i.e. quantiles; default is 0.

P

Numeric matrix, optional alternative input parameters for q and l, either of the form m:3 with m variations in the columns q.min, q.max, l or of the form m:2 with m variations in the columns q, l.

c

Numeric scalar specifying the constant sum scaling parameter, e.g. 1, 100, 1000; default is 100.

classunits

Numeric vector, optional class units (e.g. phi classes or micrometers) of the same length as columns of X.

ID

Numeric or character vector, optional sample IDs of the same length as columns of X.

rotation

Character scalar, rotation type, default is "Varimax" (cf. Dietze et al., 2012). One out of the rotations provided in GPArotation is possible (cf. rotations).

ol.rej

Numeric scalar, optional rejection threshold for overlapping criterion. All model runs with overlapping end-members greater than the specified integer will be removed.

mRt.rej

Numeric scalar, optional rejection threshold for mean total explained variance criterion. All modelled end-members below the specified value will be removed.

plot

Logical scalar, optional graphical output of the results, default is FALSE. If set to TRUE, end-member loadings and end-member scores are plotted.

pm

Logical scalar to enable pm.

…

Additional arguments passed to the plot function. Since the function returns two plots, additional graphical parameters must be specified as vector with the first element for the first plot and the second element for the second plot. If graphical parameters are natively vectors (e.g. a sequence of colours), they must be specified as matrices with each vector as a row. If colours are specified, colour should be used instead of col. ylim can only be modified for the first plot. See example section for further advice.

Value

A list with objects

q

Vector with q.

l

Vector with l.

modes

Vector with mode class.

mRt

Vector with mean total explained variance.

ol

Vector with n overlapping end-members.

loadings

Matrix with normalised rescaled end-member loadings.

Vqsn

Matrix with rescaled end-member loadings.

Vqn

Matrix with normalised factor loadings.

Details

The function value $loadings is redundant but was added for user convenience.

References

Dietze E, Hartmann K, Diekmann B, IJmker J, Lehmkuhl F, Opitz S, Stauch G, Wuennemann B, Borchers A. 2012. An end-member algorithm for deciphering modern detrital processes from lake sediments of Lake Donggi Cona, NE Tibetan Plateau, China. Sedimentary Geology 243-244: 169-180.

Examples

# NOT RUN {
## load example data set
data(X, envir = environment())

## Example 1 - perform the most simple test
q  <- 4:7
l <- seq(from = 0, to = 0.1, by = 0.02)

M1  <- test.robustness(X = X, q = q, l = l,
                       ol.rej = 1, mRt.rej = 0.8,
                       plot = TRUE,
                       colour = c(4, 7),
                       xlab = c(expression(paste("Grain size (", phi, ")",
                                                 sep = "")),
                                expression(paste("Grain size (", phi, ")",
                                                 sep = ""))))

## Example 2 -  perform the test without rejection criteria and plots
P  <- cbind(rep(q[1], length(l)),
            rep(q[3], length(l)),
            l)
M2  <- test.robustness(X = X, P = P)

## Plot 1 - end-member loadings which do not overlap and yielded mRt > 0.80.
plot(M2$Vqsn[1,], type = "l", ylim = c(0, max(M2$Vqsn, na.rm = TRUE)),
     main = "End-member loadings")
  for (i in 2:nrow(M2$Vqsn)) lines(M2$Vqsn[i,])

# Plot 2 - histogram of mode positions
hist(M2$modes,
     breaks = 1:ncol(X),
     main = "Mode positions",
     xlab = "Class")

# Plot 3 - positions of modelled end-member modes by number of end-members
# Note how scatter in end-member position decreases for the "correct" number
# of modelled end-members (6) and an appropriate weight limit (ca. 0.1).
ii <- order(M2$q, M2$modes)
modes <- t(rbind(M2$modes, M2$q))[ii,]
plot(modes[,1],
     seq(1, nrow(modes)),
     main = "Model overview",
     xlab = "Class",
     ylab = "EM number in model run",
     pch = as.character(modes[,2]),
     cex = 0.7)

# Illustrate mode positions as stem-and-leave-plot, useful as a simple
# check, which mode maxima are consistently fall into which grain-size
# class (useful to define "limits" in robust.EM).
stem(M2$modes, scale = 2)
# }