R/modelAveraging.R

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

# last changed on # Wed May 26 20:21:33 2021 ------------------------------

#' @title modelAveraging
#'
#' @description Apply Bayesian model averaging for aggregating the results of
#' multiple partial MI models. Uses STAN. Can currently deal with two-group models
#' for metric items or dichotomous items (Rasch or 2PL model).
#'
#' @param res_clusterItems Object generated by \code{\link{clusterItems}}. Delivers the
#' clustering structure as well as the model setup.
#' @param res_modelAveraging  Object generated by previous use \code{\link{modelAveraging}}.
#' An alternative to res_clusterItems.
#' Different weights can thereby be applied in a timely manner.
#' @param weights Weights in numerical order of the clusters. This represents the believe of the researcher in
#' the plausibility of the item cluster as appropriate anchors (see Schulze & Pohl, 2021).
#' With no prior knowledge or beliefs, weights should be chosen to be equally distributed across all clusters.
#' Weights have to sum to one.
#' @param CPUcores Number of CPU cores to be used. Defaults to 2.
#' @param chains Number of chains. Defaults to 2.
#' @param iter Number of total iterations. Defaults to 2'000 for metric items (blavaan)
#' and 10'000 for dichotomous items (IRT models)
#' @param burninPerc Percentage of burn-in iteration of the total iterations. Defaults to 0.5.
#' @param silent Do not print summary output? Defaults to \code{FALSE}.
#' @param anchors A list of vectors 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 ... Arguments of \code{\link{testMI}}, if \code{\link{testMI}} or
#' \code{\link{clusterItems}} has not
#' been called before to describe a measurement model.
#'
#' @references
#' Schulze, D., & Pohl, S. (2021). Measurement Invariance: Dealing with the uncertainty
#' in anchor item choice by model averaging. [Manuscript submitted]
#'
#' @return Bayesian Average of the mean difference of two groups ("muAver"). Info on
#' convergence and precision. All fitted STAN objects.
#'
#' @usage bma <- modelAveraging(res_clusterItems = NULL,
#'                              res_modelAveraging = NULL,
#'                              weights = NULL,
#'                              CPUcores = 2,
#'                              chains = 2,
#'                              iter = NULL,
#'                              burninPerc = 0.5,
#'                              silent = FALSE,
#'                              partialMIwhat,
#'                              anchors = NULL,
#'                              ...)
#'
#' @details If multiple cores are used (like in the default), STAN will maximize
#' the Viewer Tab of RStudio.
#'
#' @importFrom reshape2 "melt"
#'
#' @export

modelAveraging <- function(res_clusterItems = NULL,
                           res_modelAveraging = NULL,
                           weights = NULL,
                           CPUcores = 2,
                           chains = 2,
                           iter = NULL,
                           burninPerc = 0.5,
                           silent = FALSE,
                           # if used stand alone
                           partialMIwhat,
                           anchors = NULL,
                           # inherited arguments from clusterItems
                           items = NULL,
                           data = NULL,
                           group = NULL,
                           MIlevel = NULL,
                           stdItems = TRUE,
                           itemType = "metric",
                           dichModel = "factor") {

  res <- res_clusterItems
  partialMIwhat <- res_clusterItems$Factor$clusterWhat

  if (is.null(res_modelAveraging) && is.null(res_clusterItems)) {
    MItargetLevel <- MIlevel
    dich <- ifelse(itemType == "metric", FALSE, TRUE)

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

  if (!is.null(res_modelAveraging)) { # if only a reweighting of previously estimated models is to be done
    res <- list(data = res_modelAveraging$data,
                model = res_modelAveraging$model,
                lvs = res_modelAveraging$lvs)
    res[[res_modelAveraging$lvs]]$itemClustering$finalClustering <- res_modelAveraging[[res_modelAveraging$lvs]]$itemClustering$finalClustering
    iter <- res_modelAveraging$settings$iter
    fits <- res_modelAveraging$fittedModels
  }

  if (is.list(anchors)) { # if function is used standalone and anchor sets are given
    clus <- rep(NA, length(items))
    for (l in 1:length(anchors)) {
      clus[items %in% anchors[[l]]] <- l
    }
    res$Factor$itemClustering$finalClustering <- clus
  }

  # checks
  if ((!is.null(res_modelAveraging) | !is.null(res_clusterItems)) &&
      (!is.null(items) | !is.null(data) | !is.null(group))) {
    stop("Don't give res_clusterItems or res_modelAveraging together with items, data, or group arguments.", call. = FALSE)
  }
  if (!is.null(anchors) && !is.list(anchors)) {
    stop("anchors argument should be a list.", call. = FALSE)
  }
  if (sum(weights) != 1) {
    stop("Weights do not sum to 1.", call. = FALSE)
  }
  if (length(weights) != length(unique(res$Factor$itemClustering$finalClustering))) {
    stop("Number of weights does not fit number of clusters.", call. = FALSE)
  }
# if Bayesian estimation is actually to be done
  if (is.null(res_modelAveraging)) {
    # fit multiple partial MI models to cover item clusters

    #checks
    if (res$model$dich &&
        res$model$dichModel == "factor") {
      stop("Categorial SEM is not yet supported. Choose an IRT model instead.",
           call. = FALSE)
    }
    if (length(names(res$model$items)) > 1) {
      stop("Only unidimensional models are supported so far.", call. = FALSE)
    }

    fits <- list()
    clusters <- unique(res$Factor$itemClustering$finalClustering)
    for (currCluster in clusters[order(clusters)][!(weights %in% 0)]) { # estimate only those anchors for which the weight is != 0

      # Time estimation
      if (currCluster == clusters[order(clusters)][!(weights %in% 0)][1]) {
        cat(paste0("### Model for cluster ", currCluster,
                   ". Total time is estimated...\n"))
        allTimes <- rep(0, length(clusters))
      } else {
        allTimes[currCluster] <- mean(rowSums(time))
        finTime <- Sys.time() +
          ceiling(chains / CPUcores) *
          mean(allTimes) *
          (max(res$Factor$itemClustering$finalClustering) - currCluster + 2) # + 2 for correcting for post-processing
        finTime <- format(finTime, "%H:%M")
        cat(paste0("\n### Model for cluster ", currCluster,
                   ". Estimated time at finish: ", finTime, "\n"))
      }

  # SEM models
      if (!res$model$dich) {
        if (is.null(iter)) iter <- 2000
        partialItems <- NULL
        currNonAnchor <- which(res$Factor$itemClustering$finalClustering != currCluster)
        if ("loadings" %in% partialMIwhat) {
          partialItems <- append(partialItems,
                                 paste0("Factor =~", res$model$items$Factor[currNonAnchor]))
        }
        if ("difficulties" %in% partialMIwhat) {
          partialItems <- append(partialItems,
                                 paste0(res$model$items$Factor[currNonAnchor], "~ 1"))
        }

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

        model <- NULL
        model <- append(model, paste0("Factor =~", paste(res$model$items$Factor, collapse = " + ")))
        model <- paste0(model, collapse = "\n")

        future::plan("multicore")
        fits[[paste0("Cluster=", currCluster)]] <- bcfa(model,
                                                        data = res$data,
                                                        group = res$model$group,
                                                        std.lv = TRUE,
                                                        group.equal = groupEqual,
                                                        group.partial = partialItems,
                                                        burnin = iter*burninPerc,
                                                        sample = iter*(1-burninPerc),
                                                        n.chains = chains,
                                                        bcontrol = list(cores = CPUcores))

        # convergence and precision checks directly after single model
        conv <- rstan::summary(fits[[paste0("Cluster=", currCluster)]]@external$mcmcout)
        if (max(conv$summary[, 10]) > 1.10) {
          message(paste0("Convergence was not reached. (max PSR ",
                         round(max(conv$summary[, 10]), 2),
                         "). Try increasing iterations."))
        } else {
          if (min(conv$summary[, 9]) < 400) {
            message(paste0("Precision of the posterior distribution was below the recommended cut off of 400 effective samples. Try increasing iterations to ",
                           round((iter*(400/min(conv$summary[, 9])))/1000, 0)*1000, "."))
          }
        }
      }

  # IRT models
      if (res$model$dich) {
        if (is.null(iter)) iter <- 10000

        currAnchor <- which(res$Factor$itemClustering$finalClustering %in% currCluster)

        # preparing the data for STAN
        group <- res$data[, (res$model$group)]
        group <- as.numeric(as.factor(group)) # make sure, group is coded by positive numbers

        data <- res$data[, unlist(res$model$items)]
        data$pbn <- 1:nrow(data)
        dat0 <- reshape2::melt(data, value.name = "y",  # long data format
                               id.vars = "pbn",
                               variable.name = "item")
        gg <- rep(group, length(unlist(res$model$items)))[!is.na(dat0$y)] # response-wise grouping variable filtered for missing responses
        dat0 <- na.omit(dat0) # yields a FIML-type treatment of missing values in STAN
        dat0$item <- as.numeric(dat0$item)

        rstan_options(auto_write = TRUE)

        # STAN model input
        data2G <- list(N    = nrow(data),
                      K    = length(unlist(res$model$items)),
                      Ntot = nrow(dat0),
                      jj   = dat0$item,
                      ii   = dat0$pbn,
                      y    = dat0$y,
                      g    = group,
                      gg   = gg,
                      G    = 2,
                      NconItem = length(currAnchor),
                      conItem  = as.array(currAnchor),
                      freeItem = as.array(which(!(seq(1,length(
                        unlist(res$model$items))) %in% currAnchor)))) # as.array is needed to deal with single items for STAN

        if (res$model$dichModel == "2PL") {
          model <- "data{
                    int<lower = 1> N; // number of examinees
                    int<lower = 1> K; // number of items
                    int<lower = 1> Ntot; // number of data points
                    int<lower = 1> jj[Ntot]; // item id
                    int<lower = 1> ii[Ntot]; // person id
                    int<lower = 0> y[Ntot]; // responses
                    int<lower = 1> g[N]; // group
                    int<lower = 1> gg[Ntot]; // group
                    int<lower = 1> G; // number of groups
                    int<lower=1> NconItem; // number of constrained items
                    int<lower=1> conItem[NconItem]; // item constraint
                    int<lower=1> freeItem[K - NconItem];
                  }
                  parameters{
                    vector[N] PersPar; // person parameters
                    real Alpha_free; // freely estimated
                    real<lower=0> Psi_var; // freely estimated
                    vector[NconItem] equalB;
                    matrix[K - NconItem, 2] unequalB;
                    vector<lower=0>[NconItem] equalV;
                    matrix<lower=0>[K - NconItem, 2] unequalV;
                  }
                  transformed parameters{
                    matrix<lower=0>[K, G] v;
                    matrix[K, G] b;
                  b[conItem, 1] = equalB;
                  v[conItem, 1] = equalV;
                  b[conItem, 2] = equalB;
                  v[conItem, 2] = equalV;
                  b[freeItem, ] = unequalB;
                  v[freeItem, ] = unequalV;
                  }
                  model{
                    // prior person parameter
                    for(i in 1:N){
                      PersPar[i] ~ normal((g[i]-1)*Alpha_free, (g[i]-1)*Psi_var+((g[i]-2)*(-1)));
                     }
                      Alpha_free ~ normal(0, 5);
                      Psi_var ~ cauchy(0, 5);
                    // prior item parameter
                    to_vector(unequalB) ~ normal(0, 5); //  difficulties
                    equalB ~ normal(0, 5);
                    to_vector(unequalV) ~ normal(0, 5); //  discriminations
                    equalV ~ normal(0, 5);
                  	for(n in 1:Ntot){
                  	  target += bernoulli_logit_lpmf(y[n] | (v[jj[n], gg[n]]*PersPar[ii[n]] - b[jj[n], gg[n]]));
                  	}
                  }"
        }
        if (res$model$dichModel == "Rasch") {
          model <- "data{
                    int<lower = 1> N; // number of examinees
                    int<lower = 1> K; // number of items
                    int<lower = 1> Ntot; // number of data points
                    int<lower = 1> jj[Ntot]; // item id
                    int<lower = 1> ii[Ntot]; // person id
                    int<lower = 0> y[Ntot]; // responses
                    int<lower = 1> g[N]; // group
                    int<lower = 1> gg[Ntot]; // group
                    int<lower = 1> G; // number of groups
                    int<lower=1> NconItem; // number of constrained items
                    int<lower=1> conItem[NconItem]; // item constraint
                    int<lower=1> freeItem[K - NconItem];
                  }
                  parameters{
                    vector[N] PersPar; // person parameters
                    real Alpha_free; // freely estimated
                    vector<lower=0>[G] Psi_var; // freely estimated
                    vector[NconItem] equalB;
                    matrix[K - NconItem, 2] unequalB;
                  }
                  transformed parameters{
                    matrix[K, G] b;
                  b[conItem, 1] = equalB;
                  b[conItem, 2] = equalB;
                  b[freeItem, ] = unequalB;
                 }
                  model{
                    // prior person parameter
                    for(i in 1:N){
                      PersPar[i] ~ normal((g[i]-1)*Alpha_free, Psi_var[g[i]]);
                     }
                      Alpha_free ~ normal(0, 5);
                      Psi_var ~ cauchy(0, 5);
                    // prior item parameter
                    to_vector(unequalB) ~ normal(0, 5);
                    equalB ~ normal(0, 5);
                  	for(n in 1:Ntot){
                  	  target += bernoulli_logit_lpmf(y[n] | (PersPar[ii[n]] - b[jj[n], gg[n]]));
                  	}
                  }"
        }

        fits[[paste0("Cluster=", currCluster)]] <- stan(model_code = model,
                                                        data = data2G,
                                                        iter = iter,
                                                        warmup = iter*burninPerc,
                                                        chains = chains,
                                                        cores = CPUcores)

        # convergence and precision checks
        conv <- rstan::summary(fits[[paste0("Cluster=", currCluster)]])
        if (max(conv$summary[, 10]) > 1.10) {
          message(paste0("Convergence was not reached. (max PSR ",
                         round(max(conv$summary[, 10]), 2),
                         "). Try increasing iterations."))
        } else {
          if (min(conv$summary[, 9]) < 400) {
            message(paste0("Precision of the posterior distribution was below the recommended cut off of 400 effective samples. Try increasing iterations to ",
                           round((iter*(400/min(conv$summary[, 9])))/1000, 0)*1000, "."))
          }
        }
      }
      if (!res$model$dich) {
        fit <- fits[[paste0("Cluster=", currCluster)]]@external$mcmcout
      }
      if (res$model$dich) {
        fit <- fits[[paste0("Cluster=", currCluster)]]
      }
      time <- get_elapsed_time(fit)
    }
  }

  muAver <- NULL
  Rhat <- NULL  # set up storage for measures of convergence and efficiency
  Neff <- NULL

  # Bayesian model averaging by weighting
  for (currCluster in unique(res$Factor$itemClustering$finalClustering)) {
    if (!res$model$dich) {
      fit <- fits[[paste0("Cluster=", currCluster)]]@external$mcmcout
    }
    if (res$model$dich) {
      fit <- fits[[paste0("Cluster=", currCluster)]]
    }

    conv <- rstan::summary(fit)
    Neff <- append(Neff, min(conv$summary[, 9]))
    Rhat <- append(Rhat, max(conv$summary[, 10]))

    post <- extract(fit)
    nSims <- weights[currCluster]*(iter/2)
    draws <- sample((iter/2):iter, nSims)

    if (res$model$dichModel != "Rasch") {
      muAver <- append(muAver, post$Alpha_free[draws])
    }
    if (res$model$dichModel == "Rasch") {
      muAver <- append(muAver, post$Alpha_free[draws])
    }
  }
  resBMA <- list(muAver = muAver,
                 Rhat = Rhat,
                 Neff = Neff,
                 fittedModels = fits,
                 data = res$data,
                 lvs = "Factor",
                 model = res$model)
  resBMA$settings$iter <- iter
  resBMA$Factor$itemClustering$finalClustering <- res$Factor$itemClustering$finalClustering

  if (!silent) {
    cat(paste0("\nMaximum Rhat (aka PSR): ", round(max(Rhat), 2), " [recommended: < 1.05]"))
    cat(paste0("\nMinimum effective sample size (aka N_eff): ", round(min(Neff), 0), " [recommended: > 400]\n"))
    cat("\nMean difference in the latent variable after Bayesian model averaging:\n")
    df <- data.frame(round(c(muAver[order(muAver)][0.025*length(muAver)],
                             mean(muAver),
                             muAver[order(muAver)][0.975*length(muAver)]), 3))
    rownames(df) <- c(" 2.5% cred. int.",
                      "mean",
                      "97.5% cred. int.")
    colnames(df) <- ""
    print(df)
     }

  return(resBMA)
}