Nothing
#' Q-matrix validation
#'
#' Q-matrix validation for the (sequential) G-DINA model based on PVAF (de la Torre & Chiu, 2016; Najera, Sorrel, & Abad, 2019; Najera et al., 2020), stepwise Wald test (Ma & de la Torre, 2020) or mesa plot (de la Torre & Ma, 2016).
#' All these methods are suitable for dichotomous and ordinal response data. If too many modifications are suggested based on the default PVAF method, you are suggested to try the stepwise Wald test method, iterative procedures or predicted cutoffs.
#' You should always check the mesa plots for further examination.
#'
#' @param GDINA.obj an estimated model object of class \code{GDINA}
#' @param method which Q-matrix validation method is used? Can be either \code{"PVAF"} or \code{"wald"}.
#' @param iter implement the method iteratively? Can be \code{"none"} for non-iterative validation (by default), \code{"test"}, \code{"test.att"}, or \code{"item"} (Najera et al., 2020).
#' @param eps cutoff value for PVAF from 0 to 1. Default = 0.95. Note that it can also be -1, indicating the predicted cutoff based on Najera, Sorrel, and Abad (2019).
#' @param digits how many decimal places in each number? The default is 4.
#' @param wald.args a list of arguments for the stepwise Wald test method.
#' \describe{
#' \item{SE.type}{type of covariance matrix for the Wald test}
#' \item{alpha.level}{alpha level for the wald test}
#' \item{GDI}{it can be 0, 1 or 2; 0 means GDI is not used to choose the attribute -
#' when more than one attributes are significant, the one with the largest p-value will be selected;
#' GDI=1 means the attribute with the largest GDI will be selected; GDI=2 means the q-vector with
#' the largest GDI will be selected.}
#' \item{verbose}{print detailed information or not?}
#' \item{stepwise}{\code{TRUE} for stepwise approach and \code{FALSE} for forward approach}
#' }
#' @param iter.args a list of arguments for the iterative implementation.
#' \describe{
#' \item{empty.att}{can a Q-matrix with an empty attribute (i.e., measured by no items) be provided? Default is FALSE}
#' \item{max.iter}{maximum number of iterations. Default is 150}
#' \item{verbose}{print information after each iteration? Default is FALSE}
#' }
#' @return An object of class \code{Qval}. Elements that can be
#' extracted using \code{extract} method include:
#' \describe{
#' \item{sug.Q}{suggested Q-matrix}
#' \item{Q}{original Q-matrix}
#' \item{varsigma}{varsigma index}
#' \item{PVAF}{PVAF}
#' }
#'
#' @include GDINA.R
#' @author {Wenchao Ma, The University of Alabama, \email{wenchao.ma@@ua.edu}, \cr Miguel A. Sorrel, Universidad Autónoma de Madrid, \cr Jimmy de la Torre, The University of Hong Kong}
#' @references
#' de la Torre, J. & Chiu, C-Y. (2016). A General Method of Empirical Q-matrix Validation. \emph{Psychometrika, 81}, 253-273.
#'
#' de la Torre, J., & Ma, W. (2016, August). Cognitive diagnosis modeling: A general framework approach and its implementation in R. A Short Course at the Fourth Conference on Statistical Methods in Psychometrics, Columbia University, New York.
#'
#' Ma, W., & de la Torre, J. (2020). An empirical Q-matrix validation method for the sequential G-DINA model. \emph{British Journal of Mathematical and Statistical Psychology, 73}, 142-163.
#'
#' Najera, P., Sorrel, M. A., & Abad, F.J. (2019). Reconsidering cutoff points in the general method of empirical Q-matrix validation. \emph{Educational and Psychological Measurement, 79}, 727-753.
#'
#' Najera, P., Sorrel, M. A., de la Torre, J., & Abad, F. J. (2020). Improving robustness in Q-matrix validation using an iterative and dynamic procedure. \emph{Applied Psychological Measurement}.
#'
#' @seealso \code{\link{GDINA}}
#' @export
#' @examples
#'\dontrun{
#' ################################
#' #
#' # Binary response
#' #
#' ################################
#' dat <- sim10GDINA$simdat
#' Q <- sim10GDINA$simQ
#' Q[10,] <- c(0,1,0)
#'
#' # Fit the G-DINA model
#' mod1 <- GDINA(dat = dat, Q = Q, model = "GDINA")
#'
#' # Q-validation using de la Torre and Chiu's method
#' pvaf <- Qval(mod1,method = "PVAF",eps = 0.95)
#' pvaf
#' extract(pvaf,what = "PVAF")
#' #See also:
#' extract(pvaf,what = "varsigma")
#' extract(pvaf,what = "sug.Q")
#'
#' # Draw mesa plots using the function plot
#'
#' plot(pvaf,item=10)
#'
#' #The stepwise Wald test
#' stepwise <- Qval(mod1,method = "wald")
#' stepwise
#' extract(stepwise,what = "PVAF")
#' #See also:
#' extract(stepwise,what = "varsigma")
#' extract(stepwise,what = "sug.Q")
#'
#' #Set eps = -1 to determine the cutoff empirically
#' pvaf2 <- Qval(mod1,method = "PVAF",eps = -1)
#' pvaf2
#'
#' #Iterative procedure (test-attribute level)
#' pvaf3 <- Qval(mod1, method = "PVAF", eps = -1,
#' iter = "test.att", iter.args = list(verbose = 1))
#' pvaf3
#'
#' ################################
#' #
#' # Ordinal response
#' #
#' ################################
#' seq.est <- GDINA(sim20seqGDINA$simdat,sim20seqGDINA$simQ, sequential = TRUE)
#' stepwise <- Qval(seq.est, method = "wald")
#'}
Qval <- function(GDINA.obj, method = "PVAF", iter = "none", eps = 0.95,
digits = 4, wald.args = list(),
iter.args = list(empty.att = FALSE, max.iter = 150, verbose = FALSE)){
stopifnot(isa(GDINA.obj,"GDINA"))
if (!is.null(extract(GDINA.obj, "att.str")))
stop("Q-matrix validation is not available if attributes are structured.",
call. = FALSE)
if(method == "wald" & iter != "none"){
warning("Iterative implementation is not available for the Wald method.")
}
if(eps != -1){
if (eps > 1 || eps < 0) stop("eps must be between 0 and 1, or equal to -1.", call. = FALSE)
}
if(!(iter %in% c("none", "test", "test.att", "item"))){stop("iter must be 'none', 'test', 'test.att', or 'item'.")}
if(extract(GDINA.obj,"ngroup")>1) stop("Only available for single-group models.",call. = FALSE)
# if(any(extract(GDINA.obj,"models_numeric")<0)||any(extract(GDINA.obj,"models_numeric")>5))
# stop("Models must be GDINA, DINA, DINO, ACDM, LLM or RRUM",call. = FALSE)
if (max(extract(GDINA.obj, "Q")) > 1)
stop("Q-validation can only be used for dichotomous attribute G-DINA model.",
call. = FALSE)
updated.wald.args <- NULL
if(toupper(method)=="PVAF"){
if(is.null(iter.args$empty.att)){iter.args$empty.att <- FALSE}
if(is.null(iter.args$max.iter)){iter.args$max.iter <- 150}
if(is.null(iter.args$verbose)){iter.args$verbose <- FALSE}
iter.args <- iter.args[c("empty.att", "max.iter", "verbose")]
# if(any(extract(GDINA.obj,"models")!="GDINA"))
# warning("Saturated G-DINA model may be used to calibrate all items for better performance.",call. = FALSE)
ret <- Qval_PVAF(GDINA.obj,eps = eps, digits = digits,
iter = iter, iter.args = iter.args)
}else if (toupper(method)=="WALD"){
if(any(extract(GDINA.obj,"models")!="GDINA"))
stop("Saturated G-DINA model needs to be fitted to all items.",call. = FALSE)
if(eps == -1){
gs <- coef(GDINA.obj, what = "gs") # item parameters
eps <- plogis(-0.4045140782147791 +
4.8404850955032684E-4*extract(GDINA.obj,"nobs") +
2.8667570118638275*(1-sum(colMeans(gs))) +
-0.003315555999671906*extract(GDINA.obj,"nitem"))
}
args.default <- list(GDINA.obj = GDINA.obj, SE.type = 2,
alpha.level = 0.05, GDI = 2, PVAF = eps,
verbose = FALSE, stepwise = TRUE,digits=digits)
updated.wald.args <- modifyList(args.default,wald.args)
ret <- do.call(Qval_wald,updated.wald.args)
}
ret$method <- method
ret$iter <- iter
ret$iter.args <- iter.args
ret$wald.args <- updated.wald.args
ret$sequential <- extract(GDINA.obj,"sequential")
class(ret) <- "Qval"
return(ret)
}
Qval_wald <- function(GDINA.obj, SE.type = 2,
alpha.level = 0.05, GDI = 2, PVAF = 0.95,
verbose = FALSE, stepwise = TRUE, digits = 4,...){
seqent <- extract(GDINA.obj,"sequential")
if(seqent){
dat <- extract(GDINA.obj,"seq.dat")
}else{
dat <- extract(GDINA.obj,"dat")
}
Qc <- extract(GDINA.obj,"Qc")
Qr <- Q <- extract(GDINA.obj,"Q")
N <- extract(GDINA.obj,"nobs")
# number of categories
J <- extract(GDINA.obj,"ncat")
K <- extract(GDINA.obj,"natt")
L <- 2^K
Kj <- rowSums(attributepattern(K)[-1,])
w <- extract(GDINA.obj,"posterior.prob") #1 x L
scofun <- score(GDINA.obj,parm="prob") # a list with # of category elements
# RN <- NgRg(GDINA.obj$technicals$logposterior.i,seqdat,eta.loc(matrix(1,J,K)),1-is.na(seqdat))
expectedR <- extract(GDINA.obj,"expectedCorrect.LC")
expectedN <- extract(GDINA.obj,"expectedTotal.LC")
est.p <- t((expectedR + 1e-10)/(expectedN + 2*1e-10))
# est.p[is.nan(est.p)] <- 0.5
patt <- attributepattern(K)[-1,]
loc <- LC2LG(patt) #2^K-1 x 2^K
vsg <- varsigma(as.matrix(t(loc)),as.matrix(est.p),c(w)) # versigma
vsg0 <- vsg/vsg[,ncol(vsg)] #pvaf
vsgK <- vsg0[,1:K] # pvaf for single att
if(verbose){
cat("\nGDI\n")
print(vsg[,1:K])
cat("\nPVAF\n")
print(vsgK)
cat("First selected attribute based on GDI is",apply(vsgK, 1, which.max.randomtie))
}
#first attribute is the one with the largest GDI
first.att <- apply(vsgK, 1, which.max.randomtie)
att.monitor <- Qrr <- vector("list",J)
inichoose <- numeric(J)
fullset <- c(1:K)
# iteras <- NULL
item <- c(1:nrow(Q))
for (j in item) {
item.no <- Qc[j,1]
inichoose[j] <- currentset <- first.att[j] # initial att. ---largest GDI
att.monitor[[j]] <- c(att.monitor[[j]],currentset)
loop <- ifelse(vsgK[j,currentset]>=PVAF,FALSE,TRUE)
it <- 1
while(loop&&it<K){
difset <- setdiff(fullset,currentset)
add.a <- NULL
#********************************************************************Second round eval. Wald forward
Rm <- Rmatrix.att(length(currentset)+1)
Wp.a <- NULL
for (k in difset){ # 2nd att.
#
Qr <- extract(GDINA.obj,"Q")
Qr[j,seq_len(ncol(Qr))] <- 0
Qr[j,c(currentset,k)] <- 1
etas <- LC2LG(as.matrix(Qr))
itemparj <- extract(GDINA.obj,"catprob.parm")
itemparj[[j]] <- aggregate(expectedR[j,],by=list(etas[j,]),sum)$x/
aggregate(expectedN[j,],by=list(etas[j,]),sum)$x
index <- data.frame(Cat=rep(1:length(rowSums(Qr) ),2^rowSums(Qr)) )
index$Column <- seq_len(length(index$Cat))
sco <- scofun
sco[[j]] <- score_pj(Xj = dat[,j], # a vector of item responses to item j
parloc.j=etas[j,,drop=FALSE], # parameter locations for item j - H by 2^K matrix
catprob.j=itemparj[j], # a list with H elements giving the reduced catprob.parm for each nonzero category
logpost=indlogPost(GDINA.obj))[[1]]
if(SE.type==2){
if(extract(GDINA.obj,"att.dist")!="saturated"){
v <- inverse_crossprod(do.call(cbind,sco))
}else{
v <- inverse_crossprod(do.call(cbind,sco[-length(sco)]))
}
}else if(SE.type==3){
if(extract(GDINA.obj,"att.dist")!="saturated")
warning("structural parameters are not considered in the calculation of covariance matrix.",call. = FALSE)
v <- inverse_crossprod(do.call(cbind,sco))
}
# print(j)
# print(dim(v))
# print(index)
Varj <- v[index$Column[which(index$Cat==j)],index$Column[which(index$Cat==j)]]
# print(Varj)
Cparj <- itemparj[[j]]
# print(Cparj)
# p values for attribute in target attribute + current set
Wp.a <- rbind(Wp.a,sapply(Rm,function(x) {
# wald statistic & p values
W <-t(x %*% Cparj) %*% MASS::ginv(x %*% Varj %*% t(x)) %*% (x %*% Cparj)
p <- pchisq(W,nrow(x),lower.tail = FALSE)
return (p)
}))
}
#********************************************************************end 2nd round eval. att. Forward
remove.new <- add.new <- NULL
for (d in 1:length(difset)){
newq <- vector("numeric",K)
newq[c(currentset,difset[d])] <- 1
newq.loc <- rowMatch(patt,newq)$row.no
# p values, GDI for target att.
add.new <- rbind(add.new,c(difset[d],Wp.a[d,1+sum(difset[d]>currentset)],vsgK[j,difset[d]],vsg0[j,newq.loc]))
remove.new <- rbind(remove.new,Wp.a[d,-(1+sum(difset[d]>currentset)),drop=FALSE]) #p values for att in current set
}
# iteras <- rbind(iteras,cbind(j,it,add.new))
if(verbose){
cat("\nItem",j,"Att",k,"\n")
info <- cbind(add.new,remove.new)
colnames(info) <- c("att","p(att)","GDI-largest att","GDI-largest vec",paste0("p-A",currentset))
print(info)
}
#*********************************************************If additional elements should be added
if (any(add.new[,2]<alpha.level)){ # Yes - sig some attributes should be required
add.new <- add.new[which(add.new[,2]<alpha.level),,drop=FALSE] # which should be added
remove.new <- remove.new[which(add.new[,2]<alpha.level),,drop=FALSE]
loc.a <- which.max(add.new[,2+GDI])
add.a <- add.new[loc.a,1]
remove.a <- remove.new[loc.a,] # pvalues for all current attributes in the q-vector
if(any(remove.a>alpha.level)) {
remove.a <- currentset[which(remove.a>alpha.level)] # the att that needs to be removed
}else{
remove.a <- NULL
}
if(!stepwise) remove.a <- NULL
currentset <- sort(setdiff(c(currentset,add.a),remove.a)) # all att. chosen
att.monitor[[j]] <- c(att.monitor[[j]],add.a)
if(is.null(remove.a)){ #if no attributes will be removed
current_vec_PVAF <- add.new[which(add.new[,1,drop=FALSE]==add.a),4]
}else{ # if some need to be removed - vector PVAF need to be recalculated
newq <- vector("numeric",K)
newq[currentset] <- 1
current_vec_PVAF <- vsg0[j,rowMatch(patt,newq)$row.no]
}
if(current_vec_PVAF>=PVAF){
loop <- FALSE
}else{
loop <- TRUE
it <- it + 1
}
}else{ # No
loop <- FALSE
}
}
Q[j,] <- 0
Q[j,currentset] <- 1
j <- j+1
}
if(seqent){
ret.Q <- Qc
ret.sugQ <- cbind(Qc[,1:2],Q)
}else{
ret.Q <- extract(GDINA.obj,"Q")
ret.sugQ <- Q
}
out.vsg <- round(t(vsg), digits)
out.PVAF <- round(t(vsg0), digits)
rownames(out.vsg) <-
rownames(out.PVAF) <-
apply(patt, 1, paste, collapse = "")
qvalid <-
list(
sug.Q = ret.sugQ,
Q = ret.Q,
varsigma = out.vsg,
PVAF = out.PVAF,
eps = PVAF,
initialAtt=inichoose,
est.p = est.p
)
return(qvalid)
}
Qval_PVAF <- function(GDINA.obj,
eps = 0.95,
digits = 4,
iter = "none",
iter.args = list(empty.att = FALSE, max.iter = 150, verbose = FALSE)){
seqent <- extract(GDINA.obj,"sequential")
eps.base <- eps
fit <- GDINA.obj
Qiter <- list()
if(seqent){
Y <- extract(GDINA.obj,"seq.dat")
dat <- extract(GDINA.obj,"dat")
Qiter[["0"]] <- Qtest <- Qprov <- Q <- Qr <- Qc <- extract(GDINA.obj,"Qc")
}else{
Y <- dat <- extract(GDINA.obj,"dat")
Qiter[["0"]] <- Qtest <- Qprov <- Q <- extract(GDINA.obj,"Q")
}
model <- GDINA.obj$options$model
mono.constraint <- GDINA.obj$options$mono.constraint
if(!is.null(GDINA.obj$options$higher.order)){
ho.model <- GDINA.obj$options$higher.order$model
ho.nquad <- GDINA.obj$options$higher.order$nquad
ho.sRange <- GDINA.obj$options$higher.order$SlopeRange
ho.iRange <- GDINA.obj$options$higher.order$InterceptRange
ho.sPrior <- GDINA.obj$options$higher.order$SlopePrior
ho.iPrior <- GDINA.obj$options$higher.order$InterceptPrior
ho.prior <- GDINA.obj$options$higher.order$Prior
}
N <- extract(GDINA.obj, "nobs")
J <- extract(GDINA.obj, "nitem")
K <- extract(GDINA.obj, "natt")
if(seqent){
L <- no_LC(Q[,-c(1:2)])
Js <- nrow(Q)
} else {
L <- no_LC(Q)
}
M <- attributepattern(K)[-1,]
Kj <- rowSums(attributepattern(K)[-1, ])
w <- extract(GDINA.obj, "posterior.prob") #1 x L
YY <- Y
YY[is.na(YY)] <- 0
rc <- apply(YY, 2, function(x) {
colSums(x * exp(extract(GDINA.obj, "logposterior.i")))
})
rn <- apply(1 * (!is.na(Y)), 2, function(x) {
colSums(x * exp(extract(GDINA.obj, "logposterior.i")))
})
# est.p <- rc/c(w*N)
est.p <- (rc + 1e-10) / (rn + 2* 1e-10)
patt <- attributepattern(K)[-1, ]
loc <- eta(patt) #2^K-1 x 2^K
vsg <- varsigma(as.matrix(t(loc)), as.matrix(est.p), c(w))
PVAF <- vsg / vsg[, L - 1]
i <- 0
if(iter == "none"){
if(eps.base == -1){
gs <- coef(fit, what = "gs") # item parameters
eps <- plogis(-0.4045140782147791 +
4.8404850955032684E-4*extract(fit,"nobs") +
2.8667570118638275*(1-sum(colMeans(gs))) +
-0.003315555999671906*extract(fit,"nitem"))
}
val_q <- NULL
for (k in sort(unique(Kj))) {
tmp <- PVAF[, which(Kj == k)]
if (length(which(Kj == k)) == 1) {
tmp[which(tmp > eps)] <- 1
tmp[which(tmp <= eps)] <- 0
val_q <- cbind(val_q, tmp)
} else{
val_q <- cbind(val_q, apply(tmp, 1, function(x) {
ifelse (max(x) > eps, which.max(x), 0)
}))
}
}
if (ncol(val_q) > 1) {
for (k in 2:ncol(val_q)) {
val_q[which(val_q[, k] > 0), k] <-
val_q[which(val_q[, k] > 0), k] + sum(Kj < k)
}
}
#### modified Q-matrix and associated PVAF
loc_q <- apply(val_q, 1, function(x) {
x[which.max(x > 0)]
})
val_q <- attributepattern(K)[-1, ][loc_q, ]
out.vsg <- round(t(vsg), digits)
out.PVAF <- round(t(PVAF), digits)
rownames(out.vsg) <-
rownames(out.PVAF) <-
apply(patt, 1, paste, collapse = "")
# Q <- data.frame(Q, row.names = extract(GDINA.obj, "item.names"))
if(seqent){
ret.Q <- Qc
ret.sugQ <- cbind(Qc[,1:2],val_q)
}else{
ret.Q <- Q
ret.sugQ <- val_q
}
i <- n.iter <- convergence <- 1
qvalid <-
list(
sug.Q = ret.sugQ,
Q = ret.Q,
varsigma = out.vsg,
PVAF = out.PVAF,
eps = eps,
est.p = est.p,
n.iter = n.iter,
convergence = convergence
)
} else {
while(i < iter.args$max.iter){
i <- i + 1
if(eps.base == -1){
gs <- coef(fit, what = "gs") # item parameters
eps <- plogis(-0.4045140782147791 +
4.8404850955032684E-4*extract(fit,"nobs") +
2.8667570118638275*(1-sum(colMeans(gs))) +
-0.003315555999671906*extract(fit,"nitem"))
}
w <- extract(fit, "posterior.prob") #1 x L
rc <- apply(YY, 2, function(x) {
colSums(x * exp(extract(fit, "logposterior.i")))
})
rn <- apply(1 * (!is.na(Y)), 2, function(x) {
colSums(x * exp(extract(fit, "logposterior.i")))
})
# est.p <- rc/c(w*N)
est.p <- (rc + 1e-10) / (rn + 2* 1e-10)
loc <- eta(patt) #2^K-1 x 2^K
vsg <- varsigma(as.matrix(t(loc)), as.matrix(est.p), c(w))
PVAF <- vsg / vsg[, L - 1]
val_q <- NULL
for (k in sort(unique(Kj))) {
tmp <- PVAF[, which(Kj == k)]
if (length(which(Kj == k)) == 1) {
tmp[which(tmp > eps)] <- 1
tmp[which(tmp <= eps)] <- 0
val_q <- cbind(val_q, tmp)
} else{
val_q <- cbind(val_q, apply(tmp, 1, function(x) {
ifelse (max(x) > eps, which.max(x), 0)
}))
}
}
if (ncol(val_q) > 1) {
for (k in 2:ncol(val_q)) {
val_q[which(val_q[, k] > 0), k] <-
val_q[which(val_q[, k] > 0), k] + sum(Kj < k)
}
}
#### modified Q-matrix and associated PVAF
loc_q <- apply(val_q, 1, function(x) {
x[which.max(x > 0)]
})
val_q <- attributepattern(K)[-1, ][loc_q, ]
out.vsg <- round(t(vsg), digits)
out.PVAF <- round(t(PVAF), digits)
rownames(out.vsg) <-
rownames(out.PVAF) <-
apply(patt, 1, paste, collapse = "")
# Q <- data.frame(Q, row.names = extract(GDINA.obj, "item.names"))
if(seqent){
ret.Q <- Qc
ret.sugQ <- cbind(Qc[,1:2],val_q)
}else{
ret.Q <- Q
ret.sugQ <- val_q
}
Qsug.obj <-
list(
sug.Q = ret.sugQ,
Q = ret.Q,
varsigma = out.vsg,
PVAF = out.PVAF,
eps = eps,
est.p = est.p)
Qsug <- Qsug.obj$sug.Q
res.eps <- Qsug.obj$eps
if(seqent){
maxloc <- matrix(NA, nrow = K, ncol = Js, dimnames = list(paste0("K", 1:K), paste0("Js", 1:Js)))
} else {
maxloc <- matrix(NA, nrow = K, ncol = J, dimnames = list(paste0("K", 1:K), paste0("J", 1:J)))
}
for(kj in 1:K){maxloc[kj,] <- apply(Qsug.obj$PVAF[Kj == kj,, drop = FALSE], 2, which.max)}
best.pos <- cumsum(c(0, table(Kj)[-K])) + maxloc
if(seqent){
best.PVAF <- sapply(1:Js, function(j) Qsug.obj$PVAF[best.pos[,j],j])
rownames(best.PVAF) <- paste0("K", 1:K); colnames(best.PVAF) <- paste0("Js", 1:Js)
} else {
best.PVAF <- sapply(1:J, function(j) Qsug.obj$PVAF[best.pos[,j],j])
rownames(best.PVAF) <- paste0("K", 1:K); colnames(best.PVAF) <- paste0("J", 1:J)
}
res.PVAF <- Qsug.obj$PVAF
res.varsigma <- Qsug.obj$varsigma
if(seqent){
hit <- hit.cand <- which(sapply(1:Js, function(j) any(Qprov[j,] != Qsug[j,])))
} else {
hit <- hit.cand <- which(sapply(1:J, function(j) any(Qprov[j,] != Qsug[j,])))
}
if(length(hit.cand) == 0){
i <- i - 1
convergence <- 1 # YES convergence
break
} else {
if(seqent){
Kj.prov <- rowSums(Qprov[hit, -c(1:2), drop = F])
Kj.sug <- rowSums(Qsug[hit, -c(1:2), drop = F])
} else {
Kj.prov <- rowSums(Qprov[hit,, drop = F])
Kj.sug <- rowSums(Qsug[hit,, drop = F])
}
if(iter == "test.att"){
att.change <- ifelse(Kj.prov < Kj.sug, 1, ifelse(Kj.prov == Kj.sug, 0, -1))
new.att <- Kj.prov + att.change
for(r in 1:length(hit)){
att.r <- new.att[r]
if(seqent){
Qsug[hit[r], -c(1:2)] <- M[best.pos[att.r, hit[r]],]
} else {
Qsug[hit[r],] <- M[best.pos[att.r, hit[r]],]
}
}
}
if(iter == "item"){
if(seqent){
prov.q <- match(apply(Qprov[hit.cand, -c(1:2), drop = F], 1, paste, collapse = ""), apply(M, 1, paste, collapse = ""))
} else {
prov.q <- match(apply(Qprov[hit.cand,, drop = F], 1, paste, collapse = ""), apply(M, 1, paste, collapse = ""))
}
PVAF.diff <- c()
for(j in 1:length(hit.cand)){
PVAF.diff <- c(PVAF.diff, max(abs(best.PVAF[Kj.sug[j], hit.cand[j]] - best.PVAF[Kj.prov[j], hit.cand[j]]),
abs(best.PVAF[Kj.sug[j], hit.cand[j]] - res.PVAF[prov.q[j], hit.cand[j]])))
}
hit <- hit.cand[which.max(PVAF.diff)]
}
}
Qtest[hit,] <- Qsug[hit,]
if(any(sapply(1:i, function(x) mean(Qtest == Qiter[[x]])) == 1)){convergence <- 2; break} # NO convergence (loop)
if(any(colSums(Qtest) == 0)){convergence <- 3; break} # NO convergence (null attribute)
Qprov <- Qiter[[as.character(i)]] <- Qtest
if(iter.args$verbose){cat("Iteration =", sprintf("%03d", i), "| Item/s modified =", paste(hit, collapse = ", "), "\n")}
if(is.null(GDINA.obj$options$higher.order)){
fit <- GDINA(dat, Qprov, model, mono.constraint = mono.constraint, sequential = seqent, verbose = 0)
} else {
fit <- GDINA(dat, Qprov, model, mono.constraint = mono.constraint, att.dist = "higher.order", sequential = seqent, verbose = 0,
higher.order = list(model = ho.model, nquad = ho.nquad, SlopeRange = ho.sRange, InterceptRange = ho.iRange,
SlopePrior = ho.sPrior, InterceptPrior = ho.iPrior, Prior = ho.prior))
}
}
n.iter <- i
if(n.iter == iter.args$max.iter){convergence <- 4} # NO convergence (max.iter reached)
if(iter.args$empty.att & convergence == 3){
Qsug <- Qtest
} else {
Qsug <- Qprov
}
qvalid <-
list(
# sug.Q = ret.sugQ,
sug.Q = Qsug,
Q = ret.Q,
varsigma = out.vsg,
PVAF = out.PVAF,
eps = eps,
est.p = est.p,
n.iter = n.iter,
convergence = convergence
)
}
return(qvalid)
}
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.