R/k.aefa3.R

# Kwangwoon Automated Exploratory Factor Analysis (K.AEFA)
# Seongho Bae (seongho@kw.ac.kr)

options(mc.cores = parallel::detectCores())

# for mac
if(Sys.getenv("LANG") == "ko-Kore_KR.UTF-8"){
  Sys.setenv(LANG="ko_KR.UTF-8")
  Sys.setlocale('LC_ALL', 'ko_KR.UTF-8')
  Sys.setlocale('LC_MESSAGES', 'ko_KR.UTF-8')
  Sys.setlocale('LC_CTYPE', 'ko_KR.UTF-8')
  Sys.setlocale('LC_COLLATE', 'ko_KR.UTF-8')
  Sys.setlocale('LC_TIME', 'ko_KR.UTF-8')
  Sys.setlocale('LC_MONETARY', 'ko_KR.UTF-8')
}


##############
# aefa frontend #
##############

message("\n Check Packages: Processing......")

# check packages

# try(mirtCluster(remove=TRUE), silent=TRUE)
# if(Sys.info()["sysname"] == "Linux"){
# try(mirtCluster(), silent=F)
# } else {
#   try(mirtCluster(as.numeric(parallel::detectCores() * 2)), silent=F)
# }


message(" Loading Packages: OK")

## pre-defined function for Full-information item factor analysis
findM2 <- function(mirtModel, increase = 15000, iterations = 1000, ...){

  # for treating unstable M2() function values
  quadpts <- 0 # initial quadpts
  for(i in 1:iterations){

    quadpts <- quadpts + increase

    message('quadpts: ', paste0(quadpts), ' / iteration: ', paste0(i))
    fitStats <- M2(mirtModel, QMC = TRUE, quadpts = quadpts, ...)

    if(i > 1){
      # decision-making
      if(abs((fitStats$RMSEA_5 - fitStats.old$RMSEA_5)) < .001 &&
         abs((fitStats$TLI - fitStats.old$TLI)) < .001 &&
         abs((fitStats$CFI - fitStats.old$CFI)) < .001)
        return(fitStats)
    }

    fitStats.old <- fitStats
  }

}

## fixTYPO for likert scaling
fixTYPO <- function(cleaningData){
  if((length(which(median(psych::describe(cleaningData)$rnage) != psych::describe(cleaningData)$range)) != 0) | length(which(median(psych::describe(cleaningData)$max) != psych::describe(cleaningData)$max)) != 0){
    for(ii in which(describe(cleaningData)$max > median(describe(cleaningData)$max))){
      cleaningData[,ii] <-plyr::mapvalues(cleaningData[,ii],psych::describe(cleaningData[,which(describe(cleaningData)$max > median(describe(cleaningData)$max))])$max, NA)
    }
    for(ii in which(describe(cleaningData)$min < median(describe(cleaningData)$min))){
      cleaningData[,ii] <-plyr::mapvalues(cleaningData[,ii],psych::describe(cleaningData[,which(describe(cleaningData)$min < median(describe(cleaningData)$min))])$min, NA)
    }
  }
  return(cleaningData)
}



# reverse coding
k.recode <- function(dname, target, scale) {

  i = 0
  x <- vector() ; y <- vector()
  x <- attributes(dname)$names
  y <- grep(target,x)

  for(i in 1:length(y)){

    dname[,y[i]] <- (scale + 1 - dname[,y[i]])

  }

  return(dname)

}

k.dichotomous <- function(dframe, start, end) {
  for(i in start:end) {
    dframe[,i] <- dframe[,i] >= median(dframe[,i],na.rm=T)
    dframe[,i] <-plyr::mapvalues(dframe[,i], c(TRUE, FALSE), c(1, 0))
    dframe[,i] <- as.integer(dframe[,i])
  }
  dframe <- data.frame(dframe)
  return(dframe)
}

# faking bad & aberrant response detection
fastHMM <- function(dat = ..., ...){
  temp_col <- colnames(dat)
  temp_col2 <- paste0(temp_col, "~1")
  resp <- list()
  model <- list()
  for(i in 1:length(temp_col2)){
    resp[[i]] <- as.formula(temp_col2[i])
    model[[i]] <- multinomial(link = 'identity')
  }

  for(i in 1:100){
    if(i == 1){
      try(mod <- depmix(response = resp, data = data.frame(dat), nstates = i, family = model), silent = T)
      try(fm <- fit(mod), silent = T)
    } else {
      fm_old <- fm
      rm(mod)
      rm(fm)

      try(mod <- depmix(response = resp, data = data.frame(dat), nstates = i, family = model), silent = T)
      try(fm <- fit(mod), silent = T)

      if(BIC(fm_old) < BIC(fm)){
        return(fm_old)
      }

    }
  }
}

k.faking <- function(data = ..., covdata = NULL, formula = NULL, SE = F, SE.type = "Oakes", skipNominal = F, forceGRSM = F, assumingFake = F, masterThesis = F, forceRasch = F, unstable = F, forceMHRM = F, printFactorStructureRealtime = F, itemkeys = NULL, survey.weights = NULL, IRTonly = F, ...) { # for aberrant & faking response detection
  dataset <- data
  dname <- data

  if(IRTonly == F){
    for(j in 1:1000){
      try(HMM_result <- fastHMM(dat = dname))
      if(exists('HMM_result')){

        # HMM model check
        if(HMM_result@homogeneous == TRUE){
          message('\nThis data is homogeneous data. Please be careful for check faking response.')
        }

        stateDat <- data.frame(count(HMM_result@posterior, 'state'))
        judgementHMMnormal <- vector(length = length(HMM_result@posterior$state))
        judgementHMMnormal[which(HMM_result@posterior$state == stateDat$state[which(max(stateDat$freq) == stateDat$freq)])] <- TRUE


        # person-fit test in IRT
        if(sum(is.na(dataset)) == 0){
          dataset.mirt <- fastFIFA(x = as.data.frame(data), covdata = as.data.frame(covdata), formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, itemkeys = itemkeys, survey.weights = survey.weights, ...)
          dataset.response <- mirt::personfit(dataset.mirt, method='MAP', QMC = T)


        } else {
          dataset.mirt <- fastFIFA(x = as.data.frame(data), covdata = as.data.frame(covdata), formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, itemkeys = itemkeys, survey.weights = survey.weights, ...)
          dataset_temp <- imputeMissing(x = dataset.mirt, Theta = fscores(dataset.mirt, method = 'MAP', QMC = T), QMC = T, impute = 100)
          dataset.mirt2 <- fastFIFA(x = as.data.frame(dataset_temp), covdata = as.data.frame(covdata), formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, itemkeys = itemkeys, survey.weights = survey.weights, ...)
          dataset.response <- mirt::personfit(dataset.mirt2, method='MAP', QMC = T)
        }


        print(hist(dataset.response$Zh))
        dataset.response$Zh <- dataset.response$Zh > -2.58 #(if abnormal, dataset.response$Zh < -2 is right! : See Hyeongjun Kim (2015) @ SNU)
        IRTnormal <- data.frame(dataset.response$Zh)

        # reflect results
        if(HMM_result@homogeneous == TRUE){
          message('\nThis results may not refer detection of fake response; just a aberrant response. Like a contents non-relavant random response.')
          #output <- data.frame(dataset.response$Zh)

        }

        if(sum(is.na(dataset)) == 0){
          output <- data.frame(rowSums(data.frame(judgementHMMnormal, IRTnormal)) == 2)
        } else {
          output <- data.frame(judgementHMMnormal)
        }




        names(output) <- "normal"
        row.names(output) <- row.names(dataset)
        result <- data.frame(dataset, output)
        return(result)
      } else {

      }
    }
  } else {
    for(j in 1:1000){
      # person-fit test in IRT
      if(sum(is.na(dataset)) == 0){
        dataset.mirt <- fastFIFA(x = as.data.frame(data), covdata = as.data.frame(covdata), formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, itemkeys = itemkeys, survey.weights = survey.weights, ...)
        dataset.response <- mirt::personfit(dataset.mirt, method='MAP', QMC = T)


      } else {
        message('imputing missing values')
        dataset.mirt <- fastFIFA(x = as.data.frame(data), covdata = as.data.frame(covdata), formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, itemkeys = itemkeys, survey.weights = survey.weights, ...)
        dataset_temp <- imputeMissing(x = dataset.mirt, Theta = fscores(dataset.mirt, method = 'MAP', QMC = T), QMC = T, impute = 100)
        message('re-estimating parameters')
        dataset.mirt2 <- fastFIFA(x = as.data.frame(dataset_temp), covdata = as.data.frame(covdata), formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, itemkeys = itemkeys, survey.weights = survey.weights, ...)
        dataset.response <- mirt::personfit(dataset.mirt2, method='MAP', QMC = T)
      }

      print(hist(dataset.response$Zh))
      dataset.response$Zh <- dataset.response$Zh > -2 #(if abnormal, dataset.response$Zh < -2 is right! : See Hyeongjun Kim (2015) @ SNU)
      IRTnormal <- data.frame(dataset.response$Zh)

      # if(sum(is.na(dataset)) == 0){
      output <- data.frame(IRTnormal)
      # } else {
      # stop('Please use HMM')
      # }


      names(output) <- "normal"
      row.names(output) <- row.names(dataset)
      result <- data.frame(dataset, output)
      return(result)
    }

  }
}


aberrantZero <- function(data = ..., covdata = ..., formula = ..., ...){

  tempData <- data
  colnames(tempData) <- paste0("testVars", 1:ncol(tempData))
  for(i in 1:100000){

    if(i == 1){
      mod <- k.faking(data = tempData, ...) # TRUE search
    } else {

      mod_old <- mod
      rm(mod)

      message('current number of samples: ', (nrow(mod_old[mod_old$normal == T,1:ncol(mod_old)-1]))) # count TRUE samples
      print(apply(mod_old[mod_old$normal == T,1:ncol(mod_old)-1], 2, table)) # count TRUE sample patterns
      mod <- k.faking(data = mod_old[mod_old$normal == T,1:ncol(mod_old)-1], ...) # recalculate using TRUE SAMPLE ONLY

      if(nrow(mod) == nrow(mod_old)) { # if no difference between mod and mod_old
        # message('final normal cases: ', paste0(rownames(mod_old)))
        return(rownames(mod_old))
      }
    }

  }
}

# CTT multilevel
k.multilevel <- function(variable = ..., group = ..., scale = ..., operation = ..., x = ..., y = ...) {

  ranvar_calculate <- (scale^2-1)/12
  ranvarT_calculate <- (scale^2+2*scale-2)/24

  #########################
  # rwg family            #
  #########################

  #rwg (single measurement variable)
  if(operation == 'rwg.single') {
    RWG.result<-rwg(variable,group,ranvar=ranvar_calculate) # var, group, rectangular function - (A^2-1)/12
    RWG.T.result<-rwg(variable,group,ranvar=ranvarT_calculate) # var, group, rectangular function - (A^2+2A-2)/24


    #print(summary(RWG.RELIG))
    message('James, Demaree & Wolf (1984) rwg for single item measures\n')

    message('sorted rwg(U)')
    print(sort(RWG.result[,2],decreasing=F))

    message('\nsorted rwg(T)')
    print(sort(RWG.T.result[,2],decreasing=F))

    message('\nmean rwg(U)')
    print(mean(RWG.result[,2]))

    message('\nmean rwg(T)')
    print(mean(RWG.T.result[,2]))

    message('\nmedian rwg(U)')
    print(median(RWG.result[,2]))

    message('\nmedian rwg(T)')
    print(median(RWG.T.result[,2]))

    message('\nmean group size')
    print(mean(RWG.result[,3]))
    #max(sort(RWG.RELIG[,2]))
    #mean(sort(RWG.RELIG[,2]))
    #min(sort(RWG.RELIG[,2]))
    hist(sort(RWG.result[,2]), xlab='rwg', main='Histogram of rwg')

    message('\nANOVA Table')
    ICC.calculate <- aov(variable~as.factor(group)); print(summary(ICC.calculate))

    message('\nEta-squared')
    print(etaSquared(ICC.calculate))

    message('\nICC(1)')
    print(ICC1(ICC.calculate))

    message('\nICC(2)')
    print(ICC2(ICC.calculate))

    return(RWG.result)
  }
  else if(operation == 'rwg.multiple') {
    RWG.result<-rwg.j(variable,group,ranvar=ranvar_calculate) # var, group, rectangular function - (A^2-1)/12
    RWG.T.result<-rwg.j(variable,group,ranvar=ranvarT_calculate) # var, group, rectangular function - (A^2+2A-2)/24


    #print(summary(RWG.RELIG))
    message('James, Demaree & Wolf (1984) rwg for multi-item measures\n')
    message('sorted rwg(U)')
    print(sort(RWG.result[,2],decreasing=T))

    message('\nsorted rwg(T)')
    print(sort(RWG.T.result[,2],decreasing=T))

    message('\nmean rwg(U)')
    print(mean(RWG.result[,2]))

    message('\nmean rwg(T)')
    print(mean(RWG.T.result[,2]))

    message('\nmedian rwg(U)')
    print(median(RWG.result[,2]))

    message('\nmedian rwg(T)')
    print(median(RWG.T.result[,2]))

    hist(sort(RWG.result[,2]), xlab='rwg', main='Histogram of rwg(U)')


    message('\nANOVA Table (composite value)')
    variable2 <- rowMeans(variable)
    ICC.calculate <- aov(variable2~as.factor(group)); print(summary(ICC.calculate))

    message('\nEta-squared')
    print(etaSquared(ICC.calculate))

    message('\nICC(1)')
    print(ICC1(ICC.calculate))

    message('\nICC(2)')
    print(ICC2(ICC.calculate))

    message('\nMulti ICC(1) & ICC(2)')
    rwg.j_result <- mult.icc(x=variable, grpid=group)
    print(rwg.j_result)

    message('\nMean group size')
    print(mean(RWG.result[,3]))
    message('\nStandard deviation of group size')
    print(sd(RWG.result[,3]))

    return(RWG.result)
  }
  else if(operation == 'rwg.lindell') {
    RWG.result<-rwg.j.lindell(variable,group,ranvar=ranvar_calculate) # var, group, rectangular function - (A^2-1)/12
    RWG.T.result<-rwg.j.lindell(variable,group,ranvar=ranvarT_calculate) # var, group, rectangular function - (A^2-1)/12

    #print(summary(RWG.RELIG))
    message('Lindell, & Brandt(1997; 1999) r* wg(j) for multi-item measures\n')
    message('sorted rwg(U)')
    print(sort(RWG.result[,2],decreasing=F))

    message('\nsorted rwg(T)')
    print(sort(RWG.T.result[,2],decreasing=F))

    message('\nmean rwg(U)')
    print(mean(RWG.result[,2]))

    message('\nmean rwg(T)')
    print(mean(RWG.T.result[,2]))

    message('\nmedian rwg(U)')
    print(median(RWG.result[,2]))

    message('\nmedian rwg(T)')
    print(median(RWG.T.result[,2]))
    hist(sort(RWG.result[,2]), xlab='rwg', main='Histogram of rwg')

    message('\nANOVA Table (composite value)')
    variable2 <- rowMeans(variable)
    ICC.calculate <- aov(variable2~as.factor(group)); print(summary(ICC.calculate))

    message('\nEta-squared')
    print(etaSquared(ICC.calculate))

    message('\nICC(1)')
    print(ICC1(ICC.calculate))

    message('\nICC(2)')
    print(ICC2(ICC.calculate))

    message('\nMulti ICC(1) & ICC(2)')
    print(mult.icc(x=variable, grpid=group))

    return(RWG.result)
  }

  #########################
  # awg family            #
  #########################

  else if(operation == 'awg'){
    AWG.result<-awg(variable,group) # var, group

    #print(summary(RWG.RELIG))
    message('Brown and Hauenstein (2005) awg for multi-item measures\n')
    message('sorted awg(s)')
    print(sort(AWG.result[,2],decreasing=F))
    message('\nmean awg')
    print(mean(AWG.result[,2]))
    message('\nmean group size')
    print(mean(AWG.result[,4]))
    #max(sort(RWG.RELIG[,2]))
    #mean(sort(RWG.RELIG[,2]))
    #min(sort(RWG.RELIG[,2]))
    hist(sort(AWG.result[,2]), xlab='awg', main='Histogram of awg')
    return(AWG.result)
  }

  else if(operation == 'ad'){
    AD.result<-ad.m(variable,group) # var, group

    #print(summary(RWG.RELIG))
    message('Burke, Finkelstein and Dusig (1999) Average Deviation (AD) Agreement for multi-item measures\n')
    message('sorted ad')
    print(sort(AD.result[,2],decreasing=F))
    message('\nmean ad')
    print(mean(AD.result[,2]))
    message('\nmean group size')
    print(mean(AD.result[,3]))
    #max(sort(RWG.RELIG[,2]))
    #mean(sort(RWG.RELIG[,2]))
    #min(sort(RWG.RELIG[,2]))
    hist(sort(AD.result[,2]), xlab='ad', main='Histogram of ad')
    return(AD.result)
  }

  else if(operation == 'waba'){
    WABA.result<-waba(x, y, group) # var, group

    message('Dansereau, Alutto & Yammarino (1984) WABA II\n')
    message('Covariance Theorem')
    print(WABA.result$Cov.Theorem)
    message('\n n size of observation')
    print(WABA.result$n.obs)
    message('\n number of groups')
    print(WABA.result$n.grps)

    return(WABA.result)
  }


}

## meta-analysis Hunter & Schmidt Estimator
k.meta <- function(dframe = ..., # data frame
                   calc = ..., # Calculation mode
                   study = ..., n = ..., Rxy = ..., Rxx = ..., Ryy = ..., u = ..., moderator = ..., # correlation based meta-analysis
                   measure = ..., treatment_positive = ..., treatment_negative = ..., controlled_positive = ..., controlled_negative = ..., # treatment effect
                   reported_r = ..., selected_SD = ..., whole_SD = ..., # correcting range_restriction
                   ...) {

  if(calc == "treatment"){
    # example: meta.test <- k.meta(dframe=dat.bcg, calc="treatment", measure="RR", treatment_positive=tpos, treatment_negative=tneg, controlled_positive=cpos, controlled_negative=cneg)

    attach(dframe)
    # for 2x2 experiment
    message('Effect size meta analysis for treatment')
    message(' \n ')
    message('If you need help, please ?escalc and ?rma\n')

    dat <- escalc(measure=measure, ai=treatment_positive, bi=treatment_negative, ci=controlled_positive, di=controlled_negative, data=dframe)

    #detach(dframe)
    #detach(dframe)
    calculate <- rma(yi = dat$yi, vi = dat$vi, method="HS")

    print(calculate)
    message("r coefficient")
    print(z2r(calculate$zval))

    detach(dframe)
    return(calculate)

  }
  else if(calc == "d"){

  }
  else if(calc == "r" | calc == "correlation") {
    # example: k.meta(dframe=ABHt32, calc="r", study=ABHt32$study, n=ABHt32$n, Rxy=ABHt32$Rxy, Rxx=ABHt32$Rxx, Ryy=ABHt32$Ryy, u=NA, moderator=NA)

    data <- data.frame(study, Rxy, n, Rxx, Ryy, u, moderator)
    colnames(data) <- c("study", "Rxy", "n", "Rxx", "Ryy", "u", "moderator")
    head(data)

    N <- sum(data$n)
    k <- nrow(data)

    calculate <- MetaTable(data)

    result <- data.frame(N, k, calculate)

    print(result[1:5])
    print(result[6:10])
    print(result[11:15])
    FileDrawer(data)
    FunnelPlot(data)
    return(result)

  }
  else if(calc == "range_restriction") {
    # example: k.meta(calc='range_restriction', reported_r=0.3, selected_SD=6, whole_SD=10)
    # You can use this result where
    # See T. Yoo. & D. Kim. (2003). A Meta-Analysis for Validity Generalization: Integrating Correlation Coefficients. Digital Business Studies, 9, 61-80.

    u <- selected_SD / whole_SD
    c <- sqrt((1-u^2) * reported_r^2 + u^2)

    message("u")
    print(u)
    message("range restrict corrected rho")
    print(c)
    return(c)
  }
  else {

  }

}



## K-means++
k.kmpp <- function(X, k) {
  n <- nrow(X)
  C <- numeric(k)
  C[1] <- sample(1:n, 1)

  for (i in 2:k) {
    dm <- distmat(X, X[C, ])
    pr <- apply(dm, 1, min); pr[C] <- 0
    C[i] <- sample(1:n, 1, prob = pr)
  }

  kmeans(X, X[C, ])
}

## Polynominal Regression
k.polynomial <- function(dependent = ..., independent = ..., moderator = ..., independent_name = ..., moderator_name = ..., dependent_name = ..., rotation_axis = ..., view_axis = ...) {
  #if(centering == TRUE){
  #  Z <- dependent
  #  X <- scale(independent, center=T, scale=T)
  #  M <- scale(moderator, center=T, scale=T)
  #  X <- X[,1]
  #  M <- M[,1]
  #} else {
  Z <- dependent
  X <- independent
  M <- moderator
  #}

  #poly.regression <- lm(Z ~ X + M + I(X^2) + I(M^2) + X * M)
  poly.regression <- lm(Z ~ poly(X, M, degree=2))
  message('Polynomial Regression')
  message('
          Z = intercept + X + X^2 + M + X*M + M^2')
  message('
          ---------------------------------------
          in Terms    in Coefficients
          ---------------------------------------
          intercept:  (Intercept)
          X:          poly(X, M, degree = 2)1.0
          X^2:        poly(X, M, degree = 2)2.0
          M:          poly(X, M, degree = 2)0.1
          X*M:        poly(X, M, degree = 2)1.1
          M^2:        poly(X, M, degree = 2)0.2
          ---------------------------------------')
  print(summary(poly.regression))
  #message('VIF')
  #poly.vif <- vif(poly.regression)
  #print(poly.vif)

  #if(sum(poly.vif>10) >= 1) {
  #  par(mfrow=c(2,2))
  #  plot(poly.regression)

  #  message('\nCorrelation chart')
  #  print(cor(poly.regression$model))
  #  message(' ')

  #  stop("STOP: VIF were too high. Please reconsider about this model.")
  #} else {
  par(mfrow=c(1,1))
  try(persp(poly.regression, X ~ M, col = rainbow(50), zlab='Z', contours = "col", theta = rotation_axis, phi = view_axis)) # xlab = independent_name, ylab = moderator_name, zlab = dependent_name,
  return(poly.regression)
  #}

}

#par(mfrow=c(2,2))

k.rescale <- function(original, scalerange) {
  rescaled <- as.data.frame( lapply(original, cut, scalerange, labels=FALSE) )
  return(rescaled)
}

k.rescale.na <- function(original, scalerange, na_to_zero = ..., ...) {
  if(na_to_zero == TRUE) {
    original[is.na(original <- data.frame(original))] <- 0
    rescaled <- as.data.frame( lapply(original, cut, scalerange, labels=FALSE) )
    return(rescaled)
  } else {
    rescaled <- as.data.frame( lapply(original, cut, scalerange, labels=FALSE) )
    return(rescaled)
  }
}

splitdf <- function(dataframe, seed=NULL) {
  if (!is.null(seed)) set.seed(seed)
  index <- 1:nrow(dataframe)
  trainindex <- sample(index, trunc(length(index)/2))
  trainset <- dataframe[trainindex, ]
  testset <- dataframe[-trainindex, ]
  list(trainset=trainset,testset=testset)
}

k.split <- function(dataframe){
  split <- splitdf(dataframe, seed=12345)
  str(split)
  lapply(split,nrow)
  lapply(split,head)

  training <- split$trainset
  testing <- split$testset

  training$sample <- "developmental"
  testing$sample <- "validation"

  split <- rbind(training, testing)
  split$sample <- as.factor(split$sample)
  str(split$sample)

  return(split)
}



#################################
# Analytical tools for          #
# Likert style (ordinal) scale  #
#################################

# Multiple Imputation
# (require m > 100)
likertMI <- function(model = ..., data = ..., m = 100, fun = 'sem', estimator = 'MML', parameterization = ..., chi = 'lmrr', ...) {

  #########################
  # Multiple              #
  # Imputation for        #
  # likert style measured #
  # variables             #
  # in SEM context        #
  #########################
  # Seongho Bae           #
  # seongho@kw.ac.kr      #
  # Nov 6th 2014          #
  # Under GNU / GPL       #
  #########################

  # checking packages for multiple impuation in SEM context
  if(!require('semTools')) {
    try(install.packages('semTools', dependencies = TRUE), silent=TRUE)
  }

  if(!require('Amelia')) {
    try(install.packages('Amelia', dependencies = TRUE), silent=TRUE)
  }

  if(!require('lavaan')) {
    try(install.packages('lavaan', dependencies = TRUE), silent=TRUE)
  }

  # loading packages for multiple imputation in SEM context
  library('semTools')
  library('Amelia')
  library('lavaan')

  fit <- sem(model=model, data=data.frame(data)) # extract variable names in model syntax
  message("sample size (listwise): ", paste0(nrow(data.frame(fit@Data@X))))

  fit_MI <- runMI(model=model, data=data.frame(data[,attributes(fit)$Model@dimNames[[1]][[1]]]), m=m, fun = fun, ordered=names(data[,attributes(fit)$Model@dimNames[[1]][[1]]]), miArgs=list(ords = attributes(fit)$Model@dimNames[[1]][[1]]), estimator = estimator, chi = chi, ...)#, control=list(optim.method="L-BFGS-B"), ...)

  cat(summary(fit_MI, standardize=T))
  print(inspect(fit_MI, 'fit'))
  Sys.sleep(60)

  fit_data <- data.frame(fit_MI@Data@X)
  colnames(fit_data) <- attributes(fit_MI)$Model@dimNames[[1]][[1]]
  message("sample size (imputated): ", paste0(nrow(data.frame(fit_MI@Data@X))))

  if(fun == 'cfa' | fun == 'CFA') {
    fit_MI_theta <- cfa(model=model, data=fit_data, ordered=names(fit_data), estimator = estimator, ...) #, control=list(optim.method="L-BFGS-B")
  } else if(fun == 'sem' | fun == 'SEM') {
    fit_MI_theta <- sem(model=model, data=fit_data, ordered=names(fit_data), estimator = estimator, ...) #, control=list(optim.method="L-BFGS-B")
  } else if(fun == 'growth' | fun == 'GROWTH') {
    fit_MI_theta <- growth(model=model, data=fit_data, ordered=names(fit_data), estimator = estimator, ...) #, control=list(optim.method="L-BFGS-B")
  }

  #   return(fit_MI)
  return(fit_MI_theta)

  # fit_cfaMI <- likertMI(model=model.sem, data=rawdata, m=100, estimator='WLSMV', fun='cfa', verbose=T)
  # fit_semMI <- likertMI(model=model.sem, data=rawdata, m=100, estimator='WLSMV', fun='sem', verbose=T)
}


# Full-automatically information item factor analysis until investigating simple factor structure
likertFA <- function(data = ..., ...) {

  fa_covdata <- data.frame(data)
  # fa_covdata <- k.imputation(fa_covdata, ...) # data imputation

  message('\nStage 1 of calcuation: Discriminant coefficient based Evaluation')
  ncol_stage1 <- ncol(fa_covdata)
  message('\nCurrent number of Items: ', paste0(ncol_stage1))
  result <- k.aefa(fa_covdata, estimator='mirt', ...)

  #   test <- data.frame(coef(result))
  #   b <- abs(test[1, grep("a[0-9]$", colnames(test))]) >= .5 # F. B. Baker (2001; p. 159; .5 <= a < 2)
  #
  #   test <- test[1, colnames(b)] # replace
  #
  #   c <- vector() # saving names
  #   d <- 0 # name counter
  #   require(stringr)
  #
  #   for(i in 1:length(colnames(b))) {
  #     if((abs(test[i]) >= .5) == TRUE){
  #       d <- d+1
  #       c[d] <- str_sub(colnames(test[i]), end = -4)
  #     } else {
  #
  #     }
  #   }

  find_a_lower <- which(data.frame(MDISC(result)) >= .5)
  #   find_a_upper <- which(data.frame(MDISC(result)) < 2.00) # F. B. Baker (2001; p. 159; .5 <= a < 2)
  #   find_a <- c(find_a_lower)#, find_a_upper)

  fa_covdata_temp <- data.frame(fa_covdata[,find_a_lower])

  #   #test <- vector()
  #   test <- try(mod2values(result), silent = T)
  #   if(exists("test")){
  #     test <- data.frame(test[grep("a[0-9]+", test$name),])
  #     test <- subset(test, test$value >= .4)
  #   } else {
  #     test <- data.frame()
  #   }
  #
  #
  #   if(nrow(test)==0){
  #     message('Passing........')
  #   } else {
  #     test$item <- as.character(test$item)
  #     test$item <- as.factor(test$item)
  #     include <- as.character(test$item)
  #     myvars <- names(fa_covdata) %in% c(include)
  #     # fa_covdata_temp <- fa_covdata[myvars]
  #     fa_covdata_temp <- subset(fa_covdata, select=myvars)
  fa_covdata <- fa_covdata_temp
  #   }

  # evaluation of IRT itemfit
  message('\nStage 2 of calcuation: Item Fit Based Evaluation')


  ncol_stage2 <- ncol(fa_covdata)
  if(ncol_stage2 == 0){
    stop('All items are deleted')
  } else {
    message('\nCurrent number of Items: ', paste0(ncol_stage2))
  }


  # screening items (stage 1)
  if(ncol_stage1 == ncol_stage2){
    # evaluating model
    #result <- k.aefa(fa_covdata, estimator='mirt', ...)
    #print(summary(result))
  } else {
    # evaluating model
    result <- k.aefa(fa_covdata, estimator='mirt', ...)
    #print(summary(result))
  }

  message('estimating Theta')
  try(test2_fscores <- fscores(object = result, rotate = 'geominQ', full.scores = T, plausible.draws = 100, QMC = TRUE, method = 'MAP', MI = 100))
  message('estimating itemfit')
  try(test2 <- mirt::itemfit(result, method = 'MAP', impute = 100, QMC = TRUE), silent = T)

  if(exists("test2")) { # patch for mirt function bug...
    test2$cal <- test2$S_X2/test2$df.S_X2
    test2 <- subset(test2, test2$cal >= 3)
  } else {
    test2 <- data.frame()
    #test2 <- 0
  }

  if(nrow(test2)==0){
    message('Passing........')
  } else {
    exclude <- as.character(test2$item)
    exclude <- as.factor(exclude)
    exclude <- as.character(exclude)

    myvars <- names(fa_covdata) %in% c(exclude)
    fa_covdata <- fa_covdata[!myvars]
  }


  message('\nStage 3 of calcuation: Factor Loading Based Evaluation')

  for(i in 1:10000){

    ncol_stage3 <- ncol(fa_covdata)
    if(ncol_stage3 == 0){
      stop('All items are deleted')
    } else {
      message('\nCurrent number of Items: ', paste0(ncol_stage3))
    }

    if(ncol_stage1 == ncol_stage3 | ncol_stage2 == ncol_stage3){
      # evaluating model
      #result <- k.aefa(fa_covdata, estimator='mirt', ...)
      #print(summary(result))
    } else {
      # evaluating model
      result <- k.aefa(fa_covdata, estimator='mirt', ...)
      #print(summary(result))
    }


    if(ncol(result@Fit$F)==1){
      rotF_geomin <- data.frame(result@Fit$F)
      h2 <- data.frame(result@Fit$h2)
    } else {
      # getting factor loadings
      message('Rotating Factor Solution now')
      rotF_geomin <- GPArotation::geominQ(result@Fit$F, maxit = 100000)
      rotF_geomin <- data.frame(rotF_geomin$loadings)
      h2 <- data.frame(result@Fit$h2)
    }


    if(i > 1){
      exclude_rownum1 <- NULL
      exclude_rownum1_low <- NULL
      exclude_rownum2 <- NULL
      exclude_rownum2_low <- NULL
      exclude_rownum3 <- NULL
      exclude_rownum3_low <- NULL
      exclude <- NULL
      exclude_rownum_low <- NULL
      myvars <- NULL

      if(length(exclude_rownum1) != 0){
        try(rm(exclude_rownum1, exclude_rownum2, exclude_rownum2_low, exclude_rownum3, exclude_rownum3_low, exclude, exclude_rownum_low, myvars), silent = T)

      } else {
        try(rm(exclude_rownum1, exclude_rownum1_low, exclude_rownum2, exclude_rownum2_low, exclude_rownum3, exclude_rownum3_low, exclude, exclude_rownum_low, myvars), silent = T)

      }
    }

    # evaluation of factor structure (stage 2) -- evaluating cross loadings and very small loadngs
    exclude_rownum1 <- which(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) > 0.4) >=2) # cross loadings based delection
    if(length(exclude_rownum1) != 0){
      exclude_rownum1_low <- as.numeric(names(which(min(rowSums(abs(rotF_geomin[rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) > 0.4) >=2,]))) == rowSums(abs(rotF_geomin[rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) > 0.4) >=2,])))))
    }
    exclude_rownum2 <- which(h2 < .5) # communality based delection
    exclude_rownum2_low <- which(h2 == min(h2))
    exclude_rownum3 <- which(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) < 0.4) == ncol(rotF_geomin)) # fail to load any factors
    exclude_rownum3_low <- which(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)])) == min(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]))))

    #     exclude_rownum <- as.numeric(paste(c(exclude_rownum1, exclude_rownum2, exclude_rownum3))) # list of doing delection

    if(length(exclude_rownum3) > 0){
      message('[WARN] No loadings to any factor(s) occured!')
      exclude_rownum_low <- as.numeric(paste(c(exclude_rownum3_low))) # list of doing delection
      #       message(paste0(exclude_rownum_low))
    } else if(length(exclude_rownum1) > 0){
      message('[WARN] Cross Loadings occured!')
      exclude_rownum_low <- as.numeric(paste(c(exclude_rownum1_low))) # list of doing delection
    } else if(length(exclude_rownum2) > 0){
      message('[WARN] Communality problem occured!')
      exclude_rownum_low <- as.numeric(paste(c(exclude_rownum2_low))) # list of doing delection
    } else {
      message('[Done!]')
      exclude_rownum_low <- NULL
      #       exclude_rownum_low <- as.numeric(paste(c(exclude_rownum1_low, exclude_rownum2_low, exclude_rownum3_low))) # list of doing delection
    }

    #     exclude_rownum <- as.numeric(paste(c(exclude_rownum1, exclude_rownum3))) # list of doing delection

    # the start of variable delection
    if(length(exclude_rownum_low) > 0) { # if number of delection is non-zero
      exclude <- vector() # make vectors
      #       j <- 0 # set to zero exclude list counter

      #       for(i in c(exclude_rownum)){ # saving list of delection variable
      #       j <- j + 1
      #print(j)
      #message('Trying to extract ', paste0(i), ' factor solution') -- for debug code
      #         print(names(fa_covdata)[i])
      exclude <- names(data.frame(result@Data$data))[c(exclude_rownum_low)]
      message(paste0(exclude))
      #         exclude[j] <- names(fa_covdata)[c(exclude_rownum_low)]
      #       }

      myvars <- names(data.frame(result@Data$data)) %in% c(exclude)
      fa_covdata <- data.frame(result@Data$data)[!myvars]

    } else { # (length(exclude_rownum == 0))
      try(print(findM2(result, calcNull = T, Theta = fscores(result, full.scores = T, plausible.draws = 100, rotate = "geominQ", MI = 100, QMC = TRUE, method = 'MAP'), impute = 100)), silent = T)
      #cat('number of iteration : paste0(a))
      return(result)
    }


  }
}

k.sampling <- function(data, n){
  set.seed(1234)
  data <- data[sample(nrow(data), n), ]
  return(data)
}

# k.select <- function(model, items){
#   select <- data.frame(model@Data$data[,rownames(subset(model@Fit$F, model@Fit$F >= min(sort(model@Fit$F, decreasing = T)[1:items])))])
#   return(select)
#
#
# }

k.select <- function(model, items){
  select <- model@Data$data[,c(subset(colnames(model@Data$data), model@Fit$h2 >= min(sort(model@Fit$h2, decreasing = T)[1:items])))]
  return(select)
}

# finding problemistic data
k.fixdata <- function(data, start, end, bioend){

  data <- data[!duplicated(data[start:bioend]),] # remove duplicate cases automatically

  dropVector <- vector()
  j <- 0
  for(i in 1:nrow(data)){
    if(std(as.numeric(data[i,start:end])) == 0){
      j <- j + 1
      dropVector[j] <- i
      print(dropVector)
    } else {

    }
  }
  return(data[c(-dropVector),])
}

# surveyFA addon
fastFIFA <- function(x, covdata = NULL, formula = NULL, SE = T, SE.type = "sandwich", skipNominal = F,
                     forceGRSM = F, assumingFake = F, masterThesis = F, forceRasch = F, unstable = F,
                     forceMHRM = F, forceNormalEM = T, itemkeys = NULL, survey.weights = NULL, allowMixedResponse = T,
                     forceUIRT = F, skipIdealPoint = F, MHRM_SE_draws = 1e+4, forceNRM = F,
                     diagnosis = F, forceDefalutAccelerater = F, forceDefaultOptimizer = F, EnableFMHRM = F, testlets = NULL, plotOn = T, GenRandomPars = T,
                     fixed = ~1, random = NULL, lr.fixed = ~1, lr.random = NULL, MH_BURNIN = 1500, MH_SEMCYCLES = 1000, ...){

  if(!require('mirt')){
    install.packages('mirt')
    library('mirt')
  }
  x <- x[,psych::describe(x)$range > 0] # delete no variance items

  for(i in 1:ncol(x)){
    try(invisible(gc()), silent = T) # garbage cleaning

    # testlets
    ActualTestlets <- testlets
    if(length(testlets) != 0){
      if(!require('plyr')){
        install.packages('plyr')
        library('plyr')
      }
      TestletActivated <- T

      if(is.character(ActualTestlets)){
        ActualTestlets <- as.factor(ActualTestlets)
      }
      ActualTestlets <- as.integer(ActualTestlets)

      ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))
      # if(sum(ActualTestlets %in% 1) == 0){
      #   ActualTestlets <- ActualTestlets - min(ActualTestlets, na.rm = T) + 1
      # }
      if(sum(is.na(ActualTestlets)) == length(ActualTestlets)){
        ActualTestlets <- NULL
        TestletActivated <- F

      } else {
        ActualTestlets <- plyr::mapvalues(ActualTestlets, as.numeric(attributes(as.factor(ActualTestlets))$levels), seq(length(attributes(as.factor(c(ActualTestlets)))$levels)))

      }

    } else {
      TestletActivated <- F
    }


    if(TestletActivated){
      message('\nfactor numbers: ', paste0(i+max(na.omit(ActualTestlets))))

    } else {
      if (i == 1){
        message('\nfactor number: ', paste0(i))
      } else {
        message('\nfactor numbers: ', paste0(i))
      }
    }



    # if(sum(is.na(x)) == 0 | nrow(x) > 5000){
    NofCores <- parallel::detectCores()
    NofCores <- round(NofCores / 1.1)
    if(NofCores > 12){
      NofCores <- 12
    }
    try(invisible(mirt::mirtCluster(spec = NofCores)), silent = T)
    # }

    # optimizer config
    if(length(covdata) == 0){ # if no covariate variables
      if(TestletActivated){
        # print(ActualTestlets)
        if(sum(na.omit(ActualTestlets) > 2) != 0 | forceMHRM){
          estimationMETHOD <- 'MHRM'
          optimINPUT <- NULL
          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- 4000
        } else {
          estimationMETHOD <- 'EM'
          optimINPUT <- NULL
          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- 1000
        }

      } else if(forceMHRM == T | forceGRSM == T | assumingFake == T | masterThesis == T) {
        estimationMETHOD <- 'MHRM'
        optimINPUT <- NULL
        optimCTRL  <- NULL
        empiricalhist <- FALSE
        NCYCLES <- 4000
      } else if(length(survey.weights) != 0) {
        if (forceNormalEM == T && i < 4){
          estimationMETHOD <- 'EM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- NULL
        } else {
          estimationMETHOD <- 'QMCEM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- NULL
        }
      } else if(i < 4){ # EM will activate 1-3 factors
        if (forceNormalEM == T && unstable == T) {
          estimationMETHOD <- 'SEM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- NULL
          }

          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- 4000
        } else if(unstable == T){
          estimationMETHOD <- 'SEM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- NULL
          }

          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- NULL
        } else if (forceNormalEM == T) {
          estimationMETHOD <- 'EM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- 1000
        } else {
          estimationMETHOD <- 'EM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL  <- NULL
          empiricalhist <- TRUE
          NCYCLES <- 4000
        }
      } else {
        if(unstable == T){
          if(forceNormalEM == T && i < 4){
            estimationMETHOD <- 'EM'
          } else {
            estimationMETHOD <- 'QMCEM'
          }
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }


          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- NULL
        } else {
          estimationMETHOD <- 'MHRM'
          optimINPUT <- NULL
          optimCTRL <- NULL
          empiricalhist <- FALSE
          NCYCLES <- 4000
        }
      }
      covdataINPUT <- NULL
      formulaINPUT <- NULL

      message('estimation method: ', paste0(estimationMETHOD))

    } else { # with covariate

      covdataINPUT <- covdata
      formulaINPUT <- formula
      if(NROW(random) != 0 | NROW(lr.random) != 0){

        estimationMETHOD <- 'MHRM'
        optimINPUT <- NULL
        optimCTRL <- NULL
        empiricalhist <- FALSE
        NCYCLES <- 4000

      } else if(forceMHRM == T| forceGRSM == T | assumingFake == T | masterThesis == T | sum(na.omit(ActualTestlets) > 2) != 0){
        message('MHRM currently not supported with latent regressors')
        estimationMETHOD <- 'QMCEM'
        if(forceDefaultOptimizer){
          optimINPUT <- NULL
        } else {
          optimINPUT <- 'nlminb'
        }

        optimCTRL  <- NULL
        empiricalhist <- FALSE
        NCYCLES <- 1e+4

      } else if(length(survey.weights) != 0) {
        estimationMETHOD <- 'QMCEM'
        if(forceDefaultOptimizer){
          optimINPUT <- NULL
        } else {
          optimINPUT <- 'nlminb'
        }

        optimCTRL  <- NULL
        empiricalhist <- FALSE
        NCYCLES <- NULL
      } else if(i < 4){ # use EM 1-3 factors
        if(unstable == T){
          estimationMETHOD <- 'QMCEM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL <- NULL
          empiricalhist <- FALSE
          NCYCLES <- NULL
        } else if (forceNormalEM == T) {
          estimationMETHOD <- 'EM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL  <- NULL
          empiricalhist <- FALSE
          NCYCLES <- 1000
        } else {
          estimationMETHOD <- 'EM'
          if(forceDefaultOptimizer){
            optimINPUT <- NULL
          } else {
            optimINPUT <- 'nlminb'
          }

          optimCTRL <- NULL
          empiricalhist <- TRUE
          NCYCLES <- 4000
        }
      } else {
        estimationMETHOD <- 'QMCEM'
        if(forceDefaultOptimizer){
          optimINPUT <- NULL
        } else {
          optimINPUT <- 'nlminb'
        }

        optimCTRL <- NULL
        empiricalhist <- FALSE
        NCYCLES <- NULL
      }

      message('estimation method: ', paste0(estimationMETHOD))
      if(NROW(fixed) != 0){
        message('fixed effect formula: ', paste0(fixed))
      }
      if(NROW(random) != 0){
        message('random effect formula: ', paste0(random))
      } else {
        message('latent regression formula: ', paste0(formulaINPUT))

      }
    }
    if(estimationMETHOD == 'EM'){
      message('Empirical Histogram for find Prior distribution: ', empiricalhist)
    }
    # forcing SE estimation activate
    if((sum(is.na(data.frame(x))) != 0) && (SE.type != 'MHRM')){
      SE <- T
      if(length(covdata) == 0){
        if(estimationMETHOD == 'MHRM'){ # Richadson (BL) isn't support MHRM estimation method
          SE.type <- 'MHRM'
        } #else {
        # SE.type <- "sandwich" # Oakes
        # }
      } else {
        if(estimationMETHOD == 'MHRM'){ # Richadson (BL) isn't support MHRM estimation method
          SE.type <- 'MHRM'
        } #else {
        # SE.type <- "Oakes" # Oakes
        # }
      }
    }

    # switch SE estimation temporary
    if(length(covdata) != 0 && NROW(random) == 0 && SE.type == 'sandwich' && SE == T){
      message('sandwich estimation currently not supproted with latent regressors')
      SE.type <- "complete"
    }

    # standard error estimation config
    if((SE == T && SE.type == 'SEM') == T){
      accelerateINPUT <- 'none'
      TOLINPUT <- NULL
      SEtolINPUT <- 1e-4
      symmetricINPUT <- FALSE

      message('TOL: ', 'default', ' / SEtol: ', SEtolINPUT, ' / SE.type: ', SE.type, ' / Accelerator: ',
              accelerateINPUT, ' / Symmetric SEM: ', symmetricINPUT)
    } else if((SE == T && estimationMETHOD == 'MHRM') == T){
      if(forceDefalutAccelerater){
        accelerateINPUT <- 'Ramsay'
      } else {
        accelerateINPUT <- 'squarem'
      }

      TOLINPUT <- NULL
      SEtolINPUT <- 1e-4
      symmetricINPUT <- FALSE
      if(EnableFMHRM){
        SE.type <- 'FMHRM'
      } else {
        SE.type <- 'MHRM' # more efficient
      }

      message('TOL: ', 'default', ' / SEtol: ', SEtolINPUT, ' / SE.type: ', SE.type, ' / Accelerator: ',
              accelerateINPUT, ' / Symmetric SEM: ', symmetricINPUT)
    } else {
      if(forceDefalutAccelerater){
        accelerateINPUT <- 'Ramsay'
      } else {
        accelerateINPUT <- 'squarem'
      }


      TOLINPUT <- NULL
      SEtolINPUT <- 1e-4
      symmetricINPUT <- FALSE

      if(SE == T){
        message('Standard Error Estimation: On')
        message('TOL: ', 'default', ' / SEtol: ', SEtolINPUT, ' / SE.type: ', SE.type, ' / Accelerator: ',
                accelerateINPUT, ' / Symmetric SEM: ', symmetricINPUT)
      }
    }

    # removeEmptyRows config
    if(length(covdata) != 0){
      removeEmptyRowsConf <- FALSE
    } else {
      removeEmptyRowsConf <- TRUE
    }



    if(sum(psych::describe(x)$range == 1) == ncol(x)){ # dichotomous items

      # forceRasch (dichotomous)
      if(forceRasch == T){
        message('\nMIRT model: Rasch')
        if(NROW(random) != 0 | NROW(lr.random) != 0){
          try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
                                        covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                        GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                        accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                        removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
        }
        try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'Rasch',
                                  accelerate = accelerateINPUT, calcNull = T,
                                  technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                   removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                  formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
            rm(modTEMP)
          }
        }
        if(exists('modTEMP_MIXED')){
          if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
            rm(modTEMP_MIXED)
          }
        }
        if(exists('modTEMP') && exists('modTEMP_MIXED')){
          if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }
        if(!exists('modTEMP') && exists('modTEMP_MIXED')){
          modTEMP <- modTEMP_MIXED
          rm(modTEMP_MIXED)
        }

        # try(invisible(mirt::mirtCluster(remove = T)), silent = T)

        if(exists('modTEMP')){
          try(return(modTEMP))
        } else {
          stop('Rasch model is inappropriate with this data. turn off the forceRasch = T option.')
        }
      }

      # if(nrow(x) >= 50){

      if(length(ActualTestlets) != 0){

        if(diagnosis == F){
          message('\nMIRT model: Compensatory 4PL')

          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '4PL',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }

          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '4PL',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)

          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }


        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 3PL with upper asymptote (slip) estimated')

          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '3PLu',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }

          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PLu',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Partially compensatory 3PL')

          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'PC3PL',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 3PL')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '3PL',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PL',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # }


        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 2PL')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '2PL',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '2PL',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Partially compensatory 2PL')
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'PC2PL',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && skipIdealPoint == F){

          message('\nMIRT model: ideal point')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'ideal',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'ideal',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F){

          message('\nMIRT model: Rasch')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'Rasch',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'Rasch',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(SE.type == 'MHRM'){
          message('MHRM trials were not efficient. try to calibrate the model with EM')
          SE.type <- 'sandwich'
        } else {
          message('try to recalibrate the model with mirt::bfactor() function')
        }


        if(exists('modTEMP') == F && diagnosis == F){
          message('\nMIRT model: Compensatory 4PL')

          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '4PL',
                                       lerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }


        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 3PL with upper asymptote (slip) estimated')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '3PLu',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Partially compensatory 3PL')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'PC3PL',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 3PL')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '3PL',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }

        # }


        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 2PL')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '2PL',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Partially compensatory 2PL')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'PC2PL',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }
        if(exists('modTEMP') == F && skipIdealPoint == F){

          message('\nMIRT model: ideal point')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'ideal',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }

        if(exists('modTEMP') == F){

          message('\nMIRT model: Rasch')
          try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'Rasch',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
        }

      } else {

        if(diagnosis == F){
          message('\nMIRT model: Compensatory 4PL')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '4PL',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '4PL', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }


        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 3PL with upper asymptote (slip) estimated')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '3PLu',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PLu', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && i == 1 && diagnosis == F){

          message('\nMIRT model: Partially compensatory 3PL')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'PC3PL', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 3PL')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '3PL',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PL', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # }

        if(exists('modTEMP') == F && diagnosis == F){

          message('\nMIRT model: Compensatory 2PL')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '2PL',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '2PL', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && i == 1 && diagnosis == F){

          message('\nMIRT model: Partially compensatory 2PL')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'PC2PL', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        if(exists('modTEMP') == F && skipIdealPoint == F){

          message('\nMIRT model: ideal point')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'ideal',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'ideal', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }

        }

        if(exists('modTEMP') == F && i == 1 && diagnosis == F){

          message('\nMIRT model: Rasch')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'Rasch', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }
      }


      if(exists('modTEMP') == F && i == 1){

        message('\nMIRT model: spline response')
        if(estimationMETHOD == 'MHRM'){
          estimationMETHOD <- 'EM'
          message('switching MHRM to ', estimationMETHOD, ' ...')

        }

        try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'spline', method = estimationMETHOD,
                                  accelerate = accelerateINPUT, calcNull = T,
                                  technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                   removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                  formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
            message('Model is unstable so that try to remedy the model automatically as soon as possible')
          }
        }
      }

      if(exists('modTEMP') == F){
        if(i == 1){
          stop('Fail to find Factor solutions: Model didn\'t converge.')
        } else {
          return(modOLD)
        }
      }

    } else if(length(itemkeys) != 0){ # 2-4PLNRM

      if(length(ActualTestlets) != 0){

        message('\nMIRT model: Compensatory 4PL Nominal response')
        try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '4PLNRM',
                                  key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                  technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                   removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                  formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
        }

        if(exists('modTEMP') == F){

          message('\nMIRT model: Compensatory 3PL Nominal response')
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PLNRM',
                                    key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }


        if(exists('modTEMP') == F){

          message('\nMIRT model: Compensatory 3PL Nominal response with upper asymptote estimated')
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PLuNRM',
                                    key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }
        # }

        if(exists('modTEMP') == F){

          message('\nMIRT model: Compensatory 2PL Nominal response')
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '2PLNRM',
                                    key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

        if(exists('modTEMP') == F){
          message('\nMIRT model: Nominal response without keys')
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'nominal',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, key = NULL, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              message('Model may unstable but Trying to remedy automatically')
            }
          }
        }

      } else {

        message('\nMIRT model: Compensatory 4PL Nominal response')
        try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '4PLNRM', method = estimationMETHOD,
                                  key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                  technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                   removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                  formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
        }

        if(exists('modTEMP') == F){

          message('\nMIRT model: Compensatory 3PL Nominal response')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PLNRM', method = estimationMETHOD,
                                    key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }


        if(exists('modTEMP') == F){

          message('\nMIRT model: Compensatory 3PL Nominal response with upper asymptote estimated')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PLuNRM', method = estimationMETHOD,
                                    key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }
        # }

        if(exists('modTEMP') == F){

          message('\nMIRT model: Compensatory 2PL Nominal response')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '2PLNRM', method = estimationMETHOD,
                                    key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

        if(exists('modTEMP') == F){
          message('\nMIRT model: Nominal response without keys')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'nominal', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, key = NULL, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              message('Model may unstable but Trying to remedy automatically')
            }
          }
        }
      }


      if(exists('modTEMP') == F){
        if(i == 1){
          stop('Fail to find Factor solutions: Model didn\'t converge.')
        } else {
          return(modOLD)
        }
      }

    } else if((sum(psych::describe(x)$range == 1) != 0) && allowMixedResponse == T) { # mixed format (Construct Responses + Multiple Choices, under construction)
      if(skipNominal == F){
        if(skipIdealPoint == F){

          itemtype_mixed <- vector()
          for(i in 1:ncol(x)){
            if(psych::describe(x[,i])$range == 1){
              itemtype_mixed[i] <- 'ideal'
              dichotomous_type <- 'ideal point'
            } else {
              itemtype_mixed[i] <- 'nominal'
            }
          }
          message('\nMIRT model: nominal response + ', paste0(dichotomous_type))
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = itemtype_mixed, method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 20000000000, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = 20000), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

        if(exists('modTEMP') == F){
          itemtype_mixed <- vector()
          for(i in 1:ncol(x)){
            if(psych::describe(x[,i])$range == 1){
              if(nrow(x) >= 5000){
                itemtype_mixed[i] <- '4PL'
                dichotomous_type <- '4PL'
              } else if(nrow(x) >= 2000){
                itemtype_mixed[i] <- '3PL'
                dichotomous_type <- '3PL'
              } else {
                itemtype_mixed[i] <- '2PL'
                dichotomous_type <- '2PL'
              }
            } else {
              itemtype_mixed[i] <- 'nominal'
            }
          }
          message('\nMIRT model: nominal response + ', paste0(dichotomous_type), ' (automatically set by sample size)')
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = itemtype_mixed, method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 20000000000, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = 20000), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

      }

      # generalized partial credit model (non-sequential)
      if(exists('modTEMP') == F && forceNRM == F){
        itemtype_mixed <- vector()
        for(i in 1:ncol(x)){
          if(psych::describe(x[,i])$range == 1){
            if(nrow(x) >= 5000){
              itemtype_mixed[i] <- '4PL'
            } else if(nrow(x) >= 2000){
              itemtype_mixed[i] <- '3PL'
            } else {
              itemtype_mixed[i] <- 'ideal'
            }
          } else {
            itemtype_mixed[i] <- 'gpcm'
          }
        }
        message('\nMIRT model: Generalized partial credit + ideal or 3-4PL')
        try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = itemtype_mixed, method = estimationMETHOD,
                                  accelerate = accelerateINPUT, calcNull = T,
                                  technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                   removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                  formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
        }
      }

      # graded response model (sequential)
      if(exists('modTEMP') == F && forceNRM == F){
        itemtype_mixed <- vector()
        for(i in 1:ncol(x)){
          if(psych::describe(x[,i])$range == 1){
            if(nrow(x) >= 5000){
              itemtype_mixed[i] <- '4PL'
            } else if(nrow(x) >= 2000){
              itemtype_mixed[i] <- '3PL'
            } else {
              itemtype_mixed[i] <- 'ideal'
            }
          } else {
            itemtype_mixed[i] <- 'graded'
          }
        }
        message('\nMIRT model: Graded response + ideal or 3-4PL')
        try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = itemtype_mixed, accelerate = accelerateINPUT,
                                  calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                 SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                  TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                  optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
        } else {
          stop('retry to add forceMHRM = TRUE argument')
        }
      }

      # finally, if can not converge
      if(exists('modTEMP') == F){
        if(i == 1){
          stop('Fail to find Factor solutions: Model didn\'t converge.')
        } else {
          return(modOLD)
        }
      }
    } else { # polytomous items

      # forceRasch (PCM)
      if(forceRasch == T){

        message('\nMIRT model: Rasch')
        if(NROW(random) != 0 | NROW(lr.random) != 0){
          try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
                                        covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                        GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                        accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                        removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
        }
        try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'Rasch',
                                  accelerate = accelerateINPUT, calcNull = T,
                                  technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                   removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                  formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                  SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
        if(exists('modTEMP')){
          if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
            rm(modTEMP)
          }
        }
        if(exists('modTEMP_MIXED')){
          if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
            rm(modTEMP_MIXED)
          }
        }
        if(exists('modTEMP') && exists('modTEMP_MIXED')){
          if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }
        if(!exists('modTEMP') && exists('modTEMP_MIXED')){
          modTEMP <- modTEMP_MIXED
          rm(modTEMP_MIXED)
        }

        # try(invisible(mirt::mirtCluster(remove = T)), silent = T)

        if(exists('modTEMP')){
          try(return(modTEMP))
        } else {
          stop('Rasch model is inappropriate with this data. turn off the forceRasch = T option.')
        }
      }

      # forceGRSM (polytomous)
      # for detecting fake responses
      # if(SE == T){ # if grsm, Standard error estimation was unsuccessful.

      # } else {
      #if((max(describe(x)$range) - min(describe(x)$range)) == 0 | forceGRSM == T | assumingFake == T | masterThesis == T){
      if(forceGRSM == T | assumingFake == T | masterThesis == T){
        x <- data.frame(x)

        if(length(which(describe(x)$max > median(describe(x)$max))) != 0){ # preventing weird input (e.g: 6 in 5 point scale)
          for(ii in which(describe(x)$max > median(describe(x)$max))){
            x[,ii] <-plyr::mapvalues(x[,ii],psych::describe(x[,which(describe(x)$max > median(describe(x)$max))])$max, NA)
          }
        }

        x <- x[,which(describe(x)$range == max(describe(x)$range))]
        k <- vector()
        for(iii in 1:length(x)){
          for(j in range(na.omit(x[iii]))[1]:range(na.omit(x[iii]))[2]){
            if((sum(na.omit(x[iii]) == j) == 0) == TRUE){
              k[length(k)+1] <- colnames(x[iii])
            }
          }
        }

        x <- (x[,!colnames(x) %in% k])
        for(iiii in 1:ncol(x)){
          x[,iiii] <- as.integer(x[,iiii])
        }

        message('\nMIRT model: graded rating scale')
        if(i == 1){

          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsmIRT', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }

          if(exists('modTEMP') == F){
            if(i == 1){
              stop('Fail to find Factor solutions: Model didn\'t converge.')
            } else {
              return(modOLD)
            }
          }

        } else {

          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsm', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }

          if(exists('modTEMP') == F){
            if(i == 1){
              stop('Fail to find Factor solutions: Model didn\'t converge.')
            } else {
              return(modOLD)
            }
          }

        }

        if(modTEMP@OptimInfo$converged == 1){
          skipNominal <- T
        }
      }
      # }

      # polytomous
      # nominal model

      if(length(ActualTestlets) != 0){

        if(skipNominal == F){
          message('\nMIRT model: nominal response')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'nominal',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'nominal',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # generalized partial credit model (non-sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Generalized partial credit')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'gpcm',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'gpcm',
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # graded response model (sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Graded response')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'graded',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
                                    calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                   SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                    TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                    optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }




        if(SE.type == 'MHRM'){
          message('MHRM trials were not efficient. try to calibrate the model with EM')
          SE.type <- 'sandwich'
        } else {
          message('try to recalibrate the model with mirt::bfactor() function')
        }

        if(exists('modTEMP') == F && skipNominal == F){
          message('\nMIRT model: nominal response')
          try(modTEMP <- mirt::bfactor(data = x, model = ActualTestlets, itemtype = 'nominal',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

        # generalized partial credit model (non-sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Generalized partial credit')
          try(modTEMP <- mirt::bfactor(data = x, model = ActualTestlets, itemtype = 'gpcm',
                                       accelerate = accelerateINPUT, calcNull = T,
                                       technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                        removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                       formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

        # graded response model (sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Graded response')
          try(modTEMP <- mirt::bfactor(data = x, model = ActualTestlets, accelerate = accelerateINPUT,
                                       calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                      SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                       TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                       optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                       SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
          }
        }

      } else {
        if(skipNominal == F){
          message('\nMIRT model: nominal response')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'nominal',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'nominal', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # generalized partial credit model (non-sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Generalized partial credit')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'gpcm',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'gpcm', method = estimationMETHOD,
                                    accelerate = accelerateINPUT, calcNull = T,
                                    technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
                                                     removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
                                    formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # graded response model (sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Graded response')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'graded',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
                                    calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                   SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                    TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                    optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # graded response model (sequential)
        if(exists('modTEMP') == F && forceNRM == F){
          message('\nMIRT model: Graded rating scale')

          if(i == 1){
            if(NROW(random) != 0 | NROW(lr.random) != 0){
              try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'grsmIRT',
                                            covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                            GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                            accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                            removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
            }
            try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsmIRT', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
                                      calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                     SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                      TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                      optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                      SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
          } else {
            if(NROW(random) != 0 | NROW(lr.random) != 0){
              try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'grsm',
                                            covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                            GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                            accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                            removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
            }
            try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsm', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
                                      calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                     SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                      TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                      optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                      SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
          }

          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # graded response model (sequential)
        if(exists('modTEMP') == F && forceNRM == F && i == 1){
          message('\nMIRT model: Partial Credit')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'Rasch', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
                                    calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                   SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                    TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                    optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)

          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }

        # graded response model (sequential)
        if(exists('modTEMP') == F && forceNRM == F && i == 1){
          message('\nMIRT model: Rating Scale')
          if(NROW(random) != 0 | NROW(lr.random) != 0){
            try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'rsm',
                                          covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                                          GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
                                          accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
                                          removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
          }
          try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'rsm', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
                                    calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
                                                                   SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
                                    TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
                                    optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
                                    SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)

          if(exists('modTEMP')){
            if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP)
            }
          }
          if(exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
              rm(modTEMP_MIXED)
            }
          }
          if(exists('modTEMP') && exists('modTEMP_MIXED')){
            if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
              modTEMP <- modTEMP_MIXED
              rm(modTEMP_MIXED)
            }
          }
          if(!exists('modTEMP') && exists('modTEMP_MIXED')){
            modTEMP <- modTEMP_MIXED
            rm(modTEMP_MIXED)
          }
        }
      }




      # finally, if can not converge
      if(exists('modTEMP') == F){
        if(i == 1){
          stop('Fail to find Factor solutions: Model didn\'t converge.')
        } else {
          return(modOLD)
        }
      }
    }

    try(invisible(mirt::mirtCluster(remove = T)), silent = T)

    if(class(modTEMP)[1] == 'MixedClass'){
      MLM_rotate_formula_mod <- mirt::mirt(data = modTEMP@Data$data, model = modTEMP@Model$model, itemtype = modTEMP@Model$itemtype, pars = 'values')
      MLM_rotate_formula_mod_original <- mirt::mod2values(modTEMP)
      if(sum(MLM_rotate_formula_mod_original$name == '(Intercept)') != 0){
        MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original[!MLM_rotate_formula_mod_original$name == '(Intercept)',]

      }
      # MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original
      MLM_rotate_formula_mod$value[which(MLM_rotate_formula_mod$item %in% colnames(modTEMP@Data$data))] <- MLM_rotate_formula_mod_original$value[which(MLM_rotate_formula_mod_original$item %in% colnames(modTEMP@Data$data))]
      MLM_rotate_formula_mod$est <- F
      # print(MLM_rotate_formula_mod)

      MixedModelFlag <- T

      modTEMP_MIXED <- modTEMP

      modTEMP <- mirt::mirt(data = modTEMP@Data$data, model = modTEMP@Model$model,
                            itemtype = modTEMP@Model$itemtype, pars = MLM_rotate_formula_mod, method = 'QMCEM',
                            GenRandomPars = GenRandomPars,
                            calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws),
                            survey.weights = survey.weights)

    } else {
      MixedModelFlag <- F
    }


    if(i == 1 && length(ActualTestlets) == 0 && plotOn == T){ # ICC printing
      try(print(plot(modTEMP, type = 'infoSE')))
      try(print(plot(modTEMP, type = 'infotrace', facet_items = TRUE)))
      try(print(plot(modTEMP, type = 'trace')))
    }

    if(forceUIRT == T | TestletActivated == T){
      if(MixedModelFlag){
        return(modTEMP_MIXED)
      } else {
        return(modTEMP)
      }
    }

    if(i == 1 && modTEMP@OptimInfo$converged == F && sum(modTEMP@Model$itemtype %in% c('spline', 'nominal')) == 0){
      stop('No convergence')
    }


    if(i > 1){ # evaluation of local optimal number of factors
      if(ncol(modTEMP@Fit$F) == 1){
        rotMat <- modTEMP@Fit$F
      } else {
        rotMat <- GPArotation::geominQ(modTEMP@Fit$F, maxit = 1e+5)$loadings
      }

      if(MixedModelFlag){
        modTEMP <- modTEMP_MIXED
      }

      if(modTEMP@Fit$DIC > modOLD@Fit$DIC | modTEMP@OptimInfo$converged == F | sum(round(modTEMP@Fit$h2, 3) >= .999) != 0){ # modTEMP@Fit$AICc > modOLD@Fit$AICc |
        message('optimal factor numbers: ', paste0(i-1))
        return(modOLD)
      } #else if(sum(colSums(round(abs(rotMat), 2) > .4) < 2) != 0) {
      #message('optimal factor numbers: ', paste0(i-1))
      #return(modOLD)
      #}

    }

    if(exists('modTEMP') == T){
      if(MixedModelFlag){
        modOLD <- (modTEMP_MIXED)
        try(rm(modTEMP_MIXED), silent = T)
        try(rm(modTEMP), silent = T)

      } else {
        modOLD <- modTEMP
        try(rm(modTEMP), silent = T)
      }

    } else {
      stop('Fail to convergence')
    }
  }
}


surveyFA <- function(data = ..., covdata = NULL, formula = NULL, SE = T,
                     SE.type = "sandwich", skipNominal = F, forceGRSM = F,
                     assumingFake = F, masterThesis = F, forceRasch = F,
                     unstable = F, forceNormalEM = T, forceMHRM = F,
                     printFactorStructureRealtime = F, itemkeys = NULL,
                     survey.weights = NULL, allowMixedResponse = T, autofix = F,
                     forceUIRT = F, skipIdealPoint = F, MHRM_SE_draws = 1e+4,
                     bifactorSolution = T, skipS_X2 = F, forceNRM = F, needGlobalOptimal = T,
                     pilotTestMode = F, forceConsiderPositiveZh = F, forceDefalutAccelerater = F, forceDefaultOptimizer = F, EnableFMHRM = F, testlets = NULL,
                     minimumLeftItems = 3, plotOn = T, forceCTTmode = F, GenRandomPars = T, coefAlwaysBePositive = T,
                     fixed = ~1, random = NULL, lr.fixed = ~1, lr.random = NULL, MH_BURNIN = 1500, MH_SEMCYCLES = 1000, ...) {
  message('---------------------------------------------------------')
  message(' k.aefa: kwangwoon automated exploratory factor analysis ')
  message('---------------------------------------------------------\n')


  message('Calculating Initial Factor model')
  iteration_num <- 1
  message('Iteration: ', iteration_num, '\n')



  if(bifactorSolution) {
    rotateCriteria <- 'bifactorQ'
  } else {
    rotateCriteria <- 'geominQ'
  }

  if(!is.data.frame(data) && !is.matrix(data)){
    surveyFixMod <- data
  } else {
    if(sum(is.na(data)) != 0 && pilotTestMode == T){

      message('input data n size: ', nrow(data))
      data <- data[which(rowSums(is.na(data)) < ncol(data)*(1-3/4)),]
      message('current data n size: ', nrow(data))
    }

    surveyFixMod <- fastFIFA(x = as.data.frame(data), covdata = as.data.frame(covdata),
                             formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal,
                             forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis,
                             forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM,
                             itemkeys = itemkeys, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse,
                             autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                             forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = testlets, GenRandomPars = GenRandomPars,
                             fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
  }

  workKeys <- itemkeys
  workTestlets <- testlets
  if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
    surveyFixMod <- deepFA(surveyFixMod, survey.weights)
  }
  itemFitDone <- FALSE
  while (!itemFitDone) {
    try(invisible(gc()), silent = T) # garbage cleaning

    surveyFixModRAW <- data.frame(surveyFixMod@Data$data)
    surveyFixModCOV <- data.frame(attr(surveyFixMod@ParObjects$lrPars, "df"))

    if(ncol(surveyFixModRAW) > minimumLeftItems){

      if(class(surveyFixMod)[1] == 'MixedClass'){
        MLM_rotate_formula_mod <- mirt::mirt(data = surveyFixMod@Data$data, model = surveyFixMod@Model$model, itemtype = surveyFixMod@Model$itemtype, pars = 'values')
        MLM_rotate_formula_mod_original <- mirt::mod2values(surveyFixMod)
        if(sum(MLM_rotate_formula_mod_original$name == '(Intercept)') != 0){
          MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original[!MLM_rotate_formula_mod_original$name == '(Intercept)',]

        }
        # MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original
        MLM_rotate_formula_mod$value[which(MLM_rotate_formula_mod$item %in% colnames(surveyFixMod@Data$data))] <- MLM_rotate_formula_mod_original$value[which(MLM_rotate_formula_mod_original$item %in% colnames(surveyFixMod@Data$data))]
        MLM_rotate_formula_mod$est <- F
        # print(MLM_rotate_formula_mod)

        MixedModelFlag <- T

        surveyFixMod <- mirt::mirt(data = surveyFixMod@Data$data, model = surveyFixMod@Model$model,
                                   itemtype = surveyFixMod@Model$itemtype, pars = MLM_rotate_formula_mod, method = 'QMCEM',
                                   GenRandomPars = GenRandomPars,
                                   calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws),
                                   survey.weights = survey.weights)

      } else {
        MixedModelFlag <- F
      }

      message('\nChecking item local independence assumption')
      iteration_num <- iteration_num + 1
      message('Iteration: ', iteration_num, '\n')

      if(sum(surveyFixMod@Model$itemtype %in% 'spline') > 0){
        fscoreMethod <- 'EAP'
      } else {
        fscoreMethod <- 'MAP'
      }

      if(sum(is.na(surveyFixModRAW)) == 0 | nrow(surveyFixModRAW) > 5000){
        NofCores <- parallel::detectCores()
        NofCores <- round(NofCores / 1.1)
        if(NofCores > 8){
          NofCores <- 8
        }
        try(invisible(mirt::mirtCluster(spec = NofCores)), silent = T)
      }

      if(!forceCTTmode){
        if(sum(is.na(surveyFixMod@Data$data)) == 0){
          try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('S_X2', 'Zh', 'infit'),
                                                    method = fscoreMethod,
                                                    QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
        } else {
          try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('S_X2', 'Zh', 'infit'),
                                                    impute = 100,
                                                    method = fscoreMethod,
                                                    QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
        }

        if(!exists('surveyFixMod_itemFit')){
          if(sum(is.na(surveyFixMod@Data$data)) == 0){
            try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('Zh', 'infit'),
                                                      method = fscoreMethod,
                                                      QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
          } else {
            try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('Zh', 'infit'),
                                                      impute = 100,
                                                      method = fscoreMethod,
                                                      QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
          }
          S_X2ErrorFlag <- T
        } else {
          S_X2ErrorFlag <- F
        }
        try(invisible(mirt::mirtCluster(remove = T)), silent = T)

        if(!S_X2ErrorFlag){
          if(sum(surveyFixMod_itemFit$df.S_X2 == 0, na.rm = T) > length(surveyFixMod_itemFit$df.S_X2)/5){
            S_X2ErrorFlag <- T
          } else if(sum(is.nan(surveyFixMod_itemFit$p.S_X2)) > length(surveyFixMod_itemFit$p.S_X2)/5) {
            S_X2ErrorFlag <- T
          }
        }

        print(surveyFixMod_itemFit)

        # item evaluation

        if(S_X2ErrorFlag){
          message('S_X2 can not be calcuate normally...')
        }

      }

      # precalculation of CI for a1
      ZeroList <- vector()
      ZeroRange <- vector()

      if(surveyFixMod@Options$SE == T && (length(workTestlets) != 0 | surveyFixMod@Model$nfact == 1)){
        for(III in 1:NROW(coef(surveyFixMod))-1){
          try(vec <- data.frame(coef(surveyFixMod)[III]), silent = T)

          if(exists('vec')){
            # print(vec) # -- for debug
            if(is.na((vec[2,1] < 0 && vec[3,1] > 0))){

            } else {
              try(ZeroList[length(ZeroList)+1] <- (vec[2,1] < 0 && vec[3,1] > 0), silent = T)
              try(ZeroRange[length(ZeroRange)+1] <- psych::describe(c(vec[2,1], vec[3,1]))$range, silent = T)
            }
          }
          # print(ZeroList)# -- for debug
          # print(ZeroRange)
        }
      }

      vec2 <- coef(surveyFixMod, simplify = T)$items

      # Standardized Log-Likelihood
      if(forceCTTmode){ # force CTT-Like model
        if(sum(ZeroList) > 0){ # which item include 0 in a1
          message('\nItem discrimination include 0 / removing ', paste(colnames(surveyFixModRAW)[which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]))
          workKeys <- workKeys[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]
          workTestlets <- workTestlets[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }


        } else if(coefAlwaysBePositive && sum(vec2[,1] < 0) != 0){
          message('\nItem discrimination include negative / removing ', paste(surveyFixMod_itemFit$item[which(min(vec2[,1]) == vec2[,1])]))
          workKeys <- workKeys[-which(min(vec2[,1]) == vec2[,1])]
          workTestlets <- workTestlets[-which(min(vec2[,1]) == vec2[,1])]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(vec2[,1]) == vec2[,1])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
        } else {
          itemFitDone <- TRUE
        }
      } else if(length(which(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems] < -1.96)) != 0 && sum(is.na(surveyFixModRAW)) == 0){ # Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.
        message('\nDrasgow, F., Levine, M. V., & Williams, E. A. (1985) / removing ', paste(surveyFixMod_itemFit$item[which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]))
        workKeys <- workKeys[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
        workTestlets <- workTestlets[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
        surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                 forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                 fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
        if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
          try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
        }
        rm(surveyFixMod_itemFit)
        rm(S_X2ErrorFlag)
      } else if(length(which(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems] == -Inf)) != 0){ # Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.
        message('\nDrasgow, F., Levine, M. V., & Williams, E. A. (1985) / removing ', paste(surveyFixMod_itemFit$item[which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]))
        workKeys <- workKeys[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
        workTestlets <- workTestlets[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]

        surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                 forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                 fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
        if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
          try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
        }
        rm(surveyFixMod_itemFit)
        rm(S_X2ErrorFlag)


      } else if(sum(ZeroList) > 0){ # which item include 0 in a (currently, only a1)
        message('\nItem discrimination include 0 / removing ', paste(surveyFixMod_itemFit$item[which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]))
        workKeys <- workKeys[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]
        workTestlets <- workTestlets[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]

        surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                 forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                 fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
        if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
          try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
        }
        rm(surveyFixMod_itemFit)
        rm(S_X2ErrorFlag)


      } else if(coefAlwaysBePositive && sum(vec2[,1] < 0) != 0){
        message('\nItem discrimination include negative / removing ', paste(surveyFixMod_itemFit$item[which(min(vec2[,1]) == vec2[,1])]))
        workKeys <- workKeys[-which(min(vec2[,1]) == vec2[,1])]
        workTestlets <- workTestlets[-which(min(vec2[,1]) == vec2[,1])]

        surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(vec2[,1]) == vec2[,1])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                 forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                 fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
        if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
          try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
        }
        rm(surveyFixMod_itemFit)
        rm(S_X2ErrorFlag)
      } else if(S_X2ErrorFlag) { # if Chi-squared can not be calculate with extremely big Zh
        if(length(which(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems] > 1.96)) != 0 && sum(is.na(surveyFixModRAW)) == 0 && forceConsiderPositiveZh){ # Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.
          message('\nDrasgow, F., Levine, M. V., & Williams, E. A. (1985) / removing ', paste(surveyFixMod_itemFit$item[which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]))
          workKeys <- workKeys[-which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
          workTestlets <- workTestlets[-which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
          rm(surveyFixMod_itemFit)
          rm(S_X2ErrorFlag)

        } else {
          itemFitDone <- TRUE

        }
      } else if(S_X2ErrorFlag == F && (surveyFixMod@Model$nfact == 1 | length(workTestlets) > 0)) { # if Chi-squared can be calculate

        # chi-square exception processing
        if(sum(is.na(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) == 0 && sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) == 0){ # avoid unexpected situation

          message('all items df are 0. skipping evaluation...')
          itemFitDone <- TRUE

        } else if(sum(is.na(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) != 0 && sum(is.na(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) != surveyFixMod@Data$nitems && (sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == 0) < length(1:surveyFixMod@Data$nitems)/2)){
          message('\nremoving items df is NA / ', paste(surveyFixMod_itemFit$item[which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)]))
          workKeys <- workKeys[-which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)]
          workTestlets <- workTestlets[-which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
          rm(surveyFixMod_itemFit)
          rm(S_X2ErrorFlag)


        } else if(sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == 0) != 0 && skipS_X2 == F && (sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == 0) < length(1:surveyFixMod@Data$nitems)/2)){
          message('\nremoving items df is 0 / removing ', paste(surveyFixMod_itemFit$item[which(surveyFixMod_itemFit$df.S_X2 == 0)]))
          workKeys <- workKeys[-which(surveyFixMod_itemFit$df.S_X2 == 0)]
          workTestlets <- workTestlets[-which(surveyFixMod_itemFit$df.S_X2 == 0)]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(surveyFixMod_itemFit$df.S_X2 == 0)], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
          rm(surveyFixMod_itemFit)
          rm(S_X2ErrorFlag)


        } else if(sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) != 0 && sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) != surveyFixMod@Data$nitems){
          message('\nremoving items p.S_X2 is NA / removing ', paste(surveyFixMod_itemFit$item[which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)]))
          workKeys <- workKeys[-which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)]
          workTestlets <- workTestlets[-which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
          rm(surveyFixMod_itemFit)
          rm(S_X2ErrorFlag)


        } else if(sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) == 0 && length(which(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems] >= 3)) != 0 && skipS_X2 == F && (surveyFixMod@Model$nfact == 1 | length(workTestlets) > 0)){ # Drasgow, F., Levine, M. V., Tsien, S., Williams, B., & Mead, A. D. (1995). Fitting polytomous item response theory models to multiple-choice tests. Applied Psychological Measurement, 19(2), 143-166.
          message('\nDrasgow, F., Levine, M. V., Tsien, S., Williams, B., & Mead, A. D. (1995) / removing ', paste(surveyFixMod_itemFit$item[which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]))
          workKeys <- workKeys[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]
          workTestlets <- workTestlets[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
          rm(surveyFixMod_itemFit)
          rm(S_X2ErrorFlag)


        } else if(sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) == 0 && length(which(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems] < .05)) != 0 && skipS_X2 == F && (surveyFixMod@Model$nfact == 1 | length(workTestlets) > 0)){ # Kang, T., & Chen, T. T. (2008). Performance of the Generalized S?€륺2 Item Fit Index for Polytomous IRT Models. Journal of Educational Measurement, 45(4), 391-406.; Reise, S. P. (1990). A comparison of item- and person-fit methods of assessing model-data fit in IRT. Applied Psychological Measurement, 14, 127-137.
          message('\nKang, T., & Chen, T. T. (2008); Reise, S. P. (1990) / removing ', paste(surveyFixMod_itemFit$item[which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]))
          workKeys <- workKeys[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]
          workTestlets <- workTestlets[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]

          surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
          if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
            try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
          }
          rm(surveyFixMod_itemFit)
          rm(S_X2ErrorFlag)


        } else {
          itemFitDone <- TRUE

        }


      } else if (forceRasch == T && surveyFixMod@Model$nfact == 1) {
        if(length(which(abs(surveyFixMod_itemFit$z.infit[1:surveyFixMod@Data$nitems]) > 2)) != 0){

          message('\nRasch outfit (|z|>2)')
          surveyFixMod <- fastFIFA(surveyFixModRAW[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))],
                                   itemkeys = itemkeys[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)

        } else if(length(which(abs(surveyFixMod_itemFit$z.infit[1:surveyFixMod@Data$nitems]) > 2)) != 0){

          message('\nRasch infit (|z|>2)')
          surveyFixMod <- fastFIFA(surveyFixModRAW[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))],
                                   itemkeys = itemkeys[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                   forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM,
                                   fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)

        } else {
          itemFitDone <- TRUE
        }
      } else {
        itemFitDone <- TRUE
      }
    } else {
      itemFitDone <- TRUE
    }


    if(printFactorStructureRealtime == T){
      message('\nRealtime Factor Structure after iteration')
      if(ncol(surveyFixMod@Fit$F)==1){
        print(round(surveyFixMod@Fit$F, 2))
      } else {
        if(bifactorSolution){
          print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
        } else {
          print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
        }
      }
    }
  } # the end of while loop

  if(autofix == T){ # GET RID OF ASAP
    #
    # message('\nChecking aberrant responses')
    # iteration_num <- iteration_num + 1
    # message('Iteration: ', iteration_num, '\n')
    #
    # surveyFixModRAW <- data.frame(mirt::extract.mirt(surveyFixMod, 'data'))
    # surveyFixModCOV <- data.frame(attr(surveyFixMod@ParObjects$lrPars, "df"))
    #
    # noAberrant <- k.faking(surveyFixModRAW, IRTonly = T, itemkeys = itemkeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM, ...)
    # if(length(covdata) == 0){ # anyway, covdata is NULL
    #   surveyFixMod <- fastFIFA(surveyFixModRAW[which(noAberrant$normal==TRUE),], itemkeys = itemkeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights[which(noAberrant$normal==TRUE)], allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM, ...)
    # } else {
    #   covdata_workout <- surveyFixModCOV
    #   surveyFixMod <- fastFIFA(surveyFixModRAW[which(noAberrant$normal==TRUE),], itemkeys = itemkeys, covdata = covdata_workout[which(noAberrant$normal==TRUE),], formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights[which(noAberrant$normal==TRUE)], allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM, ...)
    # }

    if(printFactorStructureRealtime == T){
      message('\nRealtime Factor Structure after iteration')
      if(ncol(surveyFixMod@Fit$F)==1){
        print(round(surveyFixMod@Fit$F, 2))
      } else {
        if(bifactorSolution){
          print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
        } else {
          print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
        }
      }
    }

    # autofix config
    fixFactorStructure_Done <- FALSE
    surveyFixMod_Workout <- surveyFixMod

    while (!fixFactorStructure_Done) { # start of while loop
      surveyFixModRAW <- data.frame(surveyFixMod_Workout@Data$data) # update data.frame
      surveyFixModCOV <- data.frame(attr(surveyFixMod_Workout@ParObjects$lrPars, "df"))

      message('\nFixing Factor Model')
      iteration_num <- iteration_num + 1
      message('Iteration: ', iteration_num, '\n')

      tempG <- surveyFixMod_Workout@Model$itemtype
      if(sum(tempG != 'Rasch') != 0){
        LowCommunalities <- surveyFixMod_Workout@Fit$h2[which(min(surveyFixMod_Workout@Fit$h2) == surveyFixMod_Workout@Fit$h2)]
      }


      if(ncol(surveyFixMod_Workout@Fit$F) == 1){
        Fmatrix <- surveyFixMod_Workout@Fit$F
      } else {
        if(bifactorSolution){
          Fmatrix <- GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings
        } else {
          Fmatrix <- GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings
        }
      }

      NoLoadings <- surveyFixMod_Workout@Fit$h2[which(rowSums(abs(round(Fmatrix, 2)) < .4) == ncol(surveyFixMod_Workout@Fit$F))]


      # h2 have to >= .3
      if(sum(tempG != 'Rasch') == 0){
        fixFactorStructure_Done <- TRUE
      } else if(LowCommunalities < .3^2){

        surveyFixMod_New <- fastFIFA(surveyFixModRAW[,-which(min(surveyFixMod_Workout@Fit$h2) == surveyFixMod_Workout@Fit$h2)], itemkeys = itemkeys[-which(min(surveyFixMod_Workout@Fit$h2) == surveyFixMod_Workout@Fit$h2)], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                     forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, ...)
        surveyFixMod_Workout <- surveyFixMod_New

        if(printFactorStructureRealtime == T){
          message('\nRealtime Factor Structure after iteration')
          if(ncol(surveyFixMod_Workout@Fit$F)==1){
            print(round(surveyFixMod_Workout@Fit$F, 2))
          } else {
            if(bifactorSolution){
              print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
            } else {
              print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
            }
          }
        }

      } else if(length(NoLoadings) != 0){ # noloadings
        if(as.logical(length(names(which(NoLoadings == min(NoLoadings))) != 0))){
          surveyFixMod_New <- fastFIFA(surveyFixModRAW[,!colnames(surveyFixModRAW) %in% names(which(NoLoadings == min(NoLoadings)))], itemkeys = itemkeys[,!colnames(surveyFixModRAW) %in% names(which(NoLoadings == min(NoLoadings)))], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
                                       forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, ...)
          surveyFixMod_Workout <- surveyFixMod_New

          if(printFactorStructureRealtime == T){
            message('\nRealtime Factor Structure after iteration')
            if(ncol(surveyFixMod_Workout@Fit$F)==1){
              print(round(surveyFixMod_Workout@Fit$F, 2))
            } else {
              if(bifactorSolution){
                print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
              } else {
                print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
              }
            }
          }

        }
      } else {
        fixFactorStructure_Done <- TRUE
      }

    } # the end of while loop

    if(printFactorStructureRealtime == T){
      message('\nFinal Factor Structure')
      if(ncol(surveyFixMod_Workout@Fit$F)==1){
        print(round(surveyFixMod_Workout@Fit$F, 2))
      } else {
        if(bifactorSolution){
          print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
        } else {
          print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
        }
      }
    }
    return(surveyFixMod_Workout)
  } else {
    return(surveyFixMod)
  }

}

numericMI <- function(model = ..., data = ..., m = 100, fun = 'sem', estimator = 'MLMV', parameterization = 'delta', chi = 'mplus', ...) {

  #########################
  # Multiple              #
  # Imputation for        #
  # variables             #
  # in SEM context        #
  #########################
  # Seongho Bae           #
  # seongho@kw.ac.kr      #
  # May 6th 2016          #
  # Under GNU / GPL       #
  #########################

  # checking packages for multiple impuation in SEM context
  if(!require('semTools')) {
    try(install.packages('semTools', dependencies = TRUE), silent=TRUE)
  }

  if(!require('Amelia')) {
    try(install.packages('Amelia', dependencies = TRUE), silent=TRUE)
  }

  if(!require('lavaan')) {
    try(install.packages('lavaan', dependencies = TRUE), silent=TRUE)
  }

  # loading packages for multiple imputation in SEM context
  library('semTools')
  library('Amelia')
  library('lavaan')

  fit <- sem(model=model, data=data.frame(data)) # extract variable names in model syntax
  message("sample size (listwise): ", paste0(nrow(data.frame(fit@Data@X))))

  fit_MI <- runMI(model=model, data=data.frame(data[,attributes(fit)$Model@dimNames[[1]][[1]]]), m=m, fun = fun, estimator = estimator, chi = chi, ...)#, control=list(optim.method="L-BFGS-B"), ...)

  cat(summary(fit_MI, standardize=T))
  print(inspect(fit_MI, 'fit'))

  return(fit_MI)
  message('inspect(fit, "impute")')

}
cmvFA <- function(x, MHRM = F){
  initialModel <- surveyFA(x, autofix = F, forceMHRM = MHRM, bifactorSolution = T)
  workModel <- initialModel
  STOP <- FALSE

  # rotation
  while (!STOP) {
    try(invisible(gc()), silent = T)
    if(ncol(workModel@Fit$F) > 1){
      if(workModel@Model$nfact == 2){
        rotSumMat <- GPArotation::bifactorQ(workModel@Fit$F, maxit = 1e+6)$loadings[,2]
      } else {
        rotSumMat <- data.frame(GPArotation::bifactorQ(workModel@Fit$F, maxit = 1e+6)$loadings[,2:workModel@Model$nfact])
      }

      # evaluation
      evaluationMat <- abs(rotSumMat) < .1
      print(rotSumMat)
      print(evaluationMat)

      if(is.data.frame(rotSumMat)){
        evaluationMat <- rowSums(evaluationMat)
        rotSumMat <- rowSums(abs(rotSumMat))
      } else {
        evaluationMat <- as.numeric(evaluationMat)
        rotSumMat <- (abs(rotSumMat))
      }
      print(rotSumMat)
      print(evaluationMat)
      print(sum(evaluationMat))
      print(ncol(workModel@Data$data))

      if(sum(evaluationMat) == ncol(workModel@Data$data)){ # number of inappropreate items should equal with number of items in model
        message('Done')
        return(workModel)
      } else {
        message(paste0(colnames(workModel@Data$data[,which(abs(rotSumMat) == min(abs(rotSumMat)))])))
        workModel <- fastFIFA(workModel@Data$data[,-which(abs(rotSumMat) == min(abs(rotSumMat)))], forceMHRM = MHRM)
      }

    } else { # if one dimensional model
      rotSumMat <- workModel@Fit$F
      evaluationMat <- abs(rotSumMat) < .1
      print(rotSumMat)
      print(evaluationMat)
      evaluationMat <- as.numeric(evaluationMat)
      rotSumMat <- (abs(rotSumMat))
      if(sum(evaluationMat) == ncol(workModel@Data$data) | sum(evaluationMat) == 0 | ncol(workModel@Data$data) <= 5){ # number of inappropreate items should equal with number of items in model
        message('Done')
        return(workModel)
      } else {
        message(paste0(colnames(workModel@Data$data[,which(abs(rotSumMat) == min(abs(rotSumMat)))])))
        workModel <- fastFIFA(workModel@Data$data[,-which(abs(rotSumMat) == min(abs(rotSumMat)))], forceMHRM = MHRM)
      }

      # return(workModel)
    }
  }

}

bifactorFA <- function(data = ..., skipS_X2 = F, forceNormalEM = T, forceMHRM = F, covdata = NULL, formula = NULL, skipNominal = F, allowMixedResponse = T, itemkeys = NULL, needGlobalOptimal = T) {
  if(length(covdata) != 0){
    forceNormalEM <- TRUE
  } else if(forceNormalEM){
    forceNormalEM <- TRUE
  } else {
    forceNormalEM <- FALSE
  }
  mod <- surveyFA(data = data, bifactorSolution = T, skipS_X2 = skipS_X2, forceMHRM = forceMHRM, autofix = F, covdata = covdata, formula = formula, skipNominal = skipNominal, allowMixedResponse = allowMixedResponse, itemkeys = itemkeys, needGlobalOptimal = needGlobalOptimal, forceNormalEM = forceNormalEM, forceUIRT = F)
  STOP <- FALSE
  while (!STOP) {
    if(ncol(mod@Fit$F) == 1){
      rotMAT <- abs(mod@Fit$F)
    } else {
      rotMAT <- abs(GPArotation::bifactorQ(mod@Fit$F, maxit = 1e+5)$loadings)[,1]
    }

    print(rotMAT)

    if(sum(rotMAT < .9999) != ncol(mod@Data$data)){
      mod <- surveyFA(data = mod@Data$data[,-which(rotMAT == max(rotMAT))], bifactorSolution = T, skipS_X2 = skipS_X2, forceMHRM = forceMHRM, autofix = F, covdata = covdata, formula = formula, skipNominal = skipNominal, allowMixedResponse = allowMixedResponse, itemkeys = itemkeys[-which(rotMAT == max(rotMAT))], needGlobalOptimal = needGlobalOptimal, forceNormalEM = forceNormalEM, forceUIRT = F)
    } else if(sum(rotMAT > .0001) != ncol(mod@Data$data)){
      mod <- surveyFA(data = mod@Data$data[,-which(rotMAT == min(rotMAT))], bifactorSolution = T, skipS_X2 = skipS_X2, forceMHRM = forceMHRM, autofix = F, covdata = covdata, formula = formula, skipNominal = skipNominal, allowMixedResponse = allowMixedResponse, itemkeys = itemkeys[-which(rotMAT == max(rotMAT))], needGlobalOptimal = needGlobalOptimal, forceNormalEM = forceNormalEM, forceUIRT = F)
    } else {
      return(mod)
    }

  }
}


findMLCA <- function(data = ..., startN = 1, empiricalhist = F, group = NULL, empiricaloptimal = T){
  # try(invisible(gc()), silent = T)
  DICindices <- vector()
  j <- 0
  nfact <- vector()

  if(sum(psych::describe(data)$range == 0) == 0){
    workData <- data
  } else {
    workData <- data[,-which(psych::describe(data)$range == 0)]
  }
  message('starting find global optimal of latent class')
  for(i in startN:ncol(workData)){
    try(invisible(tempModel <- mirt::mdirt(data = workData, model = i, empiricalhist = empiricalhist, group = group, nruns = 100)), silent = F)
    if(exists('tempModel')){
      if(tempModel@OptimInfo$converged){
        message('class: ', i, ' / DIC: ', tempModel@Fit$DIC)
        j <- j + 1
        DICindices[j] <- tempModel@Fit$DIC
        nfact[j] <- i
        rm(tempModel)
      }
    }
  }
  bestModel <- which(min(DICindices) == DICindices)
  message('find global optimal: ', nfact[bestModel])

  if(empiricaloptimal){

    testFS <- fscores(mirt::mdirt(data = workData, model = nfact[bestModel], empiricalhist = empiricalhist, group = group, nruns = 100), QMC = T)

    membership <- vector()
    for(i in 1:nrow(testFS)){
      membership[i] <- which(testFS[i,] == max(testFS[i,]))
    }

    message('empirical optimal: ', max(membership))

    return(mirt::mdirt(data = workData, model = max(membership), empiricalhist = empiricalhist, group = group, nruns = 100))
  } else {
    return(mirt::mdirt(data = workData, model = nfact[bestModel], empiricalhist = empiricalhist, group = group, nruns = 100))
  }
}

doMLCA <- function(data = ..., startN = 1, empiricalhist = F, group = NULL){
  mirtCluster()
  if(is.data.frame(data) | is.matrix(data)){
    workModel <- findMLCA(data = data, startN = startN, empiricalhist = T, empiricaloptimal = F, group = group)
  } else {
    workModel <- data
  }
  initData <- workModel@Data$data
  workData <- workModel@Data$data

  STOP <- FALSE
  while(!STOP){

    # item fit evaluation
    workModelFit <- mirt::itemfit(workModel, impute = 100, QMC = T)
    FitSize <- workModelFit$S_X2/workModelFit$df.S_X2

    print(cbind(workModelFit, FitSize))

    if(sum(na.omit(workModelFit$S_X2[1:ncol(workData)] == "NaN")) != 0){
      workModel <- findMLCA(workData[,-which(workModelFit$S_X2[1:ncol(workData)] == "NaN")], empiricalhist = F, empiricaloptimal = T)
    } else if(sum(workModelFit$S_X2[1:ncol(workData)]/workModelFit$df.S_X2[1:ncol(workData)] > 10) != 0){
      workModel <- findMLCA(workData[,-which(max(workModelFit$S_X2[1:ncol(workData)]/workModelFit$df.S_X2[1:ncol(workData)]) == workModelFit$S_X2[1:ncol(workData)]/workModelFit$df.S_X2[1:ncol(workData)])], empiricalhist = F, empiricaloptimal = T)
    } else {
      STOP <- TRUE
    }
    rm(workData)
    workData <- initData[,colnames(workModel@Data$data)] # update work data
  }

  if(sum(is.na(initData)) != 0){
    workData <- mirt::imputeMissing(workModel, fscores(workModel, QMC = T))
    workModel <- findMLCA(workData, empiricalhist = F, empiricaloptimal = T)
  }

  print(plot(workModel, facet_items = FALSE))
  print(plot(workModel))

  fs <- fscores(workModel, QMC = T)

  class_prob <- data.frame(apply(fs, 1, function(x) sample(1:workModel@Model$nfact, 1, prob=x)))
  colnames(class_prob) <- "Class"
  mirtCluster(remove = T)

  return(class_prob)
}

deepFAengine <- function(mirtModel, survey.weights){ # for search more factors with prevent local optimal
  DICindices <- vector()
  DICindices[1] <- mirtModel@Fit$DIC


  j <- 1


  if(mirtModel@Options$method == 'EM' && length(attr(mirtModel@ParObjects$lrPars, 'formula')[[1]]) == 0 && length(survey.weights) == 0) {
    method <- 'MHRM'
    optimizer <- 'NR1'
    SE.type <- 'MHRM'
    NCYCLES <- 4000

  } else if(length(survey.weights) != 0) { # if survey weight included
    method <- 'QMCEM'
    optimizer <- mirtModel@Options$Moptim
    SE.type <- mirtModel@Options$SE.type
    NCYCLES <- mirtModel@Options$NCYCLES
  } else if(mirtModel@Options$method == 'EM' && length(attr(mirtModel@ParObjects$lrPars, 'formula')[[1]]) != 0) { # if latent regression included
    method <- 'QMCEM'
    optimizer <- mirtModel@Options$Moptim
    SE.type <- mirtModel@Options$SE.type
    NCYCLES <- mirtModel@Options$NCYCLES

  } else {
    method <- mirtModel@Options$method
    optimizer <- mirtModel@Options$Moptim
    SE.type <- mirtModel@Options$SE.type
    NCYCLES <- mirtModel@Options$NCYCLES

  }

  message('searching global optimal... / estimation method: ', method, ' / optimizer: ', optimizer)
  start <- mirtModel@Model$nfact + 1
  end <- mirtModel@Model$nfact + 3 # see http://www.tandfonline.com/doi/abs/10.1080/00273171.2012.710386

  if(start > ncol(mirtModel@Data$data)){
    return(mirtModel)
  }

  if(end > ncol(mirtModel@Data$data)){
    end <- ncol(mirtModel@Data$data)
  }

  nfact <- vector()
  nfact[1] <- mirtModel@Model$nfact
  for(i in start:end){
    try(invisible(gc()), silent = T)



    if(class(mirtModel)[1] == 'MixedClass'){
      try(invisible(tempModel <- mirt::mixedmirt(data = mirtModel@Data$data, model = i, itemtype = mirtModel@Model$itemtype, SE = F, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), fixed = mirtModel@Model$formulas$fixed, random = mirtModel@Model$formulas$random, lr.fixed = mirtModel@Model$formulas$lr.fixed, lr.random = mirtModel@Model$formulas$lr.random, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES))), silent = T)

    } else {
      try(invisible(tempModel <- mirt::mirt(data = mirtModel@Data$data, model = i, itemtype = mirtModel@Model$itemtype, SE = F, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), formula = attr(mirtModel@ParObjects$lrPars, 'formula')[[1]], method = method, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES))), silent = T)
    }


    if(exists('tempModel')){
      if(tempModel@OptimInfo$converged){ #       if(tempModel@OptimInfo$converged && tempModel@OptimInfo$secondordertest == T){
        message(i, ' factors were converged; DIC: ', tempModel@Fit$DIC)
        j <- j + 1
        DICindices[j] <- tempModel@Fit$DIC
        nfact[j] <- i
        rm(tempModel)
      }
    }
  }
  bestModel <- which(min(DICindices) == DICindices)
  message('find global optimal: ', nfact[bestModel])
  if(bestModel == 1){
    return(mirtModel)
  } else {

    if(class(mirtModel)[1] == 'MixedClass'){
      try(invisible(tempModel <- mirt::mixedmirt(data = mirtModel@Data$data, model = nfact[bestModel], itemtype = mirtModel@Model$itemtype, SE = F, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), fixed = mirtModel@Model$formulas$fixed, random = mirtModel@Model$formulas$random, lr.fixed = mirtModel@Model$formulas$lr.fixed, lr.random = mirtModel@Model$formulas$lr.random, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES))), silent = T)

    } else {
      return(mirt::mirt(data = mirtModel@Data$data, model = nfact[bestModel], itemtype = mirtModel@Model$itemtype, SE = T, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), formula = attr(mirtModel@ParObjects$lrPars, 'formula')[[1]], method = method, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES)))

    }
  }
}

deepFA <- function(mirtModel, survey.weights = NULL){
  init_nfact <- mirtModel@Model$nfact
  if(init_nfact >= ncol(mirtModel@Data$data)){
    return(mirtModel)
  }

  deepModel <- deepFAengine(mirtModel, survey.weights = survey.weights)
  if(deepModel@Model$nfact == init_nfact+3){

    deepModel <- deepFAengine(deepModel, survey.weights = survey.weights)
    if(deepModel@Model$nfact == init_nfact+6){
      deepModel <- deepFAengine(deepModel, survey.weights = survey.weights)

      if(deepModel@Model$nfact == init_nfact+9){
        deepModel <- deepFAengine(deepModel, survey.weights = survey.weights)

      } else {
        return(deepModel)
      }

    } else {
      return(deepModel)
    }


  } else {
    return(deepModel)
  }
}

jmleRaschEst <- function(data, model = 'PCM', fitCriteria = 2, as.LCA = F){
  if(!require('mixRasch')){
    install.packages('mixRasch', repos = 'http://cran.nexr.com', dependencies = T); library('mixRasch')
  }
  if(!require('psych')){
    install.packages('psych', repos = 'http://cran.nexr.com', dependencies = T); library('mixRasch')
  }
  itemList <- colnames(data)
  RaschData <- data.frame(data)
  RaschData <- RaschData[psych::describe(RaschData)$range != 0]
  try(invisible(jmleRasch <- mixRasch::mixRasch(data = RaschData, steps = max(psych::describe(RaschData)$range), max.iter = 500, conv.crit = 0.00001, model = model, maxrange = c(-6, 6), maxchange = .25, as.LCA = as.LCA)), silent = T)
  STOP <- FALSE
  i <- 0
  j <- Sys.time()
  while (!STOP) {
    i <- i+1
    if(ncol(data) > 2){
      jmleFit <- data.frame(jmleRasch$item.par$in.out)
      if(length(which(abs(jmleFit$in.Z) > fitCriteria)) != 0){
        message('infit: ', paste(colnames(RaschData)[which(max(abs(jmleFit$in.Z)) == abs(jmleFit$in.Z))]), ' // iteration: ', i)

        itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(jmleFit$in.Z)) == abs(jmleFit$in.Z))]
        RaschData <- RaschData[!itemList]

        itemList <- itemList[-c(which(max(abs(jmleFit$out.Z)) == jmleFit$out.Z))]
        try(invisible(jmleRasch <- mixRasch::mixRasch(data = RaschData, steps = max(psych::describe(RaschData)$range), max.iter = 500, conv.crit = 0.00001, model = model, maxrange = c(-6, 6), maxchange = .25, as.LCA = as.LCA)), silent = T)

      } else if(length(which(abs(jmleFit$out.Z) > fitCriteria)) != 0){
        message('outfit: ', paste(colnames(RaschData)[which(max(abs(jmleFit$out.Z)) == abs(jmleFit$out.Z))]), ' // iteration: ', i)

        itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(jmleFit$out.Z)) == abs(jmleFit$out.Z))]
        RaschData <- RaschData[!itemList]

        itemList <- itemList[-c(which(max(abs(jmleFit$out.Z)) == jmleFit$out.Z))]
        try(invisible(jmleRasch <- mixRasch::mixRasch(data = RaschData, steps = max(psych::describe(RaschData)$range), max.iter = 500, conv.crit = 0.00001, model = model, maxrange = c(-6, 6), maxchange = .25, as.LCA = as.LCA)), silent = T)

      } else {
        k <- Sys.time()
        STOP <- TRUE
      }
    } else {
      k <- Sys.time()
      STOP <- TRUE
    }
  }
  RaschJMLEresult <- new.env()
  RaschJMLEresult$itemList <- colnames(RaschData)
  RaschJMLEresult$rawData <- RaschData
  RaschJMLEresult$model <- jmleRasch
  RaschJMLEresult$TotalTime <- k-j
  return(RaschJMLEresult)
}


cmleRaschEst <- function(data, model = 'PCM'){
  if(!require('eRm')){
    install.packages('eRm', repos = 'http://cran.nexr.com', dependencies = T); library('eRm')
  }
  itemList <- colnames(data)
  RaschData <- data.frame(data)
  if(model == 'PCM'){
    cmleRasch <- eRm::PCM(RaschData)
  } else if(model == 'RSM'){
    cmleRasch <- eRm::RSM(RaschData)
  } else if(model == 'dich'){
    cmleRasch <- eRm::RM(RaschData)
  } else {
    stop('model have to define PCM, RSM, dich')
  }

  # define while loop
  STOP <- FALSE # initial is STOP = FALSE
  while (!STOP) {
    if(ncol(data) > 2){
      # cmleFit <- data.frame(cmleRasch$item.par$in.out)
      cmleFit <- data.frame(eRm::itemfit(eRm::person.parameter(cmleRasch))$i.outfitZ,
                            eRm::itemfit(eRm::person.parameter(cmleRasch))$i.outfitMSQ,
                            eRm::itemfit(eRm::person.parameter(cmleRasch))$i.infitZ,
                            eRm::itemfit(eRm::person.parameter(cmleRasch))$i.infitMSQ)
      colnames(cmleFit) <- c("out.Z", "out.MSQ", "in.Z", "in.MSQ")

      print(cmleFit)

      if(length(which(abs(cmleFit$out.Z) > 2)) != 0){
        message('outfit: ', colnames(RaschData)[which(max(abs(cmleFit$out.Z)) == abs(cmleFit$out.Z))])

        itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(cmleFit$out.Z)) == abs(cmleFit$out.Z))]
        RaschData <- RaschData[!itemList]

        itemList <- itemList[-c(which(max(abs(cmleFit$out.Z)) == cmleFit$out.Z))]

        if(model == 'PCM'){
          cmleRasch <- eRm::PCM(RaschData)
        } else if(model == 'RSM'){
          cmleRasch <- eRm::RSM(RaschData)
        } else if(model == 'dich'){
          cmleRasch <- eRm::RM(RaschData)
        } else {
          stop('model have to define PCM, RSM, dich')
        }

      } else if(length(which(abs(cmleFit$in.Z) > 2)) != 0){
        message('infit: ', colnames(RaschData)[which(max(abs(cmleFit$in.Z)) == abs(cmleFit$in.Z))])

        itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(cmleFit$in.Z)) == abs(cmleFit$in.Z))]
        RaschData <- RaschData[!itemList]

        itemList <- itemList[-c(which(max(abs(cmleFit$out.Z)) == cmleFit$out.Z))]

        if(model == 'PCM'){
          cmleRasch <- eRm::PCM(RaschData)
        } else if(model == 'RSM'){
          cmleRasch <- eRm::RSM(RaschData)
        } else if(model == 'dich'){
          cmleRasch <- eRm::RM(RaschData)
        } else {
          stop('model have to define PCM, RSM, dich')
        }
      } else {
        STOP <- TRUE
      }
    } else {
      STOP <- TRUE
    }
  }
  Raschcmleresult <- new.env()
  Raschcmleresult$itemList <- colnames(RaschData)
  Raschcmleresult$model <- cmleRasch
  return(Raschcmleresult)
}

autoMCMC2PL.ml <- function(x = NULL, group = NULL, est.b.M="h", est.b.Var="i",
                           est.a.M="h" , est.a.Var="i", burnin = 10000,
                           iter = 20000, Rhat = 1.05, autofix = T, TargetTestLength = 3, considerSigma = T, # for 2PNO Multilevel
                           testlets = NULL, survey.weights = NULL, est.slope = T, est.guess = T, verbose = F # for 1-3PNO testlet
){
  if(!require('sirt')){
    install.packages('sirt')
    library('sirt')
  }
  if(!require('plyr')){
    install.packages('plyr')
    library('plyr')
  }
  if(!require('progress')){
    install.packages('progress')
    library('progress')
  }

  if(range(x[1], na.rm = T)[2] - range(x[1], na.rm = T)[1] > 1){
    link <- 'normal'
  } else {
    link <- 'logit'
  }

  if(burnin > iter){
    iter <- burnin*2
  }

  if(verbose){
    num.progress.iter <- burnin/10
  } else {
    num.progress.iter <- burnin + iter + 1
  }

  initData <- x
  iterationTrials <- 1
  StartTime <- Sys.time()

  if(length(group) != 0){
    initData <- initData[which(!is.na(group)),]
    group <- group[which(!is.na(group))]
  } else {

  }

  if(is.character(group)){
    group <- as.factor(group)
  }

  if(is.factor(group)){
    numberOfGroups <- length(attributes(group)$levels)
  } else {
    numberOfGroups <- length(attributes(as.factor(group))$levels)
  }

  message('MCMC link: ', link)
  message('MCMC Trials: ', iterationTrials)
  message('Current number of items: ', ncol(initData))
  message('Current number of groups: ', numberOfGroups)
  message('Rhat cutoff: ', Rhat)


  group <- as.integer(group)

  if(length(testlets) != 0){
    ActualTestlets <- testlets
    TestletStatus <- (table(ActualTestlets))

    if(sum(is.na(ActualTestlets)) == 0){
      ActualTestlets <- plyr::mapvalues(ActualTestlets, as.numeric(names(TestletStatus)[which(TestletStatus == max(TestletStatus))]), NA) # CTC(M-1)
    }

    TestletStatus <- (table(ActualTestlets)) # renew
    if(length(which(TestletStatus == 1)) != 0){
      ActualTestlets <- plyr::mapvalues(ActualTestlets, names(TestletStatus)[which(TestletStatus == 1)], rep(NA, length(which(TestletStatus == 1))))
    }
    # if(sum(ActualTestlets %in% 1) == 0){
    #   ActualTestlets <- ActualTestlets - min(ActualTestlets, na.rm = T) + 1
    # }
    if(sum(is.na(ActualTestlets)) == length(ActualTestlets)){
      ActualTestlets <- NULL # disable testlet when all testlet information were NA

    } else {
      ActualTestlets <- plyr::mapvalues(ActualTestlets, as.numeric(attributes(as.factor(ActualTestlets))$levels), seq(length(attributes(as.factor(c(ActualTestlets)))$levels))) # convert text, keep same interval

    }
  } else {
    ActualTestlets <- NULL
  }

  message('initializing model...')

  if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
    if(exists('init')){
      try(rm(init))
    }
    try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
    if(!exists('init')){
      while(!exists('init')){
        try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
      }
    }
  } else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
    if(exists('init')){
      try(rm(init))
    }
    try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
    if(!exists('init')){
      while(!exists('init')){
        try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
      }
    }
  } else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'normal'){
    init <- surveyFA(data = initData, survey.weights = survey.weights, testlets = ActualTestlets)
    return(init)
  } else {
    if(exists('init')){
      try(rm(init))
    }
    try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
    if(!exists('init')){
      while(!exists('init')){
        try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
      }
    }
  }

  message('starting calibration...')
  pb <- progress_bar$new(
    format = "  processing MCMC chains [:bar] :percent (:current of :total) in :elapsed, eta: :eta",
    total = ncol(initData) - TargetTestLength, clear = FALSE, width= 120, force = T)


  STOP <- FALSE

  while(!STOP){
    ticktockUpdate <- 1-((sum(init$summary.mcmcobj$Rhat > Rhat) +
                            sum(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)] < 0) +
                            sum(cbind( init$summary.mcmcobj$Q2.5[grep("^a",init$summary.mcmcobj$parameter)] < 0 && init$summary.mcmcobj$Q97.5[grep("^a",init$summary.mcmcobj$parameter)] > 0 ))
    ) / length(init$summary.mcmcobj$Rhat))
    if(!is.na(ticktockUpdate)){
      pb$update(ticktockUpdate)
    }

    pb$tick()
    invisible(try(gc(), silent = T))
    if(sum(init$summary.mcmcobj$Rhat > Rhat) != 0){
      parmList <- init$summary.mcmcobj

      # if(length(grep("^sigma[0-9]", as.character(init$summary.mcmcobj$parameter))) != 0){
      # parmList <- parmList[-grep("^sigma[0-9]", as.character(parmList$parameter)),]
      # }
      if(!considerSigma){
        parmList <- parmList[-grep("^sigma", as.character(parmList$parameter)),]
      }
      excludeVar <- unique(na.omit(as.numeric(unlist(strsplit(unlist(as.character(parmList[which(max(parmList$Rhat) == parmList$Rhat),]$parameter)), "[^0-9]+")))))

      if(length(excludeVar) != 0 && ncol(initData) > TargetTestLength){
        if(verbose) message('Removing a item ', names(initData[excludeVar]),' / ', parmList[which(max(parmList$Rhat) == parmList$Rhat),]$parameter, ' Rhat: ', max(parmList$Rhat))

        initData <- initData[,-excludeVar]
        if(length(testlets) != 0){
          ActualTestlets <- ActualTestlets[-excludeVar]
          ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))

        }


        iterationTrials <- iterationTrials+1


        if(verbose) message('MCMC Trials: ', iterationTrials)
        if(verbose) message('Current number of items: ', ncol(initData))

        if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        } else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        } else {
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        }
      } else {
        pb$update(1)
        pb$tick()
        STOP <- TRUE
      }
    } else if(sum(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)] < 0) != 0 && autofix){
      excludeVar <- unique(na.omit(as.numeric(unlist(strsplit(unlist(as.character(init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))])), "[^0-9]+")))))
      if(length(excludeVar) != 0 && ncol(initData) > TargetTestLength){
        if(verbose) message('Removing a item ', names(initData[excludeVar]),' / ', init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))], ' value: ', min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]))

        initData <- initData[,-excludeVar]
        if(length(testlets) != 0){
          ActualTestlets <- ActualTestlets[-excludeVar]
          ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))

        }
        iterationTrials <- iterationTrials+1


        if(verbose) message('MCMC Trials: ', iterationTrials)
        if(verbose) message('Current number of items: ', ncol(initData))

        if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        } else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        } else {
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        }
      } else {
        pb$update(1)
        pb$tick()
        STOP <- TRUE
      }
    } else if(
      sum(cbind( init$summary.mcmcobj$Q2.5[grep("^a",init$summary.mcmcobj$parameter)] < 0 && init$summary.mcmcobj$Q97.5[grep("^a",init$summary.mcmcobj$parameter)] > 0 )) > 0
      && autofix){
      excludeVar <- unique(na.omit(as.numeric(unlist(strsplit(unlist(as.character(init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))])), "[^0-9]+")))))
      if(length(excludeVar) != 0 && ncol(initData) > TargetTestLength){
        if(verbose) message('Removing a item ', names(initData[excludeVar]),' / ', init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))], ' value: ', min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]))

        initData <- initData[,-excludeVar]
        if(length(testlets) != 0){
          ActualTestlets <- ActualTestlets[-excludeVar]
          ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))

        }
        iterationTrials <- iterationTrials+1


        if(verbose) message('MCMC Trials: ', iterationTrials)
        if(verbose) message('Current number of items: ', ncol(initData))

        if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        } else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        } else {
          if(exists('init')){
            try(rm(init))
          }
          try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
          if(!exists('init')){
            while(!exists('init')){
              try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
            }
          }
        }
      } else {
        pb$update(1)
        pb$tick()
        STOP <- TRUE
      }

    } else { # overall stop rule when all items normally terminated
      pb$update(1)
      pb$tick()
      STOP <- TRUE
    }
  }
  EndTime <- Sys.time()

  ReturnList <- list(Solution = init, testlet = ActualTestlets, rawData = initData, TotalTime = EndTime - StartTime)

  return(ReturnList)
}

testAssembly <- function(MIRTmodel, measurementArea, NumberOfForms = 1, meanOfdifficulty = 0, sdOfdifficulty = 2.0, numberOfItems = 16, maximumItemSelection = 1, oldFormYMIRTmodel = NULL, SCLmethod = 'Haebara', oldFormYCommonItemNumber = NULL, newFormXCommonItemNumber = NULL){
  if(!require('xxIRT')){
    install.packages('xxIRT')
    library('xxIRT')
  }
  if(!require('mirt')){
    install.packages('mirt')
    library('mirt')
  }
  if(!require('sirt')){
    install.packages('sirt')
    library('sirt')
  }
  library(xxIRT)
  library(mirt)
  library('sirt')

  if(class(MIRTmodel)[1] == 'SingleGroupClass'){
    message('mirt model provided')

    IRTpars <- data.frame(coef(MIRTmodel, IRTpars = T, simplify = T)$items[,1:3])
    colnames(IRTpars)[3] <- 'c'
  } else if(class(MIRTmodel)[1] == 'mcmc.sirt'){
    message('sirt model provided')
    if(length(grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]) == 0){
      IRTpars <- data.frame(rep(1,ncol(MIRTmodel$dat)),
                            MIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]])
    } else {
      IRTpars <- data.frame(MIRTmodel$summary.mcmcobj$MAP[grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]],
                            MIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]])
    }

    colnames(IRTpars) <- c('a', 'b')
    rownames(IRTpars) <- colnames(MIRTmodel$dat)
    if(length(MIRTmodel$summary.mcmcobj$MAP[grep("^c", MIRTmodel$summary.mcmcobj$parameter)]) != 0){
      IRTpars$c <- MIRTmodel$summary.mcmcobj$MAP[grep("^c", MIRTmodel$summary.mcmcobj$parameter)]
    } else {
      IRTpars$c <- rep(0, ncol(MIRTmodel$dat))
    }

  } else {
    stop('Did you provide right mirt or sirt mcmc model?')
  }

  if(length(oldFormYMIRTmodel) != 0){ # if SCL activated
    if(length(oldFormYCommonItemNumber) == 0 | length(newFormXCommonItemNumber) == 0 | length(newFormXCommonItemNumber) != length(oldFormYCommonItemNumber)) {
      STOP('Please provide correspond oldFormYCommonItemNumber and newFormXCommonItemNumber')
    }

    # check SCL software
    if(!require(SNSequate)){
      install.packages('SNSequate')
      library('SNSequate')
    }

    # read IRT parameters of oldform Y
    if(class(oldFormYMIRTmodel)[1] == 'SingleGroupClass'){
      message('mirt model provided')

      IRTpars2 <- data.frame(coef(oldFormYMIRTmodel, IRTpars = T, simplify = T)$items[,1:3])
      colnames(IRTpars2)[3] <- 'c'
    } else if(class(oldFormYMIRTmodel)[1] == 'mcmc.sirt'){
      message('sirt model provided')

      if(length(grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]) == 0){
        IRTpars2 <- data.frame(rep(1,ncol(oldFormYMIRTmodel$dat)),
                               oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]])
      } else {
        IRTpars2 <- data.frame(oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]],
                               oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]])
      }

      colnames(IRTpars2) <- c('a', 'b')
      rownames(IRTpars2) <- colnames(oldFormYMIRTmodel$dat)
      if(length(oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^c", oldFormYMIRTmodel$summary.mcmcobj$parameter)]) != 0){
        IRTpars2$c <- oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^c", oldFormYMIRTmodel$summary.mcmcobj$parameter)]
      } else {
        IRTpars2$c <- rep(0, ncol(oldFormYMIRTmodel$dat))
      }

    } else {
      stop('Did you provide right mirt or sirt mcmc model?')
    }

    # split common item parameters
    newCommonXIRTpars <- IRTpars[newFormXCommonItemNumber,]
    oldCommonYIRTpars <- IRTpars2[oldFormYCommonItemNumber,]

    # do calcluate linking coefficients slope A and constant B
    LinkingCoefficients <- irt.link(as.data.frame(cbind(oldCommonYIRTpars, newCommonXIRTpars)), common = 1:length(oldFormYCommonItemNumber), model = '3PL', icc = 'logistic', D = 1.702)
    names(LinkingCoefficients$Haebara) <- c('A', 'B')
    names(LinkingCoefficients$StockLord) <- c('A', 'B')

    if(SCLmethod == 'Haebara'){
      message('applying Haebara method...')
      A <- LinkingCoefficients$Haebara[1]
      B <- LinkingCoefficients$Haebara[2]
    } else if(SCLmethod == 'StockingLord'){
      message('applying Stocking-Lord method...')
      A <- LinkingCoefficients$StockLord[1]
      B <- LinkingCoefficients$StockLord[2]
    } else {
      stop('please provide "Haebara" or "StockingLord" in SCLmethod argument correctly')
    }

    # adjusting new test IRT parameters to value of old test item parameters
    IRTpars$a <- IRTpars$a/A
    IRTpars$b <- A*IRTpars$b+B
    # ignoring cx=cy, that's not require of adjustments.
  }

  items <- IRTpars
  items$content <- as.numeric(as.factor(measurementArea))

  # print(items)

  # try -2 to 2 (flat information)
  x <- xxIRT::ata(items, NumberOfForms, len = numberOfItems, maxselect = maximumItemSelection) %>%
    xxIRT::ata_obj_relative(seq(min(items$b), max(items$b), .1), "max", flatten=0.1, negative = T) %>%
    xxIRT::ata_solve()
  try(ataPlot <- plot(x), silent = T)

  if(exists('ataPlot')){
    print(plot(x))

  } else {

    # -1 to 1 (flat information)
    x <- xxIRT::ata(items, NumberOfForms, len = numberOfItems, maxselect = maximumItemSelection) %>%
      xxIRT::ata_obj_relative(seq((min(items$b)+0.5), (max(items$b)-0.5), .1), "max", flatten=0.1, negative = T) %>%
      xxIRT::ata_solve()
    try(ataPlot <- plot(x), silent = T)

    if(exists('ataPlot')){
      print(plot(x))

    } else {

      rm(x)
      x <- xxIRT::ata(items, NumberOfForms, len = numberOfItems, maxselect = maximumItemSelection)
      x <- xxIRT::ata_obj_absolute(x, "b", meanOfdifficulty * numberOfItems)
      x <- xxIRT::ata_obj_absolute(x, (x$pool$b - meanOfdifficulty)^2, sdOfdifficulty * numberOfItems)
      x <- xxIRT::ata_solve(x)

      print(plot(x))
    }


  }
  y <- xxIRT::ata_get_items(x, as.list=TRUE) ## FIXME

  # get ATA data
  ATAFormData <- list()
  ATAFormModel <- list()
  ATAFormModelValues <- list()
  for(i in 1:NROW(y)){
    message('-------------------------------------------------')
    message('mean of difficulty of Form ', i,': ', mean(y[[i]]$b))
    message('sd of difficulty of Form ', i,': ', sd(y[[i]]$b))
    message('min of difficulty of Form ', i,': ', min(y[[i]]$b))
    message('max of difficulty of Form ', i,': ', max(y[[i]]$b))
    message('-------------------------------------------------\n')
  }

  for(i in 1:NROW(y)){
    if(class(MIRTmodel)[1] == 'SingleGroupClass'){
      ATAFormData[[i]] <- data.frame(MIRTmodel@Data$data)
      ATAFormData[[i]] <- ATAFormData[[i]][,colnames(ATAFormData[[i]]) %in% rownames(y[[i]])]
    } else if(class(MIRTmodel)[1] == 'mcmc.sirt'){
      ATAFormData[[i]] <- data.frame(MIRTmodel$dat)
      ATAFormData[[i]] <- ATAFormData[[i]][colnames(ATAFormData[[i]]) %in% rownames(y[[i]])]
    }

    if(ncol(ATAFormData[[i]]) != 0 && sum(psych::describe(ATAFormData[[i]])$range == 1) == ncol(ATAFormData[[i]])){
      ATAFormModelValues[[i]] <- mirt::mirt(data = ATAFormData[[i]], model = 1, itemtype = '3PL', pars = 'values')
      ATAFormModelValues[[i]][which(ATAFormModelValues[[i]]$name == 'a1'),"value"] <- y[[i]]$a
      ATAFormModelValues[[i]][which(ATAFormModelValues[[i]]$name == 'd'),"value"] <- -1*y[[i]]$a*y[[i]]$b
      ATAFormModelValues[[i]][which(ATAFormModelValues[[i]]$name == 'g'),"value"] <- y[[i]]$c

      ATAFormModelValues[[i]]$est <- FALSE

      ATAFormModel[[i]] <- mirt::mirt(data = ATAFormData[[i]], model = 1, itemtype = '3PL', pars = ATAFormModelValues[[i]])
      names(ATAFormModel)[[i]] <- paste0('form', i)
    } else {

      # if(){
      message('warning: items were polytomous? or variable names are contain spaces?') ## FIXME: ADD 3PLNRM

      # } else {
      # message('warning: ')
      # }

      # NULL
    }

  }



  z <- list(OriginalParms = items, ATAforms = y, ATAFormModel = ATAFormModel, ATAFormSolutions = x)

  return(z)
}

fastBifactorCFA <- function(x, ga = F, itemkeys = NULL, initSolution = F, skipNRM = T, excludeUnscalableVar = F, lowerbound = .5, covdata = NULL, formula = NULL){
  if(!require('mokken')){
    install.packages('mokken')
    library('mokken')
  }
  if(!require('mirt')){
    install.packages('mirt')
    library('mirt')
  }

  x <- x[psych::describe(x)$range != 0]

  message('finding testlet structures using mokken scale analysis...')
  if(ga){
    message('mutuation probability: ',1/log(10^(log2(ncol(x)/5)) * 1000)*(log(log2(nrow(x)))+log(log(log2(nrow(x))))))
    if(length(itemkeys) != 0){
      try(modMokken <- mokken::aisp(data.frame(mirt::key2binary(data.frame(x), key = itemkeys)), lowerbound = lowerbound, verbose = T, search = 'ga', pxover = 1, pmutation = 1/log(10^(log2(ncol(x)/5)) * 1000)*(log(log2(nrow(x)))+log(log(log2(nrow(x))))), popsize = log2(nrow(x))), silent = T)

    } else {
      try(modMokken <- mokken::aisp(data.frame(x), lowerbound = lowerbound, verbose = T, search = 'ga', pxover = 1, pmutation = 1/log(10^(log2(ncol(x)/5)) * 1000)*(log(log2(nrow(x)))+log(log(log2(nrow(x))))), popsize = log2(nrow(x))), silent = T)

    }
  } else {
    if(length(itemkeys) != 0){
      try(modMokken <- mokken::aisp(data.frame(mirt::key2binary(data.frame(x), key = itemkeys)), lowerbound = lowerbound, verbose = T), silent = T)

    } else {
      try(modMokken <- mokken::aisp(data.frame(x), verbose = T, lowerbound = lowerbound), silent = T)

    }
  }

  print(modMokken)

  if(sum(modMokken == 0) == ncol(x) | sum(modMokken == 1) == ncol(x)){
    # print(modMokken)
    if(!initSolution){
      modBfactor <- surveyFA(data = data.frame(x), itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich')

    } else {
      modBfactor <- deepFA(fastFIFA(x = data.frame(x), itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich'))

    }

  } else {
    # print(modMokken)
    message('testlet structure was found!')

    if(excludeUnscalableVar){
      message('excluding unscalable varaiables...')
      testDat <- data.frame(x)[-which(modMokken == 0)]
      modMokken <- modMokken[-which(modMokken == 0)]
      modMokken[which(modMokken %in% as.numeric(names(which(table(modMokken) < 2))))] <- NA

      # print(modMokken)
    } else {
      testDat <- data.frame(x)
      modMokken[which(modMokken == 0)] <- NA
      modMokken[which(modMokken %in% as.numeric(names(which(table(modMokken) < 2))))] <- NA
      # print(modMokken)

    }
    message('number of testing variables: ', (ncol(testDat)), ' (', round(ncol(testDat)/ncol(data.frame(x))*100, digits = 2), '%)')
    if(!initSolution){
      modBfactor <- surveyFA(data = testDat, testlets = modMokken, itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich')

    } else {
      modBfactor <- fastFIFA(x = testDat, testlets = modMokken, itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich')

    }

  }

  return(modBfactor)

}

bfactor2mod <- function(model, J){ # borrow from mirt::bfactor2mod
  tmp <- tempfile('tempfile')
  unique <- sort(unique(model))
  index <- 1L:J
  tmp2 <- c()
  for(i in 1L:length(unique)){
    ind <- na.omit(index[model == unique[i]])
    comma <- rep(',', 2*length(ind))
    TF <- rep(c(TRUE,FALSE), length(ind))
    comma[TF] <- ind
    comma[length(comma)] <- ""
    tmp2 <- c(tmp2, c(paste('\nS', i, ' =', sep=''), comma))
  }
  cat(tmp2, file=tmp)
  model <- mirt::mirt.model(file=tmp, quiet = TRUE)
  unlink(tmp)
  return(model)
}

doBfactor2mod <- function(mirtDat, testlet){
  tmp <- any(sapply(colnames(mirtDat), grepl, x=paste0('G = 1-', ncol(mirtDat)))) # borrow from mirt::bfactor
  model2 <- mirt::mirt.model(paste0('G = 1-', ncol(mirtDat)), itemnames = if(tmp) colnames(mirtDat) else NULL)
  model <- bfactor2mod(testlet, length(testlet))

  model$x <- rbind(model2$x, model$x)

  return(model)
}

# item factor analysis based LCA
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

findLatentClass <- function(data = ..., nruns = 1, maxClasses = NULL, covdata = NULL, formula = NULL, SE.type = 'sandwich', checkSecondOrderTest = T, empiricalhist = T, turnOnMIRTcluster = F){
  if(!require('mirt')){
    install.packages('mirt')
    library('mirt')
  }

  if(!require('progress')){
    install.packages('progress')
    library('progress')
  }


  if(turnOnMIRTcluster){
    try(mirtCluster(remove = T))
    try(mirtCluster(spec = round(parallel::detectCores()/2)))
  }

  modelFit <- list()
  if(is.null(maxClasses)){
    maxClasses <- round(ncol(data)/5)
  }

  pb <- progress_bar$new(
    format = "(:spin) [:bar] :percent (:current of :total) in :elapsed, eta: :eta",
    total = maxClasses, clear = FALSE, width = 60)


  for(i in 1:maxClasses){
    invisible(try(testModel <- mirt::mdirt(data = data, model = i, SE = checkSecondOrderTest, verbose = F, nruns = nruns, covdata = covdata, formula = formula, SE.type = SE.type, empiricalhist = empiricalhist), silent = T))

    pb$tick()
    if(exists('testModel')){
      if(checkSecondOrderTest){

        if(testModel@OptimInfo$converged && testModel@OptimInfo$secondordertest){
          # message(round(i/maxClasses*100, 1), "% complete", ' (', i,' / ', maxClasses, ')')
          modelFit[[i]] <- testModel@Fit
        }
        rm(testModel)
      } else {

        if(testModel@OptimInfo$converged){
          # message(round(i/maxClasses*100, 1), "% complete", '(', i,' / ', maxClasses, ')')
          modelFit[[i]] <- testModel@Fit
        }
        rm(testModel)
      }

    }

  }
  # return(modelFit)

  modelFitMatrix <- (matrix(unlist(modelFit), ncol = 14, byrow = TRUE))
  modelFitNumbers <- vector()

  for(i in 1:NROW(modelFit)){
    if(is.null(modelFit[[i]]) == FALSE){
      modelFitNumbers[length(modelFitNumbers)+1] <- i
    }

  }
  rownames(modelFitMatrix) <- modelFitNumbers
  colnames(modelFitMatrix) <- names(modelFit[[1]])
  return(modelFitMatrix)
}



autoLCA <- function(data = ..., UIRT = T, nruns = 1, covdata = NULL, formula = NULL, SE.type = 'sandwich', forceMHRM = F){
  # source('https://github.com/seonghobae/k.aefa/raw/master/aFIPC.R')
  if(length(covdata) != 0 && SE.type != 'complete'){
    SE.type <- 'complete'
  }
  testMIRTmod <- surveyFA(data = data, forceMHRM = forceMHRM, forceUIRT = UIRT, covdata = covdata, formula = formula, SE.type = SE.type)
  testNumberOfClasses <- findLatentClass(data = testMIRTmod@Data$data, nruns = nruns, covdata = covdata, formula = formula, SE.type = 'complete') # sandwich not yet supported

  LCA_Judgement <- vector()
  for(j in c(7:11)){
    if(sum(is.na(testNumberOfClasses[,j])) > 0){
      try(LCA_Judgement[j] <- NA)
    } else {
      try(LCA_Judgement[j] <- (which(testNumberOfClasses[,j] == min(testNumberOfClasses[,j]))))

    }
  }

  try(FinalModel <- mirt::mdirt(testMIRTmod@Data$data, as.numeric(rownames(testNumberOfClasses)[getmode(na.omit(LCA_Judgement))]), nruns = nruns, verbose = F, SE = T, covdata = covdata, formula = formula, SE.type = 'complete'), silent = T)


  return(list(IRTmodel = testMIRTmod, LCAdecisionTable = testNumberOfClasses, FinalModel = FinalModel))
}

MLIRT <- function(data = ..., covdata = ..., model = ..., itemtype = 'Rasch', fixed = ~-1, random = NULL, lr.fixed = ~1, lr.random = NULL, SEtol = 1e-4, SE = T, TOL = .001, NCYCLES = 2000, ...){
  if(!require('mirt')){
    install.packages('mirt', repos = 'https://cran.biodisk.org')
  }

  MLM_calibration <- mirt::mixedmirt(data = data, covdata = covdata, model = model, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, itemtype = itemtype, SE = SE, SE.type = 'MHRM', TOL = TOL, technical = list(SEtol = SEtol, removeEmptyRows = TRUE, NCYCLES = NCYCLES))


  return(MLM_calibration)
  # }
}

EEMEIRT <- function(data = ..., covdata = ..., itemtype = ..., SE = T, fixed = ~-1, random = NULL, lr.fixed = ~1, lr.random = NULL, TOL = .001, SEtol = 1e-03, NCYCLES = 2000){
  testDIC <- vector()
  noConverge <- 0
  for(i in 1:ncol(data)){
    if(i > 10 && noConverge > 5){
      message(which(testDIC == min(testDIC, na.rm = T)))
      testMLM <- MLIRT(data = data, covdata = covdata, model = which(testDIC == min(testDIC, na.rm = T)),
                       itemtype = itemtype, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
                       TOL = TOL, SEtol = SEtol, NCYCLES = NCYCLES)
      return(testMLM)
    }
    try(testMLM <- MLIRT(data = data, covdata = covdata, model = i, SE = F, itemtype = itemtype, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, TOL = TOL, SEtol = SEtol, NCYCLES = NCYCLES), silent = F)
    if(exists('testMLM')){
      if(testMLM@OptimInfo$converged){
        testDIC[i] <- testMLM@Fit$DIC
        noConverge <- 0
        try(rm(testMLM))
      } else {
        testDIC[i] <- NA
        noConverge <- noConverge+1
        try(rm(testMLM))
      }
    } else {
      testDIC[i] <- NA
    }
  }
  if(sum(is.na(testDIC)) == length(testDIC)){
    stop('no solution')
  } else {
    message(which(testDIC == min(testDIC, na.rm = T)))
    testMLM <- MLIRT(data = data, covdata = covdata, model = which(testDIC == min(testDIC, na.rm = T)), itemtype = itemtype, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, TOL = TOL, SEtol = SEtol, NCYCLES = NCYCLES)
    return(testMLM)
  }

}


rotateEMEIRT <- function(EMEIRTmodel, rotate = 'bifactorQ', suppress = 0){
  if(!require('mirt')){
    install.packages('mirt')
    library('mirt')
  }
  # recursive formula
  message(EMEIRTmodel@Model$model)
  MLM_rotate_formula_mod <- mirt::mirt(data = EMEIRTmodel@Data$data, model = EMEIRTmodel@Model$model, itemtype = EMEIRTmodel@Model$itemtype, pars = 'values')
  MLM_rotate_formula_mod_original <- mirt::mod2values(EMEIRTmodel)
  if(sum(MLM_rotate_formula_mod_original$name == '(Intercept)') != 0){
    MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original[!MLM_rotate_formula_mod_original$name == '(Intercept)',]

  }
  # MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original
  MLM_rotate_formula_mod$value[which(MLM_rotate_formula_mod$item %in% colnames(EMEIRTmodel@Data$data))] <- MLM_rotate_formula_mod_original$value[which(MLM_rotate_formula_mod_original$item %in% colnames(EMEIRTmodel@Data$data))]
  MLM_rotate_formula_mod$est <- F
  # print(MLM_rotate_formula_mod)

  MLM_rotate_formula_mod_final <- mirt::mirt(data = EMEIRTmodel@Data$data, model = EMEIRTmodel@Model$model, itemtype = EMEIRTmodel@Model$itemtype, pars = MLM_rotate_formula_mod, method = 'QMCEM')
  if(ncol(MLM_rotate_formula_mod_final@Fit$F) > 1){
    if(rotate == 'bifactorQ'){
      return(GPArotation::bifactorQ(MLM_rotate_formula_mod_final@Fit$F, maxit = 1e+7))

    } else if(rotate == 'geominQ'){
      return(GPArotation::geominQ(MLM_rotate_formula_mod_final@Fit$F, maxit = 1e+7))
    } else if(rotate == 'bifactorT'){
      return(GPArotation::bifactorT(MLM_rotate_formula_mod_final@Fit$F, maxit = 1e+7))
    }else if(rotate == 'none'){
      summary(MLM_rotate_formula_mod_final, rotate = 'none', suppress = suppress)
    }

  }else{
    summary(MLM_rotate_formula_mod_final, rotate = 'none', suppress = suppress)
  }  #
}

# rotateEMEIRT(silenceMLM, rotate = 'bifactorT', suppress = .05)
doLCA <- function(data = ..., SE.type = 'Oakes', checkSecondOrderTest = T, nruns = 1, maxClasses = NULL, empiricalhist = T, covdata = NULL, formula = NULL){

  if(!require('psych')){
    install.packages('psych')
    library('psych')
  }
  if(is.data.frame(data)){
    datd <- data
  } else {
    datd <- data.frame(datd)
  }
  datd <- datd[psych::describe(datd)$range != 0] # prevent range = 0
  STOP_LCA <- FALSE
  while(!STOP_LCA){
    try(invisible(gc()), silent = T) # garbage cleaning

    preLCAoptimal <- findLatentClass(datd, SE.type = SE.type, checkSecondOrderTest = checkSecondOrderTest, nruns = nruns, maxClasses = maxClasses, empiricalhist = empiricalhist, covdata = covdata, formula = formula)
    preLCAoptimal <- data.frame(preLCAoptimal)

    optimalDecisionAIC <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$AIC == min(preLCAoptimal$AIC))]
    optimalDecisionAICc <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$AICc == min(preLCAoptimal$AICc))]
    optimalDecisionSABIC <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$SABIC == min(preLCAoptimal$SABIC))]
    optimalDecisionDIC <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$DIC == min(preLCAoptimal$DIC))]
    optimalDecisionTable <- cbind(optimalDecisionAIC, optimalDecisionAICc, optimalDecisionSABIC, optimalDecisionDIC)
    optimalDecisionTable <- table(optimalDecisionTable)

    optimalDecision <- as.numeric(names(optimalDecisionTable)[which(optimalDecisionTable == max(optimalDecisionTable))])

    preLCAmod <- mirt::mdirt(datd, optimalDecision, SE = T, SE.type = SE.type, nruns = nruns, empiricalhist = empiricalhist, covdata = covdata, formula = formula)
    preLCAmoditemFit <- data.frame(mirt::itemfit(preLCAmod))

    message('global optimal: ', optimalDecision, '')
    print(preLCAmoditemFit)
    if(sum(is.na(preLCAmoditemFit$df.S_X2)) != 0){
      datd <- datd[,!is.na(preLCAmoditemFit$df.S_X2)]
    } else if (length(which(preLCAmoditemFit$df.S_X2 == 0)) != 0 && length(which(preLCAmoditemFit$df.S_X2 == 0)) != nrow(preLCAmoditemFit)){
      datd <- datd[,-which(preLCAmoditemFit$df.S_X2 == 0)]
    } else if (length(which(preLCAmoditemFit$p.S_X2 < .05)) != 0){
      datd <- datd[,-which(preLCAmoditemFit$S_X2/preLCAmoditemFit$df.S_X2 == max(preLCAmoditemFit$S_X2/preLCAmoditemFit$df.S_X2, na.rm = T))]
    } else {
      return(list(modelFit = preLCAoptimal, LCAmodel = preLCAmod))
    }
    rm(preLCAmod)
    rm(preLCAoptimal)
    rm(preLCAmoditemFit)
  }
}


KoreanNounExtraction <- function(dat, polyReturn = F, ExtractType = 'noun'){

  if(!require('RHINO')){
    install.packages('rJava')
    install.packages('devtools')
    devtools::install_github('seonghobae/RHINO')
    library('RHINO')
  }

  if(!require('progress')){
    install.packages('progress')
    library('progress')
  }
  library('RHINO')
  invisible(.connRHINO <<- RHINO::initRhino())
  invisible(.connRHINO <- RHINO::initRhino())

  for(i in 1:ncol(dat)){
    dat[,i] <- as.character(dat[,i])
  }
  a <- vector()

  TotCells <-ncol(dat)*nrow(dat)
  k <- 0
  pb <- progress::progress_bar$new(
    format = "  extracting words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
    total = TotCells, clear = T, width= 120)

  for(i in 1:ncol(dat)){

    for(j in 1:length(dat[,i])){
      k <- k+1
      if(!is.na(dat[j,i])){

        try(pb$update(k/TotCells))
        try(pb$tick())
        dat[j,i] <- gsub("\r\n", " ", dat[j,i])
        dat[j,i] <- gsub("\r", " ", dat[j,i])
        dat[j,i] <- gsub("\n", " ", dat[j,i])
        dat[j,i] <- gsub("^\\s+|\\s+$", "", dat[j,i])
        a <- c(a, RHINO::getMorph(dat[j,i], type = ExtractType))
      } else {
        try(pb$update(k/TotCells))
        try(pb$tick())
      }


    }

  }
  a <- a[!duplicated(a)]
  a <- a[order(a)]

  b <- vector()
  for(i in 1:length(a)){
    if(is.na(a[i])){
      b[i] <- NA
    }
    else if(nchar(a[i]) > 1){
      b[i] <- a[i]
    } else {
      b[i] <- NA
    }
  }

  b <- na.omit(b)
  datTextMatrix <- data.frame(matrix(data = 0, nrow = nrow(dat), ncol = length(b)))

  colnames(datTextMatrix) <- b

  rm(pb)
  TotCells <- ncol(datTextMatrix)*ncol(dat)
  k <- 0
  pb <- progress::progress_bar$new(
    format = "  arranging words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
    total = TotCells, clear = T, width= 120)

  for(i in 1:ncol(datTextMatrix)){ # i th word

    for(j in 1:ncol(dat)){
      k <- k+1
      try(pb$update(k/TotCells))
      try(pb$tick())
      searchEngine <- grep(b[i], dat[,j]) # find a extracted word i in a raw sentense j

      if(polyReturn){
        datTextMatrix[searchEngine,i] <- datTextMatrix[searchEngine,i]+1
      } else {
        if(sum(datTextMatrix[searchEngine,i]) == 0){
          datTextMatrix[searchEngine,i] <- 1
        }
      }
    }



  }
  return (datTextMatrix)
}


KoreanExtraction <- function(dat, polyReturn = F){

  if(!require('RHINO')){
    install.packages('rJava')
    install.packages('devtools')
    devtools::install_github('seonghobae/RHINO')
    library('RHINO')
  }

  if(!require('progress')){
    install.packages('progress')
    library('progress')
  }
  library('RHINO')
  invisible(.connRHINO <<- RHINO::initRhino())
  invisible(.connRHINO <- RHINO::initRhino())

  for(i in 1:ncol(dat)){
    dat[,i] <- as.character(dat[,i])
  }
  a <- vector()

  TotCells <-ncol(dat)*nrow(dat)
  k <- 0
  pb <- progress::progress_bar$new(
    format = "  extracting words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
    total = TotCells, clear = T, width= 120)

  for(i in 1:ncol(dat)){

    for(j in 1:length(dat[,i])){
      k <- k+1
      if(!is.na(dat[j,i])){

        try(pb$update(k/TotCells))
        try(pb$tick())
        dat[j,i] <- gsub("\r\n", " ", dat[j,i])
        dat[j,i] <- gsub("\r", " ", dat[j,i])
        dat[j,i] <- gsub("\n", " ", dat[j,i])
        dat[j,i] <- gsub("^\\s+|\\s+$", "", dat[j,i])
        a <- c(a, RHINO::getMorph(dat[j,i], type = 'noun'))
        a <- c(a, RHINO::getMorph(dat[j,i], type = 'verb'))
      } else {
        try(pb$update(k/TotCells))
        try(pb$tick())
      }


    }

  }
  a <- a[!duplicated(a)]
  a <- a[order(a)]

  b <- vector()
  for(i in 1:length(a)){
    if(is.na(a[i])){
      b[i] <- NA
    }
    else if(nchar(a[i]) > 1){
      b[i] <- a[i]
    } else {
      b[i] <- NA
    }
  }

  b <- na.omit(b)
  datTextMatrix <- data.frame(matrix(data = 0, nrow = nrow(dat), ncol = length(b)))

  colnames(datTextMatrix) <- b

  rm(pb)
  TotCells <- ncol(datTextMatrix)*ncol(dat)
  k <- 0
  pb <- progress::progress_bar$new(
    format = "  arranging words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
    total = TotCells, clear = T, width= 120)

  for(i in 1:ncol(datTextMatrix)){ # i th word

    for(j in 1:ncol(dat)){
      k <- k+1
      try(pb$update(k/TotCells))
      try(pb$tick())
      searchEngine <- grep(b[i], dat[,j]) # find a extracted word i in a raw sentense j

      if(polyReturn){
        datTextMatrix[searchEngine,i] <- datTextMatrix[searchEngine,i]+1
      } else {
        if(sum(datTextMatrix[searchEngine,i]) == 0){
          datTextMatrix[searchEngine,i] <- 1
        }
      }
    }



  }
  return (datTextMatrix)
}

doLSA <- function(data, polyReturn = F, personID = NULL, SE.type = 'Oakes', checkSecondOrderTest = T, nruns = 1, maxClasses = NULL, empiricalhist = T){
  occuranceMatrix <- KoreanExtraction(data, polyReturn = polyReturn)

  if(!is.null(personID)){
    occuranceMatrix <- aggregate(occuranceMatrix, by = list(as.factor(personID)), FUN = sum)
    if(!polyReturn){
      occuranceMatrix[-1][occuranceMatrix[-1] > 0] <- 1
    }
  } else {
    occuranceMatrix <- data.frame(matrix(nrow = nrow(occuranceMatrix)), occuranceMatrix)
  }

  modLSA <- doLCA(data = occuranceMatrix[-1], SE.type = SE.type, checkSecondOrderTest = checkSecondOrderTest, nruns = nruns, maxClasses = maxClasses, empiricalhist = empiricalhist, covdata = occuranceMatrix[1])
  return(modLSA)
}

sequentialFA <- function(data, minTestLength = 3, SE.type = 'Oakes', forceUIRT = F, forceCTTmode = F, forceNormalEM = T){
  data <- data.frame(data)
  dataNames <- colnames(data)

  STOP <- FALSE
  j <- 0
  seqModels <- list()
  while(!STOP){

    if(length(dataNames) > minTestLength){

      estMod <- surveyFA(data[dataNames], SE.type = SE.type, forceUIRT = forceUIRT, forceCTTmode = forceCTTmode, minimumLeftItems = minTestLength, forceNormalEM = forceNormalEM)
      if(exists('estMod')){
        j <- j+1
        seqModels[[j]] <- estMod
        dataNames <- dataNames[!dataNames %in% colnames(estMod@Data$data)]
        rm(estMod)
      } else {
        STOP <- TRUE
      }
    } else {
      STOP <- TRUE
    }


  }
  return(seqModels)
}



fitMLIRT <- function(dat = NULL, model = 1, itemtype = 'nominal',
                     covdata, fixed = ~1, random = NULL, lr.fixed = ~1, lr.random = NULL,
                     NCYCLES = 4000, BURNIN = 1500, SEMCYCLES = 1000,
                     GenRandomPars = T, symmetric = F, accelerate = 'squarem', MHRM_SE_draws = 10000,
                     survey.weights = NULL, SEtol = 1e-4, removeEmptyRows = T){
  combine <- function (x, y) {
    combn (y, x, paste, collapse = ", ")
  }
  if(!require('progress')){
    install.packages('progress')
    library('progress')
  }

  res1 <- paste0('list(',unlist (lapply (0:NROW(random), combine, random)), ')')
  res2 <- paste0('list(',unlist (lapply (0:NROW(lr.random), combine, lr.random)), ')')
  res <- vector()
  res[1] <- length(res1)
  res[2] <- length(res2)

  if(length(res1) > length(res2)){
    res2 <- rep(res2, length(res1))
  } else {
    res1 <- rep(res1, length(res2))
  }

  DIC <- vector()

  pb <- progress_bar$new(
    format = "  estimating optimal multilevel model [:bar] :percent in :elapsed (:current of :total) eta: :eta",
    total = max(res)+2, clear = FALSE, width= 120)



  for(i in 1:max(res)){
    pb$tick()
    try(invisible(gc()))
    try(suppressMessages(invisible(modMIXED <- mirt::mixedmirt(data = dat, model = model, itemtype = itemtype,
                                                               covdata = covdata, fixed = fixed, random = eval(parse(text=res1[i])),
                                                               lr.fixed = lr.fixed, lr.random = eval(parse(text=res2[i])),
                                                               verbose = F, GenRandomPars = GenRandomPars, accelerate = accelerate,
                                                               technical = list(NCYCLES = NCYCLES, BURNIN = BURNIN, SEMCYCLES = SEMCYCLES,
                                                                                MHRM_SE_draws = MHRM_SE_draws,
                                                                                symmetric = symmetric,
                                                                                removeEmptyRows = removeEmptyRows,
                                                                                SEtol = SEtol), SE = F,
                                                               survey.weights = survey.weights))), silent = T)
    if(exists('modMIXED')){

      if(modMIXED@OptimInfo$converged){
        DIC[i] <- modMIXED@Fit$DIC
      } else {
        DIC[i] <- NA
      }

    } else {
      DIC[i] <- NA
    }
    try(invisible(rm(modMIXED)), silent = T)
  }


  if(sum(is.na(DIC)) != length(DIC)){

    pb$tick()
    try(suppressMessages(invisible(modMIXED <- mirt::mixedmirt(data = dat, model = model, itemtype = itemtype,
                                                               covdata = covdata, fixed = fixed, random = eval(parse(text=res1[which(DIC == min(DIC, na.rm = T))])),
                                                               lr.fixed = lr.fixed, lr.random = eval(parse(text=res2[which(DIC == min(DIC, na.rm = T))])),
                                                               verbose = F, GenRandomPars = GenRandomPars, accelerate = accelerate,
                                                               technical = list(NCYCLES = NCYCLES, BURNIN = BURNIN, SEMCYCLES = SEMCYCLES,
                                                                                MHRM_SE_draws = MHRM_SE_draws,
                                                                                symmetric = symmetric,
                                                                                removeEmptyRows = removeEmptyRows,
                                                                                SEtol = SEtol), SE = T,
                                                               survey.weights = survey.weights))), silent = T)
    pb$tick()
    return(modMIXED)
  } else{
    stop('No convergence')
  }

}
seonghobae/k.aefa documentation built on May 4, 2019, 1:08 p.m.