scripts/Empirical example HADS/hads_anxiety.R

In joakimwallmark/PolyOptimalIRT:

library(TestGardener)
library(mirt)

# ----------- read in data -----------
titlestr  <- "hads-7-anxiety-screen"

U <- scan("data/hads.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)

U <- U[, c(1, 3, 5, 7, 9, 11, 13)]
n <- 7

# Mirt code ---------------------------------------------------------------
colnames(U) <- paste("I", seq(1, ncol(U), 1), sep = "")

hadsAnxietyGrm = mirt(data = U, model=1, itemtype = 'graded', SE=T)
save(U, hadsAnxietyGrm, file="data/hads_anxiety_fittedmodel_mirt.RData")
for(i in 1:ncol(U)){
    ItemPlot <- itemfit(hadsAnxietyGrm,
                        group.bins=15,
                        empirical.plot = i,
                        empirical.CI = .95,
                        method = 'ML')
    print(ItemPlot)
}
itemfit(hadsAnxietyGrm) # item 5 has worst fit


# TestGardener code -------------------------------------------------------
# convert U from score format to index format by adding 1
U <- U + 1

# read item/options (csv)
sheet <- read.csv("data/hads_sheet_anxeity.csv", header = F)

key     <- NULL # NULL for scaled items

noption  <- rep(4,n)    # number of options per item

# --------- Define the option score values for each item ---------
ind       <- c(1:n)
itemVec   <- sheet[, 1] # item strings
optionMat <- sheet[,-1] # option strings

scoreList <- list() # option scores
labelList <- list() # list of just option labels
for (item in ind){
    scoreList[[item]] <- c(0:3)
    labelList[[item]] <- optionMat[item,1:noption[item]]
}

optList <- list(itemLab=itemVec, optLab=labelList, optScr=scoreList)

scrmin  <- 0
scrmax  <- 3*n
scrrng  <- c(scrmin,scrmax)

# 16 bin 9 basis is default. Much better fit with 5
dataList <- make.dataList(U, key, optList, scrrng = scrrng, nbin = 16, NumBasis = 5)
dataList <- make.dataList(U, key, optList, scrrng = scrrng)

# percentage rank density
theta     <- dataList$percntrnk
thetaQnt  <- dataList$thetaQnt
chartList <- dataList$chartList

# sum score density
ndensbasis <- 15
scoreDensity(dataList$scrvec, scrrng, ndensbasis,paste(titlestr,"- sum score"))


# ---------------  Optimal scoring: cycle of smoothing/theta estimation  ------------
#  Set number of cycles and the cell array to containing the parameter
ncycle <- 8
AnalyzeResult <- Analyze(theta, thetaQnt, dataList, ncycle=ncycle, itdisp=FALSE)
ifit <- itemFitTestGardener(AnalyzeResult, data = U-1)
ifit$SX2
ifitpar <- itemFitMirt(hadsAnxietyGrm, data = U-1)
ifitpar$SX2

parList  <- AnalyzeResult$parList
meanHvec <- AnalyzeResult$meanHvec

cycleno <- 1:ncycle
par(mfrow=c(1,1))
plot(cycleno,meanHvec, type="b", lwd=2, xlab="Cycle Number")

#  select cycle for plotting
icycle <- 6

parListi   <- parList[[icycle]]

WfdList    <- parListi$WfdList
Qvec       <- parListi$Qvec
binctr     <- parListi$binctr
theta      <- parListi$theta
arclength  <- parListi$arclength
alfine     <- parListi$alfine

#  ----------------------------------------------------------------------------
#                   Plot surprisal curves for each test question
#  ----------------------------------------------------------------------------
#  plot both the probability and surprisal curves along with data points
Wbinsmth.plot(binctr, Qvec, WfdList, dataList, Wrng=c(0,4), saveplot=F)
Wbinsmth.plot(binctr, titlestr, Qvec, WfdList, key, dataList, Wrng=c(0,4))