AUTTT: A Unified Treatment in Test Theory

Documented in mplus_montecarlo_analysis_grm

#'A Function for Analyzing Monte Carlo Data Sets Created by mplus_montecarlo_data() or simdata_grm().
#'
#' This package contains the mplus_montecarlo_analysis_grm() function,
#' which will analyze Monte Carlo data sets using the script files for various estimators. The data sets and the script files are generated by the mplus_montecarlo_data() function.
#' @param model_object A set of model specifications that can be passed from the function mplus_montecarlo_data() or simdata_grm().
#' @param estimators A list of estimators. Available estimators are: c('ML_logit',
#' 'ML_probit',
#' 'MLR_logit',
#' 'MLR_probit',
#' 'MLF_logit',
#' 'MLF_probit',
#' 'WLS_delta',
#' 'WLS_theta',
#' 'WLSM_delta',
#' 'WLSM_theta',
#' 'WLSMV_delta',
#' 'WLSMV_theta',
#' 'ULS_delta',
#' 'ULS_theta',
#' 'ULSMV_delta',
#' 'ULSMV_theta')
#' @param rep a number of replication of data set if the analysis is done on each replicated data set separately. Null if type_montecarlo = TRUE.
#' @param type_montecarlo if TRUE, the analysis is done using the type of MONTECARLO, which analyze all replications and give average parameters based on the number of replications.
#' @param run_files if TRUE, it will execute the Mplus script files based on the selected estimators.
#' @return mplus_montecarlo_analysis will analyze Monte Carlo data sets using the selected estimators. It uses
#'         MplusAutomation::runModels().
#' @export
#'
#' @examples
#' library(MplusAutomation)
#' data.2F = mplus_montecarlo_data(
#' # list your model in increasing order
#' model = list(c(1,2,3,4), #factor 1
#'               c(5,6,7,8)) # factor 2
#'eloadval = NULL,
#'vloadings = c(0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5), # a vector of factor loadings
#'thresholds = c(-2.478,	-1.818,	-0.865,	0.425), # threshold values per item for 1,2,3,4
#'factor.cor = NULL,# factor correlation
#'vfactor.cor = c(0.31), # length(vfactor.cor)
#'N = 300, # sample size
#'R = 5,# n of replications
#'seed_mplus = 4567,# seed number for Mplus
#'naming_data_files = 'data4_rep*.dat;',
#'file_dir = getwd(),
#'file_name = "generate_data4.inp",
#'run_files = TRUE
#')

#'
#' # running analysis as Monte Carlo type
#' setwd("C:/Users/Dell/mplus files")
#' mplus_montecarlo_analysis(
#' model_object = data.2F,
#' estimators = c('ULS_delta',
#'               'ULS_theta',
#'               'WLSMV_theta',
#'               'WLSMV_delta'),
#' rep = NULL,
#' type_montecarlo = TRUE,
#' run_files = TRUE)
#'
#'# running each replicate separately to get theta values
#'#' setwd("C:/Users/Dell/mplus files")
#' mplus_montecarlo_analysis(
#' model_object = data.2F,
#' estimators = c('ULS_delta',
#'               'ULS_theta',
#'               'WLSMV_theta',
#'               'WLSMV_delta'),
#' rep = 1,
#' type_montecarlo = FALSE,
#' run_files = TRUE)
#'
#' # running replicated data sets separately using loop
#'#' setwd("C:/Users/Dell/mplus files")
#' for (i in 1:10){
#' mplus_montecarlo_analysis(
#' model_object = data.2F,
#' estimators = c('ULS_delta',
#'               'ULS_theta',
#'               'WLSMV_theta',
#'               'WLSMV_delta'),
#' rep = i,
#' type_montecarlo = FALSE,
#' run_files = TRUE)
#' }
#' @references
#' \insertRef{Asparaouhov2020}{MplusMontecarlo}
#' \insertRef{Hallquist2018}{MplusMontecarlo}


mplus_montecarlo_analysis_grm <- function(
                                      model_object,
                                      estimators,
                                      rep,
                                      type_montecarlo,
                                      run_files
                                    ){

  # declare variables ####
  model <- model_object$model_spec$model
  thresholds <- model_object$model_spec$thresholds
  vfactor.cor <- model_object$model_spec$vfactor.cor
  R <- model_object$model_spec$R
  N <- model_object$model_spec$N
  vloadings <- model_object$model_spec$vloadings
  file_dir <- model_object$model_spec$file_dir
  naming_data_files <- model_object$model_spec$naming_data_files
  file_name <- model_object$model_spec$file_name



  # Check: model ####
  # Ensure that each item is unique to each factor
  # Check: one item in no more than one factor and integer ####
  items <- unlist(model)
  items_rank <- rank(items)
  for (i in 1:length(items_rank)){
    if(items_rank[i]%%1 != 0){
      stop("Currently, this function does not support having one item in more than one dimension. Or The numbers in model must integer, not fraction.")
    } else if (items[i]%%1 !=0){
      stop("Currently, this function does not support having one item in more than one dimension. Or The numbers in model must integer, not fraction.")
    }
  }


  # Check order of items in model ####
  a <- unlist(model) # original input
  b <- rank(unlist(model)) # ranks of input
  for (i in 1:length(unlist(model))){
    if (a[i] - b[i] !=0){
      stop("In specifying your model, you must input the whole numbers, which represent items, in a consecutive, increasing order from 1, 2, ..j, where j is the last item and represents the the total number of items.")
    }
  }
  # check: order of items in model ####
  for (i in 1:(length(unlist(model))-1)){
    if (a[i+1]-a[i] != 1) {
      stop("In specifying your model, you must input the whole numbers, which represent items, in a consecutive, increasing order from 1, 2, ..j, where j is the last item and represents the the total number of items.")
    }
  }

  # Check: nFs, n of factor ####
  nFs <- length(model)
  # n of factors (nFs)
  if (nFs > 0 & nFs <= 10 & nFs%%1==0 ){
    nFs = nFs
  } else {
    stop("The number of factors (nFs) must be a positive integer between 0 and 10.")
  }

  # Check: nYs, n of items ####
  nYs <- max(unlist(model))

  # n of observed indicators (nYs)
  if (nYs > 0 & nYs%%1==0 ){
    nYs = nYs
  } else {
    stop("The number of observed indicators (nYs) must be a positive integer.")
  }


  # Check: threshold values ####
  if (is.null(thresholds)){
    stop("A vector of thresholds is needed.")
  } else {
    thresholds = thresholds
    nTs <- length(thresholds)
  }

  # check: n of factor correlations ####
  # iterate pair-wise factor correlations
  Fcor.tb <- NULL
  for (i in 1:(nFs-1)){
    for (j in (i+1):nFs){
      Fcor.tb <- rbind(Fcor.tb, c(i, j))
    }
  }

  # Check: n of factor correlations given n of factors ####
  # factor correlation less than 1
  for (i in 1:length(vfactor.cor)){
    if(vfactor.cor[i]>=1){
      stop(paste("Factor correlation values must be less than 1. Check value number", i, '.', sep=' '))
    }
  }

  if (is.null(vfactor.cor)){
    stop("A vector of vfactor.cor is needed.")
  } else if (nrow(Fcor.tb)!=length(vfactor.cor)){
    stop(paste("The length of factor correlation vector (vfactor.cor) must be ",
               nrow(Fcor.tb), ",(", nFs, "*", (nFs-1), "/2)",
               ",given the number of factors:",nFs, sep=' '))
  } else {
    vfactor.cor = vfactor.cor
  }

  # Check: n of vloadings = n of items ####
  if (length(vloadings)!=length(unlist(model))){
    stop("The number of factor loadings must equal the number of items in the model.")
  }

  # Check: loadings, vloadings ####
  # On equal factor loadings and a vector of factor loadings and error variance
  # eloadval <- NULL
  # if(is.null(vloadings)&!is.null(eloadval))  { # vloadings is null and eloadval is NOT null
  #   vloadings <- rep(eloadval, nYs)
  #   errval <- 1-eloadval^2
  #   errvars <- rep(errval, nYs)
  # } else if (!is.null(vloadings)&is.null(eloadval))  {# vloadings is NOT null and eloading is null
  #   vloadings <- vloadings
  #   errvars <- 1-vloadings^2
  # } else if (is.null(vloadings)&is.null(eloadval)) {
  #   stop('A vector of factor loadings (vloadings)
  #         and an equal value for all factor loading (eloadval) are empty. Either one of them is needed.')
  # }

  # Check: N, sample size ####

  if (is.null(N)) {
    stop("A sample size that is greater than 0 is needed.")
  } else if (N <= 0) {
    stop("A sample size that is greater than 0 is needed.")
  } else {
    N <- N
  }

  # Check: R, n of replications ####

  if (is.null(R)) {
    stop("A replication number (R) must be a positive integer.")
  } else if (R <= 0) {
    stop("A replication number (R) must be a positive integer.")
  } else {
    R <- R
  }

  # # Check: seed_mplus ####
  #
  # if (is.null(seed_mplus)) {
  #   seed_mplus <- 4567
  # } else if (seed_mplus <= 0) {
  #   stop("A seed number for seed_mplus must be a positive integer.")
  # } else {
  #   R <- R
  # }

  # Check: file directory ####

  if (file.exists(file_dir)) {
    file_dir <- file_dir
  } else {
    stop("The directory is not valid. Please provide a valid directory.")
  }

  # Check: file name ####

  # if(ext <- unlist(strsplit(file_name, "\\."))[2]=="inp"){
  #   file_name <- file_name
  # } else {
  #   stop("The file name (file_name) must contain the file extension '.inp'; e.g., 'abcd.inp'.")
  # }

  # Check: Naming data files ####

  if(unlist(strsplit(naming_data_files, "\\*"))[2]== ".dat;"){
    naming_files <- naming_data_files
  } else {
    stop("Naming data files for replicated data sets must end with '*.dat;'; e.g., 'mplus_data*.dat;'")
  }

  # Check: n of Ys per factor ####

  FCTR <- factor(c(1:nFs), labels = 'F')
  Y <- factor(c(1:nYs), labels = 'y')

  if ((length(Y)/length(FCTR))<3){
    stop("There are too few observed indicators per factor.There should be 3 items per factor.")
  }
  L0 <- "! Author: Bo Klauth"
  L3 <- paste0('NAMES = ', Y[1], '-', Y[nYs], ';')
  L10 <-  paste0('CATEGORICAL= ', Y[1], '-', Y[nYs], ';')
  # # Write lines for Monte Carlo command ####
  # L0 <- "! Author: Bo Klauth"
  # # line 1: title
  # L1 <- paste0('TITLE: ', 'Generate data in an MS model, N = ', N, ', R = ', R)
  # # line 2: MONTECARLO command
  # L2 <- 'MONTECARLO:'
  # L3 <- paste0('NAMES = ', Y[1], '-', Y[nYs], ';')
  # L4 <- paste0('NOBSERVATIONS= ',N, ';')
  # L5 <- paste0('NREPS= ',R, ';')
  # L6 <- paste0('SEED= ',seed_mplus, ';')
  # L7 <- paste0('REPSAVE= All;')
  # L8 <- paste0('SAVE= ', naming_data_files)
  # L9 <- paste0('GENERATE= ',Y[1], '-', Y[nYs], '(', nTs, ')', ';')
  # L10 <-  paste0('CATEGORICAL= ', Y[1], '-', Y[nYs], ';')
  #
  # # Write MODEL POPULATION ####
  # L11 <- 'MODEL POPULATION:'
  # L12 <- '! Standardized Loadings'

  # formating into y1*loading1 y2*loading2 ...etc
  Yloadings <- paste(Y, vloadings, sep='*')


  # Write lines for factor loadings: FYlines ####
  # Convert input model into Mplus Syntax
  # unpack the factors

  D <- NULL
  F_label <- NULL
  FYlines<-NULL
  for (i in 1:length(model)){
    F_label[i] <- paste0('F', i)
    D[i] <- paste(F_label[i], 'BY', sep=' ')
    for (k in min(unlist(model[i])): max(unlist(model[i]))){
      D[i] <- paste(D[i], Yloadings[k], sep=' ') # iterate all loading lines per dimension
    }
    FYlines[i] <- paste0(D[i], ';')
  }

  FYlines

  # Write lines for factor variances: FVarlines ####
  L13 <- '! Factor variance'
  FVarlines <- NULL
  for (i in 1:length(model)){
    FVarlines[i] <- paste0('F', i, '@1;')
  }


  L14 <- '! Thresholds'

  # Write lines for thresholds ####

  Ythresholds <- paste0(
    rep(Y, each=nTs), # make Yj to match N thresholds
    rep('$',each =  nTs*nYs),
    rep(1:nTs, times= nYs),
    rep('*',each =  nTs*nYs),
    rep(thresholds, nYs))

  Ythresholds <- paste(
    paste(rep(Y, each=nTs), # make Yj to match N thresholds
          rep(1:nTs, times= nYs), sep='$'),
    rep(thresholds, nYs), sep='*')

  # turning into matrix with columns of 4
  Ythresholds <- matrix(Ythresholds, byrow = TRUE, ncol=nTs)


  # Lines of thresholds: Thlines
  Thlines <- NULL
  for (i in 1:nYs){
    Thlines[i] <-   paste0(c('[', Ythresholds[i,], '];'), collapse = ' ')
  }

  # Lines of Error Variance: Yerrvarlines ####
  L15 <- '! Error variance'

  # Yerrvars <- paste(Y, errvars, sep='*')
  # Yerrvarlines <- NULL
  # for (i in 1:length(model)){
  #   D[i] <- ""
  #   for (k in min(unlist(model[i])): max(unlist(model[i]))){
  #     D[i] <- paste(D[i], Yerrvars[k], sep=' ') # iterate all loading lines per dimension
  #   }
  #   Yerrvarlines[i] <- paste0(D[i], ';')
  # }

  # Interfactor corr lines: ####
  L16 <- '! Interfactor correlations'

  # references: vfactor.cor, Fcor.tb

  FCorlines <- NULL
  for (i in 1:nrow(Fcor.tb)){
    # use Fcor.tb to get labels of factors and match factor correlation
    FCorlines[i] <-paste0('F', Fcor.tb[i,1], ' WITH ', 'F', Fcor.tb[i,2], '*',
                          vfactor.cor[i], ';', sep=' ')
  }


  # # Write MODEL ####
  # L17 <- 'MODEL:'
  # # Acronyms ####
  # L18 <- "! MS: measurement model; N: sample size, R: number of replications."
  #
  # # consolidate lines
  # # use capture.output() to capture output using cat()
  # montecarlo_data <- capture.output(
  #   cat(L0, # Author info
  #       L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, # Monte Carlo Command
  #       L11, # MODEL POPULATION
  #       L12, FYlines, # factor loadings
  #       L13, FVarlines, # factor variance
  #       L14, Thlines, # Thresholds
  #       L15, Yerrvarlines,# Error variances
  #       L16, FCorlines, # Interfactor correlations
  #       L17, # MODEL
  #       L12, FYlines, # factor loadings
  #       L13, FVarlines, # factor variance
  #       L14, Thlines, # Thresholds
  #       # L15, Yerrvarlines,# Error variances
  #       L16, FCorlines, # Interfactor correlations
  #       sep='\n')
  # )
  #
  # # writing an mplus script to file
  # write.table(montecarlo_data,
  #             file = paste0(file_dir, "/", file_name),
  #             row.names = FALSE,
  #             col.names = FALSE,
  #             quote=FALSE) # use quote = FALSE to eliminate quotes in text
  #
  # # give message for the conditions
  # if (file.exists(paste0(file_dir, "/", file_name))) {
  #   message(paste('"', file_name,'" ', 'has been created in ', file_dir, '.',
  #                 collapse = '', sep=''))
  # } else {
  #   message(paste('Failed to create ', '"', file_name,'".', sep=''))
  # }
  #

  # Extract the model for analysis ####


  # titles
  MS1 <-  paste0('TITLE: ', 'Analyze data in an MS model, N = ', N, ', R = ', R, ';')
  MS1.2 <- paste0('TITLE: ', 'Analyze data in an MS model, N = ', N, ', R ', rep, ' of ', R, ';')
  usedata <- paste0('DATA: FILE= ', file_dir, "/",
                    unlist(strsplit(naming_data_files, "\\*"))[1], 'list', '.dat;')

  MS2 <- 'TYPE = MONTECARLO;'
  MS3 <- 'VARIABLE:'
  MS4 <- paste0('USEVARIABLES ARE ', Y[1], '-', Y[nYs], ';')
  MS6 <- 'ANALYSIS:'
  mplus_estimators <- c('ML', 'MLR', 'MLF', 'WLS','WLSM', 'WLSMV', 'ULS', 'ULSMV')
  chosen_estimators <- mplus_estimators
  parameterizations <- c('DELTA', 'THETA')
  link_functions <- c('LOGIT', 'PROBIT')

  # iterate estimator, parameterization, and link functions
  # Note: only parameterization is applied to WLSMV and ULSMV.
  est_para_linkf <- NULL
  estimator_lines <- NULL
  for (i in 1:length(chosen_estimators)){
    for (j in 1: length(parameterizations)){
      est <- paste0('ESTIMATOR = ', chosen_estimators[i], ';') # ESTIMATOR

      if(chosen_estimators[i] == 'WLSMV' | chosen_estimators[i] == 'ULSMV' |
         chosen_estimators[i] == 'ULS' | chosen_estimators[i] == 'WLSM'
         | chosen_estimators[i] == 'WLS'){
        para <- paste0('PARAMETERIZATION = ', parameterizations[j], ';') # parameterization for WLSMS/ULSMV (limited info)
        est_para_linkf[i] <- paste(est, para)
        estimator_lines <- rbind (estimator_lines, est_para_linkf[i])

      } else {
        linkf <- paste0('LINK IS ', link_functions[j], ';')
        est_para_linkf[i] <- paste(est, linkf)
        estimator_lines <- rbind (estimator_lines, est_para_linkf[i])
      }
    }

  }

  MS7 <- 'MODEL:'

  # model specifications
  D <- NULL # dimensions and loadings
  F_label <- NULL
  model_spec<-NULL

  for (i in 1:length(model)){
    F_label[i] <- paste0('F', i)
    # writing F_j BY
    D[i] <- paste(F_label[i], 'BY ', sep=' ')
    # iterate all loading lines per dimension
    D[i] <- paste0(D[i], 'Y', min(unlist(model[i])), '-', 'Y', max(unlist(model[i])),
                   # add factor loadings to be estimated
                   '*','(',
                   'Load', min(unlist(model[i])), '-', 'Load', max(unlist(model[i])),
                   ')',
                   sep=' ')

    model_spec[i] <- paste0(D[i], ';')
  }

  # model specification
  # model_spec

  # Threholds

  #nTs
  #nYs
  tau <- NULL
  for (i in 1:nYs){
    tau[i] <- paste0('[', Y[i], '$', 1, '-', Y[i], '$', nTs,'*', ']',
                     '(', Y[i], 'T',1, '-', Y[i], 'T', nTs, ')', ';')
  }
  # tau

  MS8 <- '! Factor Mean'

  # Factor mean
  factor_means <-paste0('[', 'F', 1, '-', 'F', length(model), '*]',
                        '(', 'Fmean', 1, '-', 'Fmean', length(model), ');')

  # Factor variance
  factor_variances <-paste0('F', 1, '-', 'F', length(model),
                            '*(', 'Fvariance', 1, '-', 'Fvariance', length(model), ');')

  # Model constraint

  MS9 <- 'MODEL CONSTRAINT:'

  model_constraint_m <- NULL
  # model constraint
  for (i in 1: length(model)){
    model_constraint_m[i] <- paste0('Fmean', i, '= 0;')
  }

  model_constraint_v <- NULL
  # model constraint
  for (i in 1: length(model)){
    model_constraint_v[i] <- paste0('Fvariance', i, '= 1;')
  }

  # interfactor correlation
  # Factor correlations

  FCor <- NULL
  for (i in 1:nrow(Fcor.tb)){
    # use Fcor.tb to get labels of factors and match factor correlation
    FCor[i] <-paste0('F', Fcor.tb[i,1], ' WITH ', 'F', Fcor.tb[i,2],
                     ';', sep=' ')
  }


  MS10 <- 'OUTPUT: STDYX;! Standardized solution'

  # Saving output thetas syntax
  MS11 <- 'SAVEDATA: SAVE = FSCORES; ! Saving factor scores'
  # writing file name for factor scores
  MS12 <- paste0('FILE = ',
                 gsub("\\*",
                      paste0(rep, "_fscores"),
                      naming_data_files))
  usedata2 <- gsub("list", rep, usedata)
  # need to iterate for all estimators



  # Iterate scripts through estimators

  # Estimator along with link function/parameterization
  # ML_logit <- estimator_lines[1]
  # ML_probit <- estimator_lines[2]
  # MLR_logit <- estimator_lines[3]
  # MLR_probit <- estimator_lines[4]
  # MLF_logit <- estimator_lines[5]
  # MLF_probit <- estimator_lines[6]
  # WLS_delta <- estimator_lines[7]
  # WLS_theta <- estimator_lines[8]
  # WLSM_delta <- estimator_lines[9]
  # WLSM_theta <- estimator_lines[10]
  # WLSMV_delta <- estimator_lines[11]
  # WLSMV_theta <- estimator_lines[12]
  # ULS_delta <- estimator_lines[13]
  # ULS_theta <- estimator_lines[14]
  # ULSMV_delta <- estimator_lines[15]
  # ULSMV_theta <- estimator_lines[16]

  est_names <- c('ML_logit',
                 'ML_probit',
                 'MLR_logit',
                 'MLR_probit',
                 'MLF_logit',
                 'MLF_probit',
                 'WLS_delta',
                 'WLS_theta',
                 'WLSM_delta',
                 'WLSM_theta',
                 'WLSMV_delta',
                 'WLSMV_theta',
                 'ULS_delta',
                 'ULS_theta',
                 'ULSMV_delta',
                 'ULSMV_theta')

  # add monte carlo conditions
  if (isTRUE(type_montecarlo)){
    script_container <- NULL
    for (i in 1:length(estimator_lines)){
      script_container <- capture.output(
        cat(L0, # Author info
            MS1, # title
            usedata, # put data file for analysis
            MS2, # TYPE = MONTECARLO
            MS3, # VARIABLE:
            L3, # names are y variables
            MS4, # usevariables are y variables
            L10, # categorical are
            MS6, # ANALYSIS

            # select estimator, parameterization, and link function
            estimator_lines[i],

            MS7, # MODEL:
            model_spec,
            #L14, # thresholds title
            #tau, # threshold lines
            FCor, # interactor correlations
            MS8, # Factor mean title
            factor_means, # Factor mean
            L13,# factor variance title
            factor_variances, # factor variances
            MS9, # model constraint
            model_constraint_m, # constraining factor means
            model_constraint_v, # constraining factor variance
            MS10, # output, standardized solution
            sep='\n')
      )


      # assign estimator name to script
      assign(est_names[i], script_container)
      script_container <- NULL
    }




    # creating folders
    folder <- NULL
    for (j in 1:length(estimators)){
      # create a new directory to store script of different estimator
      if(isFALSE(dir.exists(paste0(file_dir, "/", estimators[j])))){ # not exist, create
        dir.create(paste0(file_dir, "/", estimators[j]))
        folder[j] <- paste0(file_dir, "/", estimators[j])
      } else { # if exist, remove, then create
        folder[j] <- paste0(file_dir, "/", estimators[j])
      }

      # writing an mplus script data analysis to file
      # add condition from ther user: if the user's estimator, then write
      # otherwise, don't write

      for (i in 1:length(est_names)){
        if (est_names[i] %in% estimators[j]){ # matching user's estimator
          filename_obj <- paste0(folder[j], "/", est_names[i],"_replist", '.inp')
          message(paste0("creating:", filename_obj))
          # write file
          write.table(get(est_names[i]),# get the assigned object
                      file = filename_obj,
                      row.names = FALSE,
                      col.names = FALSE,
                      quote=FALSE) # use quote = FALSE to eliminate quotes in text
          if (isTRUE(run_files)){
            message(paste0("running: ", filename_obj))
            MplusAutomation::runModels(filename_obj)
            message(paste0("done: ", filename_obj))

          }
        }
      }
    }


    mplus_script <- list('ML_logit'	= ML_logit,
                         'ML_probit'	= ML_probit,
                         'MLR_logit'	= MLR_logit,
                         'MLR_probit'	= MLR_probit,
                         'MLF_logit'	= MLF_logit,
                         'MLF_probit'	= MLF_probit,
                         'WLS_delta'	= WLS_delta,
                         'WLS_theta'	= WLS_theta,
                         'WLSM_delta'	= WLSM_delta,
                         'WLSM_theta'	= WLSM_theta,
                         'WLSMV_delta'	= WLSMV_delta,
                         'WLSMV_theta'	= WLSMV_theta,
                         'ULS_delta'	= ULS_delta,
                         'ULS_theta'	= ULS_theta,
                         'ULSMV_delta'	= ULSMV_delta,
                         'ULSMV_theta'	= ULSMV_theta,
                         'script_list' = c('montecarlo_data', est_names))
  } # end type_montecarlo = TRUE











  if (isFALSE(type_montecarlo)){
    script_container <- NULL
    for (i in 1:length(estimator_lines)){
      script_container <- capture.output(
        cat(L0, # Author info
            MS1.2, # changed title from MS1 to MS1.1
            usedata2, # put data file for a specific rep for analysis
            #MS2, # TYPE = MONTECARLO
            MS3, # VARIABLE:
            L3, # names are y variables
            MS4, # usevariables are y variables
            L10, # categorical are
            MS6, # ANALYSIS

            # select estimator, parameterization, and link function
            estimator_lines[i],

            MS7, # MODEL:
            model_spec,
            #L14, # thresholds title
            #tau, # threshold lines
            FCor, # interactor correlations
            MS8, # Factor mean title
            factor_means, # Factor mean
            L13,# factor variance title
            factor_variances, # factor variances
            MS9, # model constraint
            model_constraint_m, # constraining factor means
            model_constraint_v, # constraining factor variance
            MS10, # output, standardized solution
            MS11, # SAVEDATA, FSCORES
            MS12, # naming the factor score file
            sep='\n')
      )


      # assign estimator name to script
      assign(est_names[i], script_container)
      script_container <- NULL
    }



    # writing an mplus script data analysis to file
    # adding user's estimators
    folder <- NULL
    for (j in 1:length(estimators)){
      # create a new directory to store script of different estimator
      if(isFALSE(dir.exists(paste0(file_dir, "/", estimators[j])))){ # not exist, create
        dir.create(paste0(file_dir, "/", estimators[j]))
        folder[j] <- paste0(file_dir, "/", estimators[j])
      } else { # if exist, remove, then create
        folder[j] <- paste0(file_dir, "/", estimators[j])
      }

      for (i in 1:length(est_names)){
        if (est_names[i] %in% estimators[j]){ # matching user's estimator
          filename_obj <- paste0(folder[j], "/", est_names[i],"_rep", rep, '.inp')
          message(paste0("creating:", filename_obj))
          # write file
          write.table(get(est_names[i]),# get the assigned object
                      file = filename_obj,
                      row.names = FALSE,
                      col.names = FALSE,
                      quote=FALSE) # use quote = FALSE to eliminate quotes in text

          if (isTRUE(run_files)){
            message(paste0("running: ", filename_obj))
            MplusAutomation::runModels(filename_obj)
            message(paste0("done: ", filename_obj))
          }
        }
      }
    }



    mplus_script <- list('ML_logit'	= ML_logit,
                         'ML_probit'	= ML_probit,
                         'MLR_logit'	= MLR_logit,
                         'MLR_probit'	= MLR_probit,
                         'MLF_logit'	= MLF_logit,
                         'MLF_probit'	= MLF_probit,
                         'WLS_delta'	= WLS_delta,
                         'WLS_theta'	= WLS_theta,
                         'WLSM_delta'	= WLSM_delta,
                         'WLSM_theta'	= WLSM_theta,
                         'WLSMV_delta'	= WLSMV_delta,
                         'WLSMV_theta'	= WLSMV_theta,
                         'ULS_delta'	= ULS_delta,
                         'ULS_theta'	= ULS_theta,
                         'ULSMV_delta'	= ULSMV_delta,
                         'ULSMV_theta'	= ULSMV_theta,
                         'script_list' = c('montecarlo_data', est_names))

  } # end type_montecarlo = FALSE

} # function