R/selectPCs.R

In iholzleitner/demopckg: Convenience functions around GMM

#' Select principal components using Broken Stick criterion
#'
#' Selects principal components according to Broken Stick Model (MacArthur 1957). Adapted from function \href{http://adn.biol.umontreal.ca/~numericalecology/numecolR/}{"evplot" by Francois Gillet}.
#'
#' @param PCA_output Output from geomorph::plotTangentSpace or stats::prcomp
#' @return Returns selected PCs and variance explained
#'
#' @examples
#' data(FRL_London)
#' PCA_out <- prcomp(FRL_London)
#'
#' selectPCs(PCA_out)
#'
#' @export

selectPCs <- function(PCA_output) {
  ev <- PCA_output$sdev^2
  n <- length(ev)
  bsm <- data.frame(j = seq(1:n), p = 0)
  bsm$p[1] <- 1/n
  for (i in 2:n) bsm$p[i] <- bsm$p[i - 1] + (1 / (n + 1 - i))
  bsm$p <- 100 * bsm$p / n

  test <- cbind(100 * ev / sum(ev), bsm$p[n:1])
  n_PCs <- sum(test[, 1] >= test[, 2])

  # Save number of principal components // Violates CRAN norms. Think of alternative.
  # arg_name <- deparse(substitute(PCA_output))
  # var_name <- paste("n_PCs", arg_name, sep=".")
  # assign(var_name, n_PCs, .GlobalEnv)

  # Print PCs and variance explained
  if (!is.null(PCA_output$pc.summary$importance)) {
    return(PCA_output$pc.summary$importance[, 1:n_PCs])
  } else {
    temp <- summary(PCA_output)
    return(temp$importance[, 1:n_PCs])
  }
}