Multivariate curve resolution using Alternating Least Squares
mcralls allows to resolve spectroscopic data to linear combination of individual spectra and contributions using the alternating least squares (ALS) algorithm with constraints.
mcrals( x, ncomp, cont.constraints = list(), spec.constraints = list(), spec.ini = matrix(runif(ncol(x) * ncomp), ncol(x), ncomp), cont.forced = matrix(NA, nrow(x), ncomp), spec.forced = matrix(NA, ncol(x), ncomp), cont.solver = mcrals.nnls, spec.solver = mcrals.nnls, exclrows = NULL, exclcols = NULL, verbose = FALSE, max.niter = 100, tol = 10^-6, info = "" )
x: spectra of mixtures (matrix or data frame).ncomp: number of components to calculate.cont.constraints: a list with constraints to be applied to contributions (see details).spec.constraints: a list with constraints to be applied to spectra (see details).spec.ini: a matrix with initial estimation of the pure components spectra.cont.forced: a matrix which allows to force some of the concentration values (see details).spec.forced: a matrix which allows to force some of the spectra values (see details).cont.solver: which function to use as a solver for resolving of pure components contributions (see detials).spec.solver: which function to use as a solver for resolving of pure components spectra (see detials).exclrows: rows to be excluded from calculations (numbers, names or vector with logical values).exclcols: columns to be excluded from calculations (numbers, names or vector with logical values).verbose: logical, if TRUE information about every iteration will be shown.max.niter: maximum number of iterations.tol: tolerance, when explained variance change is smaller than this value, iterations stop.info: a short text with description of the case (optional).Returns an object of mcrpure class with the following fields: - resspec: matrix with resolved spectra.
rescont: matrix with resolved contributions.
cont.constraints: list with contribution constraints provided by user.
spec.constraints: list with spectra constraints provided by user.
expvar: vector with explained variance for each component (in percent).
cumexpvar: vector with cumulative explained variance for each component (in percent).
ncomp: number of resolved components
max.niter: maximum number of iterations
info: information about the model, provided by user when build the model.
More details and examples can be found in the Bookdown tutorial.
The method implements the iterative ALS algorithm, where, at each iteration, spectra and contributions of each chemical component are estimated and then a set of constraints is applied to each. The method is well described in [1, 2].
The method assumes that the spectra (D) is a linear combination of pure components spectra (S) and pure component concentrations (C):
D = CS' + E
So the task is to get C and S by knowing D. In order to do that you need to provide:
Constraints for spectra and contributions. The constraints should be provided as a list with name of the constraint and all necessary parameters. You can see which constraints and parameters are currently supported by running constraintList(). See the code examples below or a Bookdown tutorial for more details.
Initial estimation of the pure components spectra, S. By default method uses a matrix with random numbers but you can provide a better guess (for example by running mcrpure) as a first step.
Which solver to use for resolving spectra and concentrations. There are two built in solvers: mcrals.nnls (default) and mcrals.ols. The first implements non-negative least squares method which gives non-negative (thus physically meaningful) solutions. The second is ordinary least squares and if you want to get non-negative spectra and/or contributions in this case you need to provide a non-negativity constraint.
The algorithm iteratively resolves C and S and checks how well CS' is to D. The iterations stop either when number exceeds value in max.niter or when improvements (difference between explained variance on current and previous steps) is smaller than tol value.
Parameters cont.force and spec.force allows you to force some parts of the contributions or the spectra to be equal to particular pre-defined values. In this case you need to provide the parameters (or just one of them) in form of a matrix. For example cont.force
should have as many rows as many you have in the original spectral data x and as many columns as many pure components you want to resolve. Feel all values of this matrix with NA and the values you want to force with real numbers. For example if you know that in the first measurement concentration of 2 and 3 components was zero, set the corresponding values of cont.force to zero. See also the last case in the examples section.
library(mdatools) # resolve mixture of carbonhydrates Raman spectra data(carbs) # define constraints for contributions cc <- list( constraint("nonneg") ) # define constraints for spectra cs <- list( constraint("nonneg"), constraint("norm", params = list(type = "area")) ) # because by default initial approximation is made by using random numbers # we need to seed the generator in order to get reproducable results set.seed(6) # run ALS m <- mcrals(carbs$D, ncomp = 3, cont.constraints = cc, spec.constraints = cs) summary(m) # plot cumulative and individual explained variance par(mfrow = c(1, 2)) plotVariance(m) plotCumVariance(m) # plot resolved spectra (all of them or individually) par(mfrow = c(2, 1)) plotSpectra(m) plotSpectra(m, comp = 2:3) # plot resolved contributions (all of them or individually) par(mfrow = c(2, 1)) plotContributions(m) plotContributions(m, comp = 2:3) # of course you can do this manually as well, e.g. show original # and resolved spectra par(mfrow = c(1, 1)) mdaplotg( list( "original" = prep.norm(carbs$D, "area"), "resolved" = prep.norm(mda.subset(mda.t(m$resspec), 1), "area") ), col = c("gray", "red"), type = "l" ) # in case if you have reference spectra of components you can compare them with # the resolved ones: par(mfrow = c(3, 1)) for (i in 1:3) { mdaplotg( list( "pure" = prep.norm(mda.subset(mda.t(carbs$S), 1), "area"), "resolved" = prep.norm(mda.subset(mda.t(m$resspec), 1), "area") ), col = c("gray", "red"), type = "l", lwd = c(3, 1) ) } # This example shows how to force some of the contribution values # First of all we combine the matrix with mixtures and the pure spectra, so the pure # spectra are on top of the combined matrix Dplus <- mda.rbind(mda.t(carbs$S), carbs$D) # since we know that concentration of C2 and C3 is zero in the first row (it is a pure # spectrum of first component), we can force them to be zero in the optimization procedure. # Similarly we can do this for second and third rows. cont.forced <- matrix(NA, nrow(Dplus), 3) cont.forced[1, ] <- c(NA, 0, 0) cont.forced[2, ] <- c(0, NA, 0) cont.forced[3, ] <- c(0, 0, NA) m <- mcrals(Dplus, 3, cont.forced = cont.forced, cont.constraints = cc, spec.constraints = cs) plot(m) # See bookdown tutorial for more details.
Methods for mcrals objects:
summary.mcrals | shows some statistics for the case. |
predict.mcrals | computes contributions by projection of new spectra to the resolved ones. |
Plotting methods for mcrals objects:
plotSpectra.mcr | shows plot with resolved spectra. |
plotContributions.mcr | shows plot with resolved contributions. |
plotVariance.mcr | shows plot with explained variance. |
plotCumVariance.mcr | shows plot with cumulative explained variance. |
Sergey Kucheryavskiy (svkucheryavski@gmail.com)