Nothing
# 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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.