v 1.2.0

DESCRIPTION

@@ -1,15 +1,15 @@
 Package: GDINA
 Type: Package
 Title: The Generalized DINA Model Framework
-Version: 1.1.0
-Date: 2017-1-2
+Version: 1.2.0
+Date: 2017-1-28
 Authors@R: c(person(given = "Wenchao",family = "Ma", role = c("aut", "cre"),email = "wenchao.ma@rutgers.edu"),person(given = "Jimmy", family = "de la Torre", role = "aut"))
 Description: A set of psychometric tools for cognitive diagnostic analyses for
     both dichotomous and polytomous responses. Various cognitive diagnosis models
-    (CDMs) can be estimated, include the generalized deterministic inputs, noisy
+    can be estimated, include the generalized deterministic inputs, noisy
     and gate (G-DINA) model by de la Torre (2011) <DOI:10.1007/s11336-011-9207-7>,
-    the sequential G-DINA model by Ma and de la Torre (2016)<DOI:10.1111/bmsp.
-    12070>, and many other models they subsume. Joint attribute distribution can
+    the sequential G-DINA model by Ma and de la Torre (2016) <DOI:10.1111/bmsp.12070>, 
+    and many other models they subsume. Joint attribute distribution can
     be saturated, higher-order or structured. Q-matrix validation, item and model
     fit statistics, model comparison at test and item level and differential item
     functioning can also be conducted. A graphical user interface is also provided.
@@ -18,21 +18,22 @@ LazyData: TRUE
 Depends:
     R (>= 3.1.0)
 Imports:
-    stats,
+    alabama,
+    data.table,
     graphics,
-    utils,
-    Rcpp (>= 0.12.1),
-    nloptr,
+    ggplot2,
     MASS,
+    nloptr,
     numDeriv,
-    alabama,
+    Rcpp (>= 0.12.1),
     Rsolnp,
-    ggplot2,
-    shiny,
-    shinydashboard,
-    data.table
+    stats,
+    utils
 Suggests:
-    testthat
+    CDM,
+    testthat,
+    shiny,
+    shinydashboard
 LinkingTo: Rcpp, RcppArmadillo
 URL: https://github.com/Wenchao-Ma/GDINA
 BugReports: https://github.com/Wenchao-Ma/GDINA/issues

NAMESPACE

@@ -56,13 +56,15 @@ export(monocheck)
 export(npar)
 export(personparm)
 export(plotIRF)
+export(rowCount)
+export(rowMatch)
 export(simGDINA)
 export(startGDINA)
+export(unique_only)
 import(MASS)
+import(data.table)
 import(ggplot2)
 import(numDeriv)
-import(shiny)
-import(shinydashboard)
 import(stats)
 importFrom(Rcpp,sourceCpp)
 importFrom(alabama,auglag)

NEWS.md

@@ -0,0 +1,14 @@
+# GDINA 1.2.0
+* Fixed   - infinite value issue during M-step optimization
+* Fixed   - `itemfit` function for missing data
+* Fixed   - adding monotonic constraints for the sequential models
+* Fixed   - the maximum likelihood estimation of person attribute through `personparm`
+* Changed - package imports, and suggests
+* Changed - slsqp as the first optimizer used for monotonic G-DINA 
+
+# GDINA 1.0.0
+
+* Initial release
+
+
+

R/GDINA.R

@@ -160,8 +160,8 @@
 #'    necessary. For polytomous attributes, non-zero elements indicate which level
 #'    of attributes are needed. Note that for polytomous items, the sequential G-DINA
 #'    model is used and either restricted or unrestricted category-level Q-matrix is needed.
-#'    The first column represents the item number and
-#'    the second column indicates the category number. See \code{Examples}.
+#'    The first column represents the item number, which must be numeric and match the column in the data.
+#'    The second column indicates the category number. See \code{Examples}.
 #' @param model A vector for each item/category or a scalar which will be used for all
 #'    items/categories to specify which model is fitted to each item/category. The possible options
 #'    include \code{"GDINA"},\code{"DINA"},\code{"DINO"},\code{"ACDM"},\code{"LLM"}, and \code{"RRUM"}.
@@ -220,23 +220,28 @@
 #' @param maxitr The maximum iterations of EM cycles allowed.
 #' @param higher.order.struc.parm A matrix or data frame providing higher order structural parameters. If supplied, it must be of dimension \eqn{K\times 2}.
 #'    The first column is the slope parameters and the second column is the intercept.
-#' @param diagnosis Run in advanced mode? If it is 1 or 2, some intermediate results obtained in each iteration will be saved in \code{diagnos}.
-#' @param Mstep.warning Logical; indicating whether the warning message in Mstep, if any, should be output immediately.
-#' @param optimizer Advanced option: a string indicating which optimizer should be used in M-step.
 #' @param lower.p fixed lower bound for success probabilities. If supplied as a single value, it will be assigned to all items;
 #'   It can also be specified as a numeric vector corresponding to each item. The default is 1e-4 for all items.
 #' @param upper.p fixed upper bound for success probabilities. If supplied as a single value, it will be assigned to all items;
 #'   It can also be specified as a numeric vector corresponding to each item. The default is 0.9999 for all items.
+#' @param lower.prior the lower bound for posteriori weights. The default value is -1, which means the lower bound is \eqn{1/2^K/100}.
+#' @param randomseed random seed for generating initial item parameters. The default random seed is 123456.
 #' @param smallNcorrection a numeric vector with two elements specifying the corrections applied when the expected number of
 #' examinees in some latent groups are too small. Particularly, if the expected no. of examinees is less than the second element,
 #' the first element and two times the first element will be added to the numerator and denominator of the closed-form solution of
 #' probabilities of success.
-#' @param randomseed random seed for generating initial item parameters. The default random seed is 123456.
+#' @param Mstep.warning Logical; indicating whether the warning message in Mstep, if any, should be output immediately.
+#' @param diagnosis Run in diagnostic mode? If it is 1 or 2, some intermediate results obtained in each iteration can be extracted.
+#' @param optimizer A string indicating which optimizer should be used in M-step.
+#' @param optim.control control options for optimizers in the M-step. Only available when \code{optimizer} is one specific optimization
+#' method, including \code{BFGS} from \link[stats]{optim}, \link[nloptr]{slsqp}, \link[Rsolnp]{solnp} and \link[alabama]{auglag}.
+#' For the \link[alabama]{auglag} method, \code{optim.control} specifies \code{control.outer}.
 #'
 #' @author {Wenchao Ma, Rutgers University, \email{wenchao.ma@@rutgers.edu} \cr Jimmy de la Torre, The University of Hong Kong}
 #' @seealso See \code{\link{autoGDINA}} for Q-matrix validation, item level model comparison and model calibration
 #' in one run; See \code{\link{itemfit}} for item fit analysis, \code{\link{Qval}} for Q-matrix validation,
-#' \code{\link{modelcomp}} for item level model comparison and \code{\link{simGDINA}} for data simulation
+#' \code{\link{modelcomp}} for item level model comparison and \code{\link{simGDINA}} for data simulation.
+#' Also see \code{gdina} in \pkg{CDM} package for the G-DINA model estimation.
 #'
 #' @return \code{GDINA} returns an object of class \code{GDINA}. S3 methods for \code{GDINA} objects
 #'  include \code{\link{extract}} for extracting various components, \code{\link{itemparm}}
@@ -334,11 +339,11 @@
 #' # item parameters
 #' # see ?itemparm
 #' itemparm(mod1) # item probabilities of success for each latent group
-#' itemparm(mod1, withSE = TRUE) # item probabilities of success + standard errors
+#' itemparm(mod1, withSE = TRUE) # item probabilities of success & standard errors
 #' itemparm(mod1, what = "delta") # delta parameters
 #' itemparm(mod1, what = "delta",withSE=TRUE) # delta parameters
-#' itemparm(mod1, what = "gs") # guess and slip parameters
-#' itemparm(mod1, what = "gs",withSE = TRUE) # guess and slip parameters + standard errors
+#' itemparm(mod1, what = "gs") # guessing and slip parameters
+#' itemparm(mod1, what = "gs",withSE = TRUE) # guessing and slip parameters & standard errors
 #'
 #' # person parameters
 #' # see ?personparm
@@ -635,9 +640,10 @@ GDINA <-
            mono.constraint = FALSE, group = NULL,
            empirical = !higher.order, att.prior = NULL, att.str = FALSE,
            lower.p = 0.0001,upper.p = 0.9999, smallNcorrection = c(0.0005,0.001),
-           nstarts = 1, conv.crit = 0.001, conv.type = "max.p.change",maxitr = 1000,
+           nstarts = 1, conv.crit = 0.001, lower.prior=-1,
+           conv.type = "max.p.change",maxitr = 1000,
            digits = 4,diagnosis = 0,Mstep.warning = FALSE,optimizer = "all",
-           randomseed = 123456)
+           randomseed = 123456,optim.control=list())
   {
     s1 <- Sys.time()
     GDINAcall <- match.call()
@@ -689,7 +695,7 @@ GDINA <-
     }
 
     # missing value indicator 1 - nonmissing; 0 - missing
-    missInd <- 1 - is.na(dat)
+    missInd <- 1L - is.na(dat)
     # print(dim(missInd))
     if (any(missInd == 0)) {
       # some missings individuals with one or fewer valid response are
@@ -786,6 +792,7 @@ GDINA <-
     for(j in 1:J){ #for each item
       designMlist[[j]] <- designmatrix(Kj[j],model[j])
     }
+
     while (dif > conv.crit && itr < maxitr)
     {
       item.parm_copy <- item.parm
@@ -832,7 +839,7 @@ GDINA <-
       optims <- Mstep(Kj=Kj, RN=RN, model=model, itmpar=item.parm, delta=deltas,
                       constr=mono.constraint,correction = smallNcorrection, lower.p = lower.p,
                       upper.p = upper.p, itr=itr+1,warning.immediate=Mstep.warning,
-                      optimizer = optimizer, designmatrices = designMlist)
+                      optimizer = optimizer, designmatrices = designMlist,optim.control = optim.control)
       item.parm <- optims$item.parm
       deltas <- optims$delta
 
@@ -860,9 +867,15 @@ GDINA <-
 
       if (higher.order == FALSE & empirical == TRUE)
       {
-        logprior <- likepost$logprior
+        lower.prior <- ifelse(lower.prior==-1,1/L/100,lower.prior)
+        priori <- exp(likepost$logprior)
+        priori[which(priori<lower.prior)] <- lower.prior
+        priori <- priori/sum(priori)
+        # cat("\nMin weight=",min(priori),"\n")
+        logprior <- log(priori)
       }else if (higher.order == FALSE & empirical == FALSE)
       {
+
         logprior <- log(att.prior)
       }else if(higher.order==TRUE){
         logprior <- HOlogprior

R/M.R

@@ -43,7 +43,7 @@ Mstep_optim <- function(ini,Mj,upper.pj,lower.pj,Kjj,
              llm=rep(c(qlogis(lower.pj),-1*qlogis(upper.pj)),each=2^Kjj),
              rrum=rep(c(log(lower.pj),-1*log(upper.pj)),each=2^Kjj))
 # print(model)
-  if(optimizer %in% c("Nelder-Mead","BFGS","CG","SANN","Brent")){
+  if(optimizer %in% c("Nelder-Mead","BFGS","CG")){
     ini <- initials_optim(ini=ini,Mj=Mj,upper.pj=upper.pj,
                           lower.pj=lower.pj,Kjj=Kjj,model=model)
   optims <- stats::constrOptim(ini,obj_fn,grad=gr_fn,method=optimizer,
@@ -94,7 +94,7 @@ Mstep_optim <- function(ini,Mj,upper.pj,lower.pj,Kjj,
               call. = FALSE,immediate. = warning.immediate)
     }
   }else{
-    stop(paste("Optimizer specification for item",j,"is not correct."),call. = FALSE)
+    stop(paste("Optimization method",optimizer,"is not available. Try BFGS, slsqp, solnp, auglag, Nelder-Mead, or CG."),call. = FALSE)
   }
 
     # print(optims)

R/M.mono.R

@@ -38,7 +38,7 @@ Mstep_optim_mono_gdina <- function(ini,Mj,upper.pj,lower.pj,Kjj,
                                    model,Nj,Rj,optimizer,j,itr,
                                    warning.immediate,warning.print){
 
-  # ini <- initials_optim_mono(ini=ini,Mj=Mj,bound.eps=bound.eps,Kjj=Kjj,model=model)
+   # ini <- initials_optim_mono(ini=ini,Mj=Mj,bound.eps=bound.eps,Kjj=Kjj,model=model)
 
   ineq_const_gdina_mono <- function(vdelta){
     partialorder <- partial_order2(Kjj)
@@ -53,6 +53,7 @@ Mstep_optim_mono_gdina <- function(ini,Mj,upper.pj,lower.pj,Kjj,
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [slsqp from nloptr] message:",optims$message),
               call. = FALSE,immediate. = warning.immediate)
     }
+    # cat("\n Item",j,"-slsqp\n")
   }else if(tolower(optimizer)=="solnp"){
     optims <- Rsolnp::solnp(pars=ini,fun = obj_fn,ineqfun = ineq_const_gdina_mono2,
                             ineqLB = rep(0,length(ineq_const_gdina_mono(ini))),
@@ -64,6 +65,8 @@ Mstep_optim_mono_gdina <- function(ini,Mj,upper.pj,lower.pj,Kjj,
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [solnp from Rsolnp] message: error code -",optims$convergence),
               call. = FALSE,immediate. = warning.immediate)
     }
+
+    # cat("\n Item",j,"-solnp\n")
   }else if(tolower(optimizer)=="auglag"){
 
     optims <- alabama::auglag(ini,obj_fn,gr=gr_fn,hin = ineq_const_gdina_mono2,
@@ -74,6 +77,10 @@ Mstep_optim_mono_gdina <- function(ini,Mj,upper.pj,lower.pj,Kjj,
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [auglag from alabama] message: error code -",optims$convergence),
               call. = FALSE,immediate. = warning.immediate)
     }
+    # cat("\n Item",j,"-auglag\n")
+  }else{
+    stop(paste("Optimization method",optimizer,"is not available. Try slsqp, solnp, or auglag."),call. = FALSE)
+
   }
 
   return(optims)
@@ -86,7 +93,7 @@ Mstep_optim_mono <- function(ini,Mj,upper.pj,
                              lower.pj,Kjj,model,
                              Nj,Rj,optimizer,j,itr,
                              warning.immediate,
-                             warning.print){
+                             warning.print,optim.control=list()){
 
   ini <- initials_optim_mono(ini=ini,Mj=Mj,upper.pj=upper.pj,
                              lower.pj=lower.pj,Kjj=Kjj,model=model)
@@ -113,12 +120,13 @@ Mstep_optim_mono <- function(ini,Mj,upper.pj,
                         c(qlogis(lower.pj),rep(0,length(ini)-1))),
                   rrum=c(rep(c(log(lower.pj),-1*log(upper.pj)),each=2^Kjj),
                          c(log(lower.pj),rep(0,length(ini)-1))))
-  if (optimizer %in% c("Nelder-Mead","BFGS","CG","SANN","Brent")){
+  if (optimizer %in% c("Nelder-Mead","BFGS","CG")){
 
     optims <- stats::constrOptim(ini,obj_fn,grad=gr_fn,method=optimizer,
                                  ui=ui.bfgs,ci=ci.bfgs[[model]],
                                  Nj=Nj,Rj=Rj,Mj=Mj,upper.pj=upper.pj,
-                                 lower.pj=lower.pj,Kjj=Kjj,model=model)
+                                 lower.pj=lower.pj,Kjj=Kjj,model=model,
+                                 control = optim.control)
     if (warning.print&&optims$convergence>0)  {
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [constrOptim using",optimizer,"] message:",optims$message),
               call. = FALSE,immediate. = warning.immediate)
@@ -130,7 +138,8 @@ Mstep_optim_mono <- function(ini,Mj,upper.pj,
 
     optims <- nloptr::slsqp(x0=ini,fn=obj_fn,gr=gr_fn,hin = ineq_fn_mono_1,lower=c(LB[,model]),
                             Nj=Nj,Rj=Rj,Mj=Mj,upper.pj=upper.pj,
-                            lower.pj=lower.pj,Kjj=Kjj,model=model)
+                            lower.pj=lower.pj,Kjj=Kjj,model=model,
+                            control = optim.control)
     if (warning.print&&optims$convergence<0)  {
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [slsqp from nloptr] message:",optims$message),
               call. = FALSE,immediate. = warning.immediate)
@@ -140,23 +149,24 @@ Mstep_optim_mono <- function(ini,Mj,upper.pj,
                             ineqLB = 0,LB=c(LB[,model]),ineqUB = +Inf,
                             control = list(trace=0,tol=1e-6),
                             Nj=Nj,Rj=Rj,Mj=Mj,upper.pj=upper.pj,
-                            lower.pj=lower.pj,Kjj=Kjj,model=model)
+                            lower.pj=lower.pj,Kjj=Kjj,model=model,
+                            control = optim.control)
     if (warning.print&&optims$convergence>0)  {
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [solnp from Rsolnp] message: error code -",optims$convergence),
               call. = FALSE,immediate. = warning.immediate)
     }
   }else if(tolower(optimizer)=="auglag"){
-
+    if(is.null(optim.control$trace)) optim.control$trace <- FALSE
     optims <- alabama::auglag(ini,obj_fn,gr=gr_fn,hin = ineq_fn_mono,
                               Nj=Nj,Rj=Rj,Mj=Mj,upper.pj=upper.pj,
                               lower.pj=lower.pj,Kjj=Kjj,model=model,
-                              control.outer = list(trace=FALSE))
+                              control.outer = optim.control)
     if (warning.print&&optims$convergence>0)  {
       warning(paste("Unsuccessful M-step optimization for item",j,"at iteration",itr,"\n optimizer [auglag from alabama] message: error code -",optims$convergence),
               call. = FALSE,immediate. = warning.immediate)
     }
   }else{
-    stop("Optimizer specification is not correct.",call. = FALSE)
+    stop(paste("Optimization method",optimizer,"is not available. Try BFGS, slsqp, solnp, auglag, Nelder-Mead, or CG."),call. = FALSE)
   }
   #if (tolower(optimizer)!="bfgs") print(optimizer)
   return(optims)
@@ -174,13 +184,13 @@ Moptim_set_mono <- function(ini,Mj,upper.pj,
       tryCatch({
         tryCatch({Mstep_optim_mono_gdina(ini=c(ini),Mj=Mj,upper.pj=upper.pj,
                                          lower.pj=lower.pj,Kjj=Kjj,model=model,
-                                         Nj=Nj,Rj=Rj,optimizer="solnp",
+                                         Nj=Nj,Rj=Rj,optimizer="slsqp",
                                          warning.immediate=warning.immediate,
                                          warning.print=warning.print)
         },error=function(e){
           Mstep_optim_mono_gdina(ini=c(ini),Mj=Mj,upper.pj=upper.pj,
                                  lower.pj=lower.pj,Kjj=Kjj,model=model,
-                                 Nj=Nj,Rj=Rj,optimizer="slsqp",
+                                 Nj=Nj,Rj=Rj,optimizer="solnp",
                                  warning.immediate=warning.immediate,
                                  warning.print=warning.print)
         })

R/Mstep.R

@@ -1,7 +1,7 @@
 
 Mstep <- function(Kj,RN,model,itmpar,delta,constr,
                   correction,lower.p,optimizer,
-                  upper.p,itr,warning.immediate,designmatrices){
+                  upper.p,itr,warning.immediate,designmatrices,optim.control=list()){
 
   J <- length(Kj)
 
@@ -89,7 +89,7 @@ Mstep <- function(Kj,RN,model,itmpar,delta,constr,
                                              Kjj=Kj[j],model=model[j],Nj=Nj,Rj=Rj,
                                              optimizer=optimizer,
                                              warning.immediate = warning.immediate,
-                                             warning.print=TRUE)
+                                             warning.print=TRUE,optim.control=optim.control)
           }
         }
       delta[[j]] <- optims$par

R/Msupplement.R

@@ -53,6 +53,9 @@ obj_fn <- function(vdelta,Nj,Rj,Mj,upper.pj,lower.pj,Kjj,model){ #objective func
     Pj <- Mj%*%vdelta
     Pj <- exp(Pj)
   }
+  Pj[Pj < .Machine$double.eps] <- .Machine$double.eps
+  Pj[Pj > (1 - .Machine$double.eps)] <- 1 - .Machine$double.eps
+
   -1*sum(log(Pj)*Rj+log(1-Pj)*(Nj-Rj))
 }