R/compareModelProbsSim2.R

In joakimwallmark/PolyOptimalIRT:

library(mirt)
library(TestGardener)
library(foreach)
library(doParallel)

# U is data (1-max for each item)
# folds are number of folds for CV
# R are number of iterations
# mirtModel is mirt model type
# scoreOpt is score option list for TestGardener
# cyc are cycles for TestGardener
# noBasis are number of basis functions for TestGardener
# noBins are number of bins functions for TestGardener
compareModelProbsSim2 <- function(U, folds, R, mirtModel="gpcm", thetaMethod="EAP",
                                 scoreOpt, cyc=10, noBasis, noBin) {
    # randomize folds
    foldIds <- sample(rep(1:folds, length.out=nrow(U)))

    itemCorrectMirt <- matrix(0, nrow = folds, ncol = ncol(U))
    itemCorrectTG <- matrix(0, nrow = folds, ncol = ncol(U))
    for (k in 1:folds) {
        u <- U[foldIds!=k, ] # fit with all except for given fold
        testu <- U[foldIds==k, ]

        # Fit models on u
        mirtMod <- mirt(data = data.frame(u-1), model=1, itemtype = mirtModel)
        TGMod <- fitTG(u, cyc, scoreOpt, noBasis, noBin)

        # predict thetas of testu (foldIds==k)
        mirtThetas <- fscores(mirtMod, method = thetaMethod,
                        response.pattern = testu-1)[, "F1"]
        TGThetas <- estimateThetasTG(testu, TGMod$parList[[cyc]]$WfdList, scoreOpt)

        # generate response for theta estimates using both models
        TGWfdList <- TGMod$parList[[cyc]]$WfdList
        cl <- makeCluster(detectCores(logical = F))
        registerDoParallel(cl)
        # add correct predictions from each iteration to matrix using .combine='+'
        res <- foreach(r = 1:R,
                       .packages = c('TestGardener', 'mirt'),
                       .combine = "+", .verbose = T) %dopar% {
                           sampleM <- simdata(model = mirtMod,
                                              Theta=matrix(mirtThetas, ncol = 1))
                           sampleT <- Usimulate(length(TGWfdList),
                                                TGThetas, TGWfdList)
                           # add upp number of correct for each person from each sample
                           # (TG predictions underneath mirt)
                           rbind(sampleM == testu, sampleT == testu)
                       }
        stopCluster(cl)

        # store mean percentage for each item from each person (in given fold)
        mirtMeans <- colMeans(res[1:nrow(testu), ]/R)
        TGMeans <- colMeans(res[-(1:nrow(testu)), ]/R)
        # weight by number of prediction out of total (smaller folds are weighted less)
        itemCorrectMirt <- itemCorrectMirt + mirtMeans*nrow(testu)/nrow(U)
        itemCorrectTG <- itemCorrectTG + TGMeans*nrow(testu)/nrow(U)

        if (!k%%10) { print(k/folds) } # print progress every 10th iteration
    }
    return(list(mirt = itemCorrectMirt, TG = itemCorrectTG))
}

fitTG <- function(U, cyc, scoreOpt, noBasis, noBin) {
    optList <- list(itemLab=NULL, optLab=NULL, optScr=scoreOpt)
    dataList <- make.dataList(U, NULL, optList, NumBasis = noBasis, nbin = noBin)
    initThetas  <- dataList$percntrnk
    thetaQnt  <- dataList$thetaQnt

    analyzeResult <- Analyze(initThetas, thetaQnt,
                             dataList, ncycle=cyc,
                             itdisp=FALSE)
    analyzeResult
}

estimateThetasTG <- function(U, wfdList, optscr) {
    optList <- list(itemLab=NULL, optLab=NULL, optScr=optscr)
    dataList <- make.dataList(U, NULL, optList)
    # Set initial theta values
    initThetas  <- dataList$percntrnk

    temp <- thetafun(initThetas,
                     WfdList = wfdList,
                     dataList = dataList)$theta_out
}