k.aefa3.R

# Kwangwoon Automated Exploratory Factor Analysis (K.AEFA)
# Seongho Bae ([email protected])

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(update.packages(ask = F, repos = 'http://cran.biodisk.org'))

if(!require(depmixS4)) {
  try(install.packages("depmixS4", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent=TRUE); library(depmixS4)
}

if(!require(Rsolnp)) {
  try(install.packages("Rsolnp", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent=TRUE); library(Rsolnp)
}

#if(!require(Cairo)) {
#  try(install.packages("Cairo", dependencies = TRUE, repos = 'http://cran.nexr.com'), silent=TRUE); library(Cairo)
#}
#if(!require(cairoDevice)) {
#  try(install.packages("cairoDevice", dependencies = TRUE, repos = 'http://cran.nexr.com'), silent=TRUE); library(cairoDevice)
#}
if(!require(stringr)) {
  try(install.packages("stringr", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent=TRUE); library(stringr)
}

if(!require(SQUAREM)) {
  try(install.packages("SQUAREM", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent=TRUE); library(SQUAREM)
}

if(!require(rrcovNA)) {
  try(install.packages("rrcovNA", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent=TRUE); library(rrcovNA)
}

# check packages
#if(!require(FAiR)) {
#  message("\n Installing Packages: Processing......")
#  try(install.packages("FAiR", dependencies = TRUE, repos = 'http://cran.nexr.com'), silent = T)
#}

# check packages
if(!require(bfa)) {
  try(install.packages("bfa", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

# check packages
if(!require(mirt)) {
  if(Sys.info()["sysname"] == "Linux" | Sys.info()["sysname"] == "Darwin"){
    try(install.packages("devtools", dependencies = TRUE), silent = T)
    try(library('devtools'), silent = T)
    try(install_github('philchalmers/mirt'), silent = T)
    try(install_github('philchalmers/mirtCAT'), silent = T)
  }
  else {
    try(install.packages("mirt", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
  }
  
}

if(!require(latticeExtra)) {
  try(install.packages('latticeExtra', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(plyr)) {
  try(install.packages('plyr', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(multilevel)) {
  try(install.packages('multilevel', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(nlme)) {
  try(install.packages('nlme', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(lsr)) {
  try(install.packages('lsr', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(meta)) {
  try(install.packages('meta', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(metafor)) {
  try(install.packages('metafor', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}

if(!require(rmeta)) {
  try(install.packages('rmeta', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(psychometric)) {
  try(install.packages('psychometric', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(pracma)) {
  try(install.packages('pracma', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(rsm)) {
  try(install.packages('rsm', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(car)) {
  try(install.packages('car', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(TAM)) {
  try(install.packages('TAM', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(GPArotation)) {
  try(install.packages('GPArotation', dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}  
if(!require(lavaan)) {
  try(install.packages("lavaan", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(semTools)) {
  try(install.packages("semTools", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}
if(!require(psych)) {
  try(install.packages("psych", dependencies = TRUE, repos = 'http://cran.biodisk.org'), silent = T)
}


#update.packages(ask = F, dependencies = TRUE, checkBuilt=TRUE)
message(" Checking Packages: OK")

message("\n Loading Packages: Processing......")
try(library(depmixS4), silent = T)
try(library(Rsolnp), silent = T)
#try(library(Cairo), silent = T)
#try(library(cairoDevice), silent = T)
try(library(stringr), silent = T)
try(library(SQUAREM), silent = T)
try(library(psychometric), silent = T)
try(library(psych), silent = T)
try(library(plyr), silent = T)
#try(library(FAiR), silent = T)
try(library(bfa), silent = T)
try(library(mirt), silent = T)
try(library(latticeExtra), silent = T)
try(library(pracma), silent = T)
try(library(multilevel), silent = T)
try(library(nlme), silent = T)
try(library(lsr), silent = T)
try(library(rsm), silent = T)
try(library(car), silent = T)
try(library(TAM), silent = T)
try(library(GPArotation), silent = T)
try(library(lavaan), silent = T)
try(library(semTools), silent = T)



# 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           #
  # [email protected]      #
  # 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([email protected]$data[,rownames(subset([email protected]$F, [email protected]$F >= min(sort([email protected]$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){ # [email protected]$AICc > [email protected]$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])]))</