EMMA protocol.   
This protocoll allows the repitition of the End-member modeling analysis (EMMA) of 
grain size data in R based on the grain size data set of Menges et al. (2019) contained 
in this data publication and the package EMMAgeo 
(needs to be installed, http://micha-dietze.de/pages/emmageo.html). 

This code is part of the following data publication (and should be cited like): 

Menges, J.; Hovius, N.; Andermann, C.; Dietze, M.; Swoboda, C.; Cook, K.; Adhikari, B.; 
Vieth-Hillebrand, A.; Bonnet, S.; Reimann, T.; K., Andreas; Sachse, D. (2019): 
Paleosol-derived data used for the reconstruction of environmental conditions 
during the Holocene in the upper part of the Kali Gandaki valley, Central Nepal. 
GFZ Data Services. http://doi.org/10.5880/GFZ.4.6.2019.001


## load package
library(EMMAgeo)

## read data set
data <- read.table(file = "data_16_17_mean.txt", 
                   header = TRUE, 
                   stringsAsFactors = FALSE)

## build objects from imported data set
X <- as.matrix(data[,1:ncol(data)])
ID <- rownames(data)
units <- as.numeric(x = gsub(x = colnames(data), 
                             pattern = "X",
                             replacement = ""))

## remove columns containing only zeros
units <- units[colSums(X) > 0]
X <- X[,colSums(X) > 0]

## rescale data to 100 %
X <- X / rowSums(X) * 100

## find possible weight transformation limits
l <- get.l(X = X, 
           n = 50)

## find possible numbers of end-members
q <- get.q(X = X, 
           l = l)

## model potential end-members
pdf(file = "plot_model.EM.pdf", width = 10, height = 8)
em <- model.EM(X = X,
               q = q,
               l = l,
               plot = TRUE)
dev.off()

## define limits of end-member mode clusters
limits <- cbind(c(34.5, 45.0, 52.0, 61.0, 66.5),
                c(36.5, 51.0, 54.0, 62.5, 67.5))

units[limits[,1]]
units[limits[,2]]

## model roubst end-members
pdf(file = "plot_robust.EM.pdf", width = 10, height = 10)
E_rob_1 <- robust.EM(em = em,
                     limits = limits,
                     plot = TRUE,
                     l = 0,
                     mc_n = 1000,
                     classunits= units,
                     log="x")
dev.off()