Carry out PARAFAC analysis of a Spectra2D object.
Function parafac from multiway is used.
For large data sets, computational time may be long enough that
it might desirable to run in batch mode and possibly use parallel processing.
pfacSpectra2D(spectra, parallel = FALSE, setup = FALSE, nfac = 2, ...)
Arguments
| spectra | An object of S3 class |
|---|---|
| parallel | Logical. Should parallel processing be used?
Unless you love waiting, you should use parallel processing for larger data sets.
If you are working on a shared machine and/or another process (created by you or
another user) might also try to access all or some of the cores in your CPU,
you should be careful to avoid hogging the cores.
|
| setup | Logical. If |
| nfac | Integer. The number of factors/components to compute. |
| ... | Additional parameters to be passed to function |
Value
An object of class pfac and parafac, modified to include a list
element called $method which is parafac.
Warning
To get reproducible results you will need to set.seed(). See the example.
References
R. Bro "PARAFAC. Tutorial and applications" Chemometrics and Intelligent Laboratory Systems vol. 38 pgs. 149-171 (1997).
A. Smilde, R. Bro and P. Geladi "Multi-way Analysis: Applications in the Chemical Sciences" Wiley (2004).
See also
For other data reduction methods for Spectra2D objects, see
miaSpectra2D and popSpectra2D.
Author
Bryan A. Hanson, DePauw University.
Examples
#> | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100%# plotScores uses ggplot2 graphics p1 <- plotScores(MUD1, res, leg.loc = "topright", ellipse = "cls")#> Error in (function (spectra, so, pcs = c(1, 2), ellipse = "none", tol = "none", use.sym = FALSE, leg.loc = "topright", ...) { args <- as.list(match.call())[-1] if (length(pcs) != 2) stop("You must choose exactly two PC's to plot") case <- NULL if (inherits(spectra, "Spectra")) case <- "PCA" if (inherits(spectra, "Spectra2D")) case <- "MIA" if (inherits(so, "pop")) case <- "PCA" if (is.null(case)) stop("Could not reconcile data object and scores object.") if ((case == "MIA") && (use.sym)) stop("ChemoSpec2D does not support use.sym") chkSpectra(spectra) if (case == "PCA") DF <- data.frame(so$x[, pcs], group = spectra$groups) if (case == "MIA") DF <- data.frame(so$C[, pcs], group = spectra$groups) GRPS <- dlply(DF, "group", subset, select = c(1, 2)) if ((ellipse == "cls") || (ellipse == "rob") || (ellipse == "both")) { gr <- sumGroups(spectra) for (n in 1:length(gr$group)) { if (gr$no.[n] == 1) message("Group ", gr$group[n], "\n\thas only 1 member (no ellipse possible)") if (gr$no.[n] == 2) message("Group ", gr$group[n], "\n\thas only 2 members (no ellipse possible)") if (gr$no.[n] == 3) message("Group ", gr$group[n], "\n\thas only 3 members (ellipse not drawn)") } idx <- which(gr$no. > 3) gr <- gr[idx, ] ELL <- llply(GRPS[idx], .computeEllipses) } go <- chkGraphicsOpt() if (go == "base") { if ((ellipse == "cls") || (ellipse == "rob") || (ellipse == "both")) { x.scores <- range(llply(GRPS, subset, select = 1)) y.scores <- range(llply(GRPS, subset, select = 2)) x.ell <- range(llply(ELL, function(x) { range(x[1]) })) y.ell <- range(llply(ELL, function(x) { range(x[2]) })) x.ell.r <- range(llply(ELL, function(x) { range(x[4]) })) y.ell.r <- range(llply(ELL, function(x) { range(x[5]) })) x.all <- range(x.scores, x.ell, x.ell.r) x.all <- x.all + diff(x.all) * 0.05 * c(-1, 1.15) y.all <- range(y.scores, y.ell, y.ell.r) y.all <- y.all + diff(x.all) * 0.05 * c(-1, 1.15) } if (ellipse == "none") { x.scores <- range(llply(GRPS, subset, select = 1)) y.scores <- range(llply(GRPS, subset, select = 2)) x.all <- range(x.scores) + diff(range(x.scores)) * 0.05 * c(-1, 1.15) y.all <- range(y.scores) + diff(range(y.scores)) * 0.05 * c(-1, 1.15) } dPargs <- list(PCs = DF[, 1:2], spectra = spectra, case = case, use.sym = use.sym, ... = ...) if (!"xlim" %in% names(args)) dPargs <- c(dPargs, list(xlim = x.all)) if (!"ylim" %in% names(args)) dPargs <- c(dPargs, list(ylim = y.all)) do.call(.drawPoints, dPargs) if ((ellipse == "cls") | (ellipse == "rob") | (ellipse == "both")) .drawEllipses(ELL, gr, ellipse, use.sym, ...) if (case == "PCA") { .addMethod(so) if (all(leg.loc != "none")) { leg.loc <- .prepLegendCoords(go, leg.loc, x.all[1], x.all[2], y.all[1], y.all[2]) .addLegend(spectra, leg.loc, use.sym, bty = "n") } .addEllipseInfo(ellipse) } if (case == "MIA") { if (all(leg.loc != "none")) { leg.loc <- .prepLegendCoords(go, leg.loc, x.all[1], x.all[2], y.all[1], y.all[2]) .addLegend(spectra, leg.loc, use.sym = FALSE, bty = "n") } .addEllipseInfo(ellipse) } if (tol != "none") .labelExtremes(DF[, 1:2], spectra$names, tol) } if ((go == "ggplot2") || (go == "plotly")) { args <- as.list(match.call()[-1]) xlab <- eval(args$xlab) ylab <- eval(args$ylab) .chkReqGraphicsPkgs("ggplot2") x <- y <- name <- label <- NULL if (case == "PCA") { if (!use.sym) { p <- ggplot(DF) + geom_point(aes_string(x = colnames(DF)[1], y = colnames(DF)[2]), color = spectra$colors, shape = 20, size = 3) + labs(x = xlab, y = ylab) } if (use.sym) { p <- ggplot(DF) + geom_point(aes_string(x = colnames(DF)[1], y = colnames(DF)[2]), color = "black", shape = spectra$sym) + labs(x = xlab, y = ylab) } if (go == "ggplot2") { p <- p + .ggAnnotate(so$method, x = 0.05, y = 0.98, just = "left", gp = gpar(fontsize = 8)) } } if (case == "MIA") { p <- ggplot(DF) + geom_point(aes_string(x = colnames(DF)[1], y = colnames(DF)[2]), color = spectra$colors, shape = 20, size = 3) } p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) if ((ellipse == "cls") | (ellipse == "rob") | (ellipse == "both")) { .ggPrepEllipseCoords <- function(data) { df <- data.frame(x = numeric(), y = numeric(), name = character()) for (i in 1:length(data)) { x <- c() y <- c() name <- c() total <- length(data[[i]])/2 for (j in 1:total) { x <- c(x, data[[i]][j, 1]) y <- c(y, data[[i]][j, 2]) } name <- rep(names(data)[i], total) temp <- data.frame(x, y, name) df <- rbind(df, temp) } return(df) } } if (ellipse == "cls") { cls.coords <- llply(ELL, function(x) { x[1:2] }) cls.coords <- llply(cls.coords, function(x) { do.call(cbind, x) }) df.cls <- .ggPrepEllipseCoords(cls.coords) if (!use.sym) { p <- p + geom_path(data = df.cls, aes(x = x, y = y, color = name), linetype = 2) + scale_color_manual(values = gr$color) } if (use.sym) { color.black <- rep("black", length(gr$color)) p <- p + geom_path(data = df.cls, aes(x = x, y = y, color = name), linetype = 2) + scale_color_manual(values = color.black) } if (go == "ggplot2") { p <- p + .ggAnnotate("cls") } } if (ellipse == "rob") { rob.coords <- llply(ELL, function(x) { x[4:5] }) rob.coords <- llply(rob.coords, function(x) { do.call(cbind, x) }) df.rob <- .ggPrepEllipseCoords(rob.coords) if (!use.sym) { p <- p + geom_path(data = df.rob, aes(x = x, y = y, color = name)) + scale_color_manual(values = gr$color) } if (use.sym) { color.black <- rep("black", length(gr$color)) p <- p + geom_path(data = df.rob, aes(x = x, y = y, color = name)) + scale_color_manual(values = color.black) } if (go == "ggplot2") { p <- p + .ggAnnotate("rob") } } if (ellipse == "both") { cls.coords <- llply(ELL, function(x) { x[1:2] }) cls.coords <- llply(cls.coords, function(x) { do.call(cbind, x) }) df.cls <- .ggPrepEllipseCoords(cls.coords) rob.coords <- llply(ELL, function(x) { x[4:5] }) rob.coords <- llply(rob.coords, function(x) { do.call(cbind, x) }) df.rob <- .ggPrepEllipseCoords(rob.coords) if (!use.sym) { lines <- rep(2, length(gr$color)) p <- p + geom_path(data = df.cls, aes(x = x, y = y, color = name), linetype = 2) p <- p + geom_path(data = df.rob, aes(x = x, y = y, color = name)) + scale_color_manual(values = gr$color) } if (use.sym) { color.black <- rep("black", length(gr$color)) p <- p + geom_path(data = df.cls, aes(x = x, y = y, color = name)) p <- p + geom_path(data = df.rob, aes(x = x, y = y, color = name), linetype = 2) + scale_color_manual(values = color.black) } if (go == "ggplot2") { p <- p + .ggAnnotate("both") } } if (go == "ggplot2") { if (tol != "none") { CoordList <- .getExtremeCoords(DF[, 1:2], spectra$names, tol) df <- data.frame(x = CoordList$x, y = CoordList$y, label = CoordList$l) p <- p + .ggRepel(df) } p <- p + theme(legend.position = "none") if (all(leg.loc != "none")) { p <- p + .ggAddLegend(spectra, use.sym, leg.loc) } return(p) } else { .chkReqGraphicsPkgs("plotly") p <- ggplotly(p, tooltip = c(colnames(DF[1]), colnames(DF[2]), "name")) if (tol != "none") { CoordList <- .getExtremeCoords(DF[, 1:2], spectra$names, tol) df <- data.frame(x = CoordList$x, y = CoordList$y, label = CoordList$l) p <- p %>% add_annotations(x = df$x, y = df$y, text = df$label, xref = "x", yref = "y", showarrow = TRUE, arrowhead = 4, arrowsize = 0.5, ax = 40, ay = -25, font = list(size = 12)) } p <- p %>% layout(showlegend = FALSE) return(p) } }})(spectra = MUD1, so = res, ellipse = "cls", leg.loc = "topright", xlab = "Component 1", ylab = "Component 2", use.sym = FALSE): object 'res' not found#> Error in eval(expr, envir, enclos): object 'p1' not foundp1#> Error in eval(expr, envir, enclos): object 'p1' not found# plotLoadings2D uses base graphics res1 <- plotLoadings2D(MUD1, res, load_lvls = c(1, 5, 10, 15, 25), main = "PARAFAC Comp. 1 Loadings")res2 <- plotLoadings2D(MUD1, res, load_lvls = c(1, 5, 10, 15, 25), ref = 2, ref_lvls = seq(5, 35, 5), ref_cols = rep("black", 7), main = "PARAFAC Comp. 1 Loadings + Ref. Spectrum")# Selection of loading matrix levels can be aided by the following # Use res1$names to find the index of the loadings inspectLvls(res1, which = 11, ylim = c(0, 50), main = "Histogram of Loadings Matrix")
