GDINA: The Generalized DINA Model Framework

# score function for probabilities
score_pj <- function(Xj,                   # a vector of item responses to item j
                    parloc.j,         # parameter locations for item j - H by 2^K matrix
                    catprob.j,        # a list with H elements giving the reduced catprob.parm for each nonzero category
                    logpost){          # log posterior N x 2^K

  Xj[is.na(Xj)] <- -1 # remove NA from the data
  Lj <- sapply(catprob.j,length)
  Kj <- log2(Lj)
  post <- exp(logpost) # posterior N x 2^K
  score.p <- vector("list",nrow(parloc.j)) #I(Xij>=x)/sjx - I(Xij=x-1)/[1-sjx]
  for(r in 1:length(score.p)){
    score.p[[r]] <- aggregateCol(post,parloc.j[r,])*(outer((Xj>=r),catprob.j[[r]],"/")-outer((Xj==r-1),(1-catprob.j[[r]]),"/"))
  }

  return(score.p)
}


# input is the estimation object
scorefunc2 <- function(object, item = NULL, ...){
  dat <- extract(object,"dat")

  Qc <- extract(object,"Qc")
  Q <- extract(object,"Q")

  J <- ncol(dat)
  par.loc <- eta(as.matrix(Q))
  post <- exp(indlogPost(object)) # posterior N x 2^K

  # urP <- itemparm(object,"itemprob",digits = 9) #unconditional reduced prob.
  C <- table(Qc[,1]) # number of categories
  S0 <- rep(1:J,C+1) # loc for each item including cat 0

  LC.Prob <- sequP(as.matrix(par.loc),as.matrix(object$catprob.matrix),C)
  c.LC.Prob <- LC.Prob$cPr # no cat 0
  u.LC.Prob <- LC.Prob$uPr # including cat 0
  scorej <- vector("list",J)
  if(is.null(item)) item <- 1:J
  for(j in item){
    Xj <- dat[,j]
    sj <- c.LC.Prob[which(Qc[,1]==j),,drop=FALSE] # processing function for item j
    # print(sj)
    locj <- par.loc[which(Qc[,1]==j),,drop=FALSE] # loc for item j
    Pj <- u.LC.Prob[which(S0==j),][-1,,drop=FALSE] # drop cat 0 - P(Xij=h|alpha)
    Pj0 <- u.LC.Prob[which(S0==j),][1,] # P(Xij=0|alpha)
    Dj <- list() #I(Xij=h)/P(Xij=h|alpha) - I(Xij=0)/P(Xij=0|alpha)
    for(h in seq_len(C[j])) Dj[[h]] <- outer((Xj==h),Pj[h,],"/")-outer((Xj==0),Pj0,"/")
    # print(head(Dj[[h]]))
    for(k in seq_len(C[j])) {
      scorejk <- 0
      for(h in seq_len(C[j])){
        if(k<=h) {
          tmp <- Dj[[h]]*matrix(Pj[h,]/sj[k,],nrow = nrow(Dj[[h]]),ncol = ncol(Dj[[h]]),byrow = TRUE)
        }else if(k==h+1){
          tmp <- Dj[[h]]*matrix((-1)*(Pj[h,]/(1-sj[k,])),nrow = nrow(Dj[[h]]),ncol = ncol(Dj[[h]]),byrow = TRUE)
        }else{
          tmp <- matrix(0,nrow = nrow(Dj[[h]]),ncol = ncol(Dj[[h]]))
        }
        scorejk <- scorejk + tmp
      }
      tmp <- post*scorejk
      # print(tmp[1:5,1:5])
      for(lj in sort(unique(locj[k,]))) scorej[[j]] <- cbind(scorej[[j]],rowSums(tmp[,which(locj[k,]==lj),drop=FALSE]))

    }

  }
  return(scorej)
}

# only suitable for DINA, DINO and GDINA
score_p <- function(object){
  model <- extract(object,"models_numeric")
  if(any(model>2)) stop("Score function for success probabilities only available for the G-DINA, DINA and DINO models.",call. = FALSE)
  Q <- extract(object,"Q")
  Qc <- extract(object,"Qc")
  NC <- nrow(Q)
  N <- extract(object,"nobs")
  J <- extract(object,"nitem")
  Kj <- rowSums(Q>0)
  if(extract(object,"sequential")){
    dat <- seq_coding(extract(object,"dat"),Qc)
  }else{
    dat <- extract(object,"dat")
  }

  pj <- as.matrix(extract(object,"catprob.matrix"))
  Lj <- 2^rowSums(Q>0)
  for(j in which(model %in% c(1,2))){#DINA or DINO
    pj[j,2] <- pj[j,Lj[j]]
    if (Lj[j]>2) pj[j,3:ncol(pj)] <- -1
  }

  scof <- scorefun(mX=as.matrix(dat),
                   mlogPost=as.matrix(extract(object,"logposterior.i")),
                   itmpar=pj,
                   parloc=extract(object,"eta"),
                   model=model)
  scorep <- scof$score
  ind <- scof$index + 1 # col 1: location; col 2: category
  score <- vector("list",NC)
  for (j in 1:NC){
    score[[j]] <- as.matrix(scorep[,ind[which(ind[,2]==j),1]]) #score function for P(\alpha_c) for cateogry j
  }
  names(score) <- extract(object,"item.names")
  lik <- exp(indlogLik(object))

  if(extract(object,"ngroup")>1){

    for(g in sort(unique(extract(object,"group")))){
      tmp <- ((lik-lik[,1])/colSums(c(extract(object,"posterior.prob")[g,])*t(lik)))[,-1]
      l <- length(score)+1
      # print(l)
      score[[l]] <- tmp
      score[[l]][which(extract(object,"group")!=g),] <- 0
      names(score)[l] <- paste0("G",g)
    }
  }else{
    if(any(extract(object,"att.dist")=="saturated")){
      score[[length(score)+1]] <- ((lik-lik[,1])/colSums(c(extract(object,"posterior.prob"))*t(lik)))[,-1]
    }else if(extract(object,"att.dist")=="higher.order"){
      if(extract(object,"higher.order")$model=="Rasch"){
        #To DO
      }

    }

  }
  return(score)
}


score_d <- function(object){
  Q <- extract(object,"Q")
  NC <- nrow(Q)
  N <- extract(object,"nobs")
  Kj <- rowSums(Q>0)
  pj <- extract(object,"catprob.parm")
  att <- extract(object,"attributepattern")
  linkfunc <- extract(object,"linkfunc")
  des <- extract(object,"designmatrix")
  if(extract(object,"sequential")){
    Qc <- extract(object,"Qc")
    dat <- seq_coding(extract(object,"dat"),Qc)
  }else{
    dat <- extract(object,"dat")
  }
  model <- extract(object,"models_numeric")
  scof <- scorefun(mX=as.matrix(dat),
                   mlogPost=as.matrix(extract(object,"logposterior.i")),
                   itmpar=as.matrix(extract(object,"catprob.matrix")),
                   parloc=extract(object,"eta"),
                   model=rep(0,NC))

  scorep <- scof$score
  ind <- scof$index + 1
  score <- vector("list",NC)
  for (r in 1:NC){
    scorepj <- as.matrix(scorep[,ind[which(ind[,2]==r),1]]) #score function for P(\alpha_c) for cateogry j

    if(linkfunc[r]=="identity"){
      score[[r]] <- scorepj%*%des[[r]]
    }else if(linkfunc[r]=="logit"){
      score[[r]] <- (scorepj*matrix(pj[[r]]*(1-pj[[r]]),nrow = N,ncol = length(pj[[r]]),byrow = TRUE))%*%des[[r]]
    }else if(linkfunc[r]=="log"){
      score[[r]] <- (scorepj*matrix(pj[[r]],nrow = N,ncol = length(pj[[r]]),byrow = TRUE))%*%des[[r]]
    }

  }
  names(score) <- extract(object,"item.names")
  lik <- exp(indlogLik(object))

  if(extract(object,"ngroup")>1){

    for(g in sort(unique(extract(object,"gr")))){
      tmp <- ((lik-lik[,1])/colSums(c(extract(object,"posterior.prob")[g,])*t(lik)))[,-1]
      l <- length(score)+1
      # print(l)
      score[[l]] <- tmp
      score[[l]][which(extract(object,"gr")!=g),] <- 0
      names(score)[l] <- extract(object,"group.label")[g]
    }
  }else{
    score[[length(score)+1]] <- ((lik-lik[,1])/colSums(c(extract(object,"posterior.prob"))*t(lik)))[,-1]
  }

  return(score)
}

# variance and SE of delta parameters for all models
OPG_d <- function(object,SE.type){
  Q <- extract(object,"Q")
  Qc <- extract(object,"Qc")
  J <- extract(object,"nitem")
  NC <- nrow(Q)
  NG <- extract(object,"ngroup")
  scorejh <- score_d(object) # a list of score function for delta with elements for each category
  np <- sapply(scorejh,ncol)[seq_len(NC)] # the last NG elements are for mixing proportions
# print(np)
  if(SE.type == 1){
    scorej <- vector("list",J)
    # scorejh <- scorejh[seq_len(IP.loc)]
    for(j in 1:J)  scorej[[j]] <- do.call(cbind,scorejh[which(Qc[,1]==j)])
    vars <- bdiag(lapply(scorej,inverse_crossprod))
  }else if(SE.type == 2){
    scorejh <- scorejh[seq_len(extract(object,"ncat"))] #for all categories
    vars <- inverse_crossprod(do.call(cbind,scorejh))
  }else if(SE.type==3){
    if(all(extract(object,"att.dist")=="saturated")){
      vars <- inverse_crossprod(do.call(cbind,scorejh))
      vars <- vars[1:sum(np),1:sum(np)]
    }else{
      stop("SEs base on the complete information are only available for saturated joint att. distributions.",call. = FALSE)
    }

  }
  covIndex <- data.frame(item = rep(Qc[,1],np),
                         itemcat = rep(Qc[,2],np),
                         cat = rep(1:NC,np),
                         par = unlist(lapply(np,seq_len)),
                         loc = 1:sum(np),row.names = NULL)
  se.vector <- sqrt(diag(vars))
  se <- vector("list",NC)
  for(i in 1:NC) se[[i]] <- se.vector[covIndex$loc[which(covIndex$cat==i)]]
  return(list(cov=vars,se=se,ind=covIndex))
}



# variance of prob parameters for all models
# delta method is used for additive type models
OPG_p <- function(object,SE.type){
  Q <- extract(object,"Q")
  Qc <- extract(object,"Qc")
  J <- extract(object,"nitem")
  NC <- nrow(Q)
  Kj <- rowSums(Q)
  NG <- extract(object,"ngroup")
  linkfunc <- extract(object,"linkfunc")
  des <- extract(object,"designmatrix")
  m <- extract(object,"models_numeric")



  if(all(m<=2)&all(m>=0)){

    scorejh <- score_p(object) # a list with elements for each category

    np <- sapply(scorejh,ncol)[seq_len(NC)] # the last NG elements are for mixing proportions

    if(SE.type == 1){
      # scorejh <- scorejh[seq_len(IP.loc)]
      scorej <- vector("list",J)
      for(j in 1:J)  scorej[[j]] <- do.call(cbind,scorejh[which(Qc[,1]==j)])
      vars <- bdiag(lapply(scorej,inverse_crossprod))
    }else if(SE.type == 2){
      scorejh <- scorejh[seq_len(extract(object,"ncat"))]
      vars <- inverse_crossprod(do.call(cbind,scorejh))
    }else if(SE.type==3){
      if(all(extract(object,"att.dist")=="saturated")){
        vars <- inverse_crossprod(do.call(cbind,scorejh))
        vars <- vars[1:sum(np),1:sum(np)]
      }else{
        stop("SEs base on the complete information are only available for saturated joint att. distributions.",call. = FALSE)
      }

    }
  }else{
    grad <- vector("list",NC)
    for (nc in 1:NC) {
      if (linkfunc[nc]=="identity") {
        grad[[nc]] <- des[[nc]]
      } else if (linkfunc[nc]=="logit") {
        grad[[nc]] <-diag(extract(object, "catprob.parm")[[nc]] * (1 - extract(object, "catprob.parm")[[nc]])) %*%des[[nc]]
        } else if (linkfunc[nc]=="log") {
        grad[[nc]] <- diag(extract(object, "catprob.parm")[[nc]]) %*% des[[nc]]
        }
    }
    g <- bdiag(grad)
    np <- sapply(grad,nrow)
    vars <- g %*% OPG_d(object,SE.type)$cov %*% t(g)
    }


  covIndex <- data.frame(item = rep(Qc[,1],np),
                         itemcat = rep(Qc[,2],np),
                         cat = rep(1:NC,np),
                         par = unlist(lapply(np,seq_len)),
                         loc = 1:sum(np),row.names = NULL)
  se.vector <- sqrt(diag(vars))
  se <- vector("list",NC)
  if(all(m<=2)&all(m>=0)){
  for(h in 1:NC) {

    if(m[h]==1) {
      se.pjh <- se.vector[covIndex$loc[which(covIndex$cat==h)]]
      se[[h]] <- c(rep(se.pjh[1],(2^Kj[h]-1)),se.pjh[2])
    }else if(m[h]==2){
      se.pjh <- se.vector[covIndex$loc[which(covIndex$cat==h)]]
      se[[h]] <- c(se.pjh[1],rep(se.pjh[2],(2^Kj[h]-1)))
    }else if(m[h]==0){
      se[[h]] <- se.vector[covIndex$loc[which(covIndex$cat==h)]]
    }
  }
    }else{
    for(h in 1:NC)
      se[[h]] <- se.vector[covIndex$loc[which(covIndex$cat==h)]]
  }
  return(list(cov=vars,se=se,ind=covIndex))
}