Using mirt model fit for PerFit::lz & lzstar

PerFit 패키지의 lz, lzstar는 응답의 비정상 여부를 판단하는 Person-fit statistics를 계산해 준다. 문제는 lz, lzstar는 ltm::ltm 함수를 기본으로 하는데 이 ltm은 parameter에 범위를 지정할 수가 없다. 따라서 mirt로 fit한 결과를 활용하여 lz, lzstar를 구하는 방법을 생각해보았다.

R package PerFit has functions like lz and lzstar which return person fit statistics. The problem is that these functions are based on ltm::ltm but ltm::ltm can set bounds on parameters which is easily incorporated by mirt::mirt. So I figured out the way to use mirt fit for lz and lzstar. 

한 가지 주의할 점은 mirt와 ltm은 서로 다른 방식의 parametrization을 이용한다는 점이다. 흔히 말하는 IRT 파라미터는 난이도와 변별도를 쓴다. mirt::coef(    , IRTpars=T, digits=Inf)을 써야 한다.

One thing to be cautious about is that mirt and ltm have different way of parametrization. Bare in mind. mirt::coef(   , IRTpars=T, digits=Inf)

library(mirt)

library(ltm) # install.packages(c("ltm", "PerFit"))

library(PerFit)

NCYCLES = 10000

coef2mat = function(cf) {

  stopifnot(is.list(cf))

  do.call(rbind, cf[1:(length(cf)-1)])

}

#http://www.r-bloggers.com/five-ways-to-visualize-your-pairwise-comparisons/

panel.cor <- function(x, y, digits=2, prefix="", cex.cor) 

{

  usr <- par("usr"); on.exit(par(usr)) 

  par(usr = c(0, 1, 0, 1)) 

  r <- abs(cor(x, y)) 

  txt <- format(c(r, 0.123456789), digits=digits)[1] 

  txt <- paste(prefix, txt, sep="") 

  if(missing(cex.cor)) cex <- 0.8/strwidth(txt) 

  

  test <- cor.test(x,y) 

  # borrowed from printCoefmat

  Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, 

                   cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),

                   symbols = c("***", "**", "*", ".", " ")) 

  

  text(0.5, 0.5, txt, cex = cex * r) 

  text(.8, .8, Signif, cex=cex, col=2) 

}

Pairs <- function(x, ...) {

  stopifnot(is.matrix(x) | is.data.frame(x))

  pairs(x, lower.panel=panel.smooth, upper.panel=panel.cor)

}

R <- simdata(a=c(1,0.8,1.2,1.5), d=c(-2,-1,0,1,2), 200, itemtype="dich")

# Compare the personfit from mirt and lz

# mirt::personfit(  , method="ML","EAP","MAP","WLE")$Zh

# PerFit::lz(       , method="ML",     ,"BM", "WL")$PFscores$PFscores

# Compare the personfit from mirt and lz with coefs from mirt and lz

fitMirt <- mirt(R, 1, itemtype="2PL", technical = list(NCYCLES=NCYCLES))

#fitLtm <- ltm(R ~ z1, control=list(GHk=61, iter.em=NCYCLES, iter.qN=1000))

fitLtm <- ltm(R ~ z1, control=list(GHk=61, iter.em=NCYCLES))

# Convergence?

fitMirt@OptimInfo$converged

fitLtm$convergence

# logLik

fitMirt@Fit$logLik

fitLtm$log.Lik

# Coefficients

IPMirt <- coef2mat(coef(fitMirt, IRTpars=T, digits=Inf))[,1:2]

IPMirt <- cbind(IPMirt, 0)

IPLtm <- coef(fitLtm)

IPLtm <- IPLtm[,c(2,1)]

IPLtm <- cbind(IPLtm, 0)

IPLtm

Pairs(data.frame(IPMirt[,"a"], IPLtm[,"Dscrmn"], IPMirt[,"b"], IPLtm[,"Dffclt"]))

# theta Estimates

thLtmEAP <- factor.scores(fitLtm, method="EAP", resp.patterns=R)$score.dat$z1

thLtmEB <- factor.scores(fitLtm, method="EB", resp.patterns=R)$score.dat$z1

thMirtML <- fscores(fitMirt, full.scores=T, method="ML")[,1]

thMirtEAP <- fscores(fitMirt, full.scores=T, method="EAP")[,1]

thMirtWLE <- fscores(fitMirt, full.scores=T, method="WLE")[,1]

thMirtMAP <- fscores(fitMirt, full.scores=T, method="MAP")[,1]

thIrtoysML <- irtoys::mlebme(R, IPLtm, method="ML")[,1]

thIrtoysBM <- irtoys::mlebme(R, IPLtm, method="BM")[,1]

lapply(data.frame(thLtmEAP, thMirtEAP, thLtmEB, thMirtMAP, thIrtoysBM, thMirtML,  thIrtoysML, thMirtWLE),

       function(x) c(nan = sum(is.nan(x)), na = sum(is.na(x)), inf = sum(is.infinite(x)),

                     min = min(x), max = max(x)))

plot(data.frame(thLtmEAP, thMirtEAP, thLtmEB, thMirtMAP, thIrtoysBM, thMirtML,  thIrtoysML, thMirtWLE))

Pairs(data.frame(thLtmEAP, thMirtEAP, thLtmEB, thMirtMAP, thIrtoysBM, thMirtML,  thIrtoysML, thMirtWLE))

thMirtML2 <- thMirtML

thMirtML2[is.infinite(thMirtML2)] = 6 * sign(thMirtML2[is.infinite(thMirtML2)])

Pairs(data.frame(thLtmEAP, thMirtEAP, thLtmEB, thMirtMAP, thIrtoysBM, thMirtML2,  thIrtoysML, thMirtWLE))

# -Inf vs -3.99995(for thMirtML/thIrtoysML)

# personfit from Ltm and Mirt

lzMirtML <- personfit(fitMirt, method="ML")[,"Zh"]

lzMirtEAP <- personfit(fitMirt, method="EAP")[,"Zh"]

lzMirtWLE <- personfit(fitMirt, method="WLE")[,"Zh"]

lzMirtMAP <- personfit(fitMirt, method="MAP")[,"Zh"]

lzPFML <- lz(R, Ability.PModel = "ML")$PFscores$PFscores

lzPFBM <- lz(R, Ability.PModel = "BM")$PFscores$PFscores

lzPFWL <- lz(R, Ability.PModel = "WL")$PFscores$PFscores

plot(data.frame(lzMirtML, lzPFML, lzMirtEAP, lzMirtWLE, lzPFWL, lzMirtMAP, lzPFBM))

Pairs(data.frame(lzMirtML, lzPFML, lzMirtEAP, lzMirtWLE, lzPFWL, lzMirtMAP, lzPFBM))

# personfit from Ltm using irtoys

lzPFML.ltm <- lz(R, Ability.PModel = "ML",  IP=IPLtm, Ability = thIrtoysML)$PFscores$PFscores

lzPFBM.ltm <- lz(R, Ability.PModel = "BM",  IP=IPLtm, Ability = thIrtoysBM)$PFscores$PFscores

plot(data.frame(lzPFML, lzPFML.ltm, lzPFBM, lzPFBM.ltm))

Pairs(data.frame(lzPFML, lzPFML.ltm, lzPFBM, lzPFBM.ltm))

# personfit from Mirt using Mirt

lzMirtML.mirt <- personfit(fitMirt, Theta = matrix(thMirtML, ncol=1))

lzMirtMAP.mirt <- personfit(fitMirt, Theta = matrix(thMirtMAP, ncol=1))

lzMirtWLE.mirt <- personfit(fitMirt, Theta = matrix(thMirtWLE, ncol=1))

plot(data.frame(lzMirtML, lzMirtML.mirt, lzMirtMAP, lzMirtMAP.mirt, lzMirtWLE, lzMirtWLE.mirt))

Pairs(data.frame(lzMirtML, lzMirtML.mirt, lzMirtMAP, lzMirtMAP.mirt, lzMirtWLE, lzMirtWLE.mirt))

# personfit from Ltm using Mirt

lzPFML.Mirt <- lz(R, Ability.PModel = "ML",  IP=IPMirt, Ability = thMirtML)$PFscores$PFscores

lzPFML.Mirt2 <- lz(R, Ability.PModel = "ML",  IP=IPMirt, Ability = thMirtML2)$PFscores$PFscores

lzPFBM.Mirt <- lz(R, Ability.PModel = "BM",  IP=IPMirt, Ability = thMirtMAP)$PFscores$PFscores

lzPFWL.Mirt <- lz(R, Ability.PModel = "WL",  IP=IPMirt, Ability = thMirtWLE)$PFscores$PFscores

plot(data.frame(lzPFML, lzPFML.Mirt, lzPFWL, lzPFWL.Mirt, lzPFBM, lzPFBM.Mirt))

Pairs(data.frame(lzPFML, lzPFML.Mirt, lzPFWL, lzPFWL.Mirt, lzPFBM, lzPFBM.Mirt))

Pairs(data.frame(lzPFML, lzPFML.Mirt2, lzPFWL, lzPFWL.Mirt, lzPFBM, lzPFBM.Mirt))

# personfit from Mirt using irtoys

fitMirtTemp <- mirt(R, 1, itemtype="2PL", pars="values")

fitMirtTemp[fitMirtTemp$name == "a1","value"] = IPLtm[,"Dscrmn"]

fitMirtTemp[fitMirtTemp$name == "d","value"] = -IPLtm[,"Dffclt"]*IPLtm[,"Dscrmn"]

fitMirtTemp2 <- mirt(R, 1, itemtype="2PL", pars=fitMirtTemp, TOL=NaN)

lzMirtBM.irtoys <- personfit(fitMirtTemp2, Theta = matrix(thIrtoysBM, ncol=1))

lzMirtML.irtoys <- personfit(fitMirtTemp2, Theta = matrix(thIrtoysML, ncol=1))

plot(data.frame(lzMirtMAP, lzMirtBM.irtoys, lzMirtML, lzMirtML.irtoys), bg="blue")

Pairs(data.frame(lzMirtMAP, lzMirtBM.irtoys, lzMirtML, lzMirtML.irtoys), bg="blue")

#So the correlations are at least 0.99