R/estimate_cbsem.R

In ISS-Analytics/seminr: Building and Estimating Structural Equation Models

#' seminr estimate_cbsem() function
#'
#' The \code{seminr} package provides a natural syntax for researchers to describe
#' structural equation models.
#'
#' @usage
#' estimate_cbsem(data, measurement_model = NULL,
#'                structural_model = NULL, item_associations = NULL,
#'                model = NULL, lavaan_model = NULL, estimator = "MLR", ...)
#'
#' @param data A \code{dataframe} containing the indicator measurement data.
#'
#' The entire CBSEM model can be specified in one of three ways:
#'
#' The pair of measurement and structural models, along associated items, can optionally be specified as separate model components
#'
#' @param measurement_model An optional \code{measurement_model} object representing the outer/measurement model,
#'   as generated by \code{constructs}.
#'   Note that only reflective constructs are supported for CBSEM models,
#'   though a composite measurement model can be converted into a reflective one
#'   using \code{\link{as.reflective}}.
#'
#' @param structural_model An optional \code{smMatrix} object representing the inner/structural model,
#'   as generated by \code{relationships}.
#'
#' @param item_associations An item-to-item matrix representing error
#'   covariances that are freed for estimation.
#'   This matrix is created by \code{associations()}, or defaults to NULL
#'   (no inter-item associations).
#'
#' The combination of measurement and structural models and inter-item associations can also be specified as a single \code{specified_model} object
#' Note that any given model components (measurement_model, structural_model, item_associations) will override components in the fully specified model
#'
#' @param model An optional \code{specified_model} object containing both the the outer/measurement and inner/structural models,
#'   along with any inter-item associations, as generated by \code{specify_model}.
#'
#' The entire model can also be specified in Lavaan syntax (this overrides any other specifications)
#'
#' @param lavaan_model Optionally, a single character string containing the relevant model specification in \code{lavaan} syntax.
#'
#' Any further optional parameters to alter the estimation method:
#'
#' @param estimator A character string indicating which estimation method to use
#'   in Lavaan. It defaults to "MLR" for robust estimation.
#'   See the Lavaan documentation for other supported estimators.
#'
#' @param ... Any other parameters to pass to \code{lavaan::sem} during
#'   estimation.
#'
#' @return A list of the estimated parameters for the CB-SEM model including:
#'  \item{data}{A matrix of the data upon which the model was estimated.}
#'  \item{measurement_model}{The SEMinR measurement model specification.}
#'  \item{factor_loadings}{The matrix of estimated factor loadings.}
#'  \item{associations}{A matrix of model variable associations.}
#'  \item{mmMatrix}{A Matrix of the measurement model relations.}
#'  \item{smMatrix}{A Matrix of the structural model relations.}
#'  \item{constructs}{A vector of the construct names.}
#'  \item{construct scores}{A matrix of the estimated construct scores for the CB-SEM model.}
#'  \item{item_weights}{A matrix of the estimated CFA item weights.}
#'  \item{lavaan_model}{The lavaan model syntax equivalent of the SEMinR model.}
#'  \item{lavaan_output}{The raw lavaan output generated after model estimation.}
#'
#' @references Joreskog, K. G. (1973). A general method for estimating a linear structural equation system In: Goldberger AS, Duncan OD, editors. Structural Equation Models in the Social Sciences. New York: Seminar Press.
#'
#' @seealso \code{\link{as.reflective}}
#'          \code{\link{relationships}} \code{\link{constructs}}
#'          \code{\link{paths}}
#'          \code{\link{associations}} \code{\link{item_errors}}
#'
#' @examples
#' mobi <- mobi
#'
#' #seminr syntax for creating measurement model
#' mobi_mm <- constructs(
#'   reflective("Image",        multi_items("IMAG", 1:5)),
#'   reflective("Quality",      multi_items("PERQ", 1:7)),
#'   reflective("Value",        multi_items("PERV", 1:2)),
#'   reflective("Satisfaction", multi_items("CUSA", 1:3)),
#'   reflective("Complaints",   single_item("CUSCO")),
#'   reflective("Loyalty",      multi_items("CUSL", 1:3))
#' )
#'
#' #seminr syntax for freeing up item-item covariances
#' mobi_am <- associations(
#'   item_errors(c("PERQ1", "PERQ2"), "IMAG1")
#' )
#'
#' #seminr syntax for creating structural model
#' mobi_sm <- relationships(
#'   paths(from = c("Image", "Quality"), to = c("Value", "Satisfaction")),
#'   paths(from = c("Value", "Satisfaction"), to = c("Complaints", "Loyalty")),
#'   paths(from = "Complaints",   to = "Loyalty")
#' )
#'
#' # Estimate model and get results
#' mobi_cbsem <- estimate_cbsem(mobi, mobi_mm, mobi_sm, mobi_am)
#'
#' # Use or capture the summary object for more results and metrics
#' summary(mobi_cbsem)
#'
#' cbsem_summary <- summary(mobi_cbsem)
#' cbsem_summary$descriptives$correlations$constructs
#'
#' @export
estimate_cbsem <- function(data, measurement_model=NULL, structural_model=NULL, item_associations=NULL, model=NULL, lavaan_model=NULL, estimator="MLR", ...) {
  message("Generating the seminr model for CBSEM")

  # TODO: consider higher order models (see estimate_pls() function for template)
  rawdata <- data

  if (is.null(lavaan_model)) {
    # Extract model specifications
    specified_model <- extract_models(model, measurement_model, structural_model, item_associations)
    measurement_model <- specified_model$measurement_model
    structural_model  <- specified_model$structural_model
    item_associations <- specified_model$item_associations

    # Process measurement model interactions to produce simplified mmMatrix
    post_interaction_object <- process_cbsem_interactions(measurement_model, data, structural_model, item_associations, estimator, ...)
    names(post_interaction_object$data) <- sapply(names(post_interaction_object$data), FUN=lavaanify_name, USE.NAMES = FALSE)
    mmMatrix <- post_interaction_object$mmMatrix
    data <- post_interaction_object$data

    # Rename interaction terms
    structural_model[, "source"] <- sapply(structural_model[, "source"], FUN=lavaanify_name)
    smMatrix <- structural_model

    # TODO: warning if the model is incorrectly specified
    # warnings(measurement_model, data, structural_model)

    # Create LAVAAN syntax
    measurement_syntax <- lavaan_mm_syntax(mmMatrix)
    structural_syntax <- lavaan_sm_syntax(smMatrix)
    association_syntax <- lavaan_item_associations(item_associations)

    # Put all the parts together
    lavaan_model <- paste(measurement_syntax, structural_syntax, association_syntax, sep="\n\n")
  } else {
    structural_model <- smMatrix <- lavaan2seminr(lavaan_model)$structural_model
    measurement_model <- lavaan2seminr(lavaan_model)$measurement_model

    # using process_cbsem_interactions() to produce mmMatrix
    post_interaction_object <- process_cbsem_interactions(measurement_model, data, structural_model, item_associations, estimator, ...)
    mmMatrix <- post_interaction_object$mmMatrix
  }

  # Estimate cbsem in Lavaan
  lavaan_output <- try_or_stop(
    lavaan::sem(
      model=lavaan_model, data=data, std.lv = TRUE, estimator=estimator, ...),
    "estimating CBSEM using Lavaan"
  )

  # Extract lavaan results
  constructs <- all_construct_names(measurement_model) # needed in object for reliability... move up if lavaan_model no longer supported
  lavaan_std <- lavaan::lavInspect(lavaan_output, what="std")
  HOFs <- HOCs_in_model(measurement_model, structural_model)
  if (length(HOFs) > 0) {
    loadings <- combine_first_order_second_order_loadings_cbsem(mmMatrix, rawdata, lavaan_std)
  } else {
    loadings <- lavaan_std$lambda
    class(loadings) <- "matrix"
  }

  # Arrange Coefficients Table
  estimates <- lavaan::standardizedSolution(lavaan_output)
  path_df <- estimates[estimates$op == "~",]
  all_antecedents <- all_exogenous(smMatrix)
  all_outcomes <- all_endogenous(smMatrix)
  path_matrix <- df_xtab_matrix(est.std ~ rhs + lhs, path_df,
                                all_antecedents, all_outcomes)
  rownames(path_matrix) <- gsub("_x_", "*", all_antecedents)

  # Compute results from our own methods
  tenB <- estimate_lavaan_ten_berge(lavaan_output)

  # Gather model information
  seminr_model <- list(
    data = data,
    rawdata = rawdata,
    measurement_model = measurement_model,
    factor_loadings = loadings,
    associations = item_associations,
    mmMatrix = mmMatrix,
    smMatrix = smMatrix,
    constructs = constructs,
    construct_scores = tenB$scores,
    item_weights = tenB$weights,
    path_coef = path_matrix,
    lavaan_model = lavaan_model,
    lavaan_output = lavaan_output
  )

  class(seminr_model) <- c("cbsem_model", "seminr_model")
  return(seminr_model)
}

#' seminr estimate_cfa() function
#'
#' Estimates a Confirmatory Factor Analysis (CFA) model
#'
#' @inheritParams estimate_cbsem
#'
#' @return A list of the estimated parameters for the CFA model including:
#'  \item{data}{A matrix of the data upon which the model was estimated.}
#'  \item{measurement_model}{The SEMinR measurement model specification.}
#'  \item{construct scores}{A matrix of the estimated construct scores for the CB-SEM model.}
#'  \item{item_weights}{A matrix of the estimated CFA item weights.}
#'  \item{lavaan_model}{The lavaan model syntax equivalent of the SEMinR model.}
#'  \item{lavaan_output}{The raw lavaan output generated after model estimation.}
#'
#' @usage
#' estimate_cfa(data, measurement_model = NULL, item_associations=NULL,
#'              model = NULL, lavaan_model = NULL, estimator="MLR", ...)
#'
#' @references Jöreskog, K.G. (1969) A general approach to confirmatory maximum likelihood factor analysis. Psychometrika, 34, 183-202.
#'
#' @seealso \code{\link{constructs}} \code{\link{reflective}}
#'          \code{\link{associations}} \code{\link{item_errors}}
#'          \code{\link{as.reflective}}
#'
#' #' @examples
#' mobi <- mobi
#'
#' #seminr syntax for creating measurement model
#' mobi_mm <- constructs(
#'   reflective("Image",        multi_items("IMAG", 1:5)),
#'   reflective("Expectation",  multi_items("CUEX", 1:3)),
#'   reflective("Quality",      multi_items("PERQ", 1:7))
#' )
#'
#' #seminr syntax for freeing up item-item covariances
#' mobi_am <- associations(
#'   item_errors(c("PERQ1", "PERQ2"), "CUEX3"),
#'   item_errors("IMAG1", "CUEX2")
#' )
#'
#' mobi_cfa <- estimate_cfa(mobi, mobi_mm, mobi_am)
#'
#' @export
estimate_cfa <- function(data, measurement_model=NULL, item_associations=NULL,
                         model=NULL, lavaan_model=NULL, estimator="MLR", ...) {
  message("Generating the seminr model for CFA")

  # TODO: consider higher order models (see estimate_pls() function for template)
  # TODO: warning if the model is incorrectly specified
  # warnings(measurement_model, data, structural_model)

  mmMatrix <- NULL
  rawdata <- data
  if (is.null(lavaan_model)) {
    # Extract specified models
    specified_model <- extract_models(
      model = model, measurement_model = measurement_model, item_associations = item_associations
    )
    measurement_model <- specified_model$measurement_model
    item_associations <- specified_model$item_associations
    constructs <- all_construct_names(measurement_model)

    # Create LAVAAN syntax
    mmMatrix <- mm2matrix(measurement_model)

    measurement_syntax <- lavaan_mm_syntax(mmMatrix)
    association_syntax <- lavaan_item_associations(item_associations)

    lavaan_model <- paste(measurement_syntax,
                          association_syntax,
                          sep="\n\n")
  }

  # Estimate cfa in Lavaan
  lavaan_output <- try_or_stop(
    lavaan::cfa(model=lavaan_model, data=data, std.lv = TRUE,
                estimator=estimator, ...),
    "run CFA in Lavaan"
  )

  # Extract Lavaan results
  lavaan_std <- lavaan::lavInspect(lavaan_output, what="std")
  HOFs <- HOCs_in_model(measurement_model)
  if (length(HOFs) > 0) {
    loadings <- combine_first_order_second_order_loadings_cbsem(mmMatrix, rawdata, lavaan_std)
  } else {
    loadings <- lavaan_std$lambda
    class(loadings) <- "matrix"
  }

  # Compute results from our own methods
  tenB <- estimate_lavaan_ten_berge(lavaan_output)

  # Gather model information
  seminr_model <- list(
    data = data,
    measurement_model = measurement_model,
    factor_loadings = loadings,
    constructs = constructs,
    construct_scores = tenB$scores,
    item_weights = tenB$weights,
    lavaan_model = lavaan_model,
    lavaan_output = lavaan_output
  )

  class(seminr_model) <- c("cfa_model", "seminr_model")
  return(seminr_model)
}