R/corEFA.R

In lessR: Less Code with More Comprehensive Results

corEFA <- 
function (R=mycor, n_factors, rotate=c("promax", "varimax", "none"), 
          min_loading=.2, sort=TRUE, ...) {


  # a dot in a parameter name to an underscore
  dots <- list(...)
  change <- c("n.factors", "min.loading")
  if (length(dots) > 0) {
    for (i in seq_along(dots)) {
      if (names(dots)[i] %in% change) {
        pnm <- gsub(".", "_", names(dots)[i], fixed=TRUE)
        assign(pnm, dots[[i]])
      }
    }
  }
  dots_clean <- dots[!names(dots) %in% change]

  cl <- match.call()

  rotate <- match.arg(rotate)

  if (missing(n_factors)) {
    cat("\n"); stop(call.=FALSE, "\n","------\n",
      "Specify the number of factors with:  n_factors=\n\n")
  }

  if (!("matrix" %in% class(R))) { # R is a matrix, can be called indirectly
    # cor matrix:  mycor as class out_all, mycor$R, or stand-alone matrix
    cor.nm <- deparse(substitute(R))
    .cor.exists(cor.nm)  # see if matrix exists in one of the 3 locations
    if ("out_all" %in% class(R))  # R 4.0 results in two values: matrix, array
      R <- eval(parse(text=paste(cor.nm, "$R", sep="")))  # go to $R 
  }

  title_efa <- "  EXPLORATORY FACTOR ANALYSIS"

  txer <- ""
  if (is.null(options()$knitr.in.progress)) {
    tx <- character(length = 0)
    tx[length(tx)+1] <- "Extraction: maximum likelihood"
    if (n_factors > 1)  tx[length(tx)+1] <- paste("Rotation:", rotate)
    txer <- tx
  }

  # EFA
  fa2 <- do.call(factanal,
    c(list(covmat=R, factors=n_factors, rotation="none"), dots_clean))
  if (rotate=="none" || n_factors==1) ld <- as.matrix(fa2$loadings)

  if (n_factors>1  &&  rotate!="none") {
    if (rotate == "promax") rtt <- promax(loadings(fa2))
    if (rotate == "varimax") rtt <- varimax(loadings(fa2))
    ld <- loadings(rtt)
  }

  n.ind <- nrow(ld)

  # sort option
  # mx: # factor with highest factor loading for each item
  if (sort) {
    mx <- max.col(abs(ld))  
    ld <- ld[order(mx), ]  # group items on same factor together

    fn <- max.col(abs(ld))  # each item, factor with its highest loading
    ld <- data.frame(cbind(fn, ld))

    ld.new <- matrix(nrow=0, ncol=1+n_factors)
    for (i in seq_len(n_factors)) {
      lf <- ld[ld$fn==i,]  # extract items for ith factor
      ord <- order(abs(lf[[i+1]]), decreasing=TRUE)
      lf <- lf[ord,]
      ld.new <- rbind(ld.new,lf)  # preserves row names 
    }
    ld <- ld.new[,-1]  # remove fn column
  }

  # print loadings
  # from_efa lets .prntbl adjust column widths differently 
  tx <- character(length = 0)
  tx[length(tx)+1] <- paste("Loadings (except -", min_loading, " to ",
     min_loading, ")", sep="") 
  txld <- .prntbl(ld, digits_d=3, cut=min_loading, from_efa=TRUE)
  for (i in seq_along(txld)) tx[length(tx)+1] <- txld[i]
  txld <- tx

  # print sum of squares by factor
  vx <- colSums(ld^2)
  varex <- rbind(`SS loadings` = vx)
  varex <- rbind(varex, `Proportion Var` = vx/n.ind)
  if (n_factors > 1) 
    varex <- rbind(varex, `Cumulative Var` = cumsum(vx/n.ind))
  tx <- character(length = 0)
  tx[length(tx)+1] <- "Sum of Squares"
  txss <- .prntbl(varex, 3)
  for (i in seq_along(txss)) tx[length(tx)+1] <- txss[i]
  txss <- tx


  title_cfa <- "  CONFIRMATORY FACTOR ANALYSIS CODE"

  # generate the MIMM code
  FacItems <- integer(length=n.ind)  # factor with highest loading for item
  Fac <- integer(length=n.ind)
  n.Fact <- integer(length=n_factors)

  for (i in seq_len(n.ind)) {
    max.ld <- 0
    FacItems[i] <- 0
    for (j in seq_len(n_factors)) {
      if (abs(ld[i,j]) > max.ld  &&  abs(ld[i,j]) > min_loading) {
        max.ld <- abs(ld[i,j])
        FacItems[i] <- j
      }
    }
  }


  tx <- character(length = 0)

  tx[length(tx)+1] <- "MeasModel <- "
  for (i.fact in seq_len(n_factors)) {
    n.Fact[i.fact] <- 0
    k <- 0
    for (j.item in seq_len(n.ind)) if (FacItems[j.item] == i.fact) {
      k <- k + 1
      Fac[k] <- j.item
      n.Fact[i.fact] <- n.Fact[i.fact] + 1
    }
    if (i.fact == 1)
      tx[length(tx)+1] <- "\""
    else
      tx[length(tx)+1] <- " "
    tx[length(tx)] <- paste(
      tx[length(tx)], "  F", as.character(i.fact), " =~ ", sep="")
    if (n.Fact[i.fact] > 0) {
      for (i in seq_len(n.Fact[i.fact])) {
        tx[length(tx)] <- paste(tx[length(tx)], rownames(ld)[Fac[i]], sep="")
        if (i < n.Fact[i.fact])
          tx[length(tx)] <- paste(tx[length(tx)], " + ", sep="")
      }
    }
    if (i.fact == n_factors) tx[length(tx)] <- paste(tx[length(tx)], "\n\"")
  }

  tx[length(tx)+1] <- ""
  tx[length(tx)+1] <- "fit <- lessR::cfa(MeasModel)\n"

  tx[length(tx)+1] <- "library(lavaan)"
  tx[length(tx)+1] <- "fit <- lavaan::cfa(MeasModel, data=d)"
  tx[length(tx)+1] <- "summary(fit, fit.measures=TRUE, standardized=TRUE)"
  txcfa <- tx


  # report any deleted items
  deleted <- integer(length=n.ind)
  del.count <- 0
  for (i.item in seq_len(n.ind)) if (FacItems[i.item] == 0) {
    del.count <- del.count + 1
    deleted[del.count] <- i.item
  }

  txdel <- ""
  tx <- character(length = 0)
  if (del.count > 0) {
    tx[length(tx)+1] <- paste("Deletion threshold: min_loading = ", 
        min_loading, sep="")
    tx[length(tx)+1] <- "Deleted items: "
    for (i.item in seq_len(del.count))
      tx[length(tx)] <- paste(tx[length(tx)], rownames(ld)[deleted[i.item]],
                              " ", sep="")
    txdel <- tx
  }


  class(title_efa) <- "out"
  class(txer) <- "out"
  class(txld) <- "out"
  class(txss) <- "out"
  class(title_cfa) <- "out"
  class(txcfa) <- "out"
  class(txdel) <- "out"

  output <- list(
    out_title_efa=title_efa, out_type=txer, out_loadings=txld, out_ss=txss,
    out_title_cfa=title_cfa, out_cfa_code=txcfa,
    out_deleted=txdel,

    converged=fa2$converged, n_factors=n_factors, ss_factors=vx,
    loadings=ld, call=cl)

  class(output) <- "out_all"
  return(output)

}