Nothing
#' Model all possible end-member scenarios
#'
#' This function takes a definition of weight transformation
#' limits and corresponding minimum and maximum numbers of end-members to
#' model all end-member scenarios in accordance with these parameters. Based
#' on the output the user can decide on robust end-members.
#'
#' The plot output is an overlay of several data. The coloured lines in the
#' background are end-member loadings (number noted in the plot title),
#' resulting from all possible model scenarios. If \code{col.q == TRUE} they
#' are coloured according to the number of end-members with which the model
#' was generated. This colour scheme allows to depict end-members that emerge
#' for model realisations with specific number of end-members. The thick
#' black line is a kernel density estimate curve, generated from the mode
#' positions of all end-members. The kernel bandwidth is set to 1 percent of
#' the number of grain-size classes of the input data set, which gave useful
#' results for most of our test data sets. The cumulaitve dot-line-plot is a
#' further visualisation of end-member mode positions. The function is a
#' modified wrapper function for the function \code{test.robustness()}.
#'
#'
#' @param X \code{Numeric} matrix, input data set with m samples (rows)
#' and n variables (columns).
#'
#' @param q \code{Numeric} matrix, definitions of minimum and maximum number
#' of end-members (cf. \code{get.q()}), required.
#'
#' @param l \code{Numeric} vector, weight transformation limit values,
#' corresponding to the matrix q, required.
#'
#' @param classunits \code{Numeric} vector, optional class units
#' (e.g. micrometers or phi-units) of the same length as columns of \code{X}.
#'
#' @param plot \code{Logical} scalar, option to plot the results (cf.
#' details for explanations), default is \code{TRUE}.
#'
#' @param col.q \code{Logical} scalar, option to colour end-member loadings by
#' the number of end-members which were used to create the model realisation,
#' default is \code{TRUE}.
#'
#' @param bw \code{Numeric} scalar, optional manual setting of the kde
#' bandwidth. By default, bw is calculated as 1 percent of the number of
#' grain-size classes.
#'
#' @param \dots Further arguments passed to the function.
#'
#' @return \code{List} object with all modelled end-members, each described by
#' input parameters, mode position, quality measures and value distributions.
#'
#' @author Michael Dietze, Elisabeth Dietze
#' @seealso \code{\link{EMMA}}, \code{\link{test.l.max}}
#' @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.
#' @keywords EMMA
#' @examples
#'
#' ## load example data set
#' data(example_X)
#'
#' ## define input parameters
#' l <- c(0, 0.05, 0.10)
#' q <- cbind(c(2, 2, 3), c(5, 6, 4))
#'
#' ## infer l-vector
#' em_pot <- model.EM(X = X, q = q, l = l)
#'
#' @export model.EM
model.EM <- function(
X,
q,
l,
classunits,
plot = TRUE,
col.q = TRUE,
bw,
...
) {
## check/set class units
if(missing(classunits) == TRUE) {
classunits <- seq(from = 1, to = ncol(X))
}
## check for l vs. lw
if("lw" %in% names(list(...))) {
stop('Parameter "lw" is depreciated. Use "l" instead.')
}
## check7set log option
if("log" %in% names(list(...))) {
log_plot <- list(...)$log
} else {
log_plot <- ""
}
## check/set l
if(class(q)[1] == "EMMAgeo_q") {
l <- as.numeric(rownames(q))
}
## check input data
if(length(l) != nrow(q)) {
stop("l and q are not of identical length!")
}
## check/set bw
if(missing(bw) == TRUE) {
bw <- ncol(X) / 100
}
## create P-matrix
P <- cbind(q, l)
## run test.robustness()
em <- test.robustness(X = X, P = P, ...)
## assign class units
em$modes_classunits <- classunits[em$modes]
## assign plot colour
if(col.q == TRUE) {
plot_col <- em$q - min(em$q) + 1
} else {
plot_col <- rep(x = "grey70",
times = nrow(em$loadings))
}
## optionally, plot output
if(plot == TRUE) {
## create plot object
em_plot <- em
## set negative values to zero
em_plot$loadings[em_plot$loadings < 0] <- 0
em_plot$Vqsn[em_plot$Vqsn < 0] <- 0
## get old margins
mar_old <- par()$mar
## set new margins
par(mar = c(5.1, 4.1, 4.1, 4.1))
## create empty, scaled plot
plot(NA,
xlim = range(classunits),
ylim = c(min(em_plot$loadings),
max(em_plot$loadings) * 1.1),
xlab = "Class",
ylab = "Contribution",
log = log_plot,
main = paste("Loadings (n = ",
nrow(em_plot$loadings),
")",
sep = ""))
## add all potential loadings
for(i in 1:nrow(em_plot$loadings)) {
lines(x = classunits,
y = em_plot$loadings[i,],
col = adjustcolor(col = plot_col[i],
alpha.f = 0.3))
}
## add cumulate mode position plot
par(new = TRUE)
plot(x = sort(em_plot$modes),
y = 1:length(em_plot$modes),
xlim = c(1,ncol(em_plot$loadings)),
type = "b",
lwd = 2,
col = plot_col,
ann = FALSE,
axes = FALSE)
## add mode position KDE plot
par(new = TRUE)
kde <- density(x = em_plot$modes,
from = 1,
to = ncol(X),
bw = bw)
plot(kde,
ylim = c(0, max(kde$y) * 1.2),
lwd = 2,
ann = FALSE,
axes = FALSE)
axis(side = 4)
mtext(text = "Density",
side = 4,
line = 2.5)
if(col.q == TRUE) {
legend(x = "top",
legend = paste("q = ", seq(from = min(em_plot$q),
to = max(em_plot$q))),
text.col = sort(x = unique(x = plot_col)),
horiz = TRUE,
box.lty = 0)
}
## create empty plot, scaled to initial plot
par(new = TRUE)
plot(NA,
xlim = range(classunits),
ylim = c(min(em_plot$loadings),
max(em_plot$loadings) * 1.1),
log = log_plot,
ann = FALSE,
axes = FALSE)
## restore old margins
par(mar = mar_old)
}
## set class
class(em) <- "EMMAgeo_empot"
## return output
return(em)
}
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.