R/selectPCs.R

In iholzleitner/frlgmm3D: FRL convenience functions around GMM

#'selectPCs
#'
#'Select 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, and assigns number of
#'  selected PCs to variable called "n.PCs.ARGUMENTNAME" to global environment
#'
#' @examples
#'
#'my_data<- tibble::tibble(
#'x=c(rnorm(100,1,1), rnorm(100, 1, 0.5), rnorm(100, 4, 1), rnorm(100, 1, 1.5)),
#'y=c(rnorm(100,2,1), rnorm(100, 3, 0.5), rnorm(100, 2, 1), rnorm(100, 2, 1.5)),
#'z=c(rnorm(100,1,1.5), rnorm(100, 5, 1), rnorm(100, 1, .5), rnorm(100, 2.5, 2)),
#')
#'
#'PCA.out<-prcomp(my_data)
#'
#'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
  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])
  }
}