• Elisabeth Dietze is a Post Doc at AWI Potsdam, working on biomarker based fire histories and robust interpretation of sediment proxies.

• Michael Dietze is a Post Doc at GFZ Potsdam, working on the seismic signals emitted by Earth surface processes.

• Your first task: Introduce yourself to the others
  • What is your educational background and career stage?
  • Where do you want to apply EMMA?
  • What do you expect from this course?

• A kaleidoscope of EMMA applications (5 min)
• Background and constraints (15 min)
• The R package EMMAgeo (10 min)
• Preparation is everything – getting data into R (15 min)
• Modelling your first data set with EMMAgeo (10 min)
• Including uncertainty – a protocol for robust EMMA (20 min)
• The quick way – a compact protocol for robust EMMA (10 min)
• Some thoughts on EMMA implications (10 min)

• This course focuses on
  • application cases of EMMA
  • practical experience
  • EMMA with the R package EMMAgeo
• You are looking for…
  • in depth mathematical foundations?
  • in depth insight to R function code?
  • other EMMA approaches?
    • … not here, please see respective references

• What you should be already familiar with
  • R and RStudio
  • Installing and using R packages
  • Writing scripts (and functions)
  • Structuring code and following good practice rules

• This course heavily borrows ideas and concepts from
  • Dietze E. et al. (2012). Sedimentary Geology 243–244. 169–180.
  • Dietze, E. et al. (2014). Climate of the Past 10, 91-106.
  • Dietze, E. & Dietze, M. (in revision). Robust grain-size distribution unmixing with the R package EMMAgeo. E&G Quaternary Science Journal.

The task (Time 5-10 min):

• Think about the previous examples.
• Why do we need to unmix grain-size data?
• What other approaches do you know (to interpret grain-size data)?

We will need the output of this discussion in a few moments, so do not lean back too much.

• Statistical decomposition of all data without (major) prior constraints
• Reduction of redundancy in large data sets and insight into the system
• Quantitative description of a data set by small number of components and their contribution to a sample
• Unmixing of source area, transport pathways, depositional and post-depositional processes

• Classic approaches
  • Method of moments (e.g. Folk and Ward, 1957)
  • Finite mixture modeling (parametric curve fitting)
• Non-parametric approaches:
  • Principle component analysis
  • Factor analysis
  • Cluster analysis
  • End-member modelling or unmixing (EMMA)

The task (Time 2 min):

• Think about the above and grain-size data in general
  • What might be the constraints of grain-size data that should be considered? # Which preconditions must the data fulfill?

• The data is compositional
  • no value smaller than zero (non-negativity)
  • all values of a sample add to constant sum (e.g., 100 vol.-%)
• The subpopulations preserve stable throughout time
  • no change in EM distribution (e.g., gradual fining)
  • no intense post-depositional modifications (e.g., cryoturbation)

• EMMA by Weltje 1997. Math. Geol. 29, 503-549. (providing the main principles and ideas)
• EMMA by Dietze et al., 2012 (the preceder of EMMAgeo in Matlab)
• RECA by Seidel & Hlawitschka 2015. Math. Geosc. 47, 995-1007.
• AnalySize by Paterson & Heslop 2015. G3 16, 4494-4506. (in Matlab)
• BEMMA by Yu et al. 2015. Math. Geosc., 1-19. (Bayesian, in Matlab)

EMMAgeo is the only approach to account for parametrisation uncertainties

• R packages are self contained collections of
  • functions,
  • their documentation and
  • examples/example data
• They are the default way of extending R functionality
• They guarantee that code is organised in a coherent way
• R packages are mainly written by R users that grade into developers (community driven software)

'EMMAgeo' is one out of more than 10,000 packages available on CRAN.

• EMMAgeo is hosted at the Comprehensive R Archive Network CRAN
• Installation with

install.packages("EMMAgeo")

• EMMAgeo is maintained/developed on Github
• Installation (using package devtools) with

devtools::install_github(repo = "coffeemuggler/EMMAgeo", ref = "0.9.6")

• data – examples, most important: example_X, the manually mixed data set
• inst, … – code for the EMMA GUI (use it with EMMA::GUI())
• man – documentation of each function
• R – 24 functions to perform EMMA
• DESCRIPTION – well, a description summary of the package
• NEWS – list of updates introduced with each new version

The task (Time 5 min):

• Install EMMAgeo 0.9.4 from CRAN (install.packages("EMMAgeo"))

• Install EMMAgeo 0.9.6 from Github (devtools::install_github("coffeemuggler/EMMAgeo"), maybe install package devtools in advance (install.packages("devtools")))

• Alternatively: install EMMAgeo 0.9.6 from local source

• Keep things simple!
• R works best with data frames and matrices
• Import data as
  • ASCII file (read.table()),
  • spreadsheet (xlsx::read.xlsx()) or
  • Matlab file (R.matlab::readMat())
• organise samples in rows and grain-size vol-% in columns
• optinally add sample ID, depth, age etc. as extra columns

data <- read.table(file = "dataset_gsd.txt", header = TRUE, sep = "\t", dec = ".", stringsAsFactors = FALSE)
X <- as.matrix(data[,-c(1, 2)])
s <- as.numeric(colnames(data))[-c(1, 2)]
ID <- data[,1]
depth <- data[,2]

## read an ASCII file
data <- read.table(file = "~/Desktop/X.txt", 
                header = TRUE, 
                sep = "\t", 
                stringsAsFactors = FALSE)

## read a MS Excel file (xlsx)
data <- xlsx::read.xlsx(file = "~/Desktop/X.xlsx", 
                     sheetIndex = 1, 
                     header = TRUE, 
                     stringsAsFactors = FALSE)

## extract objects of interest (function substr() is used to remove the letter X)
X <- as.matrix(data[,-1])
s <- as.numeric(substr(x = colnames(data)[-1], 
                       start = 2, 
                       stop = 10))
ID <- data$ID

• The dataset can be explored:
  • graphically (line plots of grain-size distributions, image maps)
  • in tabular form (View(X))
  • by decriptive statistics (dim(X), range(X), rowSums(X))
• The dataset should be explored to identify:
  • outliers and other "strange" samples and features
  • verify correctness of data (right core, age, whatever)
  • get an overview, envision results to be expected

The task (Time: about 5 min):

• Load the package EMMAgeo (library("EMMAgeo"))
• Load the example data set of the EMMAgeo package (data(example_X)).
• Extract the grain-size classes (colnames(X)) as a numeric vector.
• Create an ID vector (numbers from 1 to 100).
• Explore the data set visually.

library(EMMAgeo)

## Loading required package: GPArotation

## Loading required package: limSolve

## Loading required package: caTools

## Loading required package: shiny

## Loading required package: matrixStats

data("example_X")
s <- as.numeric(colnames(X))
ID <- 1:nrow(X)

par(mfcol = c(1, 2), mar = c(4, 4, 1, 1))

gsd_col <- colorRampPalette(colors = c("darkblue", "orange", "grey"))
image(x = s, y = ID, z = t(X), log = "x", col = gsd_col(100))

plot(NA, xlim = range(s), ylim = range(X), log = "x", xlab = "Size", ylab = "Vol.-%")
for(i in 1:nrow(X)) {lines(x = s, y = X[i,], col = gsd_col(100)[i])}

• Four real world end-members
  • EM_1 Overbank deposit
  • EM_2 Loess
  • EM_3 Dune sediment
  • EM_4 Alluvial fan
• Sampled, manually mixed
• Prepared and measured
• Thus, knowledge of end-member shapes and contributions

• EMMA(X, q, l, c, Vqn, classunits, ID, EM.ID, rotation, plot, ...)
  • X Data set, a numeric matrix with samples organised in rows
  • q Number of end-members to model
  • l Weight transformation limit quantile
  • c constant sum scaling factor
  • Vqn optional user-defined unscaled loadings
  • rotation factor rotation type
  • plottoggle plot output

E <- EMMA(X = X, q = 4, l = 0, c = 100, classunits = s, ID = ID, plot = TRUE, log = "x")

• X = X Our example data set, 100 samples, manually mixed
• q = 4 Four target end-members, because we know the data set
• l = 0 No weight limit transformation
• c = 100 All classes sum to 100 %
• classunits = s We provide the real µm grain-size class units
• ID = ID We provide the sample IDs, as well
• plot = TRUE We want to see the default plot output of EMMA
• log = "x" The units are in µm, so the x-axis should be in log scale

str(E)

## List of 18
##  $ loadings: num [1:4, 1:116] 0.001813 0.000891 0.000366 0.00045 0.002269 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "EM1" "EM2" "EM3" "EM4"
##   .. ..$ : chr [1:116] "0.04" "0.0439105" "0.0482033" "0.0529158" ...
##  $ scores  : num [1:100, 1:4] 0.2329 0.0748 0.507 0.2545 0.3121 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:100] "1" "2" "3" "4" ...
##   .. ..$ : NULL
##  $ Vqn     : num [1:4, 1:116] 0.858 0.331 0.179 0.213 0.847 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "EM1" "EM2" "EM3" "EM4"
##   .. ..$ : NULL
##  $ Vqsn    : num [1:4, 1:116] 0.001813 0.000891 0.000366 0.00045 0.002269 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "EM1" "EM2" "EM3" "EM4"
##   .. ..$ : chr [1:116] "0.04" "0.0439105" "0.0482033" "0.0529158" ...
##  $ Mqs     : num [1:100, 1:4] 0.2329 0.0748 0.507 0.2545 0.3121 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:100] "1" "2" "3" "4" ...
##   .. ..$ : NULL
##  $ Xm      : num [1:100, 1:116] 0.000978 0.000494 0.001354 0.000928 0.000984 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:100] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:116] "0.04" "0.0439105" "0.0482033" "0.0529158" ...
##  $ modes   : num [1:4] 14.3 39.8 339.9 493.6
##  $ Mqs.var : num [1:4] 23.5 18.7 30 27.8
##  $ Em      : Named num [1:100] 0.0741 0.0506 0.0917 0.0782 0.0893 ...
##   ..- attr(*, "names")= chr [1:100] "1" "2" "3" "4" ...
##  $ En      : Named num [1:100] 0.0741 0.0506 0.0917 0.0782 0.0893 ...
##   ..- attr(*, "names")= chr [1:100] "1" "2" "3" "4" ...
##  $ RMSEm   : num 0.135
##  $ RMSEn   : num 0.135
##  $ Rm      : Named num [1:100] 0.982 0.998 0.961 0.984 0.977 ...
##   ..- attr(*, "names")= chr [1:100] "1" "2" "3" "4" ...
##  $ Rn      : Named num [1:116] 0.912 0.912 0.923 0.932 0.931 ...
##   ..- attr(*, "names")= chr [1:116] "0.04" "0.0439105" "0.0482033" "0.0529158" ...
##  $ mRm     : num 0.967
##  $ mRn     : num 0.887
##  $ mRt     : num 0.927
##  $ ol      : num 0

• loadings (Vqsn) Normalised rescaled end-member loadings.
• scores (Mqs) Rescaled end-member scores.
• Xm Modelled data.
• Vqn Normalised end-member loadings.
• modes Mode class of end-member loadings.
• Mqs.var Explained variance of end-members
• Rm, Rn Row- and col-wise explained variance.
• Em, En Absolute row- and col-wise model error.
• ol Number of overlapping end-members.

The task (time about 5 minutes):

• Perform EMMA yourself, using the example data set.

• Try other parameter settings (see help for the function: ?EMMA).

• If you are fast, try making new plots out of the numeric output (e.g., cumulative grain-size curves, line plots of scores, ratio plots of scores)

Question: How do you know the model you use is a proper one? How to judge meaningful end-members? Which side effects do you observe?

• The big question: Which model is the right one.

• Or: how many end-members q do we need? And, what about l?

• Looks like we need more than deterministic EMMA.

• And we need a structured way to robustly model a data set.

• Define a range for the weight transformation limit
  • not less than 0 (no effect), not more than 0.5 (median)
• Define a range for the number of end-members
  • more than 1 (no variability), less than classes (redundancy)
• Run all models possible in this parameter space
• Find and extract persistent loadings, independent of parameters
• Calculate statistics of extracted robust loadings
• Model scores based on uncertain loadings

• With EMMA, we provide two generic protocols to consistenty approach this problem
  • Extended protocol: 11 steps, more profound user interaction
  • Compact protocol: 5 steps, limited user interaction
• Each protocol has its pros and cons.
• Good to start with extended protocol to understand principles
• If necessary, seamless switch between protocols is possible at many stages.

• 1 Define l_min, in most cases zero
• 2 Define l_max, as high as possible to still allow eigenspace calculation
• 3 Define l sequence
• 4 Define q_min, as many as needed to reach R² threshold in FA
• 5 Define q_max, first local EMMA R² maximum, fixed R², or fixed q
• 6 Define parameter matrix (q_min to q_max for each l)

• 7 Run all possible models
• 8 Define mode limits of robust end-members
• 9 Extract robust loadings, calculate their statistics
• 10 Optimise model for one solution of l (l_opt)
• 11 Calculate scores uncertainties

So far, so good? Now let us dive into the details…

• l_min will usually be zero, unless there is good reason.

l_min <- 0

• l_max can be up to 0.5. In reality, rescaling the data too much will violate the eigenspace calculation (eigen()). The function test.l() can be used to test if a given value is still valid. Values close to l_max yield unrealistic EMMA models. Thus, setting l_max to 95 % of the true max or lower prevents this.

l_max <- test.l(X = X, l = seq(from = 0,
                               to = 0.5,
                               length.out = 20))$l.max

l_max <- 0.6 * l_max

• Alternatively, the function test.l.max() approximates l_max in a loop
• Might be the more convenient way

l_max <- 0.6 * test.l.max(X = X, n = 20)

print(l_max)

## [1] 0.237

• The weight transformation limit l will be a sequence of values from l_minto l_max. The longer the sequence, the more models can be run (better estimate, more computation time).

l <- seq(from = l_min, to = l_max, length.out = 20)

• To exhibit any variability, a data set mus at least contain 2 end-members
• The more end-members a model uses, the better is the model result
• Testing the R² between input and output for different number of end-members provides a reasonable estimate of the minimum q
• Usually, a minimum R² value of 0.95 is used as threshold to define the minimum number of end-members to include
• For each weight transformation limit, there is a nother R². Thus, the test mus be run for all elements of land all possible q.

q_min <- test.factors(X = X, l = l, r.min = 0.95, plot = TRUE, xlim = c(2, 10), ylim = c(0.9, 1.0))$q.min

print(q_min)

##  [1] 2 3 4 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

• In EMMA, there is a tendency of decreasing model quality (R²) when more and more end-members are included - beyond a meaningful number.
• Thus, the first local maximum can be used to set the maximum number of end-members q_max. Of course, for each element of l, there is a different q-max.
• Alternatively, a fixed number of q_maxcan be used, or a given R² value of the model.

• Plot of the effect of decreasing R² for too many end-members (l = 0):

• The function test.parameters performs these tests for all l and q. It can run a loooong time.
• The output can be used to extract the maximum number of end-members q_max.

• If you have time, try the other plot types, apart from "mRt".

• See the example on the next slide, and how with increasing l also NA cases emerge.

TP <- test.parameters(X = X, q = 2:20, l = l, colour = rainbow(n = length(l)), plot = "mRt")

q_max <- TP$q.max

• Now, with l, q_min and q_max at hand, these values can be fused into a parameter matrix
• A few things must be secured:
  • q_min must be at least 2
  • q_max cannot be smaller than q_min
  • no NA values are allowed

q_min <- ifelse(test = q_min < 2, yes = 2, no = q_min)

q_max <- ifelse(test = q_max < q_min, yes = NA, no = q_max)

P <- cbind(q_min, q_max, l)

P <- na.exclude(P)

• Now, with X and P all possible models, in agreement with the parameter space can be evaulated to inspect the outcomes.
• It is also possible to set some rejection criteria (things that are not allowed to make a valid model):
  • overlapping modes
  • R² below a minimum threshold
• The function that performs this step is test.robustness(). All possible end-member loadings will be plotted.

EM.pot <- test.robustness(X = X, P = P, plot = TRUE)

• The tricky part is to define the limits for potential robust end-members.
• The most important and relevant part in the protocol.
• Use the graphical output
  • plot all loadings over each other
  • plot a histogram of mode positions
  • create a stem-and-leaf plot
• Define for each end-member the lower and upper mode class.

plot(NA, xlim = c(1, length(s)), ylim = range(EM.pot$loadings))
for(i in 1:nrow(EM.pot$loadings)) {lines(EM.pot$loadings[i,], col = adjustcolor(col = 1, alpha.f = 0.1))}

hist(x = EM.pot$modes, breaks = 90)

stem(x = EM.pot$modes)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##    6 | 222233334444444444444444444
##    6 | 
##    7 | 2222333333333333333333333333334444
##    7 | 5555555555555555555555555666
##    8 | 
##    8 | 
##    9 | 44
##    9 | 77777777777777777777777777777888888888888888888888888999999999999999
##   10 | 00000000000000000000000000000000111111112222222222222222222222222222+4
##   10 | 5555556
##   11 | 222

• define the limits as a column-wise organised matrix, first column holds lower limits, second column holds upper limits. One line for each robust end-member.

limits <- cbind(c(63, 75, 98, 102),
                c(64, 77, 99, 103))

• Now, all end-members can be extracted, combined and used to calculate their distribution statistics with robust.loadings()

loadings.robust <- robust.loadings(em = EM.pot, limits = limits, plot = TRUE)

• The weight transformation limit that yields the optimal robust model usually differs from the one of a deterministic EMMA.
• The function get.l.opt() can be used to identify the model with the best performance.
• Different performance criteria are possible (qualitiy = "mRm", "mRn", "mRt", "mEm", "mEn" and "mEt")
• Optimisation can be done with mean (loadings.robust$Vqn$mean) or median (loadings.robust$Vqn$median) loadings.

l_opt <- get.l.opt(X = X, l = l, quality = "mRt", Vqn = loadings.robust$Vqn$mean, plot = TRUE)

• The only thing that is missing are the counterpart of the robust loadings: the robust scores
• Drawing random samples of loadings and modelling the data set yields different scores
• The function robust.scores() does this job in a Monte Carlo environment.
• For tests it is good to set the argument mc_n (number om MC rounds) to a low value to increase speed.

scores_robust <- robust.scores(loadings = loadings.robust, l = l_opt, mc_n = 20, plot = TRUE)

• This was it! Let us wrap up
  • The extended protocol yields robust loadings and scores with uncertainty estimates
  • It creates a parameter space with meaningful constraints on potential parameter ranges (land q)
  • The most important step is defining the number and limits of end-members, emerging from the parameter space
  • The ten steps allow full control of all relevant modelling steps but the operation takes a lot of decisions and time

• 1 Get vector of weight transformation limits (l values)
• 2 Get range of end-members (q values)
• 3 Run all possible models
• 4 Get/set limits
• 5 Model robust end-members

• 1 Get vector of weight transformation limits (Steps 1–3)
• 2 Get range of end-members (Steps 4–6)
• 3 Run all possible models (Step 7)
• 4 Get/set limits (Step 8)
• 5 Model robust end-members (Steps 9–11)

• Assuming l_minis zero (otherwise change argument min) and l_max is 95 % of the true maximum (otherwise change argument max), and that 10 values are enough (otherwise change argument n), the function get.l() creates a vector of weight transformation limits:

l <- get.l(X = X, max = 0.6, n = 20)
print(l)

##  [1] 0.00000000 0.01247368 0.02494737 0.03742105 0.04989474 0.06236842
##  [7] 0.07484211 0.08731579 0.09978947 0.11226316 0.12473684 0.13721053
## [13] 0.14968421 0.16215789 0.17463158 0.18710526 0.19957895 0.21205263
## [19] 0.22452632 0.23700000
## attr(,"class")
## [1] "EMMAgeo_l"

• Using the l vector, and following the criteria of the extended protocol for q_minand q_max, the function get.q() generates q values for each element of l:

q <- get.q(X = X, l = l)

## [1] "1 cases of q_min > q_max removed: 13."

print(q)

##                    q_min q_max
## 0                      2     6
## 0.0124736842069489     2     6
## 0.0249473684138978     2     5
## 0.0374210526208467     3     5
## 0.0498947368277956     3     5
## 0.0623684210347445     3     5
## 0.0748421052416934     3     5
## 0.0873157894486423     3     5
## 0.0997894736555912     3     5
## 0.11226315786254       3     5
## 0.124736842069489      3     5
## 0.137210526276438      3     6
## 0.149684210483387      4     2
## 0.162157894690336      3     6
## 0.174631578897285      3     6
## 0.187105263104233      2     6
## 0.199578947311182      3     6
## 0.212052631518131      3     6
## 0.22452631572508       3     3
## 0.236999999932029      3     3
## attr(,"class")
## [1] "EMMAgeo_q"

• Using land q, all models of this parameter space are evaluated and the output is plotted.
• The plot shows all loadings, colour-coded by the number of end-members with which the model was run. It is overlain by the mode positions (dots) and a kernel density estimate of the modes.

em.pot <- model.EM(X = X, q = q, l = l, plot = TRUE)

• Like in the extended protocol, setting the limits of robust end-members is the most delicate step.
• In easy cases there is the function get.limits() that uses a kernel density estimate.
• Alternatively, the limits can/must be set manually:

limits <- cbind(c(63, 75, 98, 102),
                c(64, 77, 99, 103))

• With the limits at hand, it is possibe to evaluate both, robust loadings and scores in one step.
• The graphical output is similar to the deterministic version.
• The numeric output is more extensive, containing statistic information about loadings and scores.

em.rob <- robust.EM(em = em.pot, limits = limits, plot = TRUE, mc_n = 20)

## [1] "Parameter l missing! Set to 'mRt' by default"

Online feedback form: https://forms.gle/gsSS65xnXJGGz8LP9

The course materials are available at

www.micha-dietze.de/pages/r_courses.html

The EMMAgeo website with information, tutorials, news

www.micha-dietze.de/pages/emmageo.html

If you use EMMAgeo, please cite the following article

Dietze, E. & Dietze, M. (in revision). Robust grain-size distribution unmixing with the R package EMMAgeo. E&G Quaternary Science Journal.