R/factorAnalysis.R

In Timothy-Barry/coproanalysis: Analyze CoPRO and Related Data

#' Compute a factor analysis using the SVD method
#'
#' @param nFactors number of factors to use
#' @param corMat correlation matrix
#' @param rotation which rotation to use. Can pass a rotation matrix, or alternately a string representing a rotation (e.g., "varimax").
#' @export

computeFactorAnalysis <- function(nFactors, corMat, rotation) {
decomp <- eigen(x = corMat, symmetric = TRUE)
vectors <- decomp$vectors[,1:nFactors]
sqrtValues <- sqrt(decomp$values[1:nFactors])
lTilde <- vectors %*% diag(sqrtValues)
lLTTilde <- lTilde %*% t(lTilde)
specificVariances <- diag(corMat - lLTTilde)
communalities <- apply(lTilde, 1, function(r){sum(r^2)})
leadingEvalues <- decomp$values[1:nFactors]
if (class(rotation) == "matrix") {
  loadings <- lTilde %*% rotation
  rotMat <- rotation
} else if (rotation == "varimax") {
  rotatedLoadings <- Varimax(lTilde)
  loadings <- rotatedLoadings$loadings
  rotMat <- rotatedLoadings$Th
}
residMat <- corMat - (lLTTilde + diag(specificVariances))
list(loadings = loadings, rotMat = rotMat, communalities = communalities, specificVariances = specificVariances, leadingEvalues = leadingEvalues, leadingVectors = vectors, residMat = residMat)
}

#' getFactorScores
#'
#' @param loadings the factor loadings
#' @param dat the data
#' @export

getFactorScores <- function(loadings, dat) {
  projMat <- solve(t(loadings) %*% loadings) %*% t(loadings)
  t(projMat %*% t(as.matrix(dat))) %>% as_tibble
}

#' computeFactorAnalysisWithBootstrappedCIs
#'
#' @param dat data frame
#' @param nFactors number of factors
#' @param rotation rotation to use (either a literal matrix or string specifying analytic rotation, such as "varimax")
#' @param alpha (1-alpha) confidence interval returned
#' @param B number of bootstrap replicates
#' @param bonferoniCorrect should bonferoni correction be applied in computing CIs?
#'
#' @return A list with two entries. The first is a list containing the bootstrap samples of the factor analysis; the second is a list containing the results of the factor analysis on the entire dataset.
#' @export

getBootstrapSamples <- function(dat, nFactors, rotation, B = 5000) {
  # First, compute a "base" factor analysis
  baseFA <- computeFactorAnalysis(nFactors = nFactors, corMat = cor(dat), rotation = rotation)
  myRotation <- baseFA$rotMat
  # Which factor loading is greatest? We use that as a reference to ensure that loading is positive.
  refs <- baseFA$loadings
  # Next, compute B bootstrapped statistics
  f <- function() {
    datB <- dat[sample(1:nrow(dat), nrow(dat), TRUE),]
    corMatB <- cor(datB)
    res <- computeFactorAnalysis(nFactors, corMatB, myRotation)
    res$loadings <- pointEvectorsInSameDirection(res$loadings, refs)
    res$corMat <- corMatB
    res
  }
  bootSamples <- replicate(B, f(), FALSE)
  return(list(bootSamples = bootSamples, baseFA = baseFA))
}

#' pointEvectorsInSameDirection
#'
#' @param factorLoadings factor loadings
#' @param refs the rows to check (for each column) for sign
pointEvectorsInSameDirection <- function(factorLoadings, refs) {
  correctives <- sapply(1:ncol(factorLoadings), function(i) {ifelse(sum(refs[,i] * factorLoadings[,i]) < 0, -1, 1)})
  factorLoadings %*% diag(correctives)
}