Nothing
#' Automated Test Assembly (ATA)
#' @name ata
#' @examples
#' ## generate a pool of 100 items
#' library(Rirt)
#' n_items <- 100
#' pool <- with(model_3pl_gendata(1, n_items), data.frame(id=1:n_items, a=a, b=b, c=c))
#' pool$content <- sample(1:3, n_items, replace=TRUE)
#' pool$time <- round(rlnorm(n_items, log(60), .2))
#' pool$group <- sort(sample(1:round(n_items/3), n_items, replace=TRUE))
#' pool <- list('3pl'=pool)
#'
#' ## ex. 1: four 10-item forms, maximize b parameter
#' x <- ata(pool, 4, test_len=10, max_use=1)
#' x <- ata_relative_objective(x, "b", "max")
#' x <- ata_solve(x, time_limit=2)
#' with(x$items$'3pl', aggregate(b, by=list(form=form), mean))
#' with(x$items$'3pl', table(form))
#'
#' \donttest{
#' ## ex. 2: four 10-item forms, minimize b parameter
#' x <- ata(pool, 4, test_len=10, max_use=1)
#' x <- ata_relative_objective(x, "b", "min", negative=TRUE)
#' x <- ata_solve(x, time_limit=5)
#' with(x$items$'3pl', aggregate(b, by=list(form=form), mean))
#' with(x$items$'3pl', table(form))
#'
#' ## ex. 3: two 10-item forms, mean(b)=0, sd(b)=1
#' ## content = (3, 3, 4), avg. time = 55--65 seconds
#' constr <- data.frame(name='content',level=1:3, min=c(3,3,4), max=c(3,3,4), stringsAsFactors=FALSE)
#' constr <- rbind(constr, c('time', NA, 55*10, 65*10))
#' x <- ata(pool, 2, test_len=10, max_use=1)
#' x <- ata_absolute_objective(x, pool$'3pl'$b, target=0*10)
#' x <- ata_absolute_objective(x, (pool$'3pl'$b-0)^2, target=1*10)
#' for(i in 1:nrow(constr))
#' x <- with(constr, ata_constraint(x, name[i], min[i], max[i], level=level[i]))
#' x <- ata_solve(x)
#' with(x$items$'3pl', aggregate(b, by=list(form=form), mean))
#' with(x$items$'3pl', aggregate(b, by=list(form=form), sd))
#' with(x$items$'3pl', aggregate(time, by=list(form=form), mean))
#' with(x$items$'3pl', aggregate(content, by=list(form=form), function(x) freq(x, 1:3)$freq))
#'
#' ## ex. 4: two 10-item forms, max TIF over (-1, 1), consider item sets
#' x <- ata(pool, 2, test_len=10, max_use=1, group="group")
#' x <- ata_relative_objective(x, seq(-1, 1, .5), 'max')
#' x <- ata_solve(x, time_limit=5)
#' plot(x)
#' }
NULL
#' @rdname ata
#' @description \code{ata} creates a basic ATA model
#' @param pool the item pool(s), a list of '3pl', 'gpcm', and 'grm' items
#' @param n_forms the number of forms to be assembled
#' @param test_len test length of each form
#' @param max_use maximum use of each item
#' @param ... options, e.g. group, common_items, overlap_items
#' @return \code{ata} returns a \code{ata} object
#' @details
#' The ATA model stores the definitions of a MIP model. When \code{ata_solve}
#' is called, a real MIP object is created from the definitions.
#' @export
ata <- function(pool, n_forms=1, test_len=NULL, max_use=NULL, ...){
pool <- ata_check_item_pool(pool)
opts <- list(...)
# scaling constants
if(is.null(opts$D))
opts$D <- 1.702
# combinations of unqiue and common items in each form
form_map <- ata_form_map(n_forms, opts)
n_forms <- max(form_map)
# group item sets
groups <- ata_item_set_groups(pool, opts)
n_items <- max(unlist(groups))
# the number of items and
n_lpvar <- n_items * n_forms + 2
# LP: x's (binary) + y1 (continuous) + y2 (continuous)
obj <- c(rep(0, n_lpvar - 2), 1, 1)
names(obj) <- c(paste("f", rep(1:n_forms, each=n_items), "v", rep(1:n_items, n_forms), sep=""), "y1", "y2")
# x's are binary and y is continuous
types <- c(rep("B", n_lpvar - 2), "C", "C")
# placehoders for fixing values
bounds <- list(idx=c(), lb=c(), ub=c())
# optimization direction: TRUE to maximize, FALSE to minimize
max <- TRUE
# TRUE if the obj.value is expected to be negative
negative <- FALSE
# constraints: coefficient matrix, directions, and right-hand-side values
mat <- matrix(nrow=0, ncol=n_lpvar, dimnames=list(NULL, names(obj)))
dir <- rhs <- NULL
x <- list(n_items=n_items, n_forms=n_forms, n_lpvar=n_lpvar, pool=pool, groups=groups,
form_map=form_map, obj=obj, mat=mat, dir=dir, rhs=rhs, types=types, bounds=bounds,
max=max, negative=negative, opts=opts)
class(x) <- "ata"
# add constraint: test length
if(!is.null(test_len)) {
if(length(test_len) > 2)
stop("Invalid test_len. Expect 1 or 2 elements")
if(length(test_len) == 1)
test_len <- rep(test_len, 2)
x <- ata_constraint(x, 1, min=test_len[1], max=test_len[2])
}
# add constraint: common_items length
if(!is.null(opts$common_items)){
if(length(opts$common_items) == 1)
opts$common_items <- rep(opts$common_items, 2)
x <- ata_constraint(x, 1, min=opts$common_items[1], max=opts$common_items[2], forms=n_forms, internal_index=TRUE)
}
# add constraint: overlap_items length
if(!is.null(opts$overlap_items)){
if(length(opts$overlap_items) == 1)
opts$overlap_items <- rep(opts$overlap_items, 2)
x <- ata_constraint(x, rep(1, length(unlist(groups))), min=opts$overlap_items[1], max=opts$overlap_items[2], forms=unique(as.vector(form_map[,-1])), internal_index=TRUE)
}
# add constraint: max_use
if(!is.null(max_use))
x <- ata_item_use(x, max=max_use)
x
}
#' @rdname ata
#' @description \code{ata_relative_objective} adds a relative objective to the model
#' @param coef the coefficients of the objective function
#' @param mode optimization direction: 'max' for maximization and 'min' for minimization
#' @param tol the tolerance paraemter
#' @param negative \code{TRUE} when the objective function is expected to be negative
#' @param forms forms where objectives are added. \code{NULL} for all forms
#' @param collapse \code{TRUE} to collapse into one objective function
#' @param internal_index \code{TRUE} to use internal form indices
#' @details
#' \code{ata_obj_relative}:
#' when mode='max', maximize (y-tol), subject to y <= sum(x) <= y+tol;
#' when mode='min', minimize (y+tol), subject to y-tol <= sum(x) <= y.
#' When \code{negative} is \code{TRUE}, y < 0, tol > 0.
#' \code{coef} can be a numeric vector that has the same length with the pool,
#' or a variable name in the pool, or a numeric vector of theta points.
#' When \code{tol} is \code{NULL}, it is optimized; when it's \code{FALSE}, ignored;
#' when it's a number, fixed; when it's a range, constrained with lower and upper bounds.
#' @export
ata_relative_objective <- function(x, coef, mode=c('max', 'min'), tol=NULL, negative=FALSE, forms=NULL, collapse=FALSE, internal_index=FALSE){
if(class(x) != "ata")
stop("Not an 'ata' object")
forms <- ata_get_form_index(x, forms, collapse, internal_index)
coef <- ata_get_obj_coef(x, coef, compensate=FALSE)
x$negative <- negative
# optimization direction
x$max <- switch(match.arg(mode), "max"=TRUE, "min"=FALSE)
if(x$max){
x$obj[x$n_lpvar - 1:0] <- c(1, -1)
} else {
x$obj[x$n_lpvar - 1:0] <- c(1, 1)
}
# tolerance parameter
if(!is.null(tol))
if(length(tol) == 2) {
x$bounds$idx <- c(x$bounds$idx, x$n_lpvar)
x$bounds$lb <- c(x$bounds$lb, tol[1])
x$bounds$ub <- c(x$bounds$ub, tol[2])
} else if(is.numeric(tol)){
x$bounds$idx <- c(x$bounds$idx, x$n_lpvar)
x$bounds$lb <- c(x$bounds$lb, tol)
x$bounds$ub <- c(x$bounds$ub, tol)
} else if(!tol) {
x$obj[x$n_lpvar] <- 0
}
# objective for each form
mat <- matrix(0, nrow=nrow(forms) * nrow(coef) * 2, ncol=x$n_lpvar)
dir <- rhs <- rep(NA, nrow(forms) * nrow(coef) * 2)
for(i in 1:nrow(forms)) {
f <- forms[i, ]
ind <- outer(1:x$n_items, (f - 1) * x$n_items, "+")
ind <- as.vector(ind)
for(j in 1:nrow(coef)) {
row <- (j - 1) * 2 + (i - 1) * nrow(coef) * 2
mat[row + 1, ind] <- rep(coef[j, ], length(f))
mat[row + 2, ind] <- rep(coef[j, ], length(f))
if(x$max){
mat[row + 1, (x$n_lpvar-1):x$n_lpvar] <- c(-1, 0)
mat[row + 2, (x$n_lpvar-1):x$n_lpvar] <- c(-1, -1)
} else {
mat[row + 1, (x$n_lpvar-1):x$n_lpvar] <- c(-1, 1)
mat[row + 2, (x$n_lpvar-1):x$n_lpvar] <- c(-1, 0)
}
dir[row + 1:2] <- c(">=", "<=")
rhs[row + 1:2] <- 0
}
}
ata_append(x, mat, dir, rhs)
}
#' @rdname ata
#' @description \code{ata_absolute_objective} adds an absolute objective to the model
#' @param target the target values of the objective function
#' @param equal_tol \code{TRUE} to force upward and downward tolerance to be equal
#' @param tol_up the range of upward tolerance
#' @param tol_down the range of downward tolerance
#' @details
#' \code{ata_obj_absolute} minimizes y0+y1 subject to t-y0 <= sum(x) <= t+y1.
#' @export
ata_absolute_objective <- function(x, coef, target, equal_tol=FALSE, tol_up=NULL, tol_down=NULL, forms=NULL, collapse=FALSE, internal_index=FALSE){
if(class(x) != "ata")
stop("Not an 'ata' object")
forms <- ata_get_form_index(x, forms, collapse, internal_index)
coef <- ata_get_obj_coef(x, coef, compensate=FALSE)
if(length(target) == 1)
target <- rep(target, nrow(coef))
if(length(target) != nrow(coef))
stop("Invalid target length.")
# optimization direction
x$max <- FALSE
x$obj[x$n_lpvar - 1:0] <- 1
# objective for each form
mat <- matrix(0, nrow=nrow(forms) * nrow(coef) * 2, ncol=x$n_lpvar)
dir <- rhs <- rep(NA, nrow(forms) * nrow(coef) * 2)
for(i in 1:nrow(forms)){
f <- forms[i,]
ind <- as.vector(outer(1:x$n_items, (f - 1) * x$n_items, "+"))
for(j in 1:nrow(coef)){
row <- (j - 1) * 2 + (i - 1) * nrow(coef) * 2
mat[row + 1, ind] <- rep(coef[j, ], length(f))
mat[row + 1, x$n_lpvar - 1] <- 1
mat[row + 2, ind] <- rep(coef[j, ], length(f))
mat[row + 2, x$n_lpvar - 0] <- -1
dir[row + 1:2] <- c(">=", "<=")
rhs[row + 1:2] <- target[j]
}
}
# tolerance parameters
if(equal_tol){
mat <- rbind(mat, c(rep(0, x$n_lpvar - 2), 1, -1))
dir <- c(dir, '=')
rhs <- c(rhs, 0)
}
if(!is.null(tol_down)){
if(length(tol_down) != 2)
stop('Invalid tol_down input. Expect 2 elements.')
x$bounds$idx <- c(x$bounds$idx, x$n_lpvar - 1)
x$bounds$lb <- c(x$bounds$lb, tol_down[1])
x$bounds$ub <- c(x$bounds$ub, tol_down[2])
}
if(!is.null(tol_up)){
if(length(tol_up) != 2)
stop('Invalid tol_up input. Expect 2 elements.')
x$bounds$idx <- c(x$bounds$idx, x$n_lpvar)
x$bounds$lb <- c(x$bounds$lb, tol_up[1])
x$bounds$ub <- c(x$bounds$ub, tol_up[2])
}
ata_append(x, mat, dir, rhs)
}
#' @rdname ata
#' @description \code{ata_constraint} adds a constraint to the model
#' @param level the level of a categorical variable to be constrained
#' @param min the lower bound of the constraint
#' @param max the upper bound of the constraint
#' @details
#' When \code{level} is \code{NA}, it is assumed that the constraint is on
#' a quantitative item property; otherwise, a categorical item property.
#' \code{coef} can be a variable name, a constant, or a numeric vector that has
#' the same size as the pool.
#' @importFrom stats aggregate
#' @export
ata_constraint <- function(x, coef, min=NA, max=NA, level=NULL, forms=NULL, collapse=FALSE, internal_index=FALSE){
if(class(x) != "ata")
stop("Not an 'ata' object")
if(is.na(min) && is.na(max))
return(x)
if(!is.na(min) && !is.na(max) && min > max)
stop("min is greater than max.")
forms <- ata_get_form_index(x, forms, collapse, internal_index)
if(length(coef) == x$n_items) { # a numeric vector at the group level
coef <- coef
} else if(length(coef) == length(unlist(x$groups))) { # a numeric vector at the item level
coef <- aggregate(coef, by=list(group=unlist(x$groups)), sum, na.rm=TRUE)[, -1]
} else if(is.numeric(coef) && length(coef) == 1) { # a constant, assuming it's at the item level
coef <- rep(coef, length(unlist(x$groups)))
coef <- aggregate(coef, by=list(group=unlist(x$groups)), sum, na.rm=TRUE)[, -1]
} else if(is.character(coef) && length(coef) == 1) {
coef <- Map(function(x, g) {
if(is.null(x))
return(numeric(0))
if(!coef %in% colnames(x))
return(rep(0, length(unique(g))))
if(is.na(level) || is.null(level)) {
x <- x[, coef]
} else {
x <- as.integer(x[, coef] == level)
}
aggregate(x, by=list(group=g), sum, na.rm=TRUE)[, -1]
}, x$pool, x$groups)
coef <- unlist(coef)
} else {
stop('Invalid coefficients')
}
n <- ifelse(!is.na(min) && !is.na(max) && min != max, nrow(forms) * 2, nrow(forms))
mat <- matrix(0, nrow=n, ncol=x$n_lpvar)
dir <- rhs <- rep(NA, n)
for(i in 1:nrow(forms)) {
f <- forms[i, ]
ind <- as.vector(outer(1:x$n_items, (f - 1) * x$n_items, "+"))
if(!is.na(min) && is.na(max)) {
mat[i, ind] <- coef
dir[i] <- ">="
rhs[i] <- min
} else if(is.na(min) && !is.na(max)) {
mat[i, ind] <- coef
dir[i] <- "<="
rhs[i] <- max
} else if(min == max) {
mat[i, ind] <- coef
dir[i] <- "="
rhs[i] <- min
} else {
mat[(i - 1) * 2 + 1, ind] <- coef
mat[(i - 1) * 2 + 2, ind] <- coef
dir[(i - 1) * 2 + 1:2] <- c(">=", "<=")
rhs[(i - 1) * 2 + 1:2] <- c(min, max)
}
}
ata_append(x, mat, dir, rhs)
}
#' @rdname ata
#' @description \code{ata_item_use} limits the minimum and maximum usage for items
#' @param items a vector of item indices, \code{NULL} for all items
#' @export
ata_item_use <- function(x, min=NA, max=NA, items=NULL){
if(class(x) != "ata")
stop("not an 'ata' object")
if(is.na(min) && is.na(max))
stop('min and max are both NA')
if(is.null(items))
items <- 1:x$n_items
if(any(!items %in% 1:x$n_items))
stop("Invalid items input")
nitems <- length(items)
n <- sum(!is.na(min), !is.na(max))
mat <- matrix(0, nrow=nitems * n, ncol=x$n_lpvar)
for(i in 1:length(items)) {
ind <- items[i] + (1:x$n_forms - 1) * x$n_items
mat[(i - 1) * n + 1:n, ind] <- 1
}
if(!is.na(min) && is.na(max)){
dir <- rep(">=", nitems)
rhs <- rep(min, nitems)
} else if(is.na(min) && !is.na(max)){
dir <- rep("<=", nitems)
rhs <- rep(max, nitems)
} else {
dir <- rep(c(">=", "<="), nitems)
rhs <- rep(c(min, max), nitems)
}
ata_append(x, mat, dir, rhs)
}
#' @rdname ata
#' @description \code{ata_item_enemy} adds an enemy-item constraint to the model
#' @export
ata_item_enemy <- function(x, items){
if(class(x) != "ata")
stop("not an 'ata' object")
if(!all(items %in% 1:x$n_items))
stop("Invalid item index")
mat <- matrix(0, nrow=nrow(x$form_map), ncol=x$n_lpvar)
for(i in 1:nrow(x$form_map)){
f <- x$form_map[i, ]
ind <- items + (f - 1) * x$n_items
mat[i, ind] <- 1
}
dir <- rep("<=", x$n_forms)
rhs <- rep(1, x$n_forms)
ata_append(x, mat, dir, rhs)
}
#' @rdname ata
#' @description \code{ata_item_fix} forces an item to be selected or not selected
#' @export
ata_item_fix <- function(x, items, min=NA, max=NA, forms){
if(class(x) != "ata")
stop("Not an 'ata' object")
if(any(!items %in% 1:x$n_items))
stop("Invalid items input.")
if(length(forms) > 1)
stop('Fix values in one form at each time.')
n <- length(items)
if(length(min) == 1)
min <- rep(min, n)
if(length(max) == 1)
max <- rep(max, n)
if(length(min) != n || length(max) != n)
stop("Invalid min or max length.")
x$bounds$idx <- c(x$bounds$idx, items+(forms-1)*x$n_items)
x$bounds$lb <- c(x$bounds$lb, min)
x$bounds$ub <- c(x$bounds$ub, max)
x
}
#' @rdname ata
#' @description \code{ata_solve} solves the MIP model
#' @param solver use 'lpsolve' for lp_solve 5.5 or 'glpk' for GLPK
#' @param return_format the format of the results: use \code{'form'} to organize results in a list of forms,
#' \code{'model'} to organize results in a list of models, use \code{'simple'} to organize results in
#' data.frame after removing item paraemters.
#' @param silent \code{TRUE} to mute solution information
#' @param time_limit the time limit in seconds passed along to solvers
#' @param message \code{TRUE} to print messages from solvers
#' @return \code{ata_solve} returns a solved \code{ata} object
#' @details
#' \code{ata_solve} takes control options in \code{...}.
#' For lpsolve, see \code{lpSolveAPI::lp.control.options}.
#' For glpk, see \code{glpkAPI::glpkConstants}\cr
#' Once the model is solved, additional data are added to the model.
#' \code{status} shows the status of the solution, \code{optimum}
#' the optimal value of the objective fucntion found in the solution,
#' \code{obj_vars} the values of two critical variables in the objective
#' function, \code{result} the assembly results in a binary matrix, and
#' \code{items} the assembled items
#' @export
ata_solve <- function(x, solver=c('lpsolve', 'glpk'), return_format=c('model', 'form', 'simple'), silent=FALSE, time_limit=10, message=FALSE, ...) {
if(class(x) != "ata")
stop("Not an 'ata' object")
rs <- switch(match.arg(solver, solver),
'lpsolve'=ata_lpsolve(x, time_limit, message, ...),
'glpk'=ata_glpk(x, time_limit, message, ...))
x$code <- rs$code
x$status <- rs$status
x$optimum <- rs$optimum
x$obj_vars <- rs$obj_vars
x$start_time <- rs$start_time
x$solve_time <- rs$solve_time
if(all(rs$result == 0)) {
x$result <- x$items <- NULL
if(!silent)
warning("No solution for the LP model.\n")
} else {
if(!silent)
cat(rs$status, ', optimum: ', round(rs$optimum, 3), ' (', paste(round(rs$obj_vars, 3), collapse=', '), ')\n', sep='')
x$result <- rs$result
x$items <- ata_extract_items(x)
return_format <- match.arg(return_format, return_format)
if(return_format %in% c('model', 'simple'))
x$items <- ata_results_to_model(x$items)
if(return_format == 'simple')
x$items <- ata_results_to_dataframe(x$items)
}
x
}
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.