R/simulationIteration.R

In joakimwallmark/PolyOptimalIRT:

simulationIteration <- function(r, n, dataGenerationModel, thetas,
                               dataGenerationModelType="optimal",
                               modelsToCompare=list(c("gpcm", "EAP"), c("gpcm", "ML"), c("optimal", cycles=15)),
                               optScr=NULL,
                               usePrecalibrated=F,
                               preCalibratedOptimal=NULL,
                               preCalibratedParMod=NULL) {
    rows <- 2*length(modelsToCompare)+1
    # extra row for arc length?
    if ("optimal" %in% unlist(modelsToCompare)) rows <- rows+1

    resMat <- matrix(NA, nrow = rows, ncol = n+1)

    # generate data -------------------------------------------------------
    if (dataGenerationModelType=="optimal") {
        # get sample from model
        sample <- Usimulate(length(dataGenerationModel), thetas, dataGenerationModel)
        sample <- sample-1
    }
    else { # get sample from mirt model
        sample <- simdata(model = dataGenerationModel,
                          Theta=matrix(thetas, ncol = 1))
    }
    # get estimates -------------------------------------------------------
    resMat[1, ] <- c("sum", rowSums(sample)) # Get sum scores

    modInd <- 2
    for (model in modelsToCompare) {
        if (model[1]=="optimal") {
            optList <- list(itemLab=NULL, optLab=NULL, optScr=optScr)
            # sample needs to start with 1 instead of 0
            dataList <- make.dataList(sample+1, NULL, optList)
            # Set initial theta values
            initThetas  <- dataList$percntrnk
            if (!usePrecalibrated) { # fit optimal scores model to data
                thetaQnt  <- dataList$thetaQnt
                #  Proceed through the cycles
                cycles <- as.numeric(model[2])
                AnalyzeResult <- Analyze(initThetas, thetaQnt,
                                         dataList, ncycle=cycles,
                                         itdisp=FALSE)
                temp <- AnalyzeResult$parList[[cycles]]$theta
                arclen <- AnalyzeResult$parList[[cycles]]$alfine
                resMat[modInd, ] <- c("arcoptimal", arclen,
                                        rep(NA, ncol(resMat)-1-length(arclen)))
                modInd <- modInd+1
                wfdList <- AnalyzeResult$parList[[cycles]]$WfdList
            }
            else { # estimate theta from prev model
                wfdList <- preCalibratedOptimal
                temp <- thetafun(initThetas,
                                 WfdList = wfdList,
                                 dataList = dataList)$theta_out
            }
            resMat[modInd, ] <- c("optimal", temp)
            escores <- testscore(temp, wfdList, optList) # get expected sum score
            resMat[modInd+1, ] <- c("sumoptimal", escores)
            modInd <- modInd+2
        }
        else { #mirt with model type from modelsToCompare if not optimal
            if (!usePrecalibrated) {
                # fit mirt model to simulated data
                tryCatch({
                    preCalibratedParMod <- mirt(data = data.frame(sample),
                                                model=1, itemtype = model[1], SE=T)
                },
                warning = function(w) {
                    log(w)
                },
                error = function(e) {
                    save(sample, file = "errordata.RData");
                    NaN
                })
            }
            # estimate theta
            temp <- fscores(preCalibratedParMod,
                            method = model[2],
                            response.pattern = sample)[, "F1"]
            resMat[modInd, ] <- c(paste(model[1], model[2], sep = ""),
                                  temp)
            resMat[modInd+1, ] <- c(paste("sumparametric", model[2], sep = ""),
                                    expected.test(preCalibratedParMod, matrix(temp, ncol = 1)))
            modInd <- modInd+2
        }
    }
    resMat
}