R/module5_mst.R

In xluo11/xxIRT: Item Response Theory and Computer-Based Testing

#' Computerized Multistage Testing (MST)
#' @name mst
#' @examples
#' \dontrun{
#' ## generate item pool
#' num_item <- 300
#' pool <- with(model_3pl_gendata(1, num_item), data.frame(a=a, b=b, c=c))
#' pool$id <- 1:num_item
#' pool$content <- sample(1:3, num_item, replace=TRUE)
#' pool$time <- round(rlnorm(num_item, 4, .3))
#' pool$group <- sort(sample(1:round(num_item/3), num_item, replace=TRUE))
#' 
#' ## ex. 1: 1-2-2 MST, 2 panels, topdown
#' ## 20 items in total and 10 items in content area 1 in each route
#' ## maximize info. at -1 and 1 for easy and hard routes
#' x <- mst(pool, "1-2-2", 2, 'topdown', len=20, max_use=1)
#' x <- mst_obj(x, theta=-1, indices=1:2)
#' x <- mst_obj(x, theta=1, indices=3:4)
#' x <- mst_constraint(x, "content", 10, 10, level=1)
#' x <- mst_assemble(x, timeout=5)
#' plot(x, byroute=TRUE)
#' for(p in 1:x$num_panel)
#'   for(r in 1:x$num_route) {
#'      route <- paste(x$route[r, 1:x$num_stage], collapse='-')
#'      count <- sum(mst_get_items(x, panel_ix=p, route_ix=r)$content==1)
#'      cat('panel=', p, ', route=', route, ': ', count, ' items in content area #1\n', sep='')
#'   }
#' 
#' ## ex. 2: 1-2-3 MST, 2 panels, bottomup,
#' ## remove two routes with large theta change: 1-2-6, 1-3-4 
#' ## 10 items in total and 4 items in content area 2 in each module
#' ## maximize info. at -1, 0 and 1 for easy, medium, and hard modules
#' x <- mst(pool, "1-2-3", 2, 'bottomup', len=10, max_use=1)
#' x <- mst_route(x, c(1, 2, 6), "-")
#' x <- mst_route(x, c(1, 3, 4), "-")
#' x <- mst_obj(x, theta= 0, indices=c(1, 5))
#' x <- mst_obj(x, theta=-1, indices=c(2, 4))
#' x <- mst_obj(x, theta= 1, indices=c(3, 6))
#' x <- mst_constraint(x, "content", 4, 4, level=2)
#' x <- mst_assemble(x, timeout=10) 
#' plot(x, byroute=FALSE)
#' for(p in 1:x$num_panel)
#'   for(m in 1:x$num_module){
#'     count <- sum(mst_get_items(x, panel_ix=p, module_ix=m)$content==2)
#'     cat('panel=', p, ', module=', m, ': ', count, ' items in content area #2\n', sep='')
#'   }
#'  
#' ## ex.3: same with ex.2 (w/o content constraints), but group-based  
#' x <- mst(pool, "1-2-3", 2, 'bottomup', len=12, max_use=1, group="group")
#' x <- mst_route(x, c(1, 2, 6), "-")
#' x <- mst_route(x, c(1, 3, 4), "-")
#' x <- mst_obj(x, theta= 0, indices=c(1, 5))
#' x <- mst_obj(x, theta=-1, indices=c(2, 4))
#' x <- mst_obj(x, theta= 1, indices=c(3, 6))
#' x <- mst_assemble(x, timeout=10)
#' plot(x, byroute=FALSE)
#' for(p in 1:x$num_panel)
#'   for(m in 1:x$num_module){
#'     items <- mst_get_items(x, panel_ix=p, module_ix=m)
#'     cat('panel=', p, ', module=', m, ': ', length(unique(items$id)), ' items from ', 
#'         length(unique(items$group)), ' groups\n', sep='')
#'   }
#'   
#' ## ex.4: 2 panels of 1-2-3 top-down design
#' ## 20 total items and 10 items in content area 3
#' ## 6+ items in stage 1 & 2
#' x <- mst(pool, "1-2-3", 2, "topdown", len=20, max_use=1)
#' x <- mst_route(x, c(1, 2, 6), "-")
#' x <- mst_route(x, c(1, 3, 4), "-")
#' x <- mst_obj(x, theta=-1, indices=1)
#' x <- mst_obj(x, theta=0, indices=2:3)
#' x <- mst_obj(x, theta=1, indices=4)
#' x <- mst_constraint(x, "content", 10, 10, level=3)
#' x <- mst_stage_length(x, 1:2, min=6)
#' x <- mst_assemble(x, timeout=15)
#' head(x$items)
#' plot(x, byroute=FALSE)
#' for(p in 1:x$num_panel)
#'   for(s in 1:x$num_stage){
#'     items <- mst_get_items(x, panel_ix=p, stage_ix=s)
#'     cat('panel=', p, ', stage=', s, ': ', length(unique(items$id)), ' items\n', sep='')
#'   }
#' 
#' ## ex.5: same with ex.4, but use RDP instead of stage length to control routing errors
#' x <- mst(pool, "1-2-3", 2, "topdown", len=20, max_use=1)
#' x <- mst_route(x, c(1, 2, 6), "-")
#' x <- mst_route(x, c(1, 3, 4), "-")
#' x <- mst_obj(x, theta=-1, indices=1)
#' x <- mst_obj(x, theta=0, indices=2:3)
#' x <- mst_obj(x, theta=1, indices=4)
#' x <- mst_constraint(x, "content", 10, 10, level=3)
#' x <- mst_rdp(x, 0, 2:3, .1)
#' x <- mst_module_mininfo(x, 0, 5, 2:3)
#' x <- mst_assemble(x, timeout=15)
#' plot(x, byroute=FALSE)
#' }
NULL

#' @rdname mst
#' @description \code{mst} creates a multistage (MST) object for assembly
#' @param pool the item pool (data.frame)
#' @param design the MST design (string): e.g., "1-3", "1-2-2", "1-2-3"
#' @param num_panel the number of panels (integer)
#' @param method the design method (string): 'topdown' or 'bottomup'
#' @param len the module/route length (integer)
#' @param max_use the maximum selection of items (integer)
#' @param group the grouping variable (string or vector)
#' @details
#' There are two methods for designing a MST. The bottom-up approach adds objectives 
#' and constraints on individual modules, whereas the topdown approach adds objectives
#' and constraints directly on routes.
#' @export
mst <- function(pool, design, num_panel, method=c('topdown', 'bottomup'), len=NULL, max_use=NULL, group=NULL, ...){
  method <- match.arg(method)
  design <- as.integer(unlist(strsplit(design, split="-")))
  num_stage <- length(design)
  num_module <- sum(design)
  opts <- list(...)
  if(is.null(opts$D)) opts$D <- 1.702
  
  # module-index map
  module <- NULL
  for(s in 1:num_stage)
    for(m in 1:design[s])
      module <- rbind(module, c(stage=s, module=m))
  module <- data.frame(module, index=1:nrow(module))
  
  # route-index map
  route <- list()
  for(i in 1:num_stage)
    route[[i]] <- module[module$stage == i, "index"]
  route <- expand.grid(route)
  colnames(route) <- paste("stage", 1:num_stage, sep="")
  route$index <- 1:nrow(route)
  num_route <- nrow(route)
  
  # ata
  x <- list(pool=pool, design=design, method=method, num_item=nrow(pool), num_panel=num_panel, 
            num_stage=num_stage, num_module=num_module, num_route=num_route, module=module, route=route, 
            ata=ata(pool, num_form=num_panel*num_module, group=group), opts=opts)
  class(x) <- "mst"

  # constraint: test length
  if(!is.null(len) && length(len) == 1) x <- mst_constraint(x, 1, len, len)
  if(!is.null(len) && length(len) == 2) x <- mst_constraint(x, 1, len[1], len[2])
  if(!is.null(len) && length(len) > 2) stop("the length argument is too long.")
  # constraint: max_use
  if(!is.null(max_use)) x$ata <- ata_item_use(x$ata, max=max_use)
  # constraint: minimum stage length
  x <- mst_stage_length(x, 1:num_stage, min=1)
  
  # constraint: no item reuse in the same route
  mat <- matrix(0, nrow=x$num_item*x$num_route*x$num_panel, ncol=x$ata$num_lpvar)
  for(p in 1:x$num_panel)
    for(i in 1:x$num_item){
      ind <- as.matrix(i + (x$route[,1:x$num_stage] - 1) * x$num_item + (p - 1) * x$num_module * x$num_item)
      for(j in 1:nrow(ind))
        mat[j+(i-1)*nrow(ind), ind[j,]] <- 1
    }
  dir <- rep('<=', nrow(mat))
  rhs <- rep(1, nrow(mat))
  x$ata <- ata_append_constraints(x$ata, mat, dir, rhs)

  x
}

#' @rdname mst
#' @description \code{mst_route} adds/removes a route to/from the MST
#' @param x the MST object
#' @param route a MST route represented by a vector of module indices 
#' @param op "+" to add a route and "-" to remove a route
#' @export
mst_route <- function(x, route, op=c("+", "-")){
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  op <- match.arg(op)
  index <- apply(x$route[, 1:x$num_stage], 1, function(r) all(r == route))
  
  if(op == "+") {
    if(any(index)) stop("the route already exists")
    if(!all(route %in% 1:x$num_module)) stop("invalid route: module index is out of bound.")
    x$route <- rbind(x$route, c(route, NA))
  } else if(op == "-") {
    if(!any(index)) stop("the route hasn't been added yet")
    x$route <- x$route[!index, ]
  }
  
  # reindex routes by stages
  index <- apply(x$route[, 1:x$num_stage], 1, function(r) sum(r * 10^(x$num_stage - 1:x$num_stage)))
  x$route <- x$route[order(index), ]
  x$route$index <- 1:nrow(x$route)
  x$num_route <- nrow(x$route)
  
  x
}

#' @rdname mst
#' @description \code{mst_get_indices} maps the input indices to the actual indices
#' @keywords internal
mst_get_indices <- function(x, indices){
  if(x$method == 'topdown'){
    if(is.null(indices)) indices <- x$route[, 1:x$num_stage] else indices <- subset(x$route, x$route$index %in% indices)[, 1:x$num_stage]
  } else if(x$method == 'bottomup') {
    if(is.null(indices)) indices <- data.frame(module=1:x$num_module) else indices <- data.frame(module=indices)
  }
  indices
}

#' @rdname mst
#' @description \code{mst_obj} adds objective functions to the MST
#' @param theta a theta point or interval over which the TIF is optimized
#' @param indices the indices of the route (topdown) or the module (bottomup) where objectives are added
#' @param target the target values of the TIF objectives. \code{NULL} for maximization
#' @export
mst_obj <- function(x, theta, indices=NULL, target=NULL, ...) {
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  indices <- mst_get_indices(x, indices)
  theta <- round(theta, 2)
  
  for(i in 1:x$num_panel) {
    for(j in 1:nrow(indices)) {
      f <- unlist(indices[j, ]) + (i - 1) * x$num_module
      if(is.null(target) || is.na(target)) {
        x$ata <- ata_obj_relative(x$ata, theta, mode="max", forms=f, collapse=TRUE, ...)
      } else {
        x$ata <- ata_obj_absolute(x$ata, theta, target=target, forms=f, collapse=TRUE, ...)
      }      
    }
  }
  
  x
}

#' @rdname mst
#' @description \code{mst_constraint} adds constraints to the MST
#' @param coef the coefficients of the constraint
#' @param level the constrained level, \code{NA} for quantitative variable
#' @param min the lower bound of the constraint
#' @param max the upper bound of the constraint
#' @export
mst_constraint <- function(x, coef, min=NA, max=NA, level=NULL, indices=NULL){
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  indices <- mst_get_indices(x, indices)
  
  for(i in 1:x$num_panel){
    for(j in 1:nrow(indices)){
      f <- unlist(indices[j,] + (i - 1) * x$num_module)
      x$ata <- ata_constraint(x$ata, coef, min, max, level, forms=f, collapse=TRUE)
    }
  }
  
  x
}


#' @rdname mst
#' @description \code{mst_stage_length} sets length limits on stages
#' @param stages the stage indices
#' @export
mst_stage_length <- function(x, stages, min=NA, max=NA){
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  if(length(min) == 1) min <- rep(min, length(stages))
  if(length(max) == 1) max <- rep(max, length(stages))
  if(length(stages) != length(min) || length(stages) != length(max))
    stop("different lengths in stage, min and max")
  
  for(i in 1:length(stages)){
    if(!stages[i] %in% 1:x$num_stage) stop("invalid stage input")
    f <- subset(x$module, x$module$stage == stages[i])$index
    f <- as.vector(outer(f, (1:x$num_panel - 1) * x$num_module, "+"))
    x$ata <- ata_constraint(x$ata, 1, min[i], max[i], forms=f, collapse=FALSE)
  }
  
  x
}


#' @rdname mst
#' @description \code{mst_rdp} anchors the routing decision point (rdp) between adjacent modules
#' @param tol tolerance parameter (numeric)
#' @importFrom stats aggregate
#' @export
mst_rdp <- function(x, theta, indices, tol=0) {
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  if(length(theta) != 1) stop("rdp is not a single theta point")
  if(length(indices) != 2 || abs(indices[1] - indices[2]) != 1) stop("modules are not adjacent") 
  
  info <- round(aggregate(model_3pl_info(theta, x$pool$a, x$pool$b, x$pool$c, D=x$opts$D)[1, ], by=list(group=x$ata$group), sum)[, 2], 2)
  coef <- c(info, -1 * info)
  for(i in 1:x$num_panel)
    x$ata <- ata_constraint(x$ata, coef, -tol, tol, forms=indices + (i - 1) * x$num_module, collapse=TRUE)
  
  x
}


#' @rdname mst
#' @description \code{mst_module_mininfo} sets the minimum information for modules
#' @param thetas theta points, a vector
#' @importFrom stats aggregate
#' @export
mst_module_info <- function(x, thetas, min, max, indices) {
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  if(any(indices < 1 | indices > x$num_module)) stop("invalid module index")
  if(length(min) == 1) min <- rep(min, length(thetas))
  if(length(max) == 1) max <- rep(max, length(thetas))
  if(length(min) != length(thetas) || length(max) != length(thetas)) stop('min/max has a different length from thetas')
  
  for(i in 1:length(thetas)){
    info <- with(x$pool, model_3pl_info(thetas[i], a, b, c, D=x$opts$D))[1, ]
    coef <- aggregate(info, by=list(group=x$ata$group), sum)[, 2]
    coef <- round(coef, 2)
    for(j in 1:x$num_panel)
      x$ata <- ata_constraint(x$ata, coef, min=min[i], max=max[i], forms=indices+(j-1)*x$num_module)
  }

  x
}


#' @rdname mst
#' @description \code{mst_assemble} assembles the mst
#' @export
mst_assemble <- function(x, ...){
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  opts <- list(...)
  
  solver <- ifelse(is.null(opts$solver), 'lpsolve', opts$solver)
  x$ata <- ata_solve(x$ata, as.list=FALSE, ...)
  
  if(!is.null(x$ata$items)) {
    items <- x$ata$items
    items$module <- (items$form - 1) %% x$num_module + 1
    items$panel <- ceiling(items$form / x$num_module)
    items$stage <- x$module$stage[match(items$module, x$module$index)]
    items$form <- NULL
    x$items <- items
  }
  
  x
}


#' @rdname mst
#' @param ... further arguments
#' @export
print.mst <- function(x, ...){
  cat("The MST design has", x$num_stage, "stages,", x$num_module, "modules, and", x$num_route, "routes:\n")
  cat("route map:\n")
  print(x$route)
  if(!is.null(x$items)){
    cat("\nAssembled forms:\n")
    items <- x$items
    if(!is.data.frame(x$items)) items <- Reduce(rbind, items, NULL)
    if(nrow(items) > 10){
      print(items[1:5, ])
      cat("...\n")
      print(items[-4:0 + nrow(items),])
    } else {
      print(items)
    }
    cat("See more results in 'items'.")
  } else {
    cat("MST hasn't been assembled yet.")
  }
  invisible(x)
}


#' @rdname mst
#' @details 
#' \code{plot.mst} draws module information functions when \code{byroute=FALSE}
#' and route information functions when \code{byroute=TRUE}
#' @import ggplot2
#' @export
plot.mst <- function(x, ...){
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  if(is.null(x$items)) stop('the mst has not been assembled yet.')
  
  opts <- list(...)
  if(is.null(opts$byroute)) opts$byroute <- FALSE
  if(is.null(opts$theta)) opts$theta <- round(seq(-3, 3, .1), 1)
  
  data <- NULL
  if(opts$byroute) {
    for(i in 1:x$num_route){
      for(j in 1:x$num_panel){
        items <- mst_get_items(x, panel_ix=j, route_ix=i)
        info <- with(items, rowSums(model_3pl_info(opts$theta, a, b, c, D=x$opts$D)))
        data <- rbind(data, data.frame(t=opts$theta, info=info, panel=j, route=i))
      }
    }
    data$panel <- factor(paste("Panel", data$panel))
    data$route <- factor(data$route, levels=1:x$num_route, labels=apply(x$route[, 1:x$num_stage], 1, paste, collapse="-"))
    g <- ggplot(data, aes_string(x="t", y="info", color="route")) + 
      geom_line() + xlab(expression(theta)) + ylab("Information") + 
      theme_bw() + theme(legend.key=element_blank()) +
      guides(color=guide_legend("Routes")) +
      facet_grid(. ~ panel)
  } else {
    for(i in 1:x$num_panel){
      for(j in 1:x$num_module){
        items <- mst_get_items(x, panel_ix=i, module_ix=j)
        info <- with(items, rowSums(model_3pl_info(opts$theta, a, b, c, D=x$opts$D)))
        data <- rbind(data, data.frame(t=opts$theta, info=info, panel=items$panel[1], stage=items$stage[1], module=items$module[1]))
      }
    }
    data$panel <- factor(paste("Panel", data$panel))
    data$stage <- factor(paste("Stage", data$stage))
    data$module <- factor(paste("Module", data$module))
    g <- ggplot(data, aes_string(x="t", y="info", color="module")) + 
      geom_line() + xlab(expression(theta)) + ylab("Information") + 
      theme_bw() + theme(legend.key=element_blank()) +
      guides(color=guide_legend("Modules")) +
      facet_grid(panel ~ stage)
  }
  
  g
}


#' @rdname mst
#' @description \code{mst_get_items} extracts items from the assembly results
#' @param panel_ix the panel index, an int vector
#' @param stage_ix the stage index, an int vector
#' @param module_ix the module index, an int vector
#' @param route_ix the route index, an integer
#' @export
mst_get_items <- function(x, panel_ix=NULL, stage_ix=NULL, module_ix=NULL, route_ix=NULL){
  if(class(x) != "mst") stop("not a 'mst' object: ", class(x))
  if(is.null(x$items)) stop('the mst has not been assembled yet.')
  items <- x$items
  if(!is.null(panel_ix))
    items <- subset(items, items$panel %in% panel_ix)
  if(!is.null(stage_ix))
    items <- subset(items, items$stage %in% stage_ix)
  if(!is.null(module_ix))
    items <- subset(items, items$module %in% module_ix)
  if(!is.null(route_ix))
    items <- subset(items, items$module %in% unlist(x$route[x$route$index == route_ix, 1:x$num_stage]))
  items
}