# Kwangwoon Automated Exploratory Factor Analysis (K.AEFA)
# Seongho Bae (seongho@kw.ac.kr)
options(mc.cores = parallel::detectCores())
# for mac
if(Sys.getenv("LANG") == "ko-Kore_KR.UTF-8"){
Sys.setenv(LANG="ko_KR.UTF-8")
Sys.setlocale('LC_ALL', 'ko_KR.UTF-8')
Sys.setlocale('LC_MESSAGES', 'ko_KR.UTF-8')
Sys.setlocale('LC_CTYPE', 'ko_KR.UTF-8')
Sys.setlocale('LC_COLLATE', 'ko_KR.UTF-8')
Sys.setlocale('LC_TIME', 'ko_KR.UTF-8')
Sys.setlocale('LC_MONETARY', 'ko_KR.UTF-8')
}
##############
# aefa frontend #
##############
message("\n Check Packages: Processing......")
# check packages
try(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 #
# seongho@kw.ac.kr #
# Nov 6th 2014 #
# Under GNU / GPL #
#########################
# checking packages for multiple impuation in SEM context
if(!require('semTools')) {
try(install.packages('semTools', dependencies = TRUE), silent=TRUE)
}
if(!require('Amelia')) {
try(install.packages('Amelia', dependencies = TRUE), silent=TRUE)
}
if(!require('lavaan')) {
try(install.packages('lavaan', dependencies = TRUE), silent=TRUE)
}
# loading packages for multiple imputation in SEM context
library('semTools')
library('Amelia')
library('lavaan')
fit <- sem(model=model, data=data.frame(data)) # extract variable names in model syntax
message("sample size (listwise): ", paste0(nrow(data.frame(fit@Data@X))))
fit_MI <- runMI(model=model, data=data.frame(data[,attributes(fit)$Model@dimNames[[1]][[1]]]), m=m, fun = fun, ordered=names(data[,attributes(fit)$Model@dimNames[[1]][[1]]]), miArgs=list(ords = attributes(fit)$Model@dimNames[[1]][[1]]), estimator = estimator, chi = chi, ...)#, control=list(optim.method="L-BFGS-B"), ...)
cat(summary(fit_MI, standardize=T))
print(inspect(fit_MI, 'fit'))
Sys.sleep(60)
fit_data <- data.frame(fit_MI@Data@X)
colnames(fit_data) <- attributes(fit_MI)$Model@dimNames[[1]][[1]]
message("sample size (imputated): ", paste0(nrow(data.frame(fit_MI@Data@X))))
if(fun == 'cfa' | fun == 'CFA') {
fit_MI_theta <- cfa(model=model, data=fit_data, ordered=names(fit_data), estimator = estimator, ...) #, control=list(optim.method="L-BFGS-B")
} else if(fun == 'sem' | fun == 'SEM') {
fit_MI_theta <- sem(model=model, data=fit_data, ordered=names(fit_data), estimator = estimator, ...) #, control=list(optim.method="L-BFGS-B")
} else if(fun == 'growth' | fun == 'GROWTH') {
fit_MI_theta <- growth(model=model, data=fit_data, ordered=names(fit_data), estimator = estimator, ...) #, control=list(optim.method="L-BFGS-B")
}
# return(fit_MI)
return(fit_MI_theta)
# fit_cfaMI <- likertMI(model=model.sem, data=rawdata, m=100, estimator='WLSMV', fun='cfa', verbose=T)
# fit_semMI <- likertMI(model=model.sem, data=rawdata, m=100, estimator='WLSMV', fun='sem', verbose=T)
}
# Full-automatically information item factor analysis until investigating simple factor structure
likertFA <- function(data = ..., ...) {
fa_covdata <- data.frame(data)
# fa_covdata <- k.imputation(fa_covdata, ...) # data imputation
message('\nStage 1 of calcuation: Discriminant coefficient based Evaluation')
ncol_stage1 <- ncol(fa_covdata)
message('\nCurrent number of Items: ', paste0(ncol_stage1))
result <- k.aefa(fa_covdata, estimator='mirt', ...)
# test <- data.frame(coef(result))
# b <- abs(test[1, grep("a[0-9]$", colnames(test))]) >= .5 # F. B. Baker (2001; p. 159; .5 <= a < 2)
#
# test <- test[1, colnames(b)] # replace
#
# c <- vector() # saving names
# d <- 0 # name counter
# require(stringr)
#
# for(i in 1:length(colnames(b))) {
# if((abs(test[i]) >= .5) == TRUE){
# d <- d+1
# c[d] <- str_sub(colnames(test[i]), end = -4)
# } else {
#
# }
# }
find_a_lower <- which(data.frame(MDISC(result)) >= .5)
# find_a_upper <- which(data.frame(MDISC(result)) < 2.00) # F. B. Baker (2001; p. 159; .5 <= a < 2)
# find_a <- c(find_a_lower)#, find_a_upper)
fa_covdata_temp <- data.frame(fa_covdata[,find_a_lower])
# #test <- vector()
# test <- try(mod2values(result), silent = T)
# if(exists("test")){
# test <- data.frame(test[grep("a[0-9]+", test$name),])
# test <- subset(test, test$value >= .4)
# } else {
# test <- data.frame()
# }
#
#
# if(nrow(test)==0){
# message('Passing........')
# } else {
# test$item <- as.character(test$item)
# test$item <- as.factor(test$item)
# include <- as.character(test$item)
# myvars <- names(fa_covdata) %in% c(include)
# # fa_covdata_temp <- fa_covdata[myvars]
# fa_covdata_temp <- subset(fa_covdata, select=myvars)
fa_covdata <- fa_covdata_temp
# }
# evaluation of IRT itemfit
message('\nStage 2 of calcuation: Item Fit Based Evaluation')
ncol_stage2 <- ncol(fa_covdata)
if(ncol_stage2 == 0){
stop('All items are deleted')
} else {
message('\nCurrent number of Items: ', paste0(ncol_stage2))
}
# screening items (stage 1)
if(ncol_stage1 == ncol_stage2){
# evaluating model
#result <- k.aefa(fa_covdata, estimator='mirt', ...)
#print(summary(result))
} else {
# evaluating model
result <- k.aefa(fa_covdata, estimator='mirt', ...)
#print(summary(result))
}
message('estimating Theta')
try(test2_fscores <- fscores(object = result, rotate = 'geominQ', full.scores = T, plausible.draws = 100, QMC = TRUE, method = 'MAP', MI = 100))
message('estimating itemfit')
try(test2 <- mirt::itemfit(result, method = 'MAP', impute = 100, QMC = TRUE), silent = T)
if(exists("test2")) { # patch for mirt function bug...
test2$cal <- test2$S_X2/test2$df.S_X2
test2 <- subset(test2, test2$cal >= 3)
} else {
test2 <- data.frame()
#test2 <- 0
}
if(nrow(test2)==0){
message('Passing........')
} else {
exclude <- as.character(test2$item)
exclude <- as.factor(exclude)
exclude <- as.character(exclude)
myvars <- names(fa_covdata) %in% c(exclude)
fa_covdata <- fa_covdata[!myvars]
}
message('\nStage 3 of calcuation: Factor Loading Based Evaluation')
for(i in 1:10000){
ncol_stage3 <- ncol(fa_covdata)
if(ncol_stage3 == 0){
stop('All items are deleted')
} else {
message('\nCurrent number of Items: ', paste0(ncol_stage3))
}
if(ncol_stage1 == ncol_stage3 | ncol_stage2 == ncol_stage3){
# evaluating model
#result <- k.aefa(fa_covdata, estimator='mirt', ...)
#print(summary(result))
} else {
# evaluating model
result <- k.aefa(fa_covdata, estimator='mirt', ...)
#print(summary(result))
}
if(ncol(result@Fit$F)==1){
rotF_geomin <- data.frame(result@Fit$F)
h2 <- data.frame(result@Fit$h2)
} else {
# getting factor loadings
message('Rotating Factor Solution now')
rotF_geomin <- GPArotation::geominQ(result@Fit$F, maxit = 100000)
rotF_geomin <- data.frame(rotF_geomin$loadings)
h2 <- data.frame(result@Fit$h2)
}
if(i > 1){
exclude_rownum1 <- NULL
exclude_rownum1_low <- NULL
exclude_rownum2 <- NULL
exclude_rownum2_low <- NULL
exclude_rownum3 <- NULL
exclude_rownum3_low <- NULL
exclude <- NULL
exclude_rownum_low <- NULL
myvars <- NULL
if(length(exclude_rownum1) != 0){
try(rm(exclude_rownum1, exclude_rownum2, exclude_rownum2_low, exclude_rownum3, exclude_rownum3_low, exclude, exclude_rownum_low, myvars), silent = T)
} else {
try(rm(exclude_rownum1, exclude_rownum1_low, exclude_rownum2, exclude_rownum2_low, exclude_rownum3, exclude_rownum3_low, exclude, exclude_rownum_low, myvars), silent = T)
}
}
# evaluation of factor structure (stage 2) -- evaluating cross loadings and very small loadngs
exclude_rownum1 <- which(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) > 0.4) >=2) # cross loadings based delection
if(length(exclude_rownum1) != 0){
exclude_rownum1_low <- as.numeric(names(which(min(rowSums(abs(rotF_geomin[rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) > 0.4) >=2,]))) == rowSums(abs(rotF_geomin[rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) > 0.4) >=2,])))))
}
exclude_rownum2 <- which(h2 < .5) # communality based delection
exclude_rownum2_low <- which(h2 == min(h2))
exclude_rownum3 <- which(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]) < 0.4) == ncol(rotF_geomin)) # fail to load any factors
exclude_rownum3_low <- which(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)])) == min(rowSums(abs(rotF_geomin[1:ncol(rotF_geomin)]))))
# exclude_rownum <- as.numeric(paste(c(exclude_rownum1, exclude_rownum2, exclude_rownum3))) # list of doing delection
if(length(exclude_rownum3) > 0){
message('[WARN] No loadings to any factor(s) occured!')
exclude_rownum_low <- as.numeric(paste(c(exclude_rownum3_low))) # list of doing delection
# message(paste0(exclude_rownum_low))
} else if(length(exclude_rownum1) > 0){
message('[WARN] Cross Loadings occured!')
exclude_rownum_low <- as.numeric(paste(c(exclude_rownum1_low))) # list of doing delection
} else if(length(exclude_rownum2) > 0){
message('[WARN] Communality problem occured!')
exclude_rownum_low <- as.numeric(paste(c(exclude_rownum2_low))) # list of doing delection
} else {
message('[Done!]')
exclude_rownum_low <- NULL
# exclude_rownum_low <- as.numeric(paste(c(exclude_rownum1_low, exclude_rownum2_low, exclude_rownum3_low))) # list of doing delection
}
# exclude_rownum <- as.numeric(paste(c(exclude_rownum1, exclude_rownum3))) # list of doing delection
# the start of variable delection
if(length(exclude_rownum_low) > 0) { # if number of delection is non-zero
exclude <- vector() # make vectors
# j <- 0 # set to zero exclude list counter
# for(i in c(exclude_rownum)){ # saving list of delection variable
# j <- j + 1
#print(j)
#message('Trying to extract ', paste0(i), ' factor solution') -- for debug code
# print(names(fa_covdata)[i])
exclude <- names(data.frame(result@Data$data))[c(exclude_rownum_low)]
message(paste0(exclude))
# exclude[j] <- names(fa_covdata)[c(exclude_rownum_low)]
# }
myvars <- names(data.frame(result@Data$data)) %in% c(exclude)
fa_covdata <- data.frame(result@Data$data)[!myvars]
} else { # (length(exclude_rownum == 0))
try(print(findM2(result, calcNull = T, Theta = fscores(result, full.scores = T, plausible.draws = 100, rotate = "geominQ", MI = 100, QMC = TRUE, method = 'MAP'), impute = 100)), silent = T)
#cat('number of iteration : paste0(a))
return(result)
}
}
}
k.sampling <- function(data, n){
set.seed(1234)
data <- data[sample(nrow(data), n), ]
return(data)
}
# k.select <- function(model, items){
# select <- data.frame(model@Data$data[,rownames(subset(model@Fit$F, model@Fit$F >= min(sort(model@Fit$F, decreasing = T)[1:items])))])
# return(select)
#
#
# }
k.select <- function(model, items){
select <- model@Data$data[,c(subset(colnames(model@Data$data), model@Fit$h2 >= min(sort(model@Fit$h2, decreasing = T)[1:items])))]
return(select)
}
# finding problemistic data
k.fixdata <- function(data, start, end, bioend){
data <- data[!duplicated(data[start:bioend]),] # remove duplicate cases automatically
dropVector <- vector()
j <- 0
for(i in 1:nrow(data)){
if(std(as.numeric(data[i,start:end])) == 0){
j <- j + 1
dropVector[j] <- i
print(dropVector)
} else {
}
}
return(data[c(-dropVector),])
}
# surveyFA addon
fastFIFA <- function(x, covdata = NULL, formula = NULL, SE = T, SE.type = "sandwich", skipNominal = F,
forceGRSM = F, assumingFake = F, masterThesis = F, forceRasch = F, unstable = F,
forceMHRM = F, forceNormalEM = T, itemkeys = NULL, survey.weights = NULL, allowMixedResponse = T,
forceUIRT = F, skipIdealPoint = F, MHRM_SE_draws = 1e+4, forceNRM = F,
diagnosis = F, forceDefalutAccelerater = F, forceDefaultOptimizer = F, EnableFMHRM = F, testlets = NULL, plotOn = T, GenRandomPars = T,
fixed = ~1, random = NULL, lr.fixed = ~1, lr.random = NULL, MH_BURNIN = 1500, MH_SEMCYCLES = 1000, ...){
if(!require('mirt')){
install.packages('mirt')
library('mirt')
}
x <- x[,psych::describe(x)$range > 0] # delete no variance items
for(i in 1:ncol(x)){
try(invisible(gc()), silent = T) # garbage cleaning
# testlets
ActualTestlets <- testlets
if(length(testlets) != 0){
if(!require('plyr')){
install.packages('plyr')
library('plyr')
}
TestletActivated <- T
if(is.character(ActualTestlets)){
ActualTestlets <- as.factor(ActualTestlets)
}
ActualTestlets <- as.integer(ActualTestlets)
ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))
# if(sum(ActualTestlets %in% 1) == 0){
# ActualTestlets <- ActualTestlets - min(ActualTestlets, na.rm = T) + 1
# }
if(sum(is.na(ActualTestlets)) == length(ActualTestlets)){
ActualTestlets <- NULL
TestletActivated <- F
} else {
ActualTestlets <- plyr::mapvalues(ActualTestlets, as.numeric(attributes(as.factor(ActualTestlets))$levels), seq(length(attributes(as.factor(c(ActualTestlets)))$levels)))
}
} else {
TestletActivated <- F
}
if(TestletActivated){
message('\nfactor numbers: ', paste0(i+max(na.omit(ActualTestlets))))
} else {
if (i == 1){
message('\nfactor number: ', paste0(i))
} else {
message('\nfactor numbers: ', paste0(i))
}
}
# if(sum(is.na(x)) == 0 | nrow(x) > 5000){
NofCores <- parallel::detectCores()
NofCores <- round(NofCores / 1.1)
if(NofCores > 12){
NofCores <- 12
}
try(invisible(mirt::mirtCluster(spec = NofCores)), silent = T)
# }
# optimizer config
if(length(covdata) == 0){ # if no covariate variables
if(TestletActivated){
# print(ActualTestlets)
if(sum(na.omit(ActualTestlets) > 2) != 0 | forceMHRM){
estimationMETHOD <- 'MHRM'
optimINPUT <- NULL
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 4000
} else {
estimationMETHOD <- 'EM'
optimINPUT <- NULL
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 1000
}
} else if(forceMHRM == T | forceGRSM == T | assumingFake == T | masterThesis == T) {
estimationMETHOD <- 'MHRM'
optimINPUT <- NULL
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 4000
} else if(length(survey.weights) != 0) {
if (forceNormalEM == T && i < 4){
estimationMETHOD <- 'EM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
} else {
estimationMETHOD <- 'QMCEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
}
} else if(i < 4){ # EM will activate 1-3 factors
if (forceNormalEM == T && unstable == T) {
estimationMETHOD <- 'SEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- NULL
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 4000
} else if(unstable == T){
estimationMETHOD <- 'SEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- NULL
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
} else if (forceNormalEM == T) {
estimationMETHOD <- 'EM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 1000
} else {
estimationMETHOD <- 'EM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- TRUE
NCYCLES <- 4000
}
} else {
if(unstable == T){
if(forceNormalEM == T && i < 4){
estimationMETHOD <- 'EM'
} else {
estimationMETHOD <- 'QMCEM'
}
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
} else {
estimationMETHOD <- 'MHRM'
optimINPUT <- NULL
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 4000
}
}
covdataINPUT <- NULL
formulaINPUT <- NULL
message('estimation method: ', paste0(estimationMETHOD))
} else { # with covariate
covdataINPUT <- covdata
formulaINPUT <- formula
if(NROW(random) != 0 | NROW(lr.random) != 0){
estimationMETHOD <- 'MHRM'
optimINPUT <- NULL
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 4000
} else if(forceMHRM == T| forceGRSM == T | assumingFake == T | masterThesis == T | sum(na.omit(ActualTestlets) > 2) != 0){
message('MHRM currently not supported with latent regressors')
estimationMETHOD <- 'QMCEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 1e+4
} else if(length(survey.weights) != 0) {
estimationMETHOD <- 'QMCEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
} else if(i < 4){ # use EM 1-3 factors
if(unstable == T){
estimationMETHOD <- 'QMCEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
} else if (forceNormalEM == T) {
estimationMETHOD <- 'EM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- 1000
} else {
estimationMETHOD <- 'EM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- TRUE
NCYCLES <- 4000
}
} else {
estimationMETHOD <- 'QMCEM'
if(forceDefaultOptimizer){
optimINPUT <- NULL
} else {
optimINPUT <- 'nlminb'
}
optimCTRL <- NULL
empiricalhist <- FALSE
NCYCLES <- NULL
}
message('estimation method: ', paste0(estimationMETHOD))
if(NROW(fixed) != 0){
message('fixed effect formula: ', paste0(fixed))
}
if(NROW(random) != 0){
message('random effect formula: ', paste0(random))
} else {
message('latent regression formula: ', paste0(formulaINPUT))
}
}
if(estimationMETHOD == 'EM'){
message('Empirical Histogram for find Prior distribution: ', empiricalhist)
}
# forcing SE estimation activate
if((sum(is.na(data.frame(x))) != 0) && (SE.type != 'MHRM')){
SE <- T
if(length(covdata) == 0){
if(estimationMETHOD == 'MHRM'){ # Richadson (BL) isn't support MHRM estimation method
SE.type <- 'MHRM'
} #else {
# SE.type <- "sandwich" # Oakes
# }
} else {
if(estimationMETHOD == 'MHRM'){ # Richadson (BL) isn't support MHRM estimation method
SE.type <- 'MHRM'
} #else {
# SE.type <- "Oakes" # Oakes
# }
}
}
# switch SE estimation temporary
if(length(covdata) != 0 && NROW(random) == 0 && SE.type == 'sandwich' && SE == T){
message('sandwich estimation currently not supproted with latent regressors')
SE.type <- "complete"
}
# standard error estimation config
if((SE == T && SE.type == 'SEM') == T){
accelerateINPUT <- 'none'
TOLINPUT <- NULL
SEtolINPUT <- 1e-4
symmetricINPUT <- FALSE
message('TOL: ', 'default', ' / SEtol: ', SEtolINPUT, ' / SE.type: ', SE.type, ' / Accelerator: ',
accelerateINPUT, ' / Symmetric SEM: ', symmetricINPUT)
} else if((SE == T && estimationMETHOD == 'MHRM') == T){
if(forceDefalutAccelerater){
accelerateINPUT <- 'Ramsay'
} else {
accelerateINPUT <- 'squarem'
}
TOLINPUT <- NULL
SEtolINPUT <- 1e-4
symmetricINPUT <- FALSE
if(EnableFMHRM){
SE.type <- 'FMHRM'
} else {
SE.type <- 'MHRM' # more efficient
}
message('TOL: ', 'default', ' / SEtol: ', SEtolINPUT, ' / SE.type: ', SE.type, ' / Accelerator: ',
accelerateINPUT, ' / Symmetric SEM: ', symmetricINPUT)
} else {
if(forceDefalutAccelerater){
accelerateINPUT <- 'Ramsay'
} else {
accelerateINPUT <- 'squarem'
}
TOLINPUT <- NULL
SEtolINPUT <- 1e-4
symmetricINPUT <- FALSE
if(SE == T){
message('Standard Error Estimation: On')
message('TOL: ', 'default', ' / SEtol: ', SEtolINPUT, ' / SE.type: ', SE.type, ' / Accelerator: ',
accelerateINPUT, ' / Symmetric SEM: ', symmetricINPUT)
}
}
# removeEmptyRows config
if(length(covdata) != 0){
removeEmptyRowsConf <- FALSE
} else {
removeEmptyRowsConf <- TRUE
}
if(sum(psych::describe(x)$range == 1) == ncol(x)){ # dichotomous items
# forceRasch (dichotomous)
if(forceRasch == T){
message('\nMIRT model: Rasch')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'Rasch',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
# try(invisible(mirt::mirtCluster(remove = T)), silent = T)
if(exists('modTEMP')){
try(return(modTEMP))
} else {
stop('Rasch model is inappropriate with this data. turn off the forceRasch = T option.')
}
}
# if(nrow(x) >= 50){
if(length(ActualTestlets) != 0){
if(diagnosis == F){
message('\nMIRT model: Compensatory 4PL')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '4PL',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '4PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 3PL with upper asymptote (slip) estimated')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '3PLu',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PLu',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Partially compensatory 3PL')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'PC3PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 3PL')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '3PL',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# }
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 2PL')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = '2PL',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '2PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Partially compensatory 2PL')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'PC2PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && skipIdealPoint == F){
message('\nMIRT model: ideal point')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'ideal',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'ideal',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Rasch')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'Rasch',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'Rasch',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(SE.type == 'MHRM'){
message('MHRM trials were not efficient. try to calibrate the model with EM')
SE.type <- 'sandwich'
} else {
message('try to recalibrate the model with mirt::bfactor() function')
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 4PL')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '4PL',
lerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 3PL with upper asymptote (slip) estimated')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '3PLu',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Partially compensatory 3PL')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'PC3PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 3PL')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '3PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
# }
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 2PL')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = '2PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Partially compensatory 2PL')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'PC2PL',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
if(exists('modTEMP') == F && skipIdealPoint == F){
message('\nMIRT model: ideal point')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'ideal',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Rasch')
try(modTEMP <- mirt::bfactor(data = x, ActualTestlets, itemtype = 'Rasch',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
}
} else {
if(diagnosis == F){
message('\nMIRT model: Compensatory 4PL')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '4PL',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '4PL', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 3PL with upper asymptote (slip) estimated')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '3PLu',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PLu', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && i == 1 && diagnosis == F){
message('\nMIRT model: Partially compensatory 3PL')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'PC3PL', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 3PL')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '3PL',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PL', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# }
if(exists('modTEMP') == F && diagnosis == F){
message('\nMIRT model: Compensatory 2PL')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = '2PL',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '2PL', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && i == 1 && diagnosis == F){
message('\nMIRT model: Partially compensatory 2PL')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'PC2PL', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && skipIdealPoint == F){
message('\nMIRT model: ideal point')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'ideal',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'ideal', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') == F && i == 1 && diagnosis == F){
message('\nMIRT model: Rasch')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'Rasch', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
}
if(exists('modTEMP') == F && i == 1){
message('\nMIRT model: spline response')
if(estimationMETHOD == 'MHRM'){
estimationMETHOD <- 'EM'
message('switching MHRM to ', estimationMETHOD, ' ...')
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'spline', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
message('Model is unstable so that try to remedy the model automatically as soon as possible')
}
}
}
if(exists('modTEMP') == F){
if(i == 1){
stop('Fail to find Factor solutions: Model didn\'t converge.')
} else {
return(modOLD)
}
}
} else if(length(itemkeys) != 0){ # 2-4PLNRM
if(length(ActualTestlets) != 0){
message('\nMIRT model: Compensatory 4PL Nominal response')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '4PLNRM',
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Compensatory 3PL Nominal response')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PLNRM',
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Compensatory 3PL Nominal response with upper asymptote estimated')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '3PLuNRM',
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
# }
if(exists('modTEMP') == F){
message('\nMIRT model: Compensatory 2PL Nominal response')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = '2PLNRM',
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Nominal response without keys')
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'nominal',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, key = NULL, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
message('Model may unstable but Trying to remedy automatically')
}
}
}
} else {
message('\nMIRT model: Compensatory 4PL Nominal response')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '4PLNRM', method = estimationMETHOD,
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Compensatory 3PL Nominal response')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PLNRM', method = estimationMETHOD,
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Compensatory 3PL Nominal response with upper asymptote estimated')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '3PLuNRM', method = estimationMETHOD,
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
# }
if(exists('modTEMP') == F){
message('\nMIRT model: Compensatory 2PL Nominal response')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = '2PLNRM', method = estimationMETHOD,
key = itemkeys, accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
if(exists('modTEMP') == F){
message('\nMIRT model: Nominal response without keys')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'nominal', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, key = NULL, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
message('Model may unstable but Trying to remedy automatically')
}
}
}
}
if(exists('modTEMP') == F){
if(i == 1){
stop('Fail to find Factor solutions: Model didn\'t converge.')
} else {
return(modOLD)
}
}
} else if((sum(psych::describe(x)$range == 1) != 0) && allowMixedResponse == T) { # mixed format (Construct Responses + Multiple Choices, under construction)
if(skipNominal == F){
if(skipIdealPoint == F){
itemtype_mixed <- vector()
for(i in 1:ncol(x)){
if(psych::describe(x[,i])$range == 1){
itemtype_mixed[i] <- 'ideal'
dichotomous_type <- 'ideal point'
} else {
itemtype_mixed[i] <- 'nominal'
}
}
message('\nMIRT model: nominal response + ', paste0(dichotomous_type))
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = itemtype_mixed, method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 20000000000, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = 20000), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
if(exists('modTEMP') == F){
itemtype_mixed <- vector()
for(i in 1:ncol(x)){
if(psych::describe(x[,i])$range == 1){
if(nrow(x) >= 5000){
itemtype_mixed[i] <- '4PL'
dichotomous_type <- '4PL'
} else if(nrow(x) >= 2000){
itemtype_mixed[i] <- '3PL'
dichotomous_type <- '3PL'
} else {
itemtype_mixed[i] <- '2PL'
dichotomous_type <- '2PL'
}
} else {
itemtype_mixed[i] <- 'nominal'
}
}
message('\nMIRT model: nominal response + ', paste0(dichotomous_type), ' (automatically set by sample size)')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = itemtype_mixed, method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 20000000000, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = 20000), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
}
# generalized partial credit model (non-sequential)
if(exists('modTEMP') == F && forceNRM == F){
itemtype_mixed <- vector()
for(i in 1:ncol(x)){
if(psych::describe(x[,i])$range == 1){
if(nrow(x) >= 5000){
itemtype_mixed[i] <- '4PL'
} else if(nrow(x) >= 2000){
itemtype_mixed[i] <- '3PL'
} else {
itemtype_mixed[i] <- 'ideal'
}
} else {
itemtype_mixed[i] <- 'gpcm'
}
}
message('\nMIRT model: Generalized partial credit + ideal or 3-4PL')
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = itemtype_mixed, method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F){
itemtype_mixed <- vector()
for(i in 1:ncol(x)){
if(psych::describe(x[,i])$range == 1){
if(nrow(x) >= 5000){
itemtype_mixed[i] <- '4PL'
} else if(nrow(x) >= 2000){
itemtype_mixed[i] <- '3PL'
} else {
itemtype_mixed[i] <- 'ideal'
}
} else {
itemtype_mixed[i] <- 'graded'
}
}
message('\nMIRT model: Graded response + ideal or 3-4PL')
try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = itemtype_mixed, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
} else {
stop('retry to add forceMHRM = TRUE argument')
}
}
# finally, if can not converge
if(exists('modTEMP') == F){
if(i == 1){
stop('Fail to find Factor solutions: Model didn\'t converge.')
} else {
return(modOLD)
}
}
} else { # polytomous items
# forceRasch (PCM)
if(forceRasch == T){
message('\nMIRT model: Rasch')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'Rasch',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ... = ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
# try(invisible(mirt::mirtCluster(remove = T)), silent = T)
if(exists('modTEMP')){
try(return(modTEMP))
} else {
stop('Rasch model is inappropriate with this data. turn off the forceRasch = T option.')
}
}
# forceGRSM (polytomous)
# for detecting fake responses
# if(SE == T){ # if grsm, Standard error estimation was unsuccessful.
# } else {
#if((max(describe(x)$range) - min(describe(x)$range)) == 0 | forceGRSM == T | assumingFake == T | masterThesis == T){
if(forceGRSM == T | assumingFake == T | masterThesis == T){
x <- data.frame(x)
if(length(which(describe(x)$max > median(describe(x)$max))) != 0){ # preventing weird input (e.g: 6 in 5 point scale)
for(ii in which(describe(x)$max > median(describe(x)$max))){
x[,ii] <-plyr::mapvalues(x[,ii],psych::describe(x[,which(describe(x)$max > median(describe(x)$max))])$max, NA)
}
}
x <- x[,which(describe(x)$range == max(describe(x)$range))]
k <- vector()
for(iii in 1:length(x)){
for(j in range(na.omit(x[iii]))[1]:range(na.omit(x[iii]))[2]){
if((sum(na.omit(x[iii]) == j) == 0) == TRUE){
k[length(k)+1] <- colnames(x[iii])
}
}
}
x <- (x[,!colnames(x) %in% k])
for(iiii in 1:ncol(x)){
x[,iiii] <- as.integer(x[,iiii])
}
message('\nMIRT model: graded rating scale')
if(i == 1){
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsmIRT', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
if(exists('modTEMP') == F){
if(i == 1){
stop('Fail to find Factor solutions: Model didn\'t converge.')
} else {
return(modOLD)
}
}
} else {
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsm', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
if(exists('modTEMP') == F){
if(i == 1){
stop('Fail to find Factor solutions: Model didn\'t converge.')
} else {
return(modOLD)
}
}
}
if(modTEMP@OptimInfo$converged == 1){
skipNominal <- T
}
}
# }
# polytomous
# nominal model
if(length(ActualTestlets) != 0){
if(skipNominal == F){
message('\nMIRT model: nominal response')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'nominal',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'nominal',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# generalized partial credit model (non-sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Generalized partial credit')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'gpcm',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, itemtype = 'gpcm',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Graded response')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = doBfactor2mod(x, ActualTestlets), itemtype = 'graded',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = doBfactor2mod(x, ActualTestlets), method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(SE.type == 'MHRM'){
message('MHRM trials were not efficient. try to calibrate the model with EM')
SE.type <- 'sandwich'
} else {
message('try to recalibrate the model with mirt::bfactor() function')
}
if(exists('modTEMP') == F && skipNominal == F){
message('\nMIRT model: nominal response')
try(modTEMP <- mirt::bfactor(data = x, model = ActualTestlets, itemtype = 'nominal',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
# generalized partial credit model (non-sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Generalized partial credit')
try(modTEMP <- mirt::bfactor(data = x, model = ActualTestlets, itemtype = 'gpcm',
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Graded response')
try(modTEMP <- mirt::bfactor(data = x, model = ActualTestlets, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){rm(modTEMP)}
}
}
} else {
if(skipNominal == F){
message('\nMIRT model: nominal response')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'nominal',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'nominal', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# generalized partial credit model (non-sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Generalized partial credit')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'gpcm',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'gpcm', method = estimationMETHOD,
accelerate = accelerateINPUT, calcNull = T,
technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT, SEtol = SEtolINPUT,
removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES), TOL = TOLINPUT, covdata = covdataINPUT,
formula = formulaINPUT, optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = F)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Graded response')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'graded',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F){
message('\nMIRT model: Graded rating scale')
if(i == 1){
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'grsmIRT',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsmIRT', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
} else {
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'grsm',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'grsm', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
}
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F && i == 1){
message('\nMIRT model: Partial Credit')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'Rasch',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'Rasch', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | modTEMP@OptimInfo$secondordertest == F){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | modTEMP_MIXED@OptimInfo$secondordertest == F){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
# graded response model (sequential)
if(exists('modTEMP') == F && forceNRM == F && i == 1){
message('\nMIRT model: Rating Scale')
if(NROW(random) != 0 | NROW(lr.random) != 0){
try(modTEMP_MIXED <- fitMLIRT(dat = x, model = i, itemtype = 'rsm',
covdata = covdataINPUT, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
GenRandomPars = GenRandomPars, NCYCLES = 4000, BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, symmetric = symmetricINPUT,
accelerate = accelerateINPUT, MHRM_SE_draws = MHRM_SE_draws, survey.weights = survey.weights,
removeEmptyRows = removeEmptyRowsConf, SEtol = SEtolINPUT))
}
try(modTEMP <- mirt::mirt(data = x, model = i, itemtype = 'rsm', method = estimationMETHOD, GenRandomPars = GenRandomPars, accelerate = accelerateINPUT,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws, symmetric = symmetricINPUT,
SEtol = SEtolINPUT, removeEmptyRows = removeEmptyRowsConf, NCYCLES = NCYCLES),
TOL = TOLINPUT, covdata = covdataINPUT, formula = formulaINPUT,
optimizer = optimINPUT, solnp_args = optimCTRL, SE = SE,
SE.type = SE.type, survey.weights = survey.weights, empiricalhist = empiricalhist, ...), silent = T)
if(exists('modTEMP')){
if(modTEMP@OptimInfo$converged == F | (modTEMP@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP)
}
}
if(exists('modTEMP_MIXED')){
if(modTEMP_MIXED@OptimInfo$converged == F | (modTEMP_MIXED@OptimInfo$secondordertest == F && i == 1)){
rm(modTEMP_MIXED)
}
}
if(exists('modTEMP') && exists('modTEMP_MIXED')){
if(modTEMP_MIXED@Fit$DIC < modTEMP@Fit$DIC){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
if(!exists('modTEMP') && exists('modTEMP_MIXED')){
modTEMP <- modTEMP_MIXED
rm(modTEMP_MIXED)
}
}
}
# finally, if can not converge
if(exists('modTEMP') == F){
if(i == 1){
stop('Fail to find Factor solutions: Model didn\'t converge.')
} else {
return(modOLD)
}
}
}
try(invisible(mirt::mirtCluster(remove = T)), silent = T)
if(class(modTEMP)[1] == 'MixedClass'){
MLM_rotate_formula_mod <- mirt::mirt(data = modTEMP@Data$data, model = modTEMP@Model$model, itemtype = modTEMP@Model$itemtype, pars = 'values')
MLM_rotate_formula_mod_original <- mirt::mod2values(modTEMP)
if(sum(MLM_rotate_formula_mod_original$name == '(Intercept)') != 0){
MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original[!MLM_rotate_formula_mod_original$name == '(Intercept)',]
}
# MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original
MLM_rotate_formula_mod$value[which(MLM_rotate_formula_mod$item %in% colnames(modTEMP@Data$data))] <- MLM_rotate_formula_mod_original$value[which(MLM_rotate_formula_mod_original$item %in% colnames(modTEMP@Data$data))]
MLM_rotate_formula_mod$est <- F
# print(MLM_rotate_formula_mod)
MixedModelFlag <- T
modTEMP_MIXED <- modTEMP
modTEMP <- mirt::mirt(data = modTEMP@Data$data, model = modTEMP@Model$model,
itemtype = modTEMP@Model$itemtype, pars = MLM_rotate_formula_mod, method = 'QMCEM',
GenRandomPars = GenRandomPars,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws),
survey.weights = survey.weights)
} else {
MixedModelFlag <- F
}
if(i == 1 && length(ActualTestlets) == 0 && plotOn == T){ # ICC printing
try(print(plot(modTEMP, type = 'infoSE')))
try(print(plot(modTEMP, type = 'infotrace', facet_items = TRUE)))
try(print(plot(modTEMP, type = 'trace')))
}
if(forceUIRT == T | TestletActivated == T){
if(MixedModelFlag){
return(modTEMP_MIXED)
} else {
return(modTEMP)
}
}
if(i == 1 && modTEMP@OptimInfo$converged == F && sum(modTEMP@Model$itemtype %in% c('spline', 'nominal')) == 0){
stop('No convergence')
}
if(i > 1){ # evaluation of local optimal number of factors
if(ncol(modTEMP@Fit$F) == 1){
rotMat <- modTEMP@Fit$F
} else {
rotMat <- GPArotation::geominQ(modTEMP@Fit$F, maxit = 1e+5)$loadings
}
if(MixedModelFlag){
modTEMP <- modTEMP_MIXED
}
if(modTEMP@Fit$DIC > modOLD@Fit$DIC | modTEMP@OptimInfo$converged == F | sum(round(modTEMP@Fit$h2, 3) >= .999) != 0){ # modTEMP@Fit$AICc > modOLD@Fit$AICc |
message('optimal factor numbers: ', paste0(i-1))
return(modOLD)
} #else if(sum(colSums(round(abs(rotMat), 2) > .4) < 2) != 0) {
#message('optimal factor numbers: ', paste0(i-1))
#return(modOLD)
#}
}
if(exists('modTEMP') == T){
if(MixedModelFlag){
modOLD <- (modTEMP_MIXED)
try(rm(modTEMP_MIXED), silent = T)
try(rm(modTEMP), silent = T)
} else {
modOLD <- modTEMP
try(rm(modTEMP), silent = T)
}
} else {
stop('Fail to convergence')
}
}
}
surveyFA <- function(data = ..., covdata = NULL, formula = NULL, SE = T,
SE.type = "sandwich", skipNominal = F, forceGRSM = F,
assumingFake = F, masterThesis = F, forceRasch = F,
unstable = F, forceNormalEM = T, forceMHRM = F,
printFactorStructureRealtime = F, itemkeys = NULL,
survey.weights = NULL, allowMixedResponse = T, autofix = F,
forceUIRT = F, skipIdealPoint = F, MHRM_SE_draws = 1e+4,
bifactorSolution = T, skipS_X2 = F, forceNRM = F, needGlobalOptimal = T,
pilotTestMode = F, forceConsiderPositiveZh = F, forceDefalutAccelerater = F, forceDefaultOptimizer = F, EnableFMHRM = F, testlets = NULL,
minimumLeftItems = 3, plotOn = T, forceCTTmode = F, GenRandomPars = T, coefAlwaysBePositive = T,
fixed = ~1, random = NULL, lr.fixed = ~1, lr.random = NULL, MH_BURNIN = 1500, MH_SEMCYCLES = 1000, ...) {
message('---------------------------------------------------------')
message(' k.aefa: kwangwoon automated exploratory factor analysis ')
message('---------------------------------------------------------\n')
message('Calculating Initial Factor model')
iteration_num <- 1
message('Iteration: ', iteration_num, '\n')
if(bifactorSolution) {
rotateCriteria <- 'bifactorQ'
} else {
rotateCriteria <- 'geominQ'
}
if(!is.data.frame(data) && !is.matrix(data)){
surveyFixMod <- data
} else {
if(sum(is.na(data)) != 0 && pilotTestMode == T){
message('input data n size: ', nrow(data))
data <- data[which(rowSums(is.na(data)) < ncol(data)*(1-3/4)),]
message('current data n size: ', nrow(data))
}
surveyFixMod <- fastFIFA(x = as.data.frame(data), covdata = as.data.frame(covdata),
formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal,
forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis,
forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM,
itemkeys = itemkeys, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse,
autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = testlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
}
workKeys <- itemkeys
workTestlets <- testlets
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
surveyFixMod <- deepFA(surveyFixMod, survey.weights)
}
itemFitDone <- FALSE
while (!itemFitDone) {
try(invisible(gc()), silent = T) # garbage cleaning
surveyFixModRAW <- data.frame(surveyFixMod@Data$data)
surveyFixModCOV <- data.frame(attr(surveyFixMod@ParObjects$lrPars, "df"))
if(ncol(surveyFixModRAW) > minimumLeftItems){
if(class(surveyFixMod)[1] == 'MixedClass'){
MLM_rotate_formula_mod <- mirt::mirt(data = surveyFixMod@Data$data, model = surveyFixMod@Model$model, itemtype = surveyFixMod@Model$itemtype, pars = 'values')
MLM_rotate_formula_mod_original <- mirt::mod2values(surveyFixMod)
if(sum(MLM_rotate_formula_mod_original$name == '(Intercept)') != 0){
MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original[!MLM_rotate_formula_mod_original$name == '(Intercept)',]
}
# MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original
MLM_rotate_formula_mod$value[which(MLM_rotate_formula_mod$item %in% colnames(surveyFixMod@Data$data))] <- MLM_rotate_formula_mod_original$value[which(MLM_rotate_formula_mod_original$item %in% colnames(surveyFixMod@Data$data))]
MLM_rotate_formula_mod$est <- F
# print(MLM_rotate_formula_mod)
MixedModelFlag <- T
surveyFixMod <- mirt::mirt(data = surveyFixMod@Data$data, model = surveyFixMod@Model$model,
itemtype = surveyFixMod@Model$itemtype, pars = MLM_rotate_formula_mod, method = 'QMCEM',
GenRandomPars = GenRandomPars,
calcNull = T, technical = list(BURNIN = MH_BURNIN, SEMCYCLES = MH_SEMCYCLES, MAXQUAD = 2000000, delta = 1e-20, MHRM_SE_draws = MHRM_SE_draws),
survey.weights = survey.weights)
} else {
MixedModelFlag <- F
}
message('\nChecking item local independence assumption')
iteration_num <- iteration_num + 1
message('Iteration: ', iteration_num, '\n')
if(sum(surveyFixMod@Model$itemtype %in% 'spline') > 0){
fscoreMethod <- 'EAP'
} else {
fscoreMethod <- 'MAP'
}
if(sum(is.na(surveyFixModRAW)) == 0 | nrow(surveyFixModRAW) > 5000){
NofCores <- parallel::detectCores()
NofCores <- round(NofCores / 1.1)
if(NofCores > 8){
NofCores <- 8
}
try(invisible(mirt::mirtCluster(spec = NofCores)), silent = T)
}
if(!forceCTTmode){
if(sum(is.na(surveyFixMod@Data$data)) == 0){
try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('S_X2', 'Zh', 'infit'),
method = fscoreMethod,
QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
} else {
try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('S_X2', 'Zh', 'infit'),
impute = 100,
method = fscoreMethod,
QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
}
if(!exists('surveyFixMod_itemFit')){
if(sum(is.na(surveyFixMod@Data$data)) == 0){
try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('Zh', 'infit'),
method = fscoreMethod,
QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
} else {
try(surveyFixMod_itemFit <- mirt::itemfit(x = surveyFixMod, fit_stats = c('Zh', 'infit'),
impute = 100,
method = fscoreMethod,
QMC = T, rotate = rotateCriteria, maxit = 1e+5), silent = T)
}
S_X2ErrorFlag <- T
} else {
S_X2ErrorFlag <- F
}
try(invisible(mirt::mirtCluster(remove = T)), silent = T)
if(!S_X2ErrorFlag){
if(sum(surveyFixMod_itemFit$df.S_X2 == 0, na.rm = T) > length(surveyFixMod_itemFit$df.S_X2)/5){
S_X2ErrorFlag <- T
} else if(sum(is.nan(surveyFixMod_itemFit$p.S_X2)) > length(surveyFixMod_itemFit$p.S_X2)/5) {
S_X2ErrorFlag <- T
}
}
print(surveyFixMod_itemFit)
# item evaluation
if(S_X2ErrorFlag){
message('S_X2 can not be calcuate normally...')
}
}
# precalculation of CI for a1
ZeroList <- vector()
ZeroRange <- vector()
if(surveyFixMod@Options$SE == T && (length(workTestlets) != 0 | surveyFixMod@Model$nfact == 1)){
for(III in 1:NROW(coef(surveyFixMod))-1){
try(vec <- data.frame(coef(surveyFixMod)[III]), silent = T)
if(exists('vec')){
# print(vec) # -- for debug
if(is.na((vec[2,1] < 0 && vec[3,1] > 0))){
} else {
try(ZeroList[length(ZeroList)+1] <- (vec[2,1] < 0 && vec[3,1] > 0), silent = T)
try(ZeroRange[length(ZeroRange)+1] <- psych::describe(c(vec[2,1], vec[3,1]))$range, silent = T)
}
}
# print(ZeroList)# -- for debug
# print(ZeroRange)
}
}
vec2 <- coef(surveyFixMod, simplify = T)$items
# Standardized Log-Likelihood
if(forceCTTmode){ # force CTT-Like model
if(sum(ZeroList) > 0){ # which item include 0 in a1
message('\nItem discrimination include 0 / removing ', paste(colnames(surveyFixModRAW)[which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]))
workKeys <- workKeys[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]
workTestlets <- workTestlets[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
} else if(coefAlwaysBePositive && sum(vec2[,1] < 0) != 0){
message('\nItem discrimination include negative / removing ', paste(surveyFixMod_itemFit$item[which(min(vec2[,1]) == vec2[,1])]))
workKeys <- workKeys[-which(min(vec2[,1]) == vec2[,1])]
workTestlets <- workTestlets[-which(min(vec2[,1]) == vec2[,1])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(vec2[,1]) == vec2[,1])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
} else {
itemFitDone <- TRUE
}
} else if(length(which(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems] < -1.96)) != 0 && sum(is.na(surveyFixModRAW)) == 0){ # Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.
message('\nDrasgow, F., Levine, M. V., & Williams, E. A. (1985) / removing ', paste(surveyFixMod_itemFit$item[which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]))
workKeys <- workKeys[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
workTestlets <- workTestlets[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(length(which(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems] == -Inf)) != 0){ # Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.
message('\nDrasgow, F., Levine, M. V., & Williams, E. A. (1985) / removing ', paste(surveyFixMod_itemFit$item[which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]))
workKeys <- workKeys[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
workTestlets <- workTestlets[-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(sum(ZeroList) > 0){ # which item include 0 in a (currently, only a1)
message('\nItem discrimination include 0 / removing ', paste(surveyFixMod_itemFit$item[which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]))
workKeys <- workKeys[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]
workTestlets <- workTestlets[-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(abs(ZeroRange[ZeroList])) == abs(ZeroRange))], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(coefAlwaysBePositive && sum(vec2[,1] < 0) != 0){
message('\nItem discrimination include negative / removing ', paste(surveyFixMod_itemFit$item[which(min(vec2[,1]) == vec2[,1])]))
workKeys <- workKeys[-which(min(vec2[,1]) == vec2[,1])]
workTestlets <- workTestlets[-which(min(vec2[,1]) == vec2[,1])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(min(vec2[,1]) == vec2[,1])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(S_X2ErrorFlag) { # if Chi-squared can not be calculate with extremely big Zh
if(length(which(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems] > 1.96)) != 0 && sum(is.na(surveyFixModRAW)) == 0 && forceConsiderPositiveZh){ # Drasgow, F., Levine, M. V., & Williams, E. A. (1985). Appropriateness measurement with polychotomous item response models and standardized indices. British Journal of Mathematical and Statistical Psychology, 38(1), 67-86.
message('\nDrasgow, F., Levine, M. V., & Williams, E. A. (1985) / removing ', paste(surveyFixMod_itemFit$item[which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]))
workKeys <- workKeys[-which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
workTestlets <- workTestlets[-which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$Zh[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else {
itemFitDone <- TRUE
}
} else if(S_X2ErrorFlag == F && (surveyFixMod@Model$nfact == 1 | length(workTestlets) > 0)) { # if Chi-squared can be calculate
# chi-square exception processing
if(sum(is.na(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) == 0 && sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) == 0){ # avoid unexpected situation
message('all items df are 0. skipping evaluation...')
itemFitDone <- TRUE
} else if(sum(is.na(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) != 0 && sum(is.na(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])) != surveyFixMod@Data$nitems && (sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == 0) < length(1:surveyFixMod@Data$nitems)/2)){
message('\nremoving items df is NA / ', paste(surveyFixMod_itemFit$item[which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)]))
workKeys <- workKeys[-which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)]
workTestlets <- workTestlets[-which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(is.na(surveyFixMod_itemFit$df.S_X2) == TRUE)], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == 0) != 0 && skipS_X2 == F && (sum(na.omit(surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == 0) < length(1:surveyFixMod@Data$nitems)/2)){
message('\nremoving items df is 0 / removing ', paste(surveyFixMod_itemFit$item[which(surveyFixMod_itemFit$df.S_X2 == 0)]))
workKeys <- workKeys[-which(surveyFixMod_itemFit$df.S_X2 == 0)]
workTestlets <- workTestlets[-which(surveyFixMod_itemFit$df.S_X2 == 0)]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(surveyFixMod_itemFit$df.S_X2 == 0)], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) != 0 && sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) != surveyFixMod@Data$nitems){
message('\nremoving items p.S_X2 is NA / removing ', paste(surveyFixMod_itemFit$item[which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)]))
workKeys <- workKeys[-which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)]
workTestlets <- workTestlets[-which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(is.na(surveyFixMod_itemFit$p.S_X2) == TRUE)], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) == 0 && length(which(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems] >= 3)) != 0 && skipS_X2 == F && (surveyFixMod@Model$nfact == 1 | length(workTestlets) > 0)){ # Drasgow, F., Levine, M. V., Tsien, S., Williams, B., & Mead, A. D. (1995). Fitting polytomous item response theory models to multiple-choice tests. Applied Psychological Measurement, 19(2), 143-166.
message('\nDrasgow, F., Levine, M. V., Tsien, S., Williams, B., & Mead, A. D. (1995) / removing ', paste(surveyFixMod_itemFit$item[which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]))
workKeys <- workKeys[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]
workTestlets <- workTestlets[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else if(sum(is.na(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems])) == 0 && length(which(surveyFixMod_itemFit$p.S_X2[1:surveyFixMod@Data$nitems] < .05)) != 0 && skipS_X2 == F && (surveyFixMod@Model$nfact == 1 | length(workTestlets) > 0)){ # Kang, T., & Chen, T. T. (2008). Performance of the Generalized S?€륺2 Item Fit Index for Polytomous IRT Models. Journal of Educational Measurement, 45(4), 391-406.; Reise, S. P. (1990). A comparison of item- and person-fit methods of assessing model-data fit in IRT. Applied Psychological Measurement, 14, 127-137.
message('\nKang, T., & Chen, T. T. (2008); Reise, S. P. (1990) / removing ', paste(surveyFixMod_itemFit$item[which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]))
workKeys <- workKeys[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]
workTestlets <- workTestlets[-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])]
surveyFixMod <- fastFIFA(surveyFixModRAW[,-which(max(surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems]) == surveyFixMod_itemFit$S_X2[1:surveyFixMod@Data$nitems]/surveyFixMod_itemFit$df.S_X2[1:surveyFixMod@Data$nitems])], itemkeys = workKeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM, testlets = workTestlets, GenRandomPars = GenRandomPars,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
if(needGlobalOptimal == T && forceUIRT == F && length(testlets) == 0){
try(surveyFixMod <- deepFA(surveyFixMod, survey.weights))
}
rm(surveyFixMod_itemFit)
rm(S_X2ErrorFlag)
} else {
itemFitDone <- TRUE
}
} else if (forceRasch == T && surveyFixMod@Model$nfact == 1) {
if(length(which(abs(surveyFixMod_itemFit$z.infit[1:surveyFixMod@Data$nitems]) > 2)) != 0){
message('\nRasch outfit (|z|>2)')
surveyFixMod <- fastFIFA(surveyFixModRAW[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))],
itemkeys = itemkeys[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
} else if(length(which(abs(surveyFixMod_itemFit$z.infit[1:surveyFixMod@Data$nitems]) > 2)) != 0){
message('\nRasch infit (|z|>2)')
surveyFixMod <- fastFIFA(surveyFixModRAW[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))],
itemkeys = itemkeys[,-c(which(max(abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])) == abs(surveyFixMod_itemFit$z.outfit[1:surveyFixMod@Data$nitems])))], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, EnableFMHRM = EnableFMHRM,
fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, MH_BURNIN = MH_BURNIN, MH_SEMCYCLES = MH_SEMCYCLES, ...)
} else {
itemFitDone <- TRUE
}
} else {
itemFitDone <- TRUE
}
} else {
itemFitDone <- TRUE
}
if(printFactorStructureRealtime == T){
message('\nRealtime Factor Structure after iteration')
if(ncol(surveyFixMod@Fit$F)==1){
print(round(surveyFixMod@Fit$F, 2))
} else {
if(bifactorSolution){
print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
} else {
print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
}
}
}
} # the end of while loop
if(autofix == T){ # GET RID OF ASAP
#
# message('\nChecking aberrant responses')
# iteration_num <- iteration_num + 1
# message('Iteration: ', iteration_num, '\n')
#
# surveyFixModRAW <- data.frame(mirt::extract.mirt(surveyFixMod, 'data'))
# surveyFixModCOV <- data.frame(attr(surveyFixMod@ParObjects$lrPars, "df"))
#
# noAberrant <- k.faking(surveyFixModRAW, IRTonly = T, itemkeys = itemkeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM, ...)
# if(length(covdata) == 0){ # anyway, covdata is NULL
# surveyFixMod <- fastFIFA(surveyFixModRAW[which(noAberrant$normal==TRUE),], itemkeys = itemkeys, covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights[which(noAberrant$normal==TRUE)], allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM, ...)
# } else {
# covdata_workout <- surveyFixModCOV
# surveyFixMod <- fastFIFA(surveyFixModRAW[which(noAberrant$normal==TRUE),], itemkeys = itemkeys, covdata = covdata_workout[which(noAberrant$normal==TRUE),], formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights[which(noAberrant$normal==TRUE)], allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM, ...)
# }
if(printFactorStructureRealtime == T){
message('\nRealtime Factor Structure after iteration')
if(ncol(surveyFixMod@Fit$F)==1){
print(round(surveyFixMod@Fit$F, 2))
} else {
if(bifactorSolution){
print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
} else {
print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
}
}
}
# autofix config
fixFactorStructure_Done <- FALSE
surveyFixMod_Workout <- surveyFixMod
while (!fixFactorStructure_Done) { # start of while loop
surveyFixModRAW <- data.frame(surveyFixMod_Workout@Data$data) # update data.frame
surveyFixModCOV <- data.frame(attr(surveyFixMod_Workout@ParObjects$lrPars, "df"))
message('\nFixing Factor Model')
iteration_num <- iteration_num + 1
message('Iteration: ', iteration_num, '\n')
tempG <- surveyFixMod_Workout@Model$itemtype
if(sum(tempG != 'Rasch') != 0){
LowCommunalities <- surveyFixMod_Workout@Fit$h2[which(min(surveyFixMod_Workout@Fit$h2) == surveyFixMod_Workout@Fit$h2)]
}
if(ncol(surveyFixMod_Workout@Fit$F) == 1){
Fmatrix <- surveyFixMod_Workout@Fit$F
} else {
if(bifactorSolution){
Fmatrix <- GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings
} else {
Fmatrix <- GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings
}
}
NoLoadings <- surveyFixMod_Workout@Fit$h2[which(rowSums(abs(round(Fmatrix, 2)) < .4) == ncol(surveyFixMod_Workout@Fit$F))]
# h2 have to >= .3
if(sum(tempG != 'Rasch') == 0){
fixFactorStructure_Done <- TRUE
} else if(LowCommunalities < .3^2){
surveyFixMod_New <- fastFIFA(surveyFixModRAW[,-which(min(surveyFixMod_Workout@Fit$h2) == surveyFixMod_Workout@Fit$h2)], itemkeys = itemkeys[-which(min(surveyFixMod_Workout@Fit$h2) == surveyFixMod_Workout@Fit$h2)], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, ...)
surveyFixMod_Workout <- surveyFixMod_New
if(printFactorStructureRealtime == T){
message('\nRealtime Factor Structure after iteration')
if(ncol(surveyFixMod_Workout@Fit$F)==1){
print(round(surveyFixMod_Workout@Fit$F, 2))
} else {
if(bifactorSolution){
print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
} else {
print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
}
}
}
} else if(length(NoLoadings) != 0){ # noloadings
if(as.logical(length(names(which(NoLoadings == min(NoLoadings))) != 0))){
surveyFixMod_New <- fastFIFA(surveyFixModRAW[,!colnames(surveyFixModRAW) %in% names(which(NoLoadings == min(NoLoadings)))], itemkeys = itemkeys[,!colnames(surveyFixModRAW) %in% names(which(NoLoadings == min(NoLoadings)))], covdata = surveyFixModCOV, formula = formula, SE = SE, SE.type = SE.type, skipNominal = skipNominal, forceGRSM = forceGRSM, assumingFake = assumingFake, masterThesis = masterThesis, forceRasch = forceRasch, unstable = unstable, forceMHRM = forceMHRM, survey.weights = survey.weights, allowMixedResponse = allowMixedResponse, autofix = autofix, forceUIRT = forceUIRT, skipIdealPoint = skipIdealPoint, forceNRM = forceNRM, forceNormalEM = forceNormalEM,
forceDefalutAccelerater = forceDefalutAccelerater, forceDefaultOptimizer = forceDefaultOptimizer, ...)
surveyFixMod_Workout <- surveyFixMod_New
if(printFactorStructureRealtime == T){
message('\nRealtime Factor Structure after iteration')
if(ncol(surveyFixMod_Workout@Fit$F)==1){
print(round(surveyFixMod_Workout@Fit$F, 2))
} else {
if(bifactorSolution){
print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
} else {
print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
}
}
}
}
} else {
fixFactorStructure_Done <- TRUE
}
} # the end of while loop
if(printFactorStructureRealtime == T){
message('\nFinal Factor Structure')
if(ncol(surveyFixMod_Workout@Fit$F)==1){
print(round(surveyFixMod_Workout@Fit$F, 2))
} else {
if(bifactorSolution){
print(round(GPArotation::bifactorQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
} else {
print(round(GPArotation::geominQ(surveyFixMod_Workout@Fit$F, maxit = 10000)$loadings, 2))
}
}
}
return(surveyFixMod_Workout)
} else {
return(surveyFixMod)
}
}
numericMI <- function(model = ..., data = ..., m = 100, fun = 'sem', estimator = 'MLMV', parameterization = 'delta', chi = 'mplus', ...) {
#########################
# Multiple #
# Imputation for #
# variables #
# in SEM context #
#########################
# Seongho Bae #
# seongho@kw.ac.kr #
# May 6th 2016 #
# Under GNU / GPL #
#########################
# checking packages for multiple impuation in SEM context
if(!require('semTools')) {
try(install.packages('semTools', dependencies = TRUE), silent=TRUE)
}
if(!require('Amelia')) {
try(install.packages('Amelia', dependencies = TRUE), silent=TRUE)
}
if(!require('lavaan')) {
try(install.packages('lavaan', dependencies = TRUE), silent=TRUE)
}
# loading packages for multiple imputation in SEM context
library('semTools')
library('Amelia')
library('lavaan')
fit <- sem(model=model, data=data.frame(data)) # extract variable names in model syntax
message("sample size (listwise): ", paste0(nrow(data.frame(fit@Data@X))))
fit_MI <- runMI(model=model, data=data.frame(data[,attributes(fit)$Model@dimNames[[1]][[1]]]), m=m, fun = fun, estimator = estimator, chi = chi, ...)#, control=list(optim.method="L-BFGS-B"), ...)
cat(summary(fit_MI, standardize=T))
print(inspect(fit_MI, 'fit'))
return(fit_MI)
message('inspect(fit, "impute")')
}
cmvFA <- function(x, MHRM = F){
initialModel <- surveyFA(x, autofix = F, forceMHRM = MHRM, bifactorSolution = T)
workModel <- initialModel
STOP <- FALSE
# rotation
while (!STOP) {
try(invisible(gc()), silent = T)
if(ncol(workModel@Fit$F) > 1){
if(workModel@Model$nfact == 2){
rotSumMat <- GPArotation::bifactorQ(workModel@Fit$F, maxit = 1e+6)$loadings[,2]
} else {
rotSumMat <- data.frame(GPArotation::bifactorQ(workModel@Fit$F, maxit = 1e+6)$loadings[,2:workModel@Model$nfact])
}
# evaluation
evaluationMat <- abs(rotSumMat) < .1
print(rotSumMat)
print(evaluationMat)
if(is.data.frame(rotSumMat)){
evaluationMat <- rowSums(evaluationMat)
rotSumMat <- rowSums(abs(rotSumMat))
} else {
evaluationMat <- as.numeric(evaluationMat)
rotSumMat <- (abs(rotSumMat))
}
print(rotSumMat)
print(evaluationMat)
print(sum(evaluationMat))
print(ncol(workModel@Data$data))
if(sum(evaluationMat) == ncol(workModel@Data$data)){ # number of inappropreate items should equal with number of items in model
message('Done')
return(workModel)
} else {
message(paste0(colnames(workModel@Data$data[,which(abs(rotSumMat) == min(abs(rotSumMat)))])))
workModel <- fastFIFA(workModel@Data$data[,-which(abs(rotSumMat) == min(abs(rotSumMat)))], forceMHRM = MHRM)
}
} else { # if one dimensional model
rotSumMat <- workModel@Fit$F
evaluationMat <- abs(rotSumMat) < .1
print(rotSumMat)
print(evaluationMat)
evaluationMat <- as.numeric(evaluationMat)
rotSumMat <- (abs(rotSumMat))
if(sum(evaluationMat) == ncol(workModel@Data$data) | sum(evaluationMat) == 0 | ncol(workModel@Data$data) <= 5){ # number of inappropreate items should equal with number of items in model
message('Done')
return(workModel)
} else {
message(paste0(colnames(workModel@Data$data[,which(abs(rotSumMat) == min(abs(rotSumMat)))])))
workModel <- fastFIFA(workModel@Data$data[,-which(abs(rotSumMat) == min(abs(rotSumMat)))], forceMHRM = MHRM)
}
# return(workModel)
}
}
}
bifactorFA <- function(data = ..., skipS_X2 = F, forceNormalEM = T, forceMHRM = F, covdata = NULL, formula = NULL, skipNominal = F, allowMixedResponse = T, itemkeys = NULL, needGlobalOptimal = T) {
if(length(covdata) != 0){
forceNormalEM <- TRUE
} else if(forceNormalEM){
forceNormalEM <- TRUE
} else {
forceNormalEM <- FALSE
}
mod <- surveyFA(data = data, bifactorSolution = T, skipS_X2 = skipS_X2, forceMHRM = forceMHRM, autofix = F, covdata = covdata, formula = formula, skipNominal = skipNominal, allowMixedResponse = allowMixedResponse, itemkeys = itemkeys, needGlobalOptimal = needGlobalOptimal, forceNormalEM = forceNormalEM, forceUIRT = F)
STOP <- FALSE
while (!STOP) {
if(ncol(mod@Fit$F) == 1){
rotMAT <- abs(mod@Fit$F)
} else {
rotMAT <- abs(GPArotation::bifactorQ(mod@Fit$F, maxit = 1e+5)$loadings)[,1]
}
print(rotMAT)
if(sum(rotMAT < .9999) != ncol(mod@Data$data)){
mod <- surveyFA(data = mod@Data$data[,-which(rotMAT == max(rotMAT))], bifactorSolution = T, skipS_X2 = skipS_X2, forceMHRM = forceMHRM, autofix = F, covdata = covdata, formula = formula, skipNominal = skipNominal, allowMixedResponse = allowMixedResponse, itemkeys = itemkeys[-which(rotMAT == max(rotMAT))], needGlobalOptimal = needGlobalOptimal, forceNormalEM = forceNormalEM, forceUIRT = F)
} else if(sum(rotMAT > .0001) != ncol(mod@Data$data)){
mod <- surveyFA(data = mod@Data$data[,-which(rotMAT == min(rotMAT))], bifactorSolution = T, skipS_X2 = skipS_X2, forceMHRM = forceMHRM, autofix = F, covdata = covdata, formula = formula, skipNominal = skipNominal, allowMixedResponse = allowMixedResponse, itemkeys = itemkeys[-which(rotMAT == max(rotMAT))], needGlobalOptimal = needGlobalOptimal, forceNormalEM = forceNormalEM, forceUIRT = F)
} else {
return(mod)
}
}
}
findMLCA <- function(data = ..., startN = 1, empiricalhist = F, group = NULL, empiricaloptimal = T){
# try(invisible(gc()), silent = T)
DICindices <- vector()
j <- 0
nfact <- vector()
if(sum(psych::describe(data)$range == 0) == 0){
workData <- data
} else {
workData <- data[,-which(psych::describe(data)$range == 0)]
}
message('starting find global optimal of latent class')
for(i in startN:ncol(workData)){
try(invisible(tempModel <- mirt::mdirt(data = workData, model = i, empiricalhist = empiricalhist, group = group, nruns = 100)), silent = F)
if(exists('tempModel')){
if(tempModel@OptimInfo$converged){
message('class: ', i, ' / DIC: ', tempModel@Fit$DIC)
j <- j + 1
DICindices[j] <- tempModel@Fit$DIC
nfact[j] <- i
rm(tempModel)
}
}
}
bestModel <- which(min(DICindices) == DICindices)
message('find global optimal: ', nfact[bestModel])
if(empiricaloptimal){
testFS <- fscores(mirt::mdirt(data = workData, model = nfact[bestModel], empiricalhist = empiricalhist, group = group, nruns = 100), QMC = T)
membership <- vector()
for(i in 1:nrow(testFS)){
membership[i] <- which(testFS[i,] == max(testFS[i,]))
}
message('empirical optimal: ', max(membership))
return(mirt::mdirt(data = workData, model = max(membership), empiricalhist = empiricalhist, group = group, nruns = 100))
} else {
return(mirt::mdirt(data = workData, model = nfact[bestModel], empiricalhist = empiricalhist, group = group, nruns = 100))
}
}
doMLCA <- function(data = ..., startN = 1, empiricalhist = F, group = NULL){
mirtCluster()
if(is.data.frame(data) | is.matrix(data)){
workModel <- findMLCA(data = data, startN = startN, empiricalhist = T, empiricaloptimal = F, group = group)
} else {
workModel <- data
}
initData <- workModel@Data$data
workData <- workModel@Data$data
STOP <- FALSE
while(!STOP){
# item fit evaluation
workModelFit <- mirt::itemfit(workModel, impute = 100, QMC = T)
FitSize <- workModelFit$S_X2/workModelFit$df.S_X2
print(cbind(workModelFit, FitSize))
if(sum(na.omit(workModelFit$S_X2[1:ncol(workData)] == "NaN")) != 0){
workModel <- findMLCA(workData[,-which(workModelFit$S_X2[1:ncol(workData)] == "NaN")], empiricalhist = F, empiricaloptimal = T)
} else if(sum(workModelFit$S_X2[1:ncol(workData)]/workModelFit$df.S_X2[1:ncol(workData)] > 10) != 0){
workModel <- findMLCA(workData[,-which(max(workModelFit$S_X2[1:ncol(workData)]/workModelFit$df.S_X2[1:ncol(workData)]) == workModelFit$S_X2[1:ncol(workData)]/workModelFit$df.S_X2[1:ncol(workData)])], empiricalhist = F, empiricaloptimal = T)
} else {
STOP <- TRUE
}
rm(workData)
workData <- initData[,colnames(workModel@Data$data)] # update work data
}
if(sum(is.na(initData)) != 0){
workData <- mirt::imputeMissing(workModel, fscores(workModel, QMC = T))
workModel <- findMLCA(workData, empiricalhist = F, empiricaloptimal = T)
}
print(plot(workModel, facet_items = FALSE))
print(plot(workModel))
fs <- fscores(workModel, QMC = T)
class_prob <- data.frame(apply(fs, 1, function(x) sample(1:workModel@Model$nfact, 1, prob=x)))
colnames(class_prob) <- "Class"
mirtCluster(remove = T)
return(class_prob)
}
deepFAengine <- function(mirtModel, survey.weights){ # for search more factors with prevent local optimal
DICindices <- vector()
DICindices[1] <- mirtModel@Fit$DIC
j <- 1
if(mirtModel@Options$method == 'EM' && length(attr(mirtModel@ParObjects$lrPars, 'formula')[[1]]) == 0 && length(survey.weights) == 0) {
method <- 'MHRM'
optimizer <- 'NR1'
SE.type <- 'MHRM'
NCYCLES <- 4000
} else if(length(survey.weights) != 0) { # if survey weight included
method <- 'QMCEM'
optimizer <- mirtModel@Options$Moptim
SE.type <- mirtModel@Options$SE.type
NCYCLES <- mirtModel@Options$NCYCLES
} else if(mirtModel@Options$method == 'EM' && length(attr(mirtModel@ParObjects$lrPars, 'formula')[[1]]) != 0) { # if latent regression included
method <- 'QMCEM'
optimizer <- mirtModel@Options$Moptim
SE.type <- mirtModel@Options$SE.type
NCYCLES <- mirtModel@Options$NCYCLES
} else {
method <- mirtModel@Options$method
optimizer <- mirtModel@Options$Moptim
SE.type <- mirtModel@Options$SE.type
NCYCLES <- mirtModel@Options$NCYCLES
}
message('searching global optimal... / estimation method: ', method, ' / optimizer: ', optimizer)
start <- mirtModel@Model$nfact + 1
end <- mirtModel@Model$nfact + 3 # see http://www.tandfonline.com/doi/abs/10.1080/00273171.2012.710386
if(start > ncol(mirtModel@Data$data)){
return(mirtModel)
}
if(end > ncol(mirtModel@Data$data)){
end <- ncol(mirtModel@Data$data)
}
nfact <- vector()
nfact[1] <- mirtModel@Model$nfact
for(i in start:end){
try(invisible(gc()), silent = T)
if(class(mirtModel)[1] == 'MixedClass'){
try(invisible(tempModel <- mirt::mixedmirt(data = mirtModel@Data$data, model = i, itemtype = mirtModel@Model$itemtype, SE = F, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), fixed = mirtModel@Model$formulas$fixed, random = mirtModel@Model$formulas$random, lr.fixed = mirtModel@Model$formulas$lr.fixed, lr.random = mirtModel@Model$formulas$lr.random, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES))), silent = T)
} else {
try(invisible(tempModel <- mirt::mirt(data = mirtModel@Data$data, model = i, itemtype = mirtModel@Model$itemtype, SE = F, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), formula = attr(mirtModel@ParObjects$lrPars, 'formula')[[1]], method = method, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES))), silent = T)
}
if(exists('tempModel')){
if(tempModel@OptimInfo$converged){ # if(tempModel@OptimInfo$converged && tempModel@OptimInfo$secondordertest == T){
message(i, ' factors were converged; DIC: ', tempModel@Fit$DIC)
j <- j + 1
DICindices[j] <- tempModel@Fit$DIC
nfact[j] <- i
rm(tempModel)
}
}
}
bestModel <- which(min(DICindices) == DICindices)
message('find global optimal: ', nfact[bestModel])
if(bestModel == 1){
return(mirtModel)
} else {
if(class(mirtModel)[1] == 'MixedClass'){
try(invisible(tempModel <- mirt::mixedmirt(data = mirtModel@Data$data, model = nfact[bestModel], itemtype = mirtModel@Model$itemtype, SE = F, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), fixed = mirtModel@Model$formulas$fixed, random = mirtModel@Model$formulas$random, lr.fixed = mirtModel@Model$formulas$lr.fixed, lr.random = mirtModel@Model$formulas$lr.random, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES))), silent = T)
} else {
return(mirt::mirt(data = mirtModel@Data$data, model = nfact[bestModel], itemtype = mirtModel@Model$itemtype, SE = T, SE.type = SE.type, covdata = attr(mirtModel@ParObjects$lrPars, 'df'), formula = attr(mirtModel@ParObjects$lrPars, 'formula')[[1]], method = method, optimizer = optimizer, accelerate = mirtModel@Options$accelerate, verbose = F, technical = list(NCYCLES = NCYCLES, MAXQUAD = mirtModel@Options$MAXQUAD, SEtol = mirtModel@Options$SEtol, symmetric = mirtModel@Options$technical$symmetric, removeEmptyRows = mirtModel@Options$removeEmptyRows, MHRM_SE_draws = mirtModel@Options$MHRM_SE_draws, BURNIN = mirtModel@Options$BURNIN, SEMCYCLES = mirtModel@Options$SEMCYCLES)))
}
}
}
deepFA <- function(mirtModel, survey.weights = NULL){
init_nfact <- mirtModel@Model$nfact
if(init_nfact >= ncol(mirtModel@Data$data)){
return(mirtModel)
}
deepModel <- deepFAengine(mirtModel, survey.weights = survey.weights)
if(deepModel@Model$nfact == init_nfact+3){
deepModel <- deepFAengine(deepModel, survey.weights = survey.weights)
if(deepModel@Model$nfact == init_nfact+6){
deepModel <- deepFAengine(deepModel, survey.weights = survey.weights)
if(deepModel@Model$nfact == init_nfact+9){
deepModel <- deepFAengine(deepModel, survey.weights = survey.weights)
} else {
return(deepModel)
}
} else {
return(deepModel)
}
} else {
return(deepModel)
}
}
jmleRaschEst <- function(data, model = 'PCM', fitCriteria = 2, as.LCA = F){
if(!require('mixRasch')){
install.packages('mixRasch', repos = 'http://cran.nexr.com', dependencies = T); library('mixRasch')
}
if(!require('psych')){
install.packages('psych', repos = 'http://cran.nexr.com', dependencies = T); library('mixRasch')
}
itemList <- colnames(data)
RaschData <- data.frame(data)
RaschData <- RaschData[psych::describe(RaschData)$range != 0]
try(invisible(jmleRasch <- mixRasch::mixRasch(data = RaschData, steps = max(psych::describe(RaschData)$range), max.iter = 500, conv.crit = 0.00001, model = model, maxrange = c(-6, 6), maxchange = .25, as.LCA = as.LCA)), silent = T)
STOP <- FALSE
i <- 0
j <- Sys.time()
while (!STOP) {
i <- i+1
if(ncol(data) > 2){
jmleFit <- data.frame(jmleRasch$item.par$in.out)
if(length(which(abs(jmleFit$in.Z) > fitCriteria)) != 0){
message('infit: ', paste(colnames(RaschData)[which(max(abs(jmleFit$in.Z)) == abs(jmleFit$in.Z))]), ' // iteration: ', i)
itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(jmleFit$in.Z)) == abs(jmleFit$in.Z))]
RaschData <- RaschData[!itemList]
itemList <- itemList[-c(which(max(abs(jmleFit$out.Z)) == jmleFit$out.Z))]
try(invisible(jmleRasch <- mixRasch::mixRasch(data = RaschData, steps = max(psych::describe(RaschData)$range), max.iter = 500, conv.crit = 0.00001, model = model, maxrange = c(-6, 6), maxchange = .25, as.LCA = as.LCA)), silent = T)
} else if(length(which(abs(jmleFit$out.Z) > fitCriteria)) != 0){
message('outfit: ', paste(colnames(RaschData)[which(max(abs(jmleFit$out.Z)) == abs(jmleFit$out.Z))]), ' // iteration: ', i)
itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(jmleFit$out.Z)) == abs(jmleFit$out.Z))]
RaschData <- RaschData[!itemList]
itemList <- itemList[-c(which(max(abs(jmleFit$out.Z)) == jmleFit$out.Z))]
try(invisible(jmleRasch <- mixRasch::mixRasch(data = RaschData, steps = max(psych::describe(RaschData)$range), max.iter = 500, conv.crit = 0.00001, model = model, maxrange = c(-6, 6), maxchange = .25, as.LCA = as.LCA)), silent = T)
} else {
k <- Sys.time()
STOP <- TRUE
}
} else {
k <- Sys.time()
STOP <- TRUE
}
}
RaschJMLEresult <- new.env()
RaschJMLEresult$itemList <- colnames(RaschData)
RaschJMLEresult$rawData <- RaschData
RaschJMLEresult$model <- jmleRasch
RaschJMLEresult$TotalTime <- k-j
return(RaschJMLEresult)
}
cmleRaschEst <- function(data, model = 'PCM'){
if(!require('eRm')){
install.packages('eRm', repos = 'http://cran.nexr.com', dependencies = T); library('eRm')
}
itemList <- colnames(data)
RaschData <- data.frame(data)
if(model == 'PCM'){
cmleRasch <- eRm::PCM(RaschData)
} else if(model == 'RSM'){
cmleRasch <- eRm::RSM(RaschData)
} else if(model == 'dich'){
cmleRasch <- eRm::RM(RaschData)
} else {
stop('model have to define PCM, RSM, dich')
}
# define while loop
STOP <- FALSE # initial is STOP = FALSE
while (!STOP) {
if(ncol(data) > 2){
# cmleFit <- data.frame(cmleRasch$item.par$in.out)
cmleFit <- data.frame(eRm::itemfit(eRm::person.parameter(cmleRasch))$i.outfitZ,
eRm::itemfit(eRm::person.parameter(cmleRasch))$i.outfitMSQ,
eRm::itemfit(eRm::person.parameter(cmleRasch))$i.infitZ,
eRm::itemfit(eRm::person.parameter(cmleRasch))$i.infitMSQ)
colnames(cmleFit) <- c("out.Z", "out.MSQ", "in.Z", "in.MSQ")
print(cmleFit)
if(length(which(abs(cmleFit$out.Z) > 2)) != 0){
message('outfit: ', colnames(RaschData)[which(max(abs(cmleFit$out.Z)) == abs(cmleFit$out.Z))])
itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(cmleFit$out.Z)) == abs(cmleFit$out.Z))]
RaschData <- RaschData[!itemList]
itemList <- itemList[-c(which(max(abs(cmleFit$out.Z)) == cmleFit$out.Z))]
if(model == 'PCM'){
cmleRasch <- eRm::PCM(RaschData)
} else if(model == 'RSM'){
cmleRasch <- eRm::RSM(RaschData)
} else if(model == 'dich'){
cmleRasch <- eRm::RM(RaschData)
} else {
stop('model have to define PCM, RSM, dich')
}
} else if(length(which(abs(cmleFit$in.Z) > 2)) != 0){
message('infit: ', colnames(RaschData)[which(max(abs(cmleFit$in.Z)) == abs(cmleFit$in.Z))])
itemList <- names(RaschData) %in% colnames(RaschData)[which(max(abs(cmleFit$in.Z)) == abs(cmleFit$in.Z))]
RaschData <- RaschData[!itemList]
itemList <- itemList[-c(which(max(abs(cmleFit$out.Z)) == cmleFit$out.Z))]
if(model == 'PCM'){
cmleRasch <- eRm::PCM(RaschData)
} else if(model == 'RSM'){
cmleRasch <- eRm::RSM(RaschData)
} else if(model == 'dich'){
cmleRasch <- eRm::RM(RaschData)
} else {
stop('model have to define PCM, RSM, dich')
}
} else {
STOP <- TRUE
}
} else {
STOP <- TRUE
}
}
Raschcmleresult <- new.env()
Raschcmleresult$itemList <- colnames(RaschData)
Raschcmleresult$model <- cmleRasch
return(Raschcmleresult)
}
autoMCMC2PL.ml <- function(x = NULL, group = NULL, est.b.M="h", est.b.Var="i",
est.a.M="h" , est.a.Var="i", burnin = 10000,
iter = 20000, Rhat = 1.05, autofix = T, TargetTestLength = 3, considerSigma = T, # for 2PNO Multilevel
testlets = NULL, survey.weights = NULL, est.slope = T, est.guess = T, verbose = F # for 1-3PNO testlet
){
if(!require('sirt')){
install.packages('sirt')
library('sirt')
}
if(!require('plyr')){
install.packages('plyr')
library('plyr')
}
if(!require('progress')){
install.packages('progress')
library('progress')
}
if(range(x[1], na.rm = T)[2] - range(x[1], na.rm = T)[1] > 1){
link <- 'normal'
} else {
link <- 'logit'
}
if(burnin > iter){
iter <- burnin*2
}
if(verbose){
num.progress.iter <- burnin/10
} else {
num.progress.iter <- burnin + iter + 1
}
initData <- x
iterationTrials <- 1
StartTime <- Sys.time()
if(length(group) != 0){
initData <- initData[which(!is.na(group)),]
group <- group[which(!is.na(group))]
} else {
}
if(is.character(group)){
group <- as.factor(group)
}
if(is.factor(group)){
numberOfGroups <- length(attributes(group)$levels)
} else {
numberOfGroups <- length(attributes(as.factor(group))$levels)
}
message('MCMC link: ', link)
message('MCMC Trials: ', iterationTrials)
message('Current number of items: ', ncol(initData))
message('Current number of groups: ', numberOfGroups)
message('Rhat cutoff: ', Rhat)
group <- as.integer(group)
if(length(testlets) != 0){
ActualTestlets <- testlets
TestletStatus <- (table(ActualTestlets))
if(sum(is.na(ActualTestlets)) == 0){
ActualTestlets <- plyr::mapvalues(ActualTestlets, as.numeric(names(TestletStatus)[which(TestletStatus == max(TestletStatus))]), NA) # CTC(M-1)
}
TestletStatus <- (table(ActualTestlets)) # renew
if(length(which(TestletStatus == 1)) != 0){
ActualTestlets <- plyr::mapvalues(ActualTestlets, names(TestletStatus)[which(TestletStatus == 1)], rep(NA, length(which(TestletStatus == 1))))
}
# if(sum(ActualTestlets %in% 1) == 0){
# ActualTestlets <- ActualTestlets - min(ActualTestlets, na.rm = T) + 1
# }
if(sum(is.na(ActualTestlets)) == length(ActualTestlets)){
ActualTestlets <- NULL # disable testlet when all testlet information were NA
} else {
ActualTestlets <- plyr::mapvalues(ActualTestlets, as.numeric(attributes(as.factor(ActualTestlets))$levels), seq(length(attributes(as.factor(c(ActualTestlets)))$levels))) # convert text, keep same interval
}
} else {
ActualTestlets <- NULL
}
message('initializing model...')
if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'normal'){
init <- surveyFA(data = initData, survey.weights = survey.weights, testlets = ActualTestlets)
return(init)
} else {
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
}
message('starting calibration...')
pb <- progress_bar$new(
format = " processing MCMC chains [:bar] :percent (:current of :total) in :elapsed, eta: :eta",
total = ncol(initData) - TargetTestLength, clear = FALSE, width= 120, force = T)
STOP <- FALSE
while(!STOP){
ticktockUpdate <- 1-((sum(init$summary.mcmcobj$Rhat > Rhat) +
sum(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)] < 0) +
sum(cbind( init$summary.mcmcobj$Q2.5[grep("^a",init$summary.mcmcobj$parameter)] < 0 && init$summary.mcmcobj$Q97.5[grep("^a",init$summary.mcmcobj$parameter)] > 0 ))
) / length(init$summary.mcmcobj$Rhat))
if(!is.na(ticktockUpdate)){
pb$update(ticktockUpdate)
}
pb$tick()
invisible(try(gc(), silent = T))
if(sum(init$summary.mcmcobj$Rhat > Rhat) != 0){
parmList <- init$summary.mcmcobj
# if(length(grep("^sigma[0-9]", as.character(init$summary.mcmcobj$parameter))) != 0){
# parmList <- parmList[-grep("^sigma[0-9]", as.character(parmList$parameter)),]
# }
if(!considerSigma){
parmList <- parmList[-grep("^sigma", as.character(parmList$parameter)),]
}
excludeVar <- unique(na.omit(as.numeric(unlist(strsplit(unlist(as.character(parmList[which(max(parmList$Rhat) == parmList$Rhat),]$parameter)), "[^0-9]+")))))
if(length(excludeVar) != 0 && ncol(initData) > TargetTestLength){
if(verbose) message('Removing a item ', names(initData[excludeVar]),' / ', parmList[which(max(parmList$Rhat) == parmList$Rhat),]$parameter, ' Rhat: ', max(parmList$Rhat))
initData <- initData[,-excludeVar]
if(length(testlets) != 0){
ActualTestlets <- ActualTestlets[-excludeVar]
ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))
}
iterationTrials <- iterationTrials+1
if(verbose) message('MCMC Trials: ', iterationTrials)
if(verbose) message('Current number of items: ', ncol(initData))
if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else {
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
}
} else {
pb$update(1)
pb$tick()
STOP <- TRUE
}
} else if(sum(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)] < 0) != 0 && autofix){
excludeVar <- unique(na.omit(as.numeric(unlist(strsplit(unlist(as.character(init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))])), "[^0-9]+")))))
if(length(excludeVar) != 0 && ncol(initData) > TargetTestLength){
if(verbose) message('Removing a item ', names(initData[excludeVar]),' / ', init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))], ' value: ', min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]))
initData <- initData[,-excludeVar]
if(length(testlets) != 0){
ActualTestlets <- ActualTestlets[-excludeVar]
ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))
}
iterationTrials <- iterationTrials+1
if(verbose) message('MCMC Trials: ', iterationTrials)
if(verbose) message('Current number of items: ', ncol(initData))
if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else {
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
}
} else {
pb$update(1)
pb$tick()
STOP <- TRUE
}
} else if(
sum(cbind( init$summary.mcmcobj$Q2.5[grep("^a",init$summary.mcmcobj$parameter)] < 0 && init$summary.mcmcobj$Q97.5[grep("^a",init$summary.mcmcobj$parameter)] > 0 )) > 0
&& autofix){
excludeVar <- unique(na.omit(as.numeric(unlist(strsplit(unlist(as.character(init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))])), "[^0-9]+")))))
if(length(excludeVar) != 0 && ncol(initData) > TargetTestLength){
if(verbose) message('Removing a item ', names(initData[excludeVar]),' / ', init$summary.mcmcobj$parameter[which(min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]) == (init$summary.mcmcobj$MAP))], ' value: ', min(init$summary.mcmcobj$MAP[grep("^a",init$summary.mcmcobj$parameter)]))
initData <- initData[,-excludeVar]
if(length(testlets) != 0){
ActualTestlets <- ActualTestlets[-excludeVar]
ActualTestlets <- plyr::mapvalues(ActualTestlets, names(which(table(ActualTestlets) == 1)), rep(NA, length(names(which(table(ActualTestlets) == 1)))))
}
iterationTrials <- iterationTrials+1
if(verbose) message('MCMC Trials: ', iterationTrials)
if(verbose) message('Current number of items: ', ncol(initData))
if(length(group) == 0 && length(ActualTestlets) == 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, est.slope = est.slope, weights = survey.weights, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else if(length(group) == 0 && length(ActualTestlets) != 0 && link == 'logit'){
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.3pno.testlet(dat = initData, testlets = ActualTestlets, weights = survey.weights, est.slope = est.slope, est.guess = est.guess, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
} else {
if(exists('init')){
try(rm(init))
}
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
if(!exists('init')){
while(!exists('init')){
try(init <- sirt::mcmc.2pno.ml(dat = initData, group = group, link = link, est.b.M=est.b.M, est.b.Var=est.b.Var , est.a.M=est.a.M, est.a.Var=est.a.Var, burnin = burnin, iter = iter, N.sampvalues = iter, progress.iter = num.progress.iter), silent = T)
}
}
}
} else {
pb$update(1)
pb$tick()
STOP <- TRUE
}
} else { # overall stop rule when all items normally terminated
pb$update(1)
pb$tick()
STOP <- TRUE
}
}
EndTime <- Sys.time()
ReturnList <- list(Solution = init, testlet = ActualTestlets, rawData = initData, TotalTime = EndTime - StartTime)
return(ReturnList)
}
testAssembly <- function(MIRTmodel, measurementArea, NumberOfForms = 1, meanOfdifficulty = 0, sdOfdifficulty = 2.0, numberOfItems = 16, maximumItemSelection = 1, oldFormYMIRTmodel = NULL, SCLmethod = 'Haebara', oldFormYCommonItemNumber = NULL, newFormXCommonItemNumber = NULL){
if(!require('xxIRT')){
install.packages('xxIRT')
library('xxIRT')
}
if(!require('mirt')){
install.packages('mirt')
library('mirt')
}
if(!require('sirt')){
install.packages('sirt')
library('sirt')
}
library(xxIRT)
library(mirt)
library('sirt')
if(class(MIRTmodel)[1] == 'SingleGroupClass'){
message('mirt model provided')
IRTpars <- data.frame(coef(MIRTmodel, IRTpars = T, simplify = T)$items[,1:3])
colnames(IRTpars)[3] <- 'c'
} else if(class(MIRTmodel)[1] == 'mcmc.sirt'){
message('sirt model provided')
if(length(grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]) == 0){
IRTpars <- data.frame(rep(1,ncol(MIRTmodel$dat)),
MIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]])
} else {
IRTpars <- data.frame(MIRTmodel$summary.mcmcobj$MAP[grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]],
MIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(MIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(MIRTmodel$summary.mcmcobj$parameter))]])
}
colnames(IRTpars) <- c('a', 'b')
rownames(IRTpars) <- colnames(MIRTmodel$dat)
if(length(MIRTmodel$summary.mcmcobj$MAP[grep("^c", MIRTmodel$summary.mcmcobj$parameter)]) != 0){
IRTpars$c <- MIRTmodel$summary.mcmcobj$MAP[grep("^c", MIRTmodel$summary.mcmcobj$parameter)]
} else {
IRTpars$c <- rep(0, ncol(MIRTmodel$dat))
}
} else {
stop('Did you provide right mirt or sirt mcmc model?')
}
if(length(oldFormYMIRTmodel) != 0){ # if SCL activated
if(length(oldFormYCommonItemNumber) == 0 | length(newFormXCommonItemNumber) == 0 | length(newFormXCommonItemNumber) != length(oldFormYCommonItemNumber)) {
STOP('Please provide correspond oldFormYCommonItemNumber and newFormXCommonItemNumber')
}
# check SCL software
if(!require(SNSequate)){
install.packages('SNSequate')
library('SNSequate')
}
# read IRT parameters of oldform Y
if(class(oldFormYMIRTmodel)[1] == 'SingleGroupClass'){
message('mirt model provided')
IRTpars2 <- data.frame(coef(oldFormYMIRTmodel, IRTpars = T, simplify = T)$items[,1:3])
colnames(IRTpars2)[3] <- 'c'
} else if(class(oldFormYMIRTmodel)[1] == 'mcmc.sirt'){
message('sirt model provided')
if(length(grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]) == 0){
IRTpars2 <- data.frame(rep(1,ncol(oldFormYMIRTmodel$dat)),
oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]])
} else {
IRTpars2 <- data.frame(oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^a",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^a_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]],
oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))[!grep("^b",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter)) %in% grep("^b_marg",as.character(oldFormYMIRTmodel$summary.mcmcobj$parameter))]])
}
colnames(IRTpars2) <- c('a', 'b')
rownames(IRTpars2) <- colnames(oldFormYMIRTmodel$dat)
if(length(oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^c", oldFormYMIRTmodel$summary.mcmcobj$parameter)]) != 0){
IRTpars2$c <- oldFormYMIRTmodel$summary.mcmcobj$MAP[grep("^c", oldFormYMIRTmodel$summary.mcmcobj$parameter)]
} else {
IRTpars2$c <- rep(0, ncol(oldFormYMIRTmodel$dat))
}
} else {
stop('Did you provide right mirt or sirt mcmc model?')
}
# split common item parameters
newCommonXIRTpars <- IRTpars[newFormXCommonItemNumber,]
oldCommonYIRTpars <- IRTpars2[oldFormYCommonItemNumber,]
# do calcluate linking coefficients slope A and constant B
LinkingCoefficients <- irt.link(as.data.frame(cbind(oldCommonYIRTpars, newCommonXIRTpars)), common = 1:length(oldFormYCommonItemNumber), model = '3PL', icc = 'logistic', D = 1.702)
names(LinkingCoefficients$Haebara) <- c('A', 'B')
names(LinkingCoefficients$StockLord) <- c('A', 'B')
if(SCLmethod == 'Haebara'){
message('applying Haebara method...')
A <- LinkingCoefficients$Haebara[1]
B <- LinkingCoefficients$Haebara[2]
} else if(SCLmethod == 'StockingLord'){
message('applying Stocking-Lord method...')
A <- LinkingCoefficients$StockLord[1]
B <- LinkingCoefficients$StockLord[2]
} else {
stop('please provide "Haebara" or "StockingLord" in SCLmethod argument correctly')
}
# adjusting new test IRT parameters to value of old test item parameters
IRTpars$a <- IRTpars$a/A
IRTpars$b <- A*IRTpars$b+B
# ignoring cx=cy, that's not require of adjustments.
}
items <- IRTpars
items$content <- as.numeric(as.factor(measurementArea))
# print(items)
# try -2 to 2 (flat information)
x <- xxIRT::ata(items, NumberOfForms, len = numberOfItems, maxselect = maximumItemSelection) %>%
xxIRT::ata_obj_relative(seq(min(items$b), max(items$b), .1), "max", flatten=0.1, negative = T) %>%
xxIRT::ata_solve()
try(ataPlot <- plot(x), silent = T)
if(exists('ataPlot')){
print(plot(x))
} else {
# -1 to 1 (flat information)
x <- xxIRT::ata(items, NumberOfForms, len = numberOfItems, maxselect = maximumItemSelection) %>%
xxIRT::ata_obj_relative(seq((min(items$b)+0.5), (max(items$b)-0.5), .1), "max", flatten=0.1, negative = T) %>%
xxIRT::ata_solve()
try(ataPlot <- plot(x), silent = T)
if(exists('ataPlot')){
print(plot(x))
} else {
rm(x)
x <- xxIRT::ata(items, NumberOfForms, len = numberOfItems, maxselect = maximumItemSelection)
x <- xxIRT::ata_obj_absolute(x, "b", meanOfdifficulty * numberOfItems)
x <- xxIRT::ata_obj_absolute(x, (x$pool$b - meanOfdifficulty)^2, sdOfdifficulty * numberOfItems)
x <- xxIRT::ata_solve(x)
print(plot(x))
}
}
y <- xxIRT::ata_get_items(x, as.list=TRUE) ## FIXME
# get ATA data
ATAFormData <- list()
ATAFormModel <- list()
ATAFormModelValues <- list()
for(i in 1:NROW(y)){
message('-------------------------------------------------')
message('mean of difficulty of Form ', i,': ', mean(y[[i]]$b))
message('sd of difficulty of Form ', i,': ', sd(y[[i]]$b))
message('min of difficulty of Form ', i,': ', min(y[[i]]$b))
message('max of difficulty of Form ', i,': ', max(y[[i]]$b))
message('-------------------------------------------------\n')
}
for(i in 1:NROW(y)){
if(class(MIRTmodel)[1] == 'SingleGroupClass'){
ATAFormData[[i]] <- data.frame(MIRTmodel@Data$data)
ATAFormData[[i]] <- ATAFormData[[i]][,colnames(ATAFormData[[i]]) %in% rownames(y[[i]])]
} else if(class(MIRTmodel)[1] == 'mcmc.sirt'){
ATAFormData[[i]] <- data.frame(MIRTmodel$dat)
ATAFormData[[i]] <- ATAFormData[[i]][colnames(ATAFormData[[i]]) %in% rownames(y[[i]])]
}
if(ncol(ATAFormData[[i]]) != 0 && sum(psych::describe(ATAFormData[[i]])$range == 1) == ncol(ATAFormData[[i]])){
ATAFormModelValues[[i]] <- mirt::mirt(data = ATAFormData[[i]], model = 1, itemtype = '3PL', pars = 'values')
ATAFormModelValues[[i]][which(ATAFormModelValues[[i]]$name == 'a1'),"value"] <- y[[i]]$a
ATAFormModelValues[[i]][which(ATAFormModelValues[[i]]$name == 'd'),"value"] <- -1*y[[i]]$a*y[[i]]$b
ATAFormModelValues[[i]][which(ATAFormModelValues[[i]]$name == 'g'),"value"] <- y[[i]]$c
ATAFormModelValues[[i]]$est <- FALSE
ATAFormModel[[i]] <- mirt::mirt(data = ATAFormData[[i]], model = 1, itemtype = '3PL', pars = ATAFormModelValues[[i]])
names(ATAFormModel)[[i]] <- paste0('form', i)
} else {
# if(){
message('warning: items were polytomous? or variable names are contain spaces?') ## FIXME: ADD 3PLNRM
# } else {
# message('warning: ')
# }
# NULL
}
}
z <- list(OriginalParms = items, ATAforms = y, ATAFormModel = ATAFormModel, ATAFormSolutions = x)
return(z)
}
fastBifactorCFA <- function(x, ga = F, itemkeys = NULL, initSolution = F, skipNRM = T, excludeUnscalableVar = F, lowerbound = .5, covdata = NULL, formula = NULL){
if(!require('mokken')){
install.packages('mokken')
library('mokken')
}
if(!require('mirt')){
install.packages('mirt')
library('mirt')
}
x <- x[psych::describe(x)$range != 0]
message('finding testlet structures using mokken scale analysis...')
if(ga){
message('mutuation probability: ',1/log(10^(log2(ncol(x)/5)) * 1000)*(log(log2(nrow(x)))+log(log(log2(nrow(x))))))
if(length(itemkeys) != 0){
try(modMokken <- mokken::aisp(data.frame(mirt::key2binary(data.frame(x), key = itemkeys)), lowerbound = lowerbound, verbose = T, search = 'ga', pxover = 1, pmutation = 1/log(10^(log2(ncol(x)/5)) * 1000)*(log(log2(nrow(x)))+log(log(log2(nrow(x))))), popsize = log2(nrow(x))), silent = T)
} else {
try(modMokken <- mokken::aisp(data.frame(x), lowerbound = lowerbound, verbose = T, search = 'ga', pxover = 1, pmutation = 1/log(10^(log2(ncol(x)/5)) * 1000)*(log(log2(nrow(x)))+log(log(log2(nrow(x))))), popsize = log2(nrow(x))), silent = T)
}
} else {
if(length(itemkeys) != 0){
try(modMokken <- mokken::aisp(data.frame(mirt::key2binary(data.frame(x), key = itemkeys)), lowerbound = lowerbound, verbose = T), silent = T)
} else {
try(modMokken <- mokken::aisp(data.frame(x), verbose = T, lowerbound = lowerbound), silent = T)
}
}
print(modMokken)
if(sum(modMokken == 0) == ncol(x) | sum(modMokken == 1) == ncol(x)){
# print(modMokken)
if(!initSolution){
modBfactor <- surveyFA(data = data.frame(x), itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich')
} else {
modBfactor <- deepFA(fastFIFA(x = data.frame(x), itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich'))
}
} else {
# print(modMokken)
message('testlet structure was found!')
if(excludeUnscalableVar){
message('excluding unscalable varaiables...')
testDat <- data.frame(x)[-which(modMokken == 0)]
modMokken <- modMokken[-which(modMokken == 0)]
modMokken[which(modMokken %in% as.numeric(names(which(table(modMokken) < 2))))] <- NA
# print(modMokken)
} else {
testDat <- data.frame(x)
modMokken[which(modMokken == 0)] <- NA
modMokken[which(modMokken %in% as.numeric(names(which(table(modMokken) < 2))))] <- NA
# print(modMokken)
}
message('number of testing variables: ', (ncol(testDat)), ' (', round(ncol(testDat)/ncol(data.frame(x))*100, digits = 2), '%)')
if(!initSolution){
modBfactor <- surveyFA(data = testDat, testlets = modMokken, itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich')
} else {
modBfactor <- fastFIFA(x = testDat, testlets = modMokken, itemkeys = itemkeys, skipNominal = skipNRM, covdata = covdata, formula = formula, SE.type = 'sandwich')
}
}
return(modBfactor)
}
bfactor2mod <- function(model, J){ # borrow from mirt::bfactor2mod
tmp <- tempfile('tempfile')
unique <- sort(unique(model))
index <- 1L:J
tmp2 <- c()
for(i in 1L:length(unique)){
ind <- na.omit(index[model == unique[i]])
comma <- rep(',', 2*length(ind))
TF <- rep(c(TRUE,FALSE), length(ind))
comma[TF] <- ind
comma[length(comma)] <- ""
tmp2 <- c(tmp2, c(paste('\nS', i, ' =', sep=''), comma))
}
cat(tmp2, file=tmp)
model <- mirt::mirt.model(file=tmp, quiet = TRUE)
unlink(tmp)
return(model)
}
doBfactor2mod <- function(mirtDat, testlet){
tmp <- any(sapply(colnames(mirtDat), grepl, x=paste0('G = 1-', ncol(mirtDat)))) # borrow from mirt::bfactor
model2 <- mirt::mirt.model(paste0('G = 1-', ncol(mirtDat)), itemnames = if(tmp) colnames(mirtDat) else NULL)
model <- bfactor2mod(testlet, length(testlet))
model$x <- rbind(model2$x, model$x)
return(model)
}
# item factor analysis based LCA
getmode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
findLatentClass <- function(data = ..., nruns = 1, maxClasses = NULL, covdata = NULL, formula = NULL, SE.type = 'sandwich', checkSecondOrderTest = T, empiricalhist = T, turnOnMIRTcluster = F){
if(!require('mirt')){
install.packages('mirt')
library('mirt')
}
if(!require('progress')){
install.packages('progress')
library('progress')
}
if(turnOnMIRTcluster){
try(mirtCluster(remove = T))
try(mirtCluster(spec = round(parallel::detectCores()/2)))
}
modelFit <- list()
if(is.null(maxClasses)){
maxClasses <- round(ncol(data)/5)
}
pb <- progress_bar$new(
format = "(:spin) [:bar] :percent (:current of :total) in :elapsed, eta: :eta",
total = maxClasses, clear = FALSE, width = 60)
for(i in 1:maxClasses){
invisible(try(testModel <- mirt::mdirt(data = data, model = i, SE = checkSecondOrderTest, verbose = F, nruns = nruns, covdata = covdata, formula = formula, SE.type = SE.type, empiricalhist = empiricalhist), silent = T))
pb$tick()
if(exists('testModel')){
if(checkSecondOrderTest){
if(testModel@OptimInfo$converged && testModel@OptimInfo$secondordertest){
# message(round(i/maxClasses*100, 1), "% complete", ' (', i,' / ', maxClasses, ')')
modelFit[[i]] <- testModel@Fit
}
rm(testModel)
} else {
if(testModel@OptimInfo$converged){
# message(round(i/maxClasses*100, 1), "% complete", '(', i,' / ', maxClasses, ')')
modelFit[[i]] <- testModel@Fit
}
rm(testModel)
}
}
}
# return(modelFit)
modelFitMatrix <- (matrix(unlist(modelFit), ncol = 14, byrow = TRUE))
modelFitNumbers <- vector()
for(i in 1:NROW(modelFit)){
if(is.null(modelFit[[i]]) == FALSE){
modelFitNumbers[length(modelFitNumbers)+1] <- i
}
}
rownames(modelFitMatrix) <- modelFitNumbers
colnames(modelFitMatrix) <- names(modelFit[[1]])
return(modelFitMatrix)
}
autoLCA <- function(data = ..., UIRT = T, nruns = 1, covdata = NULL, formula = NULL, SE.type = 'sandwich', forceMHRM = F){
# source('https://github.com/seonghobae/k.aefa/raw/master/aFIPC.R')
if(length(covdata) != 0 && SE.type != 'complete'){
SE.type <- 'complete'
}
testMIRTmod <- surveyFA(data = data, forceMHRM = forceMHRM, forceUIRT = UIRT, covdata = covdata, formula = formula, SE.type = SE.type)
testNumberOfClasses <- findLatentClass(data = testMIRTmod@Data$data, nruns = nruns, covdata = covdata, formula = formula, SE.type = 'complete') # sandwich not yet supported
LCA_Judgement <- vector()
for(j in c(7:11)){
if(sum(is.na(testNumberOfClasses[,j])) > 0){
try(LCA_Judgement[j] <- NA)
} else {
try(LCA_Judgement[j] <- (which(testNumberOfClasses[,j] == min(testNumberOfClasses[,j]))))
}
}
try(FinalModel <- mirt::mdirt(testMIRTmod@Data$data, as.numeric(rownames(testNumberOfClasses)[getmode(na.omit(LCA_Judgement))]), nruns = nruns, verbose = F, SE = T, covdata = covdata, formula = formula, SE.type = 'complete'), silent = T)
return(list(IRTmodel = testMIRTmod, LCAdecisionTable = testNumberOfClasses, FinalModel = FinalModel))
}
MLIRT <- function(data = ..., covdata = ..., model = ..., itemtype = 'Rasch', fixed = ~-1, random = NULL, lr.fixed = ~1, lr.random = NULL, SEtol = 1e-4, SE = T, TOL = .001, NCYCLES = 2000, ...){
if(!require('mirt')){
install.packages('mirt', repos = 'https://cran.biodisk.org')
}
MLM_calibration <- mirt::mixedmirt(data = data, covdata = covdata, model = model, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, itemtype = itemtype, SE = SE, SE.type = 'MHRM', TOL = TOL, technical = list(SEtol = SEtol, removeEmptyRows = TRUE, NCYCLES = NCYCLES))
return(MLM_calibration)
# }
}
EEMEIRT <- function(data = ..., covdata = ..., itemtype = ..., SE = T, fixed = ~-1, random = NULL, lr.fixed = ~1, lr.random = NULL, TOL = .001, SEtol = 1e-03, NCYCLES = 2000){
testDIC <- vector()
noConverge <- 0
for(i in 1:ncol(data)){
if(i > 10 && noConverge > 5){
message(which(testDIC == min(testDIC, na.rm = T)))
testMLM <- MLIRT(data = data, covdata = covdata, model = which(testDIC == min(testDIC, na.rm = T)),
itemtype = itemtype, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random,
TOL = TOL, SEtol = SEtol, NCYCLES = NCYCLES)
return(testMLM)
}
try(testMLM <- MLIRT(data = data, covdata = covdata, model = i, SE = F, itemtype = itemtype, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, TOL = TOL, SEtol = SEtol, NCYCLES = NCYCLES), silent = F)
if(exists('testMLM')){
if(testMLM@OptimInfo$converged){
testDIC[i] <- testMLM@Fit$DIC
noConverge <- 0
try(rm(testMLM))
} else {
testDIC[i] <- NA
noConverge <- noConverge+1
try(rm(testMLM))
}
} else {
testDIC[i] <- NA
}
}
if(sum(is.na(testDIC)) == length(testDIC)){
stop('no solution')
} else {
message(which(testDIC == min(testDIC, na.rm = T)))
testMLM <- MLIRT(data = data, covdata = covdata, model = which(testDIC == min(testDIC, na.rm = T)), itemtype = itemtype, fixed = fixed, random = random, lr.fixed = lr.fixed, lr.random = lr.random, TOL = TOL, SEtol = SEtol, NCYCLES = NCYCLES)
return(testMLM)
}
}
rotateEMEIRT <- function(EMEIRTmodel, rotate = 'bifactorQ', suppress = 0){
if(!require('mirt')){
install.packages('mirt')
library('mirt')
}
# recursive formula
message(EMEIRTmodel@Model$model)
MLM_rotate_formula_mod <- mirt::mirt(data = EMEIRTmodel@Data$data, model = EMEIRTmodel@Model$model, itemtype = EMEIRTmodel@Model$itemtype, pars = 'values')
MLM_rotate_formula_mod_original <- mirt::mod2values(EMEIRTmodel)
if(sum(MLM_rotate_formula_mod_original$name == '(Intercept)') != 0){
MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original[!MLM_rotate_formula_mod_original$name == '(Intercept)',]
}
# MLM_rotate_formula_mod_original <- MLM_rotate_formula_mod_original
MLM_rotate_formula_mod$value[which(MLM_rotate_formula_mod$item %in% colnames(EMEIRTmodel@Data$data))] <- MLM_rotate_formula_mod_original$value[which(MLM_rotate_formula_mod_original$item %in% colnames(EMEIRTmodel@Data$data))]
MLM_rotate_formula_mod$est <- F
# print(MLM_rotate_formula_mod)
MLM_rotate_formula_mod_final <- mirt::mirt(data = EMEIRTmodel@Data$data, model = EMEIRTmodel@Model$model, itemtype = EMEIRTmodel@Model$itemtype, pars = MLM_rotate_formula_mod, method = 'QMCEM')
if(ncol(MLM_rotate_formula_mod_final@Fit$F) > 1){
if(rotate == 'bifactorQ'){
return(GPArotation::bifactorQ(MLM_rotate_formula_mod_final@Fit$F, maxit = 1e+7))
} else if(rotate == 'geominQ'){
return(GPArotation::geominQ(MLM_rotate_formula_mod_final@Fit$F, maxit = 1e+7))
} else if(rotate == 'bifactorT'){
return(GPArotation::bifactorT(MLM_rotate_formula_mod_final@Fit$F, maxit = 1e+7))
}else if(rotate == 'none'){
summary(MLM_rotate_formula_mod_final, rotate = 'none', suppress = suppress)
}
}else{
summary(MLM_rotate_formula_mod_final, rotate = 'none', suppress = suppress)
} #
}
# rotateEMEIRT(silenceMLM, rotate = 'bifactorT', suppress = .05)
doLCA <- function(data = ..., SE.type = 'Oakes', checkSecondOrderTest = T, nruns = 1, maxClasses = NULL, empiricalhist = T, covdata = NULL, formula = NULL){
if(!require('psych')){
install.packages('psych')
library('psych')
}
if(is.data.frame(data)){
datd <- data
} else {
datd <- data.frame(datd)
}
datd <- datd[psych::describe(datd)$range != 0] # prevent range = 0
STOP_LCA <- FALSE
while(!STOP_LCA){
try(invisible(gc()), silent = T) # garbage cleaning
preLCAoptimal <- findLatentClass(datd, SE.type = SE.type, checkSecondOrderTest = checkSecondOrderTest, nruns = nruns, maxClasses = maxClasses, empiricalhist = empiricalhist, covdata = covdata, formula = formula)
preLCAoptimal <- data.frame(preLCAoptimal)
optimalDecisionAIC <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$AIC == min(preLCAoptimal$AIC))]
optimalDecisionAICc <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$AICc == min(preLCAoptimal$AICc))]
optimalDecisionSABIC <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$SABIC == min(preLCAoptimal$SABIC))]
optimalDecisionDIC <- as.numeric(rownames(preLCAoptimal))[which(preLCAoptimal$DIC == min(preLCAoptimal$DIC))]
optimalDecisionTable <- cbind(optimalDecisionAIC, optimalDecisionAICc, optimalDecisionSABIC, optimalDecisionDIC)
optimalDecisionTable <- table(optimalDecisionTable)
optimalDecision <- as.numeric(names(optimalDecisionTable)[which(optimalDecisionTable == max(optimalDecisionTable))])
preLCAmod <- mirt::mdirt(datd, optimalDecision, SE = T, SE.type = SE.type, nruns = nruns, empiricalhist = empiricalhist, covdata = covdata, formula = formula)
preLCAmoditemFit <- data.frame(mirt::itemfit(preLCAmod))
message('global optimal: ', optimalDecision, '')
print(preLCAmoditemFit)
if(sum(is.na(preLCAmoditemFit$df.S_X2)) != 0){
datd <- datd[,!is.na(preLCAmoditemFit$df.S_X2)]
} else if (length(which(preLCAmoditemFit$df.S_X2 == 0)) != 0 && length(which(preLCAmoditemFit$df.S_X2 == 0)) != nrow(preLCAmoditemFit)){
datd <- datd[,-which(preLCAmoditemFit$df.S_X2 == 0)]
} else if (length(which(preLCAmoditemFit$p.S_X2 < .05)) != 0){
datd <- datd[,-which(preLCAmoditemFit$S_X2/preLCAmoditemFit$df.S_X2 == max(preLCAmoditemFit$S_X2/preLCAmoditemFit$df.S_X2, na.rm = T))]
} else {
return(list(modelFit = preLCAoptimal, LCAmodel = preLCAmod))
}
rm(preLCAmod)
rm(preLCAoptimal)
rm(preLCAmoditemFit)
}
}
KoreanNounExtraction <- function(dat, polyReturn = F, ExtractType = 'noun'){
if(!require('RHINO')){
install.packages('rJava')
install.packages('devtools')
devtools::install_github('seonghobae/RHINO')
library('RHINO')
}
if(!require('progress')){
install.packages('progress')
library('progress')
}
library('RHINO')
invisible(.connRHINO <<- RHINO::initRhino())
invisible(.connRHINO <- RHINO::initRhino())
for(i in 1:ncol(dat)){
dat[,i] <- as.character(dat[,i])
}
a <- vector()
TotCells <-ncol(dat)*nrow(dat)
k <- 0
pb <- progress::progress_bar$new(
format = " extracting words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
total = TotCells, clear = T, width= 120)
for(i in 1:ncol(dat)){
for(j in 1:length(dat[,i])){
k <- k+1
if(!is.na(dat[j,i])){
try(pb$update(k/TotCells))
try(pb$tick())
dat[j,i] <- gsub("\r\n", " ", dat[j,i])
dat[j,i] <- gsub("\r", " ", dat[j,i])
dat[j,i] <- gsub("\n", " ", dat[j,i])
dat[j,i] <- gsub("^\\s+|\\s+$", "", dat[j,i])
a <- c(a, RHINO::getMorph(dat[j,i], type = ExtractType))
} else {
try(pb$update(k/TotCells))
try(pb$tick())
}
}
}
a <- a[!duplicated(a)]
a <- a[order(a)]
b <- vector()
for(i in 1:length(a)){
if(is.na(a[i])){
b[i] <- NA
}
else if(nchar(a[i]) > 1){
b[i] <- a[i]
} else {
b[i] <- NA
}
}
b <- na.omit(b)
datTextMatrix <- data.frame(matrix(data = 0, nrow = nrow(dat), ncol = length(b)))
colnames(datTextMatrix) <- b
rm(pb)
TotCells <- ncol(datTextMatrix)*ncol(dat)
k <- 0
pb <- progress::progress_bar$new(
format = " arranging words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
total = TotCells, clear = T, width= 120)
for(i in 1:ncol(datTextMatrix)){ # i th word
for(j in 1:ncol(dat)){
k <- k+1
try(pb$update(k/TotCells))
try(pb$tick())
searchEngine <- grep(b[i], dat[,j]) # find a extracted word i in a raw sentense j
if(polyReturn){
datTextMatrix[searchEngine,i] <- datTextMatrix[searchEngine,i]+1
} else {
if(sum(datTextMatrix[searchEngine,i]) == 0){
datTextMatrix[searchEngine,i] <- 1
}
}
}
}
return (datTextMatrix)
}
KoreanExtraction <- function(dat, polyReturn = F){
if(!require('RHINO')){
install.packages('rJava')
install.packages('devtools')
devtools::install_github('seonghobae/RHINO')
library('RHINO')
}
if(!require('progress')){
install.packages('progress')
library('progress')
}
library('RHINO')
invisible(.connRHINO <<- RHINO::initRhino())
invisible(.connRHINO <- RHINO::initRhino())
for(i in 1:ncol(dat)){
dat[,i] <- as.character(dat[,i])
}
a <- vector()
TotCells <-ncol(dat)*nrow(dat)
k <- 0
pb <- progress::progress_bar$new(
format = " extracting words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
total = TotCells, clear = T, width= 120)
for(i in 1:ncol(dat)){
for(j in 1:length(dat[,i])){
k <- k+1
if(!is.na(dat[j,i])){
try(pb$update(k/TotCells))
try(pb$tick())
dat[j,i] <- gsub("\r\n", " ", dat[j,i])
dat[j,i] <- gsub("\r", " ", dat[j,i])
dat[j,i] <- gsub("\n", " ", dat[j,i])
dat[j,i] <- gsub("^\\s+|\\s+$", "", dat[j,i])
a <- c(a, RHINO::getMorph(dat[j,i], type = 'noun'))
a <- c(a, RHINO::getMorph(dat[j,i], type = 'verb'))
} else {
try(pb$update(k/TotCells))
try(pb$tick())
}
}
}
a <- a[!duplicated(a)]
a <- a[order(a)]
b <- vector()
for(i in 1:length(a)){
if(is.na(a[i])){
b[i] <- NA
}
else if(nchar(a[i]) > 1){
b[i] <- a[i]
} else {
b[i] <- NA
}
}
b <- na.omit(b)
datTextMatrix <- data.frame(matrix(data = 0, nrow = nrow(dat), ncol = length(b)))
colnames(datTextMatrix) <- b
rm(pb)
TotCells <- ncol(datTextMatrix)*ncol(dat)
k <- 0
pb <- progress::progress_bar$new(
format = " arranging words [:bar] :percent in :elapsed (:current of :total cells, ETA: :eta)",
total = TotCells, clear = T, width= 120)
for(i in 1:ncol(datTextMatrix)){ # i th word
for(j in 1:ncol(dat)){
k <- k+1
try(pb$update(k/TotCells))
try(pb$tick())
searchEngine <- grep(b[i], dat[,j]) # find a extracted word i in a raw sentense j
if(polyReturn){
datTextMatrix[searchEngine,i] <- datTextMatrix[searchEngine,i]+1
} else {
if(sum(datTextMatrix[searchEngine,i]) == 0){
datTextMatrix[searchEngine,i] <- 1
}
}
}
}
return (datTextMatrix)
}
doLSA <- function(data, polyReturn = F, personID = NULL, SE.type = 'Oakes', checkSecondOrderTest = T, nruns = 1, maxClasses = NULL, empiricalhist = T){
occuranceMatrix <- KoreanExtraction(data, polyReturn = polyReturn)
if(!is.null(personID)){
occuranceMatrix <- aggregate(occuranceMatrix, by = list(as.factor(personID)), FUN = sum)
if(!polyReturn){
occuranceMatrix[-1][occuranceMatrix[-1] > 0] <- 1
}
} else {
occuranceMatrix <- data.frame(matrix(nrow = nrow(occuranceMatrix)), occuranceMatrix)
}
modLSA <- doLCA(data = occuranceMatrix[-1], SE.type = SE.type, checkSecondOrderTest = checkSecondOrderTest, nruns = nruns, maxClasses = maxClasses, empiricalhist = empiricalhist, covdata = occuranceMatrix[1])
return(modLSA)
}
sequentialFA <- function(data, minTestLength = 3, SE.type = 'Oakes', forceUIRT = F, forceCTTmode = F, forceNormalEM = T){
data <- data.frame(data)
dataNames <- colnames(data)
STOP <- FALSE
j <- 0
seqModels <- list()
while(!STOP){
if(length(dataNames) > minTestLength){
estMod <- surveyFA(data[dataNames], SE.type = SE.type, forceUIRT = forceUIRT, forceCTTmode = forceCTTmode, minimumLeftItems = minTestLength, forceNormalEM = forceNormalEM)
if(exists('estMod')){
j <- j+1
seqModels[[j]] <- estMod
dataNames <- dataNames[!dataNames %in% colnames(estMod@Data$data)]
rm(estMod)
} else {
STOP <- TRUE
}
} else {
STOP <- TRUE
}
}
return(seqModels)
}
fitMLIRT <- function(dat = NULL, model = 1, itemtype = 'nominal',
covdata, fixed = ~1, random = NULL, lr.fixed = ~1, lr.random = NULL,
NCYCLES = 4000, BURNIN = 1500, SEMCYCLES = 1000,
GenRandomPars = T, symmetric = F, accelerate = 'squarem', MHRM_SE_draws = 10000,
survey.weights = NULL, SEtol = 1e-4, removeEmptyRows = T){
combine <- function (x, y) {
combn (y, x, paste, collapse = ", ")
}
if(!require('progress')){
install.packages('progress')
library('progress')
}
res1 <- paste0('list(',unlist (lapply (0:NROW(random), combine, random)), ')')
res2 <- paste0('list(',unlist (lapply (0:NROW(lr.random), combine, lr.random)), ')')
res <- vector()
res[1] <- length(res1)
res[2] <- length(res2)
if(length(res1) > length(res2)){
res2 <- rep(res2, length(res1))
} else {
res1 <- rep(res1, length(res2))
}
DIC <- vector()
pb <- progress_bar$new(
format = " estimating optimal multilevel model [:bar] :percent in :elapsed (:current of :total) eta: :eta",
total = max(res)+2, clear = FALSE, width= 120)
for(i in 1:max(res)){
pb$tick()
try(invisible(gc()))
try(suppressMessages(invisible(modMIXED <- mirt::mixedmirt(data = dat, model = model, itemtype = itemtype,
covdata = covdata, fixed = fixed, random = eval(parse(text=res1[i])),
lr.fixed = lr.fixed, lr.random = eval(parse(text=res2[i])),
verbose = F, GenRandomPars = GenRandomPars, accelerate = accelerate,
technical = list(NCYCLES = NCYCLES, BURNIN = BURNIN, SEMCYCLES = SEMCYCLES,
MHRM_SE_draws = MHRM_SE_draws,
symmetric = symmetric,
removeEmptyRows = removeEmptyRows,
SEtol = SEtol), SE = F,
survey.weights = survey.weights))), silent = T)
if(exists('modMIXED')){
if(modMIXED@OptimInfo$converged){
DIC[i] <- modMIXED@Fit$DIC
} else {
DIC[i] <- NA
}
} else {
DIC[i] <- NA
}
try(invisible(rm(modMIXED)), silent = T)
}
if(sum(is.na(DIC)) != length(DIC)){
pb$tick()
try(suppressMessages(invisible(modMIXED <- mirt::mixedmirt(data = dat, model = model, itemtype = itemtype,
covdata = covdata, fixed = fixed, random = eval(parse(text=res1[which(DIC == min(DIC, na.rm = T))])),
lr.fixed = lr.fixed, lr.random = eval(parse(text=res2[which(DIC == min(DIC, na.rm = T))])),
verbose = F, GenRandomPars = GenRandomPars, accelerate = accelerate,
technical = list(NCYCLES = NCYCLES, BURNIN = BURNIN, SEMCYCLES = SEMCYCLES,
MHRM_SE_draws = MHRM_SE_draws,
symmetric = symmetric,
removeEmptyRows = removeEmptyRows,
SEtol = SEtol), SE = T,
survey.weights = survey.weights))), silent = T)
pb$tick()
return(modMIXED)
} else{
stop('No convergence')
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.