Nothing
###############################################################################
# JMLE: #
# #
# Joint maximum likelihood estimation of item parameters and attribute #
# profiles in cognitive diagnostic models. #
# #
# Inputs: #
# (1) Y: a matrix of binary responses (1=correct, 0=incorrect). Rows #
# represent persons, and columns represent items. #
# (2) Q: the Q-matrix of the test. Rows represent items, and colums represent #
# attributes. #
# (3) model: currently has three options, "DINA", "DINO", "NIDA", "GNIDA", and#
# "RRUM". #
# (4) NP.method: "Hamming", the plain hamming distance method; #
# "Weighted", the hamming distance weighted by inversed variance #
# "Penalized", the hamming distance weighted by inversed variance #
# and specified penalizing weights for guess and slip#
# Additional input for the "penalized" NP.method: #
# (5) wg = weight assigned to guess #
# (6) ws = weight assigned to slip #
# (7) conv.crit.par: the critical value for the change in item parameter #
# values to determine convergence #
# (8) conv.crit.att: the critical value for the percentage of attributes that #
# are changed to determine convergence #
# #
# Outputs: #
# (1) item parameters estimates and standard errors for each item #
# (2) attribute profiles for each examinee #
# #
###############################################################################
JMLE <- function(Y, Q, model=c("DINA", "DINO", "NIDA", "GNIDA", "RRUM"),
NP.method=c("Weighted", "Hamming", "Penalized"),
wg=1, ws=1, conv.crit.par=0.001, conv.crit.att=0.01, max.ite=100)
{
#####
# 1 #
##### Check dimension consistency and convert data to the right formats
Y <- as.matrix(Y)
Q <- as.matrix(Q)
nitem <- nrow(Q)
natt <- ncol(Q)
nperson <- nrow(Y)
check <- NULL
check <- CheckInput(Y, Q)
if (!is.null(check)) return(warning(check))
model <- match.arg(model)
NP.method <- match.arg(NP.method)
#####
# 2 #
##### Step 1: Initial model-free estimation of the ability pattern
if (model %in% c("DINA", "NIDA", "GNIDA", "RRUM"))
{
gate="AND"
} else if (model == "DINO")
{
gate="OR"
} else
{
return(warning("Model specification is not valid."))
}
alpha.est.NP <- AlphaNP(Y, Q, gate=gate, method=NP.method, wg=1, ws=1)
alpha.est <- alpha.est.NP$alpha.est
loss.matrix <- alpha.est.NP$loss.matrix
#####
# 3 #
##### Step 2: Iterative MLE estimation for item parameters and ability pattern
d.par <- d.att <- d.undefined <- 1
ite <- 0
conv <- "Convergence criteria met."
while (((max(d.par) > conv.crit.par & d.att > 0) || d.undefined > 0 || (d.att > conv.crit.att & max(d.par) > 0)) & ite < max.ite)
{
ite <- ite + 1
# MLE of item parameters
if (model == "DINA" || model == "DINO")
{
temp <- ParMLE(Y, Q, alpha.est, model)
par.est <- list(slip=temp$slip, guess=temp$guess, se.slip=temp$se.slip, se.guess=temp$se.guess)
undefined.flag <- 1 - ((1 - is.infinite(par.est$slip)) * (1 - is.infinite(par.est$guess))
* (1 - is.nan(par.est$slip)) * (1 - is.nan(par.est$guess)))
} else if (model == "NIDA" || model == "GNIDA")
{
temp <- ParMLE(Y, Q, alpha.est, model)
par.est <- list(slip=temp$slip, guess=temp$guess, se.slip=temp$se.slip, se.guess=temp$se.guess)
undefined.flag <- rep(0, nitem)
} else if (model == "RRUM")
{
temp <- ParMLE(Y, Q, alpha.est, model)
par.est <- list(pi=temp$pi, r=temp$r, se.pi=temp$se.pi, se.r=temp$se.r)
undefined.flag <- rep(0, nitem)
} else
{
return(warning("Model specification is not valid."))
}
# MLE of alpha
alpha.out.MLE <- AlphaMLE(Y, Q, par.est, model, undefined.flag)
alpha.est <- alpha.out.MLE$alpha.est
loglike.matrix <- alpha.out.MLE$loglike.matrix
# Compute loglikelihood
loglike <- 0
for (i in 1:ncol(loglike.matrix)){
loglike <- loglike + loglike.matrix[alpha.out.MLE$est.class[i],i]
}
# Compute changes
if (ite > 1)
{
d.par <- abs(unlist(par.est) - unlist(par.est.old))
d.par <- d.par[(1 - is.infinite(d.par)) * (1 - is.na(d.par))]
d.undefined <- sum(abs(undefined.flag - undefined.flag.old))
if (all(alpha.out.MLE.old$class.tie == alpha.out.MLE$class.tie))
{
d.att <- 0
} else {
d.att <- sum(apply((alpha.out.MLE.old$class.tie != alpha.out.MLE$class.tie), 1, any) &
apply((alpha.est != alpha.est.old), 1, any)) / nrow(Y)
}
}
par.est.old <- par.est
alpha.out.MLE.old <- alpha.out.MLE
alpha.est.old <- alpha.est
undefined.flag.old <- undefined.flag
cat(sprintf("Iteration: %d, loglike = %.5f, diff.par = %.5f, diff.att = %.5f, diff.undefined = %.5f", ite, loglike, max(d.par), d.att, d.undefined))
cat("\n")
}
if (ite == max.ite) conv <- "Maximum iteration reached."
if (max(d.par) > conv.crit.par & d.att == 0) conv <- "Examinee profile has stabilized except for randomness from ties."
# AIC and BIC model fit statistics
if (model %in% c("DINA", "DINO")) npar <- 2 * nitem + 2 ^ natt - 1
if (model == "NIDA") npar <- 2 * natt + 2 ^ natt - 1
if (model == "GNIDA") npar <- 2 * natt * nitem + 2 ^ natt - 1
if (model == "RRUM") npar <- nitem * (natt + 1) + 2 ^ natt - 1
AIC <- -2 * loglike + 2 * npar
BIC <- -2 * loglike + npar * log(nperson)
output <- list(alpha.est=alpha.est, par.est=par.est, n.tie=alpha.out.MLE$n.tie, undefined.flag=undefined.flag,
loglike=loglike, convergence=conv, n.ite=ite, AIC=AIC, BIC=BIC, loglike.matrix=loglike.matrix,
est.class=alpha.out.MLE$est.class, NP.loss.matrix=loss.matrix,
NP.alpha.est=alpha.est.NP$alpha.est, NP.method=NP.method,
NP.est.class=alpha.est.NP$est.class, pattern=alpha.est.NP$pattern, model=model, Q=Q, Y=Y)
class(output) <- "JMLE"
return(output)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.