R/internal.functions.R

In NPCDTools: The Nonparametric Classification Methods for Cognitive Diagnosis

#=====Estimation of Proficiency Class via Hamming distance=====
hamming = function(Ideal,Y){
  M=nrow(Ideal)
  J=ncol(Ideal)
  N=nrow(Y)

  H.class=NULL
  h.ntie=0
  for (i in 1:N)
  {
    ham=apply(abs(matrix(rep(Y[i,], M), M, J, byrow=TRUE)-Ideal), 1, sum)
    min.ham=which(ham==min(ham))
    if (length(min.ham)!=1)
    {
      h.ntie=h.ntie+1
      min.ham=sample(min.ham,1,prob=rep(1/length(min.ham),length(min.ham))) ## temporarily fix ties
    }else h.ntie=h.ntie
    H.class=c(H.class, min.ham)
  }
  return(H.class)
}

#=====Estimation of Proficiency Class via weighted Hamming distance=====
whamming=function(Ideal,Y){
  M=nrow(Ideal)
  J=ncol(Ideal)
  N=nrow(Y)

  p.bar=apply(Y,2,mean)
  weight=1/(p.bar*(1-p.bar))
  WH.class=NULL
  ntie=0
  H=NULL
  for (i in 1:N)
  {
    ham=apply(matrix(rep(weight, M), M, J, byrow=TRUE)*abs(matrix(rep(Y[i,], M), M, J, byrow=TRUE)-Ideal), 1, sum)
    H=rbind(H,ham)
    min.ham=which(ham==min(ham))
    if (length(min.ham)!=1)
    {
      ntie=ntie+1
      min.ham=sample(min.ham,1,prob=rep(1/length(min.ham),length(min.ham))) ## temporarily fix ties
    }else ntie=ntie
    WH.class=c(WH.class, min.ham)
  }
  return(WH.class)
}

class.generate = function(K){
  M <- diag(K)
  for (i in 2:K){
    M <- rbind(M,t(apply(utils::combn(K,i),2,function(x){apply(M[x,],2,sum)})))
  }
  M <- rbind(0,M)
  return(M)
}

aepc.generate = function(mcq,O,J,K,key,Class){
  Q <- mcq[, -c(1:2)]
  Item.info <- mcq[,1:2]
  item.no <- mcq[,1]
  coded.op <- mcq[,2]
  num.coded <- tabulate(item.no)   # number of coded option  save.m
  no.options <- rep(O, J)

  # "scored" option
  eta.class <- matrix(0,J,2^K)
  for(j in 1:J){

    Qj <- Q[which(item.no==j),,drop=FALSE]  # won't change data type
    Kj <- rowSums(Qj)
    kj.order <- order(Kj,decreasing = TRUE)  #get location;
    coded.op.j <- coded.op[which(item.no==j)]

    if (num.coded[j]>1){

      key.loc <- which(coded.op.j==key[j])
      if(key.loc!=kj.order[1]){
        kj.order <- c(key.loc,setdiff(kj.order,key.loc))  # make the answer goes first
      }


      Qj <- Qj[kj.order,]
      #et.label[[j]] <- c(apply(Qj,1,paste,collapse=""),paste(rep("*",ncol(Qj)),collapse = ""))
      A <- Class%*%t(Qj)
      B <- rbind(0,t(1*(A==(matrix(1,2^K)%*%t(rowSums(Qj))))))
      #matrix(unlist(lapply(apply(B,2,function(x){which(x==max(x))}),function(y){if (length(y)<O){c(y,rep(0,O-length(y)))}else{y=y}})),nrow=2^K,byrow=TRUE)

      #tmp <- eta(Qj)
      eta.j <- apply(B,2,which.max)
      max.j <- max(eta.j)
      eta.j <- eta.j - 1
      eta.j[eta.j==0] <- max.j
      eta.j <- num.coded[j]+1-eta.j


    }else{

      A <- Class%*%t(Qj)
      B <- rbind(0,t(1*(A==(matrix(1,2^K)%*%t(rowSums(Qj))))))

      eta.j <- apply(B,2,which.max)
      eta.j <- eta.j-1
    }

    eta.class[j,] <- eta.j

  }
  return(eta.class)
}

score.option = function(mcq,O){


  J = max(mcq[,1])
  Q <- mcq[, -c(1:2)]

  Item.info <- mcq[,1:2]
  item.no <- mcq[,1]
  coded.op <- mcq[,2]
  num.coded <- tabulate(item.no)   # number of coded option  save.m

  op <- vector(mode="list",length=J)

  mc.q = as.data.frame(mcq)
  save.m <- c()
  key <- c()
  for (i in 1:J){
    key[i] <- mc.q$Option[min(which(mc.q$Item==i))]
    save.m <- c(save.m,length(which(mc.q$Item==i)))
  }


  for (i in 1:J){

    Qj <- Q[which(item.no==i),,drop=FALSE]  # won't change data type
    Kj <- rowSums(Qj)
    kj.order <- order(Kj,decreasing = TRUE)  #get location;
    coded.op.j <- coded.op[which(item.no==i)]

    sub.info <- Item.info[which(Item.info[,1]==i),,drop=FALSE]

    if (num.coded[i]>1){

      key.loc <- which(coded.op.j==key[i])
      if(key.loc!=kj.order[1]){
        kj.order <- c(key.loc,setdiff(kj.order,key.loc))  # make the answer goes first
      }

      if (num.coded[i]==O){
        g = c(O:1)
      } else {g = c(c(nrow(sub.info):1),rep(0,O-nrow(sub.info)))}


      ans = c(sub.info[,2][kj.order],setdiff(c(1:O),sub.info[,2][kj.order]))
      score <- cbind(ans,g)


    } else {
      ans =c(unname(sub.info[,2]),setdiff(c(1:O),sub.info[,2]))
      g = c(1,rep(0,(O-1)))
      score <- cbind(ans,g)
    }

    op[[i]] <- score

  }

  return(op)
}

AlphaPermute <- function(dim) {
  
  alpha <- matrix(c(0, 1), 2, 1)
  
  for (i in 1:(dim - 1)) {
    alpha <- rbind(alpha, alpha)
    alpha <- cbind(alpha, c(rep(0, 2 ^ i), rep(1, 2 ^ i)))
  }
  
  return(alpha)
}


#--------- MCNPC Internal Function-----------


# Higher-order attribute structure:
# Student higher-order ability delta_i is generated from N(0, 1).
# Attribute difficulty lambda0_k is generated from N(0, 1),
# and attribute discrimination lambda1_k is generated from a lognormal distribution.
att.generate <- function(K, N,
                         att.structure = c("unif", "mvn", "higher_order"),
                         cor.att = 0.5,
                         lambda0 = NULL,
                         lambda1 = NULL,
                         sdlog.lambda1 = 0.5) {
  
  att.structure <- match.arg(att.structure)
  Class <- class.generate(K)
  
  if (att.structure == "unif") {
    
    true.att <- Class[sample(1:nrow(Class), N, replace = TRUE), , drop = FALSE]
    rownames(true.att) <- NULL
    
  } else if (att.structure == "mvn") {
    
    mu <- rep(0, K)
    sigma <- matrix(cor.att, nrow = K, ncol = K)
    diag(sigma) <- 1
    
    if (any(eigen(sigma)$values <= 0)) {
      stop("The correlation matrix is not positive definite. Try a different cor.att.")
    }
    
    prob.att <- MASS::mvrnorm(N, mu = mu, Sigma = sigma)
    true.att <- ifelse(prob.att > 0, 1, 0)
    
  } else if (att.structure == "higher_order") {
    
    delta <- stats::rnorm(N, mean = 0, sd = 1)
    
    if (is.null(lambda0)) {
      lambda0 <- stats::rnorm(K, mean = 0, sd = 1)
    }
    
    if (is.null(lambda1)) {
      lambda1 <- stats::rlnorm(K, meanlog = 0, sdlog = sdlog.lambda1)
    }
    
    if (length(lambda0) != K | length(lambda1) != K) {
      stop("lambda0 and lambda1 must both have length K.")
    }
    
    eta <- outer(delta, lambda1, "*") -
      matrix(lambda1 * lambda0, nrow = N, ncol = K, byrow = TRUE)
    
    p.att <- stats::plogis(eta)
    
    true.att <- matrix(
      stats::rbinom(N * K, size = 1, prob = as.vector(p.att)),
      nrow = N,
      ncol = K
    )
  }
  
  colnames(true.att) <- paste0("A", 1:K)
  return(true.att)
}


class.generate = function(K) {
  M <- diag(K)
  for (i in 2:K) {
    M <- rbind(M, t(apply(utils::combn(K, i), 2, function(x) {
      apply(M[x, ], 2, sum)
    })))
  }
  M <- rbind(0, M)
  return(M)
}

epc.generate = function(mcQ,H,key,LS = NULL){
  # for ideal responses
  
  Q <- mcQ[, -c(1:2),drop = FALSE]
  Item.info <- mcQ[,1:2]
  item.no <- mcQ[,1]
  coded.op <- mcQ[,2]
  num.coded <- tabulate(item.no)   # number of coded option  save.m
  
  J = length(unique(item.no))
  K = ncol(Q)
  if (is.null(LS)){
    Class <- class.generate(K)
    M = 2^K
  } else {
    Class <- LS
    M = nrow(LS)
  }
  
  if(length(H) == 1){
    no.options <- rep(H, J)
  } else {
    no.options <- H
  }
  
  
  # "scored" option
  eta.class <- matrix(0,J,M)
  # B.list = vector(length = J,mode = "list")
  for(j in 1:J){
    
    j.id = unique(item.no)[j]
    
    Qj <- Q[which(item.no==j.id),,drop=FALSE]  # won't change data type
    row.names(Qj) = NULL
    coded.op.j <- coded.op[which(item.no==j.id)]
    
    
    if (num.coded[j.id]>1){
      
      kj.order <- l.q(Qj)
      
      Qj <- Qj[order(kj.order),]
      A <- Class%*%t(Qj)
      B <- t(1*(A==(matrix(1,M)%*%t(rowSums(Qj)))))
      
      eta.j = apply(B,2,function(x){
        l = which(x==1)
        ifelse(length(l)==0,0,max(l))})
      
    }else{
      
      A <- Class%*%t(Qj)
      B <- rbind(0,t(1*(A==(matrix(1,M)%*%t(rowSums(Qj))))))
      
      eta.j <- apply(B,2,which.max)
      eta.j <- eta.j-1
    }
    
    eta.class[j,] <- eta.j
    
  }
  return(eta.class)
}


l.q <- function(Qj){
  ### requirement of input Qj
  # a matrix
  # the first row must be the key
  # no IDs
  
  # rules to order coded options
  # 1. The key is at the highest level
  # 2. coded distractors that require more skills are at the higher level
  # 3. coded distractors that require the same number of skills, 1 at the inconsistent entry......
  
  Qj_nokey <- Qj[-1,,drop = FALSE]
  Kj <- rowSums(Qj_nokey)
  Kj.order <- rank.Kj <- rank(Kj)
  if(any(duplicated(rank.Kj))){
    # currently, only one set of ties is considered
    # at most three coded distractors are considered
    loc.tie = c(which(duplicated(rank.Kj))[1]-1,which(duplicated(rank.Kj)))
    
    Qj_tie <- Qj_nokey[loc.tie,]
    tie.order <- order(apply(Qj_tie, 1, function(x) {
      as.integer(paste(x, collapse = ""))
    }))
    
    
    Kj.order[loc.tie] = mean(rank.Kj[loc.tie]-tie.order) + tie.order
    
  } else{
    # all coded distractors require different numbers of skills
    Kj.order <- rank(Kj)
  }
  
  
  Kj.order <- c(nrow(Qj),Kj.order)
  
  return(Kj.order)
  
}


prob.generate = function(mcQ, H, att, lambda) {
  # att must be a matrix
  
  item.no <- unique(mcQ[, 1])
  J = length(item.no)
  K = ncol(mcQ) - 2
  N = nrow(att)
  
  LatentClass = class.generate(K)
  
  Ideal <- LatentClass %*% t(mcQ[, 3:(2 + K)])
  Ideal.met <- 1 * (Ideal == (matrix(1, 2^K) %*% t(rowSums(mcQ[, 3:(2 + K)]))))
  
  prob <- vector(length = 2^K, mode = "list")
  for (i in 1:2^K) {
    # 0.06; 0.82; 0.25: same as Jimmy
    #op <- mc.q[,1:2][which(Ideal.conj[i,]>0),,drop=FALSE]
    sub.r <- matrix(1 / H, nrow = J, ncol = H)
    i.eta <- Ideal.met[i, ]
    
    for (j in 1:J) {
      q = lambda[j]
      
      work.q <- mcQ[which(mcQ[, 1] == item.no[j]), , drop = FALSE]
      work.eta <- i.eta[which(mcQ[, 1] == item.no[j])]
      
      if (sum(work.eta) == 1 || work.eta[1] == 1) {
        sub.r[j, ] = (1 - q) / (H - 1)
        sub.r[j, work.q[which(work.eta == 1)[1], 2]] = q
      } else if (sum(work.eta) > 1) {
        possible.improper.Q = work.q[work.eta == 1, -(1:2)]
        
        ori.order <-
          order(rowSums(possible.improper.Q), decreasing = TRUE)
        all.include = all(apply(possible.improper.Q, 1, function(row)
          all(row <= possible.improper.Q[ori.order[1], ])))
        
        if (all.include) {
          work.eta[which(work.eta == 1)[-ori.order[1]]] = 0
        }
        
        sub.r[j, ] = (1 - q) / (H - 1)
        sub.r[j, work.q[which(work.eta == 1), 2]] = (q + (1 - q) * (length(which(
          work.eta == 1
        )) - 1) / (H - 1)) / length(which(work.eta == 1))
      }
    }
    prob[[i]] <- sub.r
  }
  prob.each <- vector(length = N, mode = "list")
  for (i in 1:N){
    prob.each[[i]] = prob[[which(apply(LatentClass,1,function(y){
      all.equal(att[i,],y)
    })==TRUE)]]
  }
  return(prob.each)
}


score.option = function(mcQ, H) {
  J = length(unique(mcQ[, 1]))
  Q <- mcQ[,-c(1:2), drop = FALSE]
  K = ncol(Q)
  
  Item.info <- mcQ[, 1:2, drop = FALSE]
  item.no <- mcQ[, 1]
  coded.op <- mcQ[, 2]
  num.coded <- tabulate(item.no)   # number of coded option  save.m
  
  op <- vector(mode = "list", length = J)
  
  mc.q = as.data.frame(mcQ)
  colnames(mc.q) = c("Item", "Option", paste0("K", 1:K))
  
  save.m <- c()
  key <- c()
  for (i in 1:J) {
    j.id = unique(item.no)[i]
    
    key[i] <- mc.q$Option[min(which(mc.q$Item == j.id))]
    save.m <- c(save.m, length(which(mc.q$Item == j.id)))
  }
  
  if(length(H) == 1){
    no.options <- rep(H, J)
  } else {
    no.options <- H
  }
  
  for (i in 1:J) {
    
    O = no.options[i]
    
    j.id = unique(item.no)[i]
    
    Qj <- Q[which(item.no == j.id), , drop = FALSE]  # won't change data type
    
    ### wrong!!! the new index was updated
    # Kj <- rowSums(Qj)
    # kj.order <- order(Kj, decreasing = TRUE)  #get location;
    
    kj.order <- l.q(Qj)  #get location;
    coded.op.j <- coded.op[which(item.no == j.id)]
    
    sub.info <- Item.info[which(Item.info[, 1] == j.id), , drop = FALSE]
    
    if (num.coded[j.id] > 1) {
      key.loc <- which(coded.op.j == key[i])
      # if (key.loc != kj.order[1]) {
      #   kj.order <-
      #     c(key.loc, setdiff(kj.order, key.loc))  # make the answer goes first
      # }
      
      if (num.coded[j.id] == O) {
        g = c(O:1)
      } else {
        g = c(c(nrow(sub.info):1), rep(0, O - nrow(sub.info)))
      }
      
      
      ans = c(sub.info[, 2][order(kj.order,decreasing = TRUE)], setdiff(c(1:O), sub.info[, 2][order(kj.order,decreasing = TRUE)]))
      score <- cbind(ans, g)
      
      
    } else {
      ans = c(unname(sub.info[, 2]), setdiff(c(1:O), sub.info[, 2]))
      g = c(1, rep(0, (O - 1)))
      score <- cbind(ans, g)
    }
    
    op[[i]] <- score
    
  }
  
  return(op)
}