# PCA on Coe objects
# 1. PCA calculation and builder --------------------------
#' Principal component analysis on Coe objects
#'
#' Performs a PCA on \link{Coe} objects, using \link{prcomp}.
#'
#' By default, methods on \link{Coe} object do not scale the input data but center them.
#' There is also a generic method (eg for traditional morphometrics) that centers and scales data.
#' @aliases PCA
#' @rdname PCA
#' @param x a \link{Coe} object or an appropriate object (eg \link{prcomp}) for \code{as_PCA}
#' @param fac any factor or data.frame to be passed to \code{as_PCA} and for use with \link{plot.PCA}
#' @param scale. logical whether to scale the input data
#' @param center logical whether to center the input data
#' @return a 'PCA' object on which to apply \link{plot.PCA}, among others. This list has several
#' components, most of them inherited from the \code{prcomp} object:
#' \enumerate{
#' \item \code{sdev} the standard deviations of the principal components
#' (i.e., the square roots of the eigenvalues of the
#' covariance/correlation matrix, though the calculation
#' is actually done with the singular values of the data matrix)
#' \item \code{eig} the cumulated proportion of variance along the PC axes
#' \item \code{rotation} the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors).
#' The function princomp returns this in the element loadings.
#' \item \code{center}, scale the centering and scaling used
#' \item \code{x} PCA scores (the value of the rotated data (the centred (and scaled if requested)
#' data multiplied by the rotation matrix))
#' \item other components are inherited from the \code{Coe} object passed to \code{PCA},
#' eg \code{fac}, \code{mshape}, \code{method}, \code{baseline1} and \code{baseline2}, etc. They
#' are documented in the corresponding \code{*Coe} file.
#' }
#' @family multivariate
#' @examples
#' bot.f <- efourier(bot, 12)
#' bot.p <- PCA(bot.f)
#' bot.p
#' plot(bot.p, morpho=FALSE)
#' plot(bot.p, 'type')
#'
#' op <- npoly(olea, 5)
#' op.p <- PCA(op)
#' op.p
#' plot(op.p, 1, morpho=TRUE)
#'
#' wp <- fgProcrustes(wings, tol=1e-4)
#' wpp <- PCA(wp)
#' wpp
#' plot(wpp, 1)
#'
#' # "foreign prcomp"
#' head(iris)
#' iris.p <- prcomp(iris[, 1:4])
#' iris.p <- as_PCA(iris.p, iris[, 5])
#' class(iris.p)
#' plot(iris.p, 1)
#' @export
PCA <- function(x, scale., center, fac) {
UseMethod("PCA")
}
#' @rdname PCA
#' @export
PCA.OutCoe <- function(x, scale. = FALSE, center = TRUE, fac) {
OutCoe <- x
PCA <- prcomp(OutCoe$coe, scale. = scale., center = center)
eig <- (PCA$sdev^2)
PCA$eig <- eig/sum(eig)
PCA$fac <- OutCoe$fac
PCA$mshape <- apply(OutCoe$coe, 2, mean)
PCA$method <- OutCoe$method
if (!is.null(OutCoe$baseline1)){
PCA$baseline1 <- OutCoe$baseline1
PCA$baseline2 <- OutCoe$baseline2}
PCA$cuts <- OutCoe$cuts
class(PCA) <- c("PCA", class(PCA))
return(PCA)
}
#' @rdname PCA
#' @export
PCA.OpnCoe <- function(x, scale. = FALSE, center = TRUE, fac) {
OpnCoe <- x
PCA <- prcomp(OpnCoe$coe, scale. = scale., center = center)
eig <- (PCA$sdev^2)
PCA$eig <- eig/sum(eig)
PCA$fac <- OpnCoe$fac
PCA$mshape <- apply(OpnCoe$coe, 2, mean)
PCA$method <- OpnCoe$method
PCA$mod <- OpnCoe$mod #the only diff so far
PCA$baseline1 <- OpnCoe$baseline1
PCA$baseline2 <- OpnCoe$baseline2
PCA$cuts <- OpnCoe$cuts
class(PCA) <- c("PCA", class(PCA))
return(PCA)
}
#' @rdname PCA
#' @export
PCA.LdkCoe <- function(x, scale. = FALSE, center = TRUE, fac) {
LdkCoe <- x
# LdkCoe$coe <- a2m(l2a(Coe$coo))
PCA <- prcomp(LdkCoe$coe, scale. = scale., center = center)
eig <- (PCA$sdev^2)
PCA$eig <- eig/sum(eig)
PCA$fac <- LdkCoe$fac
PCA$mshape <- apply(LdkCoe$coe, 2, mean)
PCA$method <- "procrustes"
PCA$cuts <- LdkCoe$cuts
PCA$links <- LdkCoe$links
# PCA$mod <- OpnCoe$mod #the only diff so far PCA$baseline1
# <- OpnCoe$baseline1 PCA$baseline2 <- OpnCoe$baseline2
class(PCA) <- c("PCA", class(PCA))
return(PCA)
}
#' @rdname PCA
#' @export
PCA.TraCoe <- function(x, scale. = TRUE, center = TRUE, fac) {
TraCoe <- x
# LdkCoe$coe <- a2m(l2a(Coe$coo))
PCA <- prcomp(TraCoe$coe, scale. = scale., center = center)
eig <- (PCA$sdev^2)
PCA$eig <- eig/sum(eig)
PCA$fac <- TraCoe$fac
PCA$mshape <- NULL
PCA$method <- NULL
class(PCA) <- c("PCA", class(PCA))
return(PCA)
}
#' @rdname PCA
#' @export
PCA.default <- function(x, scale. = TRUE, center = TRUE, fac=dplyr::tibble()) {
PCA <- prcomp(x, scale. = scale., center = center)
eig <- (PCA$sdev^2)
PCA$eig <- eig/sum(eig)
if (!is.null(fac)) fac <- tibble::as_tibble()(fac)
PCA$fac <- fac
PCA$method <- NULL
# PCA$baseline2 <- OpnCoe$baseline2
class(PCA) <- c("PCA", class(PCA))
return(PCA)
}
# 2. PCA Bridges ------------------------------------------
#' @rdname PCA
#' @export
as_PCA <- function(x, fac){
UseMethod("as_PCA")
}
#' @export
as_PCA.default <- function(x, fac){
if (class(x)[1] != "PCA"){
class(x) <- c("PCA", class(x))
if (!missing(fac)) x$fac <- tibble::as_tibble(fac)
return(x)}}
#' @export
print.PCA <- function(x, ...){
cat("A PCA object\n")
cat(rep("-", 20), "\n", sep = "")
cat(" -", ifelse(is.null(nrow(x$x)), 1, nrow(x$x)), "shapes \n")
# Method printer
if (length(x$method)>1) {
cat(" - $method: [", paste0(x$method, collapse=" + "), "analyses ]\n")
} else {
cat(" - $method: [", x$method, "analysis ]\n")}
# we print the fac
.print_fac(x$fac)
cat(" - All components: ", paste(names(x), collapse=", "), ".\n", sep="")
}
# 3. PCA methods ------------------------------------------
#' Get paired individual on a Coe, PCA or LDA objects
#'
#' If you have paired individuals, i.e. before and after a treatment or for repeated measures,
#' and if you have coded coded it into \code{$fac}, this methods allows you to retrieve the corresponding PC/LD scores,
#' or coefficients for \link{Coe} objects.
#' @param x any \link{Coe}, \link{PCA} of \link{LDA} object.
#' @param fac factor or column name or id corresponding to the pairing factor.
#' @param range numeric the range of coefficients for \code{Coe}, or PC (LD) axes on which to return scores.
#' @return a list with components \code{x1} all coefficients/scores corresponding to the
#' first level of the \code{fac} provided; \code{x2} same thing for the second level;
#' \code{fac} the corresponding \code{fac}.
#' @examples
#' bot2 <- bot1 <- coo_scale(coo_center(coo_sample(bot, 60)))
#' bot1$fac$session <- factor(rep("session1", 40))
#' # we simulate an measurement error
#' bot2 <- coo_jitter(bot1, amount=0.01)
#' bot2$fac$session <- factor(rep("session2", 40))
#' botc <- combine(bot1, bot2)
#' botcf <- efourier(botc, 12)
#'
#' # we gonna plot the PCA with the two measurement sessions and the two types
#' botcp <- PCA(botcf)
#' plot(botcp, "type", col=col_summer(2), pch=rep(c(1, 20), each=40), eigen=FALSE)
#' bot.pairs <- get_pairs(botcp, fac = "session", range=1:2)
# # with bot.pairs we can add segments between the two sessions
#' segments(bot.pairs$session1[, 1], bot.pairs$session1[, 2],
#' bot.pairs$session2[, 1], bot.pairs$session2[, 2],
#' col=col_summer(2)[bot.pairs$fac$type])
#'
#'
#' @export
get_pairs <- function(x, fac, range){UseMethod("get_pairs")}
#' @export
get_pairs.Coe <- function(x, fac, range){
# we check and prepare
fac <- fac_dispatcher(x, fac)
if (nlevels(fac) != 2) stop("more than two levels for the 'fac' provided")
tab <- table(fac)
if (length(unique(tab))!=1) stop("some mismatches between pairs")
# we get paired individuals
if (missing(range)) range <- 1:ncol(x$coe)
fl <- levels(fac)
x1 <- x$coe[fac==fl[1], range]
x2 <- x$coe[fac==fl[2], range]
res <- list(x1, x2, fac=x$fac[fac==fl[1],])
names(res)[1:2] <- fl
return(res)
}
#' @export
get_pairs.PCA <- function(x, fac, range){
# we check and prepare
fac <- fac_dispatcher(x, fac)
if (nlevels(fac) != 2) stop("more than two levels for the 'fac' provided")
tab <- table(fac)
if (length(unique(tab))!=1) stop("some mismatches between pairs")
# we get paired individuals
if (missing(range)) range <- 1:ncol(x$x)
fl <- levels(fac)
x1 <- x$x[fac==fl[1], range]
x2 <- x$x[fac==fl[2], range]
res <- list(x1, x2, fac=x$fac[fac==fl[1],])
names(res)[1:2] <- fl
return(res)
}
#' @export
get_pairs.LDA <- get_pairs.PCA
# 4. redo PCA ---------
#' "Redo" a PCA on a new Coe
#'
#' Basically reapply rotation to a new Coe object.
#' @param PCA a \link{PCA} object
#' @param Coe a \link{Coe} object
#' @note Quite experimental. Dimensions of the matrices and methods must match.
#' @examples
#' b <- filter(bot, type=="beer")
#' w <- filter(bot, type=="whisky")
#'
#' bf <- efourier(b, 8)
#' bp <- PCA(bf)
#'
#' wf <- efourier(w, 8)
#'
#' # and we use the "beer" PCA on the whisky coefficients
#' wp <- rePCA(bp, wf)
#'
#' plot(wp)
#'
#' plot(bp, eig=FALSE)
#' points(wp$x[, 1:2], col="red", pch=4)
#'
#'@export
rePCA <- function(PCA, Coe){
UseMethod("rePCA")
}
#'@export
rePCA.default <- function(PCA, Coe){
stop("method only defined for PCA objetcs")
}
#'@export
rePCA.PCA <- function(PCA, Coe){
if (Coe$method != PCA$method)
warning("methods differ between Coe and PCA")
scores <- PCA$x
rot <- PCA$rotation
coe <- Coe$coe
if (any(colnames(coe) != rownames(rot))) warning("matrices coefficients must match")
# we prepare a new PCA object
PCA2 <- PCA
PCA2$x <- matrix(NA, nrow(coe), ncol(rot), dimnames = list(rownames(coe), colnames(rot)))
PCA2$fac <- Coe$fac
# we recenter
coe <- apply(coe, 2, function(x) x - mean(x))
# learn matrix calculus bitch
for (PC in 1:ncol(rot)){
for (ind in 1:nrow(coe)){
PCA2$x[ind, PC] <- sum(coe[ind, ] * rot[, PC])
}
}
return(PCA2)
}
#' Calculates convex hull area/volume of PCA scores
#'
#' May be useful to compare shape diversity. Expressed in PCA units that should
#' only be compared within the same PCA.
#'
#' @param x a PCA object
#' @param fac (optionnal) column name or ID from the $fac slot.
#' @param xax the first PC axis to use (1 by default)
#' @param yax the second PC axis (2 by default)
#' @param zax the third PC axis (3 by default only for volume)
#'
#' @return
#' If fac is not provided global area/volume is returned; otherwise a named
#' list for every level of fac
#'
#' @details get_chull_area is calculated using \link{coo_chull} followed by \link{coo_area};
#' get_chull_volume is calculated using geometry::convexhulln
#'
#' @examples
#' bp <- PCA(efourier(bot, 12))
#' get_chull_area(bp)
#' get_chull_area(bp, 1)
#'
#' get_chull_volume(bp)
#' get_chull_volume(bp, 1)
#' @export
get_chull_area <- function (x, fac, xax = 1, yax = 2) {
if (!is_PCA(x)) stop("'x' must be a PCA object")
# no fac provided
if (missing(fac)){
xy <- x$x[, c(xax, yax)]
return(coo_area(coo_chull(xy)))
}
# else... if a fac is provided
fac <- fac_dispatcher(x, fac)
x <- x$x[, c(xax, yax)]
# we prepare the list
res <- list()
# we loop over to extract subset of coordinates
for (i in seq(along = levels(fac))) {
xy.i <- x[fac == levels(fac)[i], ]
# boring but prevents the numeric/2rows matrices
if (is.matrix(xy.i)) {
if (nrow(xy.i) > 2) {
res[[i]] <- coo_chull(xy.i)
} else {
res[[i]] <- NULL
}
} else {
res[[i]] <- NULL
}
}
# we calculate the area and return the results
names(res) <- levels(fac)
res <- lapply(res, coo_area)
return(res)
}
#' @rdname get_chull_area
#' @export
get_chull_volume <- function (x, fac, xax = 1, yax = 2, zax = 3) {
if (!is_PCA(x)) stop("'x' must be a PCA object")
# no fac provided
if (missing(fac)){
xy <- x$x[, c(xax, yax, zax)]
res <- geometry::convhulln(xy, options="FA")$vol
return(res)
}
# else...fac provided
fac <- fac_dispatcher(x, fac)
# we prepare the list
x <- x$x[, c(xax, yax, zax)]
res <- list()
# we loop over to extract subset of coordinates
for (i in seq(along = levels(fac))) {
xy.i <- x[fac == levels(fac)[i], ]
# boring but prevents the numeric/2rows matrices
if (is.matrix(xy.i)) {
if (nrow(xy.i) >= 4) {
res[[i]] <- geometry::convhulln(xy.i, options="FA")$vol
} else {
res[[i]] <- NULL
}
} else {
res[[i]] <- NULL
}
}
# we calculate the area and return the results
names(res) <- levels(fac)
return(res)
}
# 5. Flip PCA axes -------------
#' Flips PCA axes
#'
#' Simply multiply by -1, corresponding scores and rotation vectors for PCA objects.
#' PC orientation being arbitrary, this may help to have a better display.
#' @param x a PCA object
#' @param axs numeric which PC(s) to flip
#' @examples
#' bp <- bot %>% efourier(6) %>% PCA
#' bp %>% plot
#' bp %>% flip_PCaxes(1) %>% plot()
#' @export
flip_PCaxes <- function(x, axs){
UseMethod("flip_PCaxes")
}
#' @export
flip_PCaxes.default <- function(x, axs){
message("only defined on PCA objects")
}
#' @export
flip_PCaxes.PCA <- function(x, axs){
x$x[, axs] %<>% `*`(-1)
x$rotation[, axs] %<>% `*`(-1)
x
}
#### end PCA
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.