R/partialMI.R

In Dani-Schulze/measurementInvariance: (Partial) measurement invariance analysis

# Fri Sep 17 13:29:30 2021 ------------------------------

#' @title partial MI model
#'
#' @description Estimate a single unidimensional partial MI model with a given anchor item set (aka cluster).
#' Can currently deal with two-group models
#' for metric items or dichotomous items (Rasch, 2PL, or probit model).
#'
#' @param res_clusterItems Optional object generated by \code{\link{clusterItems}}.
#' Alternatively, a MI model
#' can be described using the arguments of \code{\link{testMI}} (see ...).
#' @param anchor Either a number indicating which cluster from the results of \code{\link{clusterItems}}
#' to use for anchoring, or a vector of item names which shall serve as anchors.
#' @param partialMIwhat String, either \code{"loadings"}, \code{"difficulties"}, or
#' \code{c("loadings", "difficulties")}. Has not to be specified when using res_clusterItems but when setting up
#' a new MI model. This is an equivalent to clusterWhat in \code{\link{clusterMI}}.
#' @param silent Do not print summary output? Defaults to \code{FALSE}.
#' @param ... Arguments of \code{\link{testMI}}, if \code{\link{testMI}} or
#' \code{\link{clusterItems}} has not
#' been called before to describe a measurement model.
#'
#' @return Parameter estimates of the model as defined by the model argument.
#' Results are directly printed to the console but can be further inspected by l
#' avaan's or mirt's functions.
#'
#' @usage res <- partialMI(res_clusterItems = NULL,
#'                         anchor = NULL,
#'                         partialMIwhat = NULL,
#'                         silent = FALSE,
#'                         ...)
#'
#' @export


partialMI <- function(res_clusterItems = NULL,
                      anchor,
                      partialMIwhat,
                      silent = FALSE,
                      # inherited arguments from clusterItems
                      items,
                      data,
                      group,
                      MIlevel,
                      stdItems = TRUE,
                      itemType = "metric",
                      dichModel = "factor") {

  if (!is.null(res_clusterItems)) {
    res <- res_clusterItems
    partialMIwhat <- res_clusterItems$Factor$clusterWhat
  } else { # if clusterItems has not been called before
     #pars0 <- lavParTable(model) # test if model fits res_clusterItems
     #lvsModel <- unique(pars0[pars0$op == "=~", "lhs"])
     #if (any(lvs[order(lvs)] != lvsModel[order(lvsModel)])) {
     #   stop("Factors in model do not match factors from previous clusterItems/testMI model.",
     #         call. = FALSE)
   MItargetLevel <- MIlevel
   dich <- ifelse(itemType == "metric", FALSE, TRUE)

   if (is.numeric(anchor)) {
     stop("No input from clusterItem is given, thus give items as anchor, not cluster numbers.",
          call. = FALSE)
   }

   res <- setUpModel(items = items,
                      data = data,
                      group = group,
                      MItargetLevel = MItargetLevel,
                      stdItems = stdItems,
                      dich = dich,
                      dichModel = dichModel)
   }
  items <- res$model$items
  lvs <- names(items)

    #  itemsModel <- list()
    #  for (lv in lvs) { # get item assignments in model
    #    itemsModel[[lv]] <- unique(pars0[pars0$op == "=~" & pars0$lhs == lv, "rhs"])
    #  }
    #  if (any(!(unlist(items) %in% unlist(itemsModel))) |
    #      any(!(unlist(itemsModel) %in% unlist(items)))) {
    #    stop(paste0("Items in model do not match items from previous clusterItems/testMI model.",
    #                " Problem involving item ",
    #                unlist(items)[!(unlist(items) %in% unlist(itemsModel))],
    #                unlist(itemsModel)[!(unlist(itemsModel) %in% unlist(items))], ".\n"),
    #         call. = FALSE)
    #  }
    #}
  if (is.null(items) && is.null(res_clusterItems)) {
    stop("Input is missing. Give either res_clusterItems or items, data, and so forth.",
         call. = FALSE)
  }

  clusters <- list()  # internal representation
  clusters$Factor <- anchor
  ## old switch for getting desired clusters per factor. A rather badly written section.
  #lvsPlaces <- lapply(res$lvs, FUN = regexpr, anchorClusters) # parse input string
  #lvsPlaces <- append(lvsPlaces, nchar(anchorClusters) + 1) # add final stop in order to make loop below work
  #if (any(lvsPlaces < 0) &&
  #    nchar(anchorClusters) > (sum(nchar(res$lvs)) + 2*length(res$lvs))) {
  #  stop("Error while reading anchor. Couldn't retrieve factor names.", call. = FALSE)   # check for misspelled factor names
  #}
  #for (lv in res$lvs) {
  #  clusters[[lv]] <- as.numeric(substring(anchorClusters, lvsPlaces[[which(
  #    res$lvs %in% lv)]] + attr(lvsPlaces[[which(res$lvs %in% lv)]], "match.length"),
  #    lvsPlaces[[which(res$lvs %in% lv) + 1]] - 1))
  #}

  #if (any(lapply(clusters, FUN = length) == 0)) {
  #  stop("Error while reading anchor. Check your statements.")
  #}

# lavaan
  if (res$model$dichModel == "factor") {
    partialItems <- NULL
    for (lv in lvs) {
      if (is.numeric(anchor)) { # check whether anchor gives a cluster from clusterItems or gives items
        currNonAnchor <- which(res[[lv]]$itemClustering$finalClustering != clusters[[lv]])
      } else {
        currNonAnchor <- which(!(res$model$items[[lv]] %in% anchor))
      }

      if ("loadings" %in% partialMIwhat) {
        partialItems <- append(partialItems,
                              paste0(lv, "=~", res$model$items[[lv]][currNonAnchor]))
      }
      if ("difficulties" %in% partialMIwhat) {
        partialItems <- append(partialItems,
                               paste0(res$model$items[[lv]][currNonAnchor], "~ 1"))
      }
    }

    int <- "intercepts"
    if (res$model$dich) {
      int <- "thresholds"
    }
    catItems <- NULL
    if (res$model$dich) {
      catItems <- unlist(res$items)
    }

    if (res$Factor$MItargetLevel == "weak") {
      groupEqual <-  "loadings"
    }
    if (res$Factor$MItargetLevel == "strong") {
      groupEqual <-  c("loadings", int)
    }

    model <- NULL
    model <- append(model, paste0(lvs, "=~", paste(items[[lvs]], collapse = " + ")))
    model <- paste0(model, collapse = "\n")

    output <- cfa(model,
                  data = res$data,
                  estimator = res$model$estim,
                  group = res$model$group,
                  meanstructure = TRUE,
                  std.lv = TRUE,
                  ordered = catItems,
                  group.equal = groupEqual,
                  group.partial = partialItems,
                  missing = res$model$missings)
  } else {
# mirt
    if (length(lvs) > 2) {
      print("This will take a while..." )
    }
    dat0 <- res$data[, -which(colnames(res$data) %in% res$model$group)] # prepare data for mirt
    dat0 <- dat0[, unlist(res$model$items)]
    partialItems <- NULL
    mod <- NULL
    invariance <- NULL
    for (lv in lvs) {
      if (is.numeric(anchor)) { # check whether anchor gives a cluster from clusterItems or gives items
        partialItems <- append(partialItems,
                               res$model$items[[lv]][which(res[[lv]]$itemClustering$finalClustering %in% clusters[[lv]])])
      } else {
        partialItems <- append(partialItems,
                               res$model$items[[lv]][which(res$model$items[[lv]] %in% anchor)])
      }

      mod <- append(mod, paste0(lv, " = ", paste0(res$model$items[[lv]], collapse = ",")))
    }

    if ("loadings" %in% partialMIwhat) {
      mod <- append(mod, paste0("CONSTRAINB = (", paste0(partialItems, collapse = ","), ", a1)"))
    }
    if ("difficulties" %in% partialMIwhat) {
      mod <- append(mod, paste0("CONSTRAINB = (", paste0(partialItems, collapse = ","), ", d)"))
    }
    mod <- paste0(mod, collapse = "\n")

    if (res$Factor$MItargetLevel == "weak") {
      invariance <-  append(invariance, "free_variances")
    }
    if (res$Factor$MItargetLevel == "strong") {
      invariance <-  append(invariance,c("free_variances", "free_means"))
      if ("difficulties" == partialMIwhat) {
        invariance <-  append(invariance, "slopes")
      }
    }

    if ("difficulties" %in% partialMIwhat && res$Factor$MItargetLevel == "weak") {
      stop("MI Level statement (weak) and partial MI statement (difficulties) makes no sense in combination",
           call. = FALSE)
    }

    output <- multipleGroup(data = dat0,
                            model = mod,
                            itemtype = res$model$dichModel,
                            method = res$model$estim,
                            group = as.factor(res$data[, res$model$group]),
                            invariance = invariance,
                            SE = TRUE,
                            verbose = FALSE)
  }

  if (!silent) {
    if (res$model$dichModel == "factor") {      # lavaan
      lavaan::summary(output, fit.measures = TRUE)
    } else {                                       # mirt
      cat("mirt model estimates\n\n")
      header <- capture.output(show(output))
      header <- header[grep("Full-info", header):(length(header) - 2)]
      cat(paste(header, collapse = "\n"))
      cat("\n")
      if (res$model$missings == "none") {
        print(t(round(M2(output), 3)))
      } else {
        cat("Further fit statistics could not be computed due to missing values.")
      }
      cat("\n\nParameter Estimates:\n\n")
      print(coef(output, simplify = TRUE))
    }
  }

  return(output)
}