R/pls_fit.R

In plssem: Complex Partial Least Squares Structural Equation Modeling

getFitPLSModel <- function(model, consistent = TRUE, quick = model@status$quick) {
  if (!consistent && model@info$path.estimator == "ols")
    return(getFitPLSModelUncorrected(model, quick = quick))

  lambda    <- model@matrices$lambda
  gamma     <- model@matrices$gamma
  preds     <- model@matrices$preds
  etas      <- model@info$etas
  xis       <- model@info$xis
  lvs       <- model@info$lvs
  lvs.lin   <- model@info$lvs.linear
  inds      <- model@info$allInds
  inds.a    <- model@info$inds.a
  inds.b    <- model@info$inds.b
  indsLvs   <- model@info$indsLvs
  modes     <- model@info$modes
  mode.a    <- model@info$mode.a
  mode.b    <- model@info$mode.b
  ptl       <- model@parTableInput
  SC        <- model@matrices$SC
  estimator <- model@info$path.estimator

  fitMeasurement <- fitLambda <- fitWeights <- lambda
  fitMeasurement[TRUE] <- fitLambda[TRUE] <- fitWeights[TRUE] <- 0

  C <- model@matrices$C

  for (lv in lvs.lin) {
    inds.lv <- indsLvs[[lv]]
    mode.lv <- modes[[lv]]

    wq <- lambda[inds.lv, lv]
    lq <- SC[inds.lv, lv]
    pq <- switch(mode.lv, A = lq, B = wq, NA_real_)

    fitMeasurement[inds.lv, lv] <- pq
    fitWeights[inds.lv, lv]     <- wq
    fitLambda[inds.lv, lv]      <- lq
  }

  if (consistent) {
    Q                   <- getConstructQualities(model)
    fitMeasurement      <- getConsistentLoadings(model, Q = Q)
    fitLambda[, mode.a] <- fitMeasurement[, mode.a]
    C                   <- getConsistentCorrMat(model, Q = Q)

  } else {
    Q <- numeric(0)
    attr(Q, "admissible") <- TRUE

  }

  fitStructural       <- gamma
  fitStructural[TRUE] <- 0

  switch(estimator,
    ols = {

      # paths
      for (lv in lvs) {
        predsLv <- lvs[preds[, lv, drop = TRUE]]

        if (length(predsLv))
          fitStructural[predsLv, lv] <- getOlsPathCoefs(lv, predsLv, C)
      }

      # (residual) covariances
      fitCov     <- C
      fitCovProj <- t(fitStructural) %*% C %*% fitStructural
      fitCovRes  <- diag2(fitCov) - diag2(fitCovProj)
      fitCov[etas, etas]  <- fitCovRes[etas, etas]
      fitCov[etas, xis]   <- fitCov[xis, etas] <- 0
    },

    gls = {
      gmod <- model@glsPathModel

      success <- TRUE
      tryCatch({
        glsModelCovMatrix(gmod) <- C # update input
        gfit <- glsEstimateParameters(gmod) # fit model

      }, error = function(e) {
        success <<- FALSE
        pls_msg_warn(
          "Estimation of the structural model using GLS failed!",
          "Attempting to use OLS instead!",
          "Message:", conditionMessage(e)
        )
      })

      if (!success) {
        # switch to ols and mark as inadmissible
        model@info$path.estimator  <- "ols"
        model@status$is.admissible <- FALSE

        return( # this is not computationally efficient, but it's simple
          getFitPLSModel(model = model, consistent = consistent)
        )
      }

      # paths
      gamma <- gfit@matrices$gamma

      for (lv in lvs) {
        predsLv <- lvs[preds[, lv, drop = TRUE]]

        if (length(predsLv))
          fitStructural[predsLv, lv] <- gamma[lv, predsLv]
      }

      fitCov <- gfit@matrices$psi[rownames(C), colnames(C)]
    },

    # Shouldn't happen
    pls_msg_stop(
      "Unrecognized path estimator! Estimator:", estimator
    )
  )

  k        <- length(inds)
  fitTheta <- matrix(0, nrow = k, ncol = k, dimnames = list(inds, inds))

  if (!quick) {
    crossLoaded <- apply(
      X      = fitMeasurement,
      MARGIN = 1L,
      FUN    = \(x) sum(abs(x) > .Machine$double.xmin) > 1L
    )

    pls_warnif(any(crossLoaded),
      "Did not expect any cross loaded indicators,\n",
      "when calculating indicator residuals!"
    )
  }

  fitThetaFull <- model@matrices$SC[inds, inds]

  # keep formative blocks
  for (b in mode.b) {
    idx <- indsLvs[[b]]
    fitTheta[idx, idx] <- fitThetaFull[idx, idx]
  }

  for (ind in inds.a) {                                  # Guard for NaN in fitMeasurement
    j   <- max(which.max(abs(fitMeasurement[ind, ])), 1) # max(numeric(0), 1) = 1
    r   <- fitMeasurement[ind, j]
    v   <- SC[ind, ind]

    fitTheta[ind, ind] <- v - r^2
  }

  if (quick) {
    return(list(
      fitMeasurement    = fitMeasurement,
      fitStructural     = fitStructural,
      fitCov            = fitCov,
      fitTheta          = fitTheta,
      fitWeights        = fitWeights,
      fitLambda         = fitLambda,
      fitC              = C,
      Q                 = Q,
      status.admissible = model@status$is.admissible
    ))
  }

  list(
    fitMeasurement    = plssemMatrix(fitMeasurement, symmetric = FALSE),
    fitStructural     = plssemMatrix(fitStructural,  symmetric = FALSE),
    fitCov            = plssemMatrix(fitCov,         symmetric = TRUE),
    fitTheta          = plssemMatrix(fitTheta,       symmetric = TRUE),
    fitWeights        = plssemMatrix(fitWeights,     symmetric = FALSE),
    fitLambda         = plssemMatrix(fitLambda,      symmetric = FALSE),
    fitC              = plssemMatrix(C,              symmetric = FALSE),
    Q                 = plssemVector(Q),
    status.admissible = model@status$is.admissible
  )
}


getFitPLSModelUncorrected <- function(model, quick = model@status$quick) {
  matrices <- model@matrices
  info     <- model@info

  lambda  <- matrices$lambda
  C       <- matrices$C
  S       <- matrices$S
  inds    <- info$allInds
  lvs     <- info$lvs
  etas    <- info$etas
  xis     <- info$xis
  mode.b  <- info$mode.b
  inds.a  <- info$inds.a
  indsLvs <- info$indsLvs

  # Measurement model -------------------------------------------------------
  fitWeights <- lambda

  fitLambda   <- lambda
  fitLambda[] <- 0

  selected <- matrices$select$lambda
  loadings <- matrices$SC[inds, lvs, drop = FALSE]
  fitLambda[selected] <- loadings[selected]

  fitMeasurement <- fitLambda
  fitMeasurement[, mode.b] <- fitWeights[, mode.b, drop = FALSE]

  # Structural model --------------------------------------------------------
  fitStructural   <- matrices$gamma
  fitStructural[] <- 0
  fitCov          <- C

  if (length(etas)) {
    fitCov[etas, etas] <- 0
    fitCov[etas, xis]  <- 0
    fitCov[xis, etas]  <- 0
  }

  for (lv in etas) {
    pred.idx <- which(matrices$preds[, lv, drop = TRUE])
    if (!length(pred.idx)) next

    beta <- solve(
      C[pred.idx, pred.idx, drop = FALSE],
      C[pred.idx, lv, drop = FALSE]
    )

    fitStructural[pred.idx, lv] <- beta
    fitCov[lv, lv] <- C[lv, lv] - sum(beta * C[pred.idx, lv])
  }

  # Indicator residuals -----------------------------------------------------
  fitTheta <- matrix(
    0,
    nrow = length(inds), ncol = length(inds),
    dimnames = list(inds, inds)
  )

  for (b in mode.b) {
    idx <- indsLvs[[b]]
    fitTheta[idx, idx] <- S[idx, idx, drop = FALSE]
  }

  if (!quick) {
    crossLoaded <- rowSums(
      abs(fitMeasurement) > .Machine$double.xmin
    ) > 1L

    pls_warnif(any(crossLoaded),
               "Did not expect any cross loaded indicators,\n",
               "when calculating indicator residuals!")
  }

  if (length(inds.a)) {
    measurement.a <- fitMeasurement[inds.a, , drop = FALSE]
    construct.idx <- max.col(abs(measurement.a), ties.method = "first")
    loading <- measurement.a[cbind(seq_along(inds.a), construct.idx)]
    residual <- diag(S[inds.a, inds.a, drop = FALSE]) - loading^2
    theta.idx <- match(inds.a, inds)

    fitTheta[cbind(theta.idx, theta.idx)] <- residual
  }

  Q <- numeric(0)
  attr(Q, "admissible") <- TRUE

  if (quick) {
    return(list(
      fitMeasurement    = fitMeasurement,
      fitStructural     = fitStructural,
      fitCov            = fitCov,
      fitTheta          = fitTheta,
      fitWeights        = fitWeights,
      fitLambda         = fitLambda,
      fitC              = C,
      Q                 = Q,
      status.admissible = model@status$is.admissible
    ))
  }

  list(
    fitMeasurement    = plssemMatrix(fitMeasurement, symmetric = FALSE),
    fitStructural     = plssemMatrix(fitStructural,  symmetric = FALSE),
    fitCov            = plssemMatrix(fitCov,         symmetric = TRUE),
    fitTheta          = plssemMatrix(fitTheta,       symmetric = TRUE),
    fitWeights        = plssemMatrix(fitWeights,     symmetric = FALSE),
    fitLambda         = plssemMatrix(fitLambda,      symmetric = FALSE),
    fitC              = plssemMatrix(C,              symmetric = FALSE),
    Q                 = plssemVector(Q),
    status.admissible = model@status$is.admissible
  )
}


modelFitIsAdmissible <- function(fit, tol = 1e-12) {
  # Simple check to see if model fit is (in)admissible
  atol <- abs(tol)
  ltol <- 1 + atol # loadings

  Q.admissible <- (
    is.null(attr(fit$Q, "admissible")) ||
    isTRUE(attr(fit$Q, "admissible"))
  )

  (
    !anyNA(fit$fitWeights)                               &&
    !anyNA(fit$fitLambda)                                &&
    !anyNA(fit$fitStructural)                            &&
    !anyNA(fit$fitTheta)                                 &&
    !anyNA(fit$fitCov)                                   &&
    !anyNA(fit$fitC)                                     &&
    isPositiveDefinite(fit$fitC)                         &&
    all(diag(fit$fitTheta) >= -atol)                     &&
    all(diag(fit$fitCov) >= -atol)                       &&
    all(fit$fitLambda  >= -ltol & fit$fitLambda <= ltol) && # weights can exceed +/- 1, but not loadings
    Q.admissible                                         &&
    fit$status.admissible # check flag from the input model
  )
}


getParamVecNames <- function(model) {
  selectLambda <- model@matrices$select$lambda
  modes        <- model@info$modes
  lvs.linear   <- model@info$lvs.linear
  lambda       <- selectLambda

  for (j in lvs.linear) {
    op <- switch(modes[[j]], A = "=~", B = "<~", "=~")
    for (i in rownames(lambda))
      lambda[i, j] <- paste0(j, op, i)
  }

  selectGamma <- model@matrices$select$gamma
  gamma       <- selectGamma
  for (j in colnames(gamma)) for (i in rownames(gamma))
    gamma[i, j] <- paste0(j, "~", i)

  selectCov <- model@matrices$select$cov
  psi       <- selectCov
  for (j in colnames(psi)) for (i in rownames(psi))
    psi[i, j] <- paste0(j, "~~", i)

  selectTheta <- model@matrices$select$theta
  theta       <- selectTheta
  for (j in colnames(theta)) for (i in rownames(theta))
    theta[i, j] <- paste0(j, "~~", i)

  thresholds <- model@thresholdStruct@thresholds
  customParams <- names(model@matrices$customExpressions)

  c(
    lambda[selectLambda],
    gamma[selectGamma],
    psi[selectCov],
    theta[selectTheta],
    names(thresholds),
    customParams
  )
}


getParamVecLabels <- function(model) {
  parTable <- addReverseCovariancesToParTable(
    model@parTableInput
  )

  nm <- getParNamesFromParTable(parTable)
  lab <- parTable$label
  keep <- lab != ""

  stats::setNames(lab[keep], nm = nm[keep])
}


extractCoefs <- function(model) {
  fit <- model@fit
  thresholdStruct <- model@thresholdStruct

  lambda       <- fit$fitMeasurement
  selectLambda <- model@matrices$select$lambda

  gamma       <- fit$fitStructural
  selectGamma <- model@matrices$select$gamma

  fitCov    <- fit$fitCov
  selectCov <- model@matrices$select$cov

  fitTheta    <- fit$fitTheta
  selectTheta <- model@matrices$select$theta

  thr  <- thresholdStruct@thresholds
  pars <- c(
    lambda[selectLambda],
    gamma[selectGamma],
    fitCov[selectCov],
    fitTheta[selectTheta]
  )

  names(pars) <- model@params$names[seq_along(pars)]

  custom <- evalCustomExpressions(
    pars = c(pars, thr), labels = model@params$labels,
    expressions = model@matrices$customExpressions
  )

  plssemVector(c(pars, thr, custom))
}


computeFactorScores <- function(model) {
  W <- model@matrices$lambda
  X <- model@data

  if (model@info$is.probit) {
    ordered <- model@info$ordered

    for (ord in ordered)
      X[,ord] <- plsMapOrderedToExpectations(X[,ord])
  }

  F <- X %*% W

  if (!model@info$standardized || model@info$is.probit)
    F <- Rfast::standardise(F)

  F
}


getEstimatorFromInfo <- function(info) {
  consistent <- info$consistent
  is.mcpls   <- info$is.mcpls
  is.mlm     <- info$is.mlm
  is.ord     <- info$is.probit || (info$is.mcpls && length(info$ordered))

  estimator <- "PLS"
  if (consistent || is.mcpls) estimator <- paste0(estimator, "c")
  if (is.mlm)                 estimator <- paste0(estimator, "-MLM")
  if (is.ord)                 estimator <- paste0("Ord", estimator)
  if (is.mcpls)               estimator <- paste0("MC", estimator)

  estimator
}


refreshModelParams <- function(model, update.names = TRUE) {
  # Should we update names?
  if (update.names) {
    model@params$names <- getParamVecNames(model)
    model@params$labels <- getParamVecLabels(model)
  }

  # Single level params
  model@params$values <- extractCoefs(model)
  model@params$se     <- rep(NA_real_, length(model@params$values))

  # Multilevel/Mixed-Effect params
  if (isMLM(model))
    model <- refreshLmerParams(model)

  model
}


refreshLmerParams <- function(model) {
  lmerFit <- modelFitLmer(model)

  if (!isMLM(model) || is.null(lmerFit))
    return(model)

  coefs.x <- model@params$values
  coefs.y <- lmerFit$values

  common  <- intersect(names(coefs.x), names(coefs.y))
  new     <- setdiff(names(coefs.y),   names(coefs.x))

  coefs.x[common] <- coefs.y[common]
  coefs.all       <- c(coefs.x, coefs.y[new])

  model@params$values <- plssemVector(coefs.all)
  model@params$se     <- rep(NA_real_, length(coefs.all))

  model
}


plsMatricesLavRep <- function(object) {
  combined <- combinedModel(object)
  fit      <- modelFit(combined)

  list(
    lambda = fit$fitLambda,
    wmat   = fit$fitWeights,
    theta  = fit$fitTheta,
    C      = fit$fitC,
    psi    = fit$fitCov,
    gamma  = fit$fitStructural
  )
}