R/CEA.R
In traets/idefix: Efficient Designs for Discrete Choice Experiments
Defines functions
SeqCEA
CEAcore_ucpp
CEA
Documented in
CEA
SeqCEA
#' Coordinate Exchange algorithm for MNL models.
#'
#' The algorithm improves an initial start design by considering changes on an
#' attribute-by-attribute basis. By doing this, it tries to minimize the
#' D(B)-error based on a multinomial logit model. This routine is repeated for
#' multiple starting designs.
#'
#' Each iteration will loop through all profiles from the initial design,
#' evaluating the change in D(B)-error for every level in each attribute.
#' The algorithm stops when an iteration occured without replacing a profile or
#' when \code{max.iter} is reached.
#'
#' By specifying a numeric vector in \code{par.draws}, the D-error will be
#' calculated and the design will be optimised locally. By specifying a matrix,
#' in which each row is a draw from a multivariate distribution, the DB-error
#' will be calculated, and the design will be optimised globally. Whenever there
#' are alternative specific constants, \code{par.draws} should be a list
#' containing two matrices: The first matrix containing the parameter draws for
#' the alternative specific constant parameters. The second matrix containing
#' the draws for the rest of the parameters.
#'
#' The DB-error is calculated by taking the mean over D-errors. It could be that
#' for some draws the design results in an infinite D-error. The percentage of
#' draws for which this was true for the final design can be found in the output
#' \code{inf.error}.
#'
#' Alternative specific constants can be specified in \code{alt.cte}. The length
#' of this binary vector should equal \code{n.alts}, were \code{0} indicates the
#' absence of an alternative specific constant and \code{1} the opposite.
#'
#' \code{start.des} is a list with one or several matrices corresponding to
#' initial start design(s). In each matrix each row is a profile. The number of rows
#' equals \code{n.sets * n.alts}, and the
#' number of columns equals the number of columns of the design matrix + the
#' number of non-zero elements in \code{alt.cte}. Consider that for a
#' categorical attribute with *p* levels, there are *p - 1* columns in the design
#' matrix, whereas for a continuous attribute there is only one column. If
#' \code{start.des = NULL}, \code{n.start} random initial designs will be
#' generated. If start designs are provided, \code{n.start} is ignored.
#'
#' If \code{no.choice} is \code{TRUE}, in each choice set an alternative with
#' one alternative specific constant is added. The return value of the
#' D(B)-error is however based on the design without the no choice option.
#'
#' When \code{parallel} is \code{TRUE}, \code{\link[parallel]{detectCores}} will
#' be used to decide upon the number of available cores. That number minus 1
#' cores will be used to search for efficient designs. The computation time will
#' decrease significantly when \code{parallel = TRUE}.
#'
#' @param lvls A numeric vector which contains for each attribute the number
#' of levels.
#' @param coding Type of coding that needs to be used for each attribute.
#' @param c.lvls A list containing numeric vectors with the attribute levels for
#' each continuous attribute. The default is \code{NULL}.
#' @param n.sets Numeric value indicating the number of choice sets.
#' @param n.alts Numeric value indicating the number of alternatives per choice
#' set.
#' @param par.draws A matrix or a list, depending on \code{alt.cte}.
#' @param alt.cte A binary vector indicating for each alternative whether an
#' alternative specific constant is desired. The default is \code{NULL}.
#' @param no.choice A logical value indicating whether a no choice alternative
#' should be added to each choice set. The default is \code{FALSE}.
#' @param start.des A list containing one or more matrices corresponding to
#' initial start design(s). The default is \code{NULL}.
#' @param parallel Logical value indicating whether computations should be done
#' over multiple cores. The default is \code{TRUE}.
#' @param max.iter A numeric value indicating the maximum number allowed
#' iterations. The default is \code{Inf}.
#' @param n.start A numeric value indicating the number of random start designs
#' to use. The default is 12.
#' @param best A logical value indicating whether only the best design should be
#' returned. The default is \code{TRUE}.
#' @return
#' If \code{best = TRUE} the design with the lowest D(B)-error is returned.
#' If \code{best = FALSE}, the results of all (provided) start designs are
#' returned. \item{design}{A
#' numeric matrix wich contains an efficient design.} \item{error}{Numeric
#' value indicating the D(B)-error of the design.} \item{inf.error}{Numeric
#' value indicating the percentage of draws for which the D-error was
#' \code{Inf}.} \item{probs}{Numeric matrix containing the probabilities of
#' each alternative in each choice set. If a sample matrix was provided in
#' \code{par.draws}, this is the average over all draws.}
#' @examples
#' \donttest{
#' # DB-efficient designs
#' # 3 Attributes, all dummy coded. 1 alternative specific constant = 7 parameters
#' mu <- c(1.2, 0.8, 0.2, -0.3, -1.2, 1.6, 2.2) # Prior parameter vector
#' v <- diag(length(mu)) # Prior variance.
#' set.seed(123)
#' pd <- MASS::mvrnorm(n = 10, mu = mu, Sigma = v) # 10 draws.
#' p.d <- list(matrix(pd[,1], ncol = 1), pd[,2:7])
#' CEA(lvls = c(3, 3, 3), coding = c("D", "D", "D"), par.draws = p.d,
#' n.alts = 2, n.sets = 8, parallel = FALSE, alt.cte = c(0, 1))
#'
#' # DB-efficient design with categorical and continuous factors
#' # 2 categorical attributes with 4 and 2 levels (effect coded) and 1
#' # continuous attribute (= 5 parameters)
#' mu <- c(0.5, 0.8, 0.2, 0.4, 0.3)
#' v <- diag(length(mu)) # Prior variance.
#' set.seed(123)
#' pd <- MASS::mvrnorm(n = 3, mu = mu, Sigma = v) # 10 draws.
#' CEA(lvls = c(4, 2, 3), coding = c("E", "E", "C"), par.draws = pd,
#' c.lvls = list(c(2, 4, 6)), n.alts = 2, n.sets = 6, parallel = FALSE)
#'
#' # DB-efficient design with start design provided.
#' # 3 Attributes with 3 levels, all dummy coded (= 6 parameters).
#' mu <- c(0.8, 0.2, -0.3, -0.2, 0.7, 0.4)
#' v <- diag(length(mu)) # Prior variance.
#' sd <- list(example_design)
#' set.seed(123)
#' ps <- MASS::mvrnorm(n = 10, mu = mu, Sigma = v) # 10 draws.
#' CEA(lvls = c(3, 3, 3), coding = c("D", "D", "D"), par.draws = ps,
#' n.alts = 2, n.sets = 8, parallel = FALSE, start.des = sd)
#'}
#' @importFrom Rdpack reprompt
#' @export
CEA <- function(lvls, coding, c.lvls = NULL, n.sets, n.alts, par.draws,
alt.cte = NULL, no.choice = FALSE, start.des = NULL,
parallel = TRUE, max.iter = Inf, n.start = 12, best = TRUE) {
### Error handling for creating initial random design
# Error lvls vector. At least 2 attributes
if (length(lvls) < 2 || (!(is.numeric(lvls)))) {
stop("lvls argument is incorrect.")
}
# Error correct coding types.
codings.types <- c("E", "D", "C")
if (!all(coding %in% codings.types) || (length(coding) != length(lvls))) {
stop("coding argument is incorrect.")
}
# Continuous attributes.
contins <- which(coding == "C")
n.contins <- length(contins)
# error continuous levels (Not a list)
if (!is.null(c.lvls) && !is.list(c.lvls)) {
stop('c.lvls should be a list.')
}
# Error continuous specified and NULL.
if (n.contins > 0 && is.null(c.lvls)) {
stop("when 'coding' contains C, 'c.lvls' should be specified")
}
# Error continuous levels specification.
if (!is.null(c.lvls)) {
if (length(c.lvls) != n.contins) {
stop("length of 'c.lvls' does not match number of specified continuous attributes in 'coding'")
}
# Error c.lvls same number of levels.
if (!isTRUE(all.equal(lvls[contins], lengths(c.lvls)))) {
stop("the number of levels provided in 'c.lvls' does not match the expected based on 'lvls'")
}
}
### Error handling for design specifications
# If no alternative constant is given, create the variable as a vector of 0s
if (is.null(alt.cte)) {
alt.cte <- rep(0L, n.alts)
}
#init
n.cte <- length(which(alt.cte == 1))
# If only one draw is given, transform it to a matrix
if (!is.list(par.draws)) {
if (is.vector(par.draws)) {
par.draws <- matrix(par.draws, nrow = 1)
}
}
# Alternative constant errors
if (length(alt.cte) != n.alts) {
stop("'n.alts' does not match the 'alt.cte' vector")
}
if (!all(alt.cte %in% c(0, 1))) {
stop("'alt.cte' should only contain zero or ones.")
}
# No choice errors
if (!is.logical(no.choice)) {
stop("'no.choice' should be TRUE or FALSE")
}
if (no.choice) {
if (!isTRUE(all.equal(alt.cte[n.alts], 1))) {
stop("if 'no.choice' is TRUE, the last alternative constant should equal 1.")
}
ncsek <- seq(n.alts, (n.sets * n.alts), n.alts) # Rows in the design matrix
# to be assigned as zeros in all attributes
} else {
ncsek <- NULL
}
# Handling par.draws with alternative specific contstants.
# This conditional is when there is only one alternative constant
if (isTRUE(all.equal(n.cte, 1))) {
if (!(is.list(par.draws))) {
stop("par.draws should be a list")
}
if (!isTRUE(all.equal(length(par.draws), 2))) {
stop("'par.draws' should contain two components")
}
# If only draws for the constant are given in a vector, transform it to a
# matrix
if (is.vector(par.draws[[1]])) {
par.draws[[1]] <- matrix(par.draws[[1]], ncol = 1)
}
if (!(all(unlist(lapply(par.draws, is.matrix))))) {
stop("'par.draws' should contain two matrices")
}
dims <- as.data.frame(lapply(par.draws, dim))
if (!isTRUE(all.equal(dims[1, 1], dims[1, 2]))) {
stop("the number of rows in the components of 'par.draws' should be equal")
}
par.draws <- do.call("cbind", par.draws) # Transform draws to a matrix
} else {
if (n.cte > 1.2) {
if (!(is.list(par.draws))) {
stop("par.draws should be a list")
}
if (!isTRUE(all.equal(length(par.draws), 2))) {
stop("'par.draws' should contain two components")
}
if (!(all(unlist(lapply(par.draws, is.matrix))))) {
stop("'par.draws' should contain two matrices")
}
if (!isTRUE(all.equal(ncol(par.draws[[1]]), n.cte))) {
stop("the first component of 'par.draws' should contain the same number of columns as there are non zero elements in 'alt.cte'")
}
dims <- as.data.frame(lapply(par.draws, dim))
if (!isTRUE(all.equal(dims[1, 1], dims[1, 2]))) {
stop("the number of rows in the components of 'par.draws' should be equal")
}
par.draws <- do.call("cbind", par.draws) # Transform draws to a matrix
}
}
# Error identifying model.
if (n.sets < ncol(par.draws)) {
stop("Model is unidentified. Increase the number of choice sets or decrease parameters to estimate.")
}
### Create all levels for each attribute
# Change into correct coding.
coding <- dplyr::recode(coding, D = "contr.treatment", E = "contr.sum")
# Create all combinations of attribute levels.
levels.list <- lapply(X = as.list(lvls), function(x) (1:x))
# Replace continuous.
levels.list[contins] <- c.lvls
# Transform to matrix
levels.list[contins] <- lapply(X = levels.list[contins], as.matrix)
# Transform categorical attributes
categ <- which(coding %in% c("contr.treatment", "contr.sum"))
n.categ <- length(categ)
if (n.categ > 0) {
levels.list[categ] <- lapply(X = levels.list[categ], factor)
# Apply coding
if (n.contins > 0) {
for (i in 1:length(lvls)) {
if (!(i %in% contins)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
}else {
for (i in 1:length(lvls)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
# Compute all possible values for each categorical attribute
levels.list[categ] <- lapply(X = levels.list[categ], stats::contrasts)
}
# Set colnames for the design matrix
c.nam = list()
for (i in 1:length(lvls)) {
if (coding[i] == "contr.treatment") {
c.nam[[i]] <- paste("Var", i, 2:lvls[i], sep = "")
} else {
if (coding[i] == "contr.sum") {
c.nam[[i]] <- paste("Var", i, 1:(lvls[i] - 1), sep = "")
} else {
c.nam[[i]] <- paste("Var", i, sep = "")
}
}
}
c.nam = unlist(c.nam)
# Create alternative specific design.
cte.des <- Altspec(alt.cte = alt.cte, n.sets = n.sets)
# Count the number of colums of the design matrix without constants
ncol.des.noconst <- sum(unlist(lapply(levels.list,ncol)))
#ncol.des.noconst <- sum(unlist(lapply(levels.list,function(x){
# if (is.matrix(x)) { ncol(x) }
# else{if (is.numeric(x)) { 1 }}
#})))
if (!identical(as.integer(ncol(par.draws)),
as.integer(n.cte + ncol.des.noconst))) {
stop("The sum of the number of columns in the components of 'par.draws' should equal the number of columns of design matrix (including alternative specific constants)")
}
### Random initial design.
if (is.null(start.des)) {
#create start designs
nr.starts <- n.start
start.des <- vector(mode = 'list', length = nr.starts)
okstart <- FALSE
while (okstart == FALSE) {
for (i in 1:nr.starts) {
# r is to know which levels to take in each attribute
r <- NULL
start <- NULL
for (j in 1:length(lvls)) {
r <- round(stats::runif((n.sets * n.alts), 1, lvls[j]))
start <- cbind(start, levels.list[[j]][r,])
}
colnames(start) <- c.nam
start.des[[i]] <- cbind(cte.des, start)
if (no.choice) {
start.des[[i]][ncsek, (ncol(cte.des) + 1):(ncol(cte.des) + ncol.des.noconst)] <- c(rep(0, ncol.des.noconst))
}
}
# Compute D-optimality for each design and each draw
d.start <- lapply(start.des, StartDB, par.draws, n.alts)
# If the DB-optimality of any starting desing is finite, continue
if (any(is.finite(unlist(lapply(d.start, mean, na.rm = TRUE))))) {
okstart <- TRUE
}
}
} else {
if (!is.list(start.des)) {
stop("'start.des' should be a list")
}
if (!(all(unlist(lapply(start.des, is.matrix))))) {
stop("'start.des' should contain matrices as components")
}
# Save the dimension of each starting design
dimstart <- as.matrix(lapply(start.des, dim))
# Save the number of random starts given
nr.starts <- length(dimstart)
if (nr.starts > 1.5) {
if (!isTRUE(all.equal(length(unique(unlist(dimstart))), 2))) {
stop("start designs have different dimensions")
}
}
if (!isTRUE(all.equal(n.alts * n.sets, unique(unlist(dimstart))[1]))) {
stop("number of rows of start design(s) does not match with 'n.alts' * 'n.sets'")
}
if (!isTRUE(all.equal(as.integer(n.cte + ncol.des.noconst),
unique(unlist(dimstart))[2]))) {
stop("number of columns of start design(s) does not match with the number of columns in the design matrix")
}
# Compute D-optimality for each design and each draw
d.start <- lapply(start.des, StartDB, par.draws, n.alts)
if (!any(is.finite(unlist(lapply(d.start, mean, na.rm = TRUE))))) {
stop("One or more of the provided start designs resulted in an unvalid db-error.")
}
}
### Improving the initial design
if (parallel) {
no_cores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(no_cores)
parallel::clusterExport(cl, c("n.sets", "par.draws", "n.alts", "n.cte",
"alt.cte", "no.choice", "max.iter","ncsek"),
envir = environment())
deslist <- parallel::parLapply(cl, start.des, CEAcore_ucpp, par.draws,
levels.list, n.alts, n.sets, n.cte, alt.cte,
no.choice, max.iter, ncsek)
parallel::stopCluster(cl)
} else {
deslist <- lapply(start.des, CEAcore_ucpp, par.draws, levels.list,
n.alts, n.sets, n.cte, alt.cte, no.choice,
max.iter = max.iter, ncsek)
}
bestdes <- deslist[[which.min(unlist(lapply(deslist, function(x) (x$error))))]]
ifelse(best, return(bestdes), return(deslist))
}
# Core of Coordinate Exchange
CEAcore_ucpp <- function(des, par.draws, levels.list, n.alts, n.sets, n.cte,
alt.cte, no.choice, max.iter, ncsek) {
converge <- FALSE # Boolean for convergence
change <- FALSE # Boolean for a change in the attribute
it <- 1 # Indicator of number of iterations
n.samples <- nrow(par.draws) # Number of samples from distribution of betas
n.par <- ncol(des) # Number of parameters
# Beginning of algorithm
while (!converge & it <= max.iter) {
# Compute DB-optimality of initial design
db.start <- mean(apply(par.draws, 1, Derr_ucpp, des = des,
n.alts = n.alts), na.rm = TRUE)
it <- it + 1 # Increase iteration
# Save design before iteration.
iter.des <- des
# For every row in the design.
sek <- 1:nrow(des)
if (no.choice) {
sek <- sek[-ncsek] # If no.choice is given, then it is not improved
}
# Loop for each row of the initial design
for (i in sek) {
# Loop for each attribute in the row
for (j in 1:length(levels.list)) {
#%%%%%%%
# This can be improved by removing mods object and just replace the
# initial design with the new attributes to compute the db.error and
# then if one of those is better than the initial db error. Replace
# finally the initial design (as it is in modfed)
# Initialize mods object to track changes of the design in each attribute
mods <- vector("list", nrow(levels.list[[j]]))
mods <- lapply(mods,function(x){des})
# Initialize db object for db error of each level
db <- numeric(nrow(levels.list[[j]]))
# Indicator of columns modified
ncol.ok <- sum(unlist(lapply(levels.list[(1:j - 1)],ncol)))
# Indicator of columns to change
ncol.mod <- ifelse( j == 1, (n.cte + j),
(n.cte + 1 + ncol.ok))
# ncol.left tracks which columns are left to be improved in the row
ncol.left <- sum(unlist(lapply(levels.list[-(1:j)],ncol)))
# Loop for each level in the attribute
for (k in 1:nrow(levels.list[[j]])) {
# Update design with new attribute
#mods[[k]][i, (n.cte + j):(n.par - ncol.left) ] <- levels.list[[j]][k,]
mods[[k]][i, ncol.mod:(n.par - ncol.left) ] <- levels.list[[j]][k,]
# Calculate D-optimality for each draw
d.errors <- apply(par.draws, 1, Derr_ucpp, des = mods[[k]],
n.alts = n.alts)
# Compute DB-optimality
db[k] <- mean(d.errors, na.rm = TRUE)
}
pr <- which.min(db)
db <- min(db)
# Change if lower db error.
if (!is.na(db) && !is.na(db.start)) {
if (db < db.start) {
des <- mods[[pr]]
db.start <- db
}
}
} # End loop for each attribute in the row
} # End loop for each row of the initial design
converge <- isTRUE(all.equal(des, iter.des)) # Convergence if no profile is swapped this iteration.
} # End while (after convergence)
#return(it, des)
# calculate percentage NA values.
na.percentage <- 0
d.errors <- apply(par.draws, 1, Derr_ucpp, des = des, n.alts = n.alts)
if (any(is.na(d.errors))) {
na.percentage <- scales::percent(sum(is.na(d.errors)) / n.samples)
}
# Utility balance.
# ub <- apply(par.draws, 1, Utbal, des = des, n.alts = n.alts)
# ub2 <- apply(par.draws, 1, InfoDes2, des = des, n.alts = n.alts, utbal = TRUE)
# ub_ucpp <- apply(par.draws, 1, InfoDes_cpp, des = des, n_alts = n.alts,
# utbal = TRUE)
# Utility balance using c++ function
ub <- apply(par.draws, 1, InfoDes_cpp, des = des, n_alts = n.alts,
utbal = TRUE)
pmat <- matrix(rowMeans(ub), ncol = n.alts, byrow = TRUE)
rownames(pmat) <- paste("set", 1:n.sets, sep = "")
colnames(pmat) <- paste(paste("Pr(", paste("alt", 1:n.alts, sep = ""),
sep = ""), ")", sep = "")
if (no.choice) {
colnames(pmat)[n.alts] <- "Pr(no choice)"
}
# Rownames design.
des.names <- Rcnames(n.sets = n.sets, n.alts = n.alts, alt.cte = alt.cte,
no.choice = no.choice)
rownames(des) <- des.names[[1]]
# Colnames alternative specific constants.
if (n.cte != 0 && !is.null(colnames(des))) {
colnames(des)[1:n.cte] <- des.names[[2]]
}
# Return design, D(B)error, percentage NA's, utility balance.
return(list("design" = des, "error" = db.start, "inf.error" = na.percentage,
"probs" = pmat))
}
#' Sequential Coordinate Exchange algorithm for MNL model.
#'
#' Selects the choice set that minimizes the DB-error when added to an initial
#' design, given (updated) parameter values.
#'
#' This algorithm is ideally used in an adaptive context. The algorithm will
#' select the next DB-efficient choice set given parameter values and possible
#' previously generated choice sets. In an adaptive context these parameter
#' values are updated after each observed response.
#'
#' Previously generated choice sets, which together form an initial design, can
#' be provided in \code{des}. When no design is provided, the algorithm will
#' select the most efficient choice set based on the fisher information of the
#' prior covariance matrix \code{prior.covar}.
#'
#' If \code{alt.cte = NULL}, \code{par.draws} should be a matrix in which each
#' row is a sample from the multivariate parameter distribution. In case that
#' \code{alt.cte} is not \code{NULL}, a list containing two matrices should be
#' provided to \code{par.draws}. The first matrix containing the parameter draws
#' for the alternative specific parameters. The second matrix containing the
#' draws for the rest of the parameters.
#'
#' The list of potential choice sets is created by selecting randomly a level for
#' each attribute in an alternative/profile. \code{n.cs} controls the number of
#' potential choice sets to consider. The default is \code{
#' NULL}, which means that the number of possible choice sets is the product of
#' attribute levels considered in the experiment. For instance, an experiment
#' with 3 attribute and 3 levels each will consider 3^3 = 27 possible choice sets.
#'
#' The \code{weights} argument can be used when the \code{par.draws} have
#' weights. This is for example the case when parameter values are updated using
#' \code{\link{ImpsampMNL}}.
#'
#' When \code{parallel} is \code{TRUE}, \code{\link[parallel]{detectCores}} will
#' be used to decide upon the number of available cores. That number minus 1
#' cores will be used to search for the optimal choice set. For small problems
#' (6 parameters), \code{parallel = TRUE} can be slower. For larger problems the
#' computation time will decrease significantly.
#'
#' *Note:* this function is faster than \code{\link[idefix]{SeqMOD}}, but
#' the output is not as stable. This happens because this function
#' makes a random search to get the choice set, whereas
#' \code{\link[idefix]{SeqMOD}} makes an exhaustive search.
#' @inheritParams CEA
#' @param par.draws A matrix or a list, depending on \code{alt.cte}.
#' @param des A design matrix in which each row is a profile. If alternative
#' specific constants are present, those should be included as the first
#' column(s) of the design. Can be generated with \code{\link{Modfed}} or
#' \code{\link{CEA}}
#' @param n.cs An integer indicating the number of possible random choice sets to
#' consider in the search for the next best choice set possible. The default is
#' \code{NULL}.
#' @param prior.covar Covariance matrix of the prior distribution.
#' @param weights A vector containing the weights of the draws. Default is
#' \code{NULL}. See also \code{\link{ImpsampMNL}}.
#' @param parallel Logical value indicating whether computations should be done
#' over multiple cores.
#' @param no.choice An integer indicating the no choice alternative. The default
#' is \code{NULL}.
#' @param reduce Logical value indicating whether the candidate set should be
#' reduced or not.
#' @return \item{set}{A matrix representing a DB efficient choice set.}
#' \item{error}{A numeric value indicating the DB-error of the whole
#' design.}
#' @importFrom Rdpack reprompt
#' @references \insertRef{idefix}{idefix}
#' @references \insertRef{ju}{idefix}
#' @references \insertRef{cea}{idefix}
#' @references \insertRef{cea_discrete}{idefix}
#' @examples
#' # DB efficient choice set, given a design and parameter draws.
#' # 3 attributes with 3 levels each
#' m <- c(0.3, 0.2, -0.3, -0.2, 1.1, 2.4) # mean (total = 6 parameters).
#' pc <- diag(length(m)) # covariance matrix
#' set.seed(123)
#' sample <- MASS::mvrnorm(n = 10, mu = m, Sigma = pc)
#' # Initial design.
#' des <- example_design
#' # Efficient choice set to add.
#' SeqCEA(des = des, lvls = c(3, 3, 3), coding = c("D", "D", "D"), n.alts = 2,
#' par.draws = sample, prior.covar = pc, parallel = FALSE)
#'
#' # DB efficient choice set, given parameter draws.
#' # with alternative specific constants
#' des <- example_design2
#' ac <- c(1, 1, 0) # Alternative specific constants.
#' m <- c(0.3, 0.2, -0.3, -0.2, 1.1, 2.4, 1.8, 1.2) # mean
#' pc <- diag(length(m)) # covariance matrix
#' pos <- MASS::mvrnorm(n = 10, mu = m, Sigma = pc)
#' sample <- list(pos[ , 1:2], pos[ , 3:8])
#' # Efficient choice set.
#' SeqCEA(des = des, lvls = c(3, 3, 3), coding = c("D", "D", "D"), n.alts = 3,
#' par.draws = sample, alt.cte = ac, prior.covar = pc, parallel = FALSE)
#' @export
SeqCEA <- function(des = NULL, lvls, coding, c.lvls = NULL, n.alts, par.draws,
prior.covar, alt.cte = NULL, no.choice = NULL,
weights = NULL, parallel = TRUE, reduce = TRUE, n.cs = NULL) {
# Error handling initial design
if (is.null(des)) {
n.sets <- 1L
} else {
if (!is.matrix(des)) {
stop("'des' should be a matrix or NULL")
}
if (!isTRUE(nrow(des) %% n.alts == 0)) {
stop("'n.alts' does not seem to match with the number of rows in 'des'")
}
n.sets <- nrow(des) / n.alts
}
### Error handling for design specifications
# No choice errors
if (!is.null(no.choice)) {
if (!is.wholenumber(no.choice)) {
stop("'no.choice' should be an integer or NULL")
}
if (any(isTRUE(no.choice > (n.alts + 0.2)), isTRUE(no.choice < 0.2))) {
stop("'no.choice' does not indicate one of the alternatives")
}
if (is.null(alt.cte)) {
stop("if there is a no choice alternative, 'alt.cte' should be specified")
}
if (!isTRUE(all.equal(alt.cte[no.choice], 1))) {
stop("the no choice alternative should correspond with a 1 in 'alt.cte'")
}
}
# Alternative constant errors
if (!is.null(alt.cte)) {
if (length(alt.cte) != n.alts) {
stop("'n.alts' does not match the 'alt.cte' vector")
}
if (!all(alt.cte %in% c(0, 1))) {
stop("'alt.cte' should only contain zero or ones.")
}
# alternative specific constants
n.cte <- length(which(alt.cte == 1L))
if (isTRUE(all.equal(n.cte, 0L))) {
alt.cte <- NULL
cte.des <- NULL
}
# Handling errors when there are alternative constants
if (!is.null(alt.cte)) {
if (n.cte > 0.2) {
if (!is.list(par.draws)) {
stop("'par.draws' should be a list when 'alt.cte' is not NULL")
}
if (!isTRUE(all.equal(length(par.draws), 2))) {
stop("'par.draws' should contain two components")
}
# If there is only one specific constant and is a vector, then it is
# transformed to a matrix
if (isTRUE(all.equal(n.cte, 1))) {
if (is.vector(par.draws[[1]])) {
par.draws[[1]] <- matrix(par.draws[[1]], ncol = 1)
}
}
if (!(all(unlist(lapply(par.draws, is.matrix))))) {
stop("'par.draws' should contain two matrices")
}
if (!isTRUE(all.equal(ncol(par.draws[[1]]), n.cte))) {
stop("the first component of 'par.draws' should contain the same number
of columns as there are non zero elements in 'alt.cte'")
}
dims <- as.data.frame(lapply(par.draws, dim))
if (!isTRUE(all.equal(dims[1, 1], dims[1, 2]))) {
stop("the number of rows in the components of 'par.draws' should be equal")
}
par.draws <- do.call("cbind", par.draws) # Transform par.draws to a matrix
}
# Create alternative specific design.
cte.des <- Altspec(alt.cte = alt.cte, n.sets = n.sets)
cte.set <- matrix(cte.des[1:n.alts, ], ncol = n.cte, byrow = FALSE)
}
} else {
cte.des <- NULL
n.cte <- 0
# if no alternative constants
if (!is.matrix(par.draws)) {
stop("'par.draws'should be a matrix when 'alt.cte' = NULL")
}
}
# Weights errors
n.par <- ncol(par.draws)
if (!is.null(weights)) {
if (!isTRUE(all.equal(length(weights), nrow(par.draws)))) {
stop("length of 'weights' does not match number total number of rows in 'par.draws'")
}
} else {
weights <- rep(1L, nrow(par.draws))
}
### Create all levels for each attribute
# Continuous attributes.
contins <- which(coding == "C")
n.contins <- length(contins)
# Change into correct coding.
coding <- dplyr::recode(coding, D = "contr.treatment", E = "contr.sum")
# Create all combinations of attribute levels.
levels.list <- lapply(X = as.list(lvls), function(x) (1:x))
# Replace continuous.
levels.list[contins] <- c.lvls
# Transform to matrix
levels.list[contins] <- lapply(X = levels.list[contins], as.matrix)
# Transform categorical attributes
categ <- which(coding %in% c("contr.treatment", "contr.sum"))
n.categ <- length(categ)
if (n.categ > 0) {
levels.list[categ] <- lapply(X = levels.list[categ], factor)
# Apply coding
if (n.contins > 0) {
for (i in 1:length(lvls)) {
if (!(i %in% contins)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
}else {
for (i in 1:length(lvls)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
# Compute all possible values for each categorical attribute
levels.list[categ] <- lapply(X = levels.list[categ], stats::contrasts)
}
# Set colnames for the design matrix
c.nam = list()
for (i in 1:length(lvls)) {
if (coding[i] == "contr.treatment") {
c.nam[[i]] <- paste("Var", i, 2:lvls[i], sep = "")
} else {
if (coding[i] == "contr.sum") {
c.nam[[i]] <- paste("Var", i, 1:(lvls[i] - 1), sep = "")
} else {
c.nam[[i]] <- paste("Var", i, sep = "")
}
}
}
c.nam = unlist(c.nam)
# Count the number of colums of the design matrix without constants
ncol.des.noconst <- sum(unlist(lapply(levels.list,ncol)))
# if (!identical(as.integer(n.par), as.integer(n.cte + ncol.des.noconst))) {
# stop("The sum of the number of columns in the components of 'par.draws' should equal the number of columns of design matrix (including alternative specific constants)")
# }
if (!identical(as.integer(ncol(prior.covar)),
as.integer(n.cte + ncol.des.noconst))) {
stop("number of columns of 'prior.covar' does not equal the number of columns of design matrix (including alternative specific constants)")
}
## When a design is supplied
if (!is.null(des)) {
# Error par.draws
if (!isTRUE(all.equal(ncol(des), n.par))) {
stop("number of columns in 'par.draws' does not match the number of columns in 'des'")
}
# Error dimension of initial design and design matrix
if (!identical(as.integer(ncol(des)), as.integer(n.cte + ncol.des.noconst))) {
stop("number of columns in 'des' does not match the number of columns of design matrix (including alternative specific constants)")
}
# Starting and initializing values.
i.cov <- solve(prior.covar)
d.start <- apply(par.draws, 1, DerrC_ucpp, des = des, n.alts = n.alts,
i.cov = i.cov)
db.start <- mean(d.start, na.rm = TRUE)
full.comb <- Newsets_ucpp(levels.list = levels.list, n.alts = n.alts,
no.choice = no.choice, reduce = reduce,
c.names = c.nam, n.cs = n.cs)
# Adding alternative constants
if (!is.null(cte.des)) {
full.comb <- lapply(full.comb, function(x) cbind(cte.set, x))
}
# Adding these new choice sets to the initial design
full.des <- lapply(full.comb, function(x) rbind(des, x))
# Select the next best choice set
if (parallel) {
no_cores <- parallel::detectCores() - 1L
cl <- parallel::makeCluster(no_cores)
db.errors <- parallel::parLapply(cl, full.des, DBerrS.P_ucpp, par.draws,
n.alts, i.cov, weights)
parallel::stopCluster(cl)
} else {
# For each potential set, select best.
db.errors <- lapply(full.des, DBerrS.P_ucpp, par.draws, n.alts, i.cov,
weights)
}
dbs <- unlist(db.errors, use.names = FALSE)
set <- full.comb[[which.min(dbs)]]
row.names(set) <- NULL
colnames(set) <- colnames(des)
db <- min(dbs)
#return best set and db error design.
return(list("set" = set, "error" = db))
} # End if when a design is supplied
else {
# Starting and initializing values.
i.cov <- solve(prior.covar)
full.comb <- Newsets_ucpp(levels.list = levels.list, n.alts = n.alts,
no.choice = no.choice, reduce = reduce,
c.names = c.nam, n.cs = n.cs)
# Adding alternative constants
if (!is.null(cte.des)) {
full.comb <- lapply(full.comb, function(x) cbind(cte.set, x))
}
# Select the next best choice set
if (parallel) {
no_cores <- parallel::detectCores() - 1L
cl <- parallel::makeCluster(no_cores)
db.errors <- parallel::parLapply(cl, full.comb, DBerrS.P_ucpp, par.draws,
n.alts, i.cov, weights)
parallel::stopCluster(cl)
} else {
# For each potential set, select best.
db.errors <- lapply(full.comb, DBerrS.P_ucpp, par.draws, n.alts, i.cov,
weights)
}
dbs <- unlist(db.errors, use.names = FALSE)
set <- full.comb[[which.min(dbs)]]
row.names(set) <- NULL
if (n.cte > 0) {
colnames(set) <- c(paste("alt",1:n.cte,".cte", sep = ""), c.nam)
}
db <- min(dbs)
#return best set and db error design.
return(list("set" = set, "error" = db))
}
}
traets/idefix documentation built on April 7, 2022, 1:56 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions SeqCEA CEAcore_ucpp CEA
Documented in CEA SeqCEA
#' Coordinate Exchange algorithm for MNL models.
#'
#' The algorithm improves an initial start design by considering changes on an
#' attribute-by-attribute basis. By doing this, it tries to minimize the
#' D(B)-error based on a multinomial logit model. This routine is repeated for
#' multiple starting designs.
#'
#' Each iteration will loop through all profiles from the initial design,
#' evaluating the change in D(B)-error for every level in each attribute.
#' The algorithm stops when an iteration occured without replacing a profile or
#' when \code{max.iter} is reached.
#'
#' By specifying a numeric vector in \code{par.draws}, the D-error will be
#' calculated and the design will be optimised locally. By specifying a matrix,
#' in which each row is a draw from a multivariate distribution, the DB-error
#' will be calculated, and the design will be optimised globally. Whenever there
#' are alternative specific constants, \code{par.draws} should be a list
#' containing two matrices: The first matrix containing the parameter draws for
#' the alternative specific constant parameters. The second matrix containing
#' the draws for the rest of the parameters.
#'
#' The DB-error is calculated by taking the mean over D-errors. It could be that
#' for some draws the design results in an infinite D-error. The percentage of
#' draws for which this was true for the final design can be found in the output
#' \code{inf.error}.
#'
#' Alternative specific constants can be specified in \code{alt.cte}. The length
#' of this binary vector should equal \code{n.alts}, were \code{0} indicates the
#' absence of an alternative specific constant and \code{1} the opposite.
#'
#' \code{start.des} is a list with one or several matrices corresponding to
#' initial start design(s). In each matrix each row is a profile. The number of rows
#' equals \code{n.sets * n.alts}, and the
#' number of columns equals the number of columns of the design matrix + the
#' number of non-zero elements in \code{alt.cte}. Consider that for a
#' categorical attribute with *p* levels, there are *p - 1* columns in the design
#' matrix, whereas for a continuous attribute there is only one column. If
#' \code{start.des = NULL}, \code{n.start} random initial designs will be
#' generated. If start designs are provided, \code{n.start} is ignored.
#'
#' If \code{no.choice} is \code{TRUE}, in each choice set an alternative with
#' one alternative specific constant is added. The return value of the
#' D(B)-error is however based on the design without the no choice option.
#'
#' When \code{parallel} is \code{TRUE}, \code{\link[parallel]{detectCores}} will
#' be used to decide upon the number of available cores. That number minus 1
#' cores will be used to search for efficient designs. The computation time will
#' decrease significantly when \code{parallel = TRUE}.
#'
#' @param lvls A numeric vector which contains for each attribute the number
#' of levels.
#' @param coding Type of coding that needs to be used for each attribute.
#' @param c.lvls A list containing numeric vectors with the attribute levels for
#' each continuous attribute. The default is \code{NULL}.
#' @param n.sets Numeric value indicating the number of choice sets.
#' @param n.alts Numeric value indicating the number of alternatives per choice
#' set.
#' @param par.draws A matrix or a list, depending on \code{alt.cte}.
#' @param alt.cte A binary vector indicating for each alternative whether an
#' alternative specific constant is desired. The default is \code{NULL}.
#' @param no.choice A logical value indicating whether a no choice alternative
#' should be added to each choice set. The default is \code{FALSE}.
#' @param start.des A list containing one or more matrices corresponding to
#' initial start design(s). The default is \code{NULL}.
#' @param parallel Logical value indicating whether computations should be done
#' over multiple cores. The default is \code{TRUE}.
#' @param max.iter A numeric value indicating the maximum number allowed
#' iterations. The default is \code{Inf}.
#' @param n.start A numeric value indicating the number of random start designs
#' to use. The default is 12.
#' @param best A logical value indicating whether only the best design should be
#' returned. The default is \code{TRUE}.
#' @return
#' If \code{best = TRUE} the design with the lowest D(B)-error is returned.
#' If \code{best = FALSE}, the results of all (provided) start designs are
#' returned. \item{design}{A
#' numeric matrix wich contains an efficient design.} \item{error}{Numeric
#' value indicating the D(B)-error of the design.} \item{inf.error}{Numeric
#' value indicating the percentage of draws for which the D-error was
#' \code{Inf}.} \item{probs}{Numeric matrix containing the probabilities of
#' each alternative in each choice set. If a sample matrix was provided in
#' \code{par.draws}, this is the average over all draws.}
#' @examples
#' \donttest{
#' # DB-efficient designs
#' # 3 Attributes, all dummy coded. 1 alternative specific constant = 7 parameters
#' mu <- c(1.2, 0.8, 0.2, -0.3, -1.2, 1.6, 2.2) # Prior parameter vector
#' v <- diag(length(mu)) # Prior variance.
#' set.seed(123)
#' pd <- MASS::mvrnorm(n = 10, mu = mu, Sigma = v) # 10 draws.
#' p.d <- list(matrix(pd[,1], ncol = 1), pd[,2:7])
#' CEA(lvls = c(3, 3, 3), coding = c("D", "D", "D"), par.draws = p.d,
#' n.alts = 2, n.sets = 8, parallel = FALSE, alt.cte = c(0, 1))
#'
#' # DB-efficient design with categorical and continuous factors
#' # 2 categorical attributes with 4 and 2 levels (effect coded) and 1
#' # continuous attribute (= 5 parameters)
#' mu <- c(0.5, 0.8, 0.2, 0.4, 0.3)
#' v <- diag(length(mu)) # Prior variance.
#' set.seed(123)
#' pd <- MASS::mvrnorm(n = 3, mu = mu, Sigma = v) # 10 draws.
#' CEA(lvls = c(4, 2, 3), coding = c("E", "E", "C"), par.draws = pd,
#' c.lvls = list(c(2, 4, 6)), n.alts = 2, n.sets = 6, parallel = FALSE)
#'
#' # DB-efficient design with start design provided.
#' # 3 Attributes with 3 levels, all dummy coded (= 6 parameters).
#' mu <- c(0.8, 0.2, -0.3, -0.2, 0.7, 0.4)
#' v <- diag(length(mu)) # Prior variance.
#' sd <- list(example_design)
#' set.seed(123)
#' ps <- MASS::mvrnorm(n = 10, mu = mu, Sigma = v) # 10 draws.
#' CEA(lvls = c(3, 3, 3), coding = c("D", "D", "D"), par.draws = ps,
#' n.alts = 2, n.sets = 8, parallel = FALSE, start.des = sd)
#'}
#' @importFrom Rdpack reprompt
#' @export
CEA <- function(lvls, coding, c.lvls = NULL, n.sets, n.alts, par.draws,
alt.cte = NULL, no.choice = FALSE, start.des = NULL,
parallel = TRUE, max.iter = Inf, n.start = 12, best = TRUE) {
### Error handling for creating initial random design
# Error lvls vector. At least 2 attributes
if (length(lvls) < 2 || (!(is.numeric(lvls)))) {
stop("lvls argument is incorrect.")
}
# Error correct coding types.
codings.types <- c("E", "D", "C")
if (!all(coding %in% codings.types) || (length(coding) != length(lvls))) {
stop("coding argument is incorrect.")
}
# Continuous attributes.
contins <- which(coding == "C")
n.contins <- length(contins)
# error continuous levels (Not a list)
if (!is.null(c.lvls) && !is.list(c.lvls)) {
stop('c.lvls should be a list.')
}
# Error continuous specified and NULL.
if (n.contins > 0 && is.null(c.lvls)) {
stop("when 'coding' contains C, 'c.lvls' should be specified")
}
# Error continuous levels specification.
if (!is.null(c.lvls)) {
if (length(c.lvls) != n.contins) {
stop("length of 'c.lvls' does not match number of specified continuous attributes in 'coding'")
}
# Error c.lvls same number of levels.
if (!isTRUE(all.equal(lvls[contins], lengths(c.lvls)))) {
stop("the number of levels provided in 'c.lvls' does not match the expected based on 'lvls'")
}
}
### Error handling for design specifications
# If no alternative constant is given, create the variable as a vector of 0s
if (is.null(alt.cte)) {
alt.cte <- rep(0L, n.alts)
}
#init
n.cte <- length(which(alt.cte == 1))
# If only one draw is given, transform it to a matrix
if (!is.list(par.draws)) {
if (is.vector(par.draws)) {
par.draws <- matrix(par.draws, nrow = 1)
}
}
# Alternative constant errors
if (length(alt.cte) != n.alts) {
stop("'n.alts' does not match the 'alt.cte' vector")
}
if (!all(alt.cte %in% c(0, 1))) {
stop("'alt.cte' should only contain zero or ones.")
}
# No choice errors
if (!is.logical(no.choice)) {
stop("'no.choice' should be TRUE or FALSE")
}
if (no.choice) {
if (!isTRUE(all.equal(alt.cte[n.alts], 1))) {
stop("if 'no.choice' is TRUE, the last alternative constant should equal 1.")
}
ncsek <- seq(n.alts, (n.sets * n.alts), n.alts) # Rows in the design matrix
# to be assigned as zeros in all attributes
} else {
ncsek <- NULL
}
# Handling par.draws with alternative specific contstants.
# This conditional is when there is only one alternative constant
if (isTRUE(all.equal(n.cte, 1))) {
if (!(is.list(par.draws))) {
stop("par.draws should be a list")
}
if (!isTRUE(all.equal(length(par.draws), 2))) {
stop("'par.draws' should contain two components")
}
# If only draws for the constant are given in a vector, transform it to a
# matrix
if (is.vector(par.draws[[1]])) {
par.draws[[1]] <- matrix(par.draws[[1]], ncol = 1)
}
if (!(all(unlist(lapply(par.draws, is.matrix))))) {
stop("'par.draws' should contain two matrices")
}
dims <- as.data.frame(lapply(par.draws, dim))
if (!isTRUE(all.equal(dims[1, 1], dims[1, 2]))) {
stop("the number of rows in the components of 'par.draws' should be equal")
}
par.draws <- do.call("cbind", par.draws) # Transform draws to a matrix
} else {
if (n.cte > 1.2) {
if (!(is.list(par.draws))) {
stop("par.draws should be a list")
}
if (!isTRUE(all.equal(length(par.draws), 2))) {
stop("'par.draws' should contain two components")
}
if (!(all(unlist(lapply(par.draws, is.matrix))))) {
stop("'par.draws' should contain two matrices")
}
if (!isTRUE(all.equal(ncol(par.draws[[1]]), n.cte))) {
stop("the first component of 'par.draws' should contain the same number of columns as there are non zero elements in 'alt.cte'")
}
dims <- as.data.frame(lapply(par.draws, dim))
if (!isTRUE(all.equal(dims[1, 1], dims[1, 2]))) {
stop("the number of rows in the components of 'par.draws' should be equal")
}
par.draws <- do.call("cbind", par.draws) # Transform draws to a matrix
}
}
# Error identifying model.
if (n.sets < ncol(par.draws)) {
stop("Model is unidentified. Increase the number of choice sets or decrease parameters to estimate.")
}
### Create all levels for each attribute
# Change into correct coding.
coding <- dplyr::recode(coding, D = "contr.treatment", E = "contr.sum")
# Create all combinations of attribute levels.
levels.list <- lapply(X = as.list(lvls), function(x) (1:x))
# Replace continuous.
levels.list[contins] <- c.lvls
# Transform to matrix
levels.list[contins] <- lapply(X = levels.list[contins], as.matrix)
# Transform categorical attributes
categ <- which(coding %in% c("contr.treatment", "contr.sum"))
n.categ <- length(categ)
if (n.categ > 0) {
levels.list[categ] <- lapply(X = levels.list[categ], factor)
# Apply coding
if (n.contins > 0) {
for (i in 1:length(lvls)) {
if (!(i %in% contins)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
}else {
for (i in 1:length(lvls)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
# Compute all possible values for each categorical attribute
levels.list[categ] <- lapply(X = levels.list[categ], stats::contrasts)
}
# Set colnames for the design matrix
c.nam = list()
for (i in 1:length(lvls)) {
if (coding[i] == "contr.treatment") {
c.nam[[i]] <- paste("Var", i, 2:lvls[i], sep = "")
} else {
if (coding[i] == "contr.sum") {
c.nam[[i]] <- paste("Var", i, 1:(lvls[i] - 1), sep = "")
} else {
c.nam[[i]] <- paste("Var", i, sep = "")
}
}
}
c.nam = unlist(c.nam)
# Create alternative specific design.
cte.des <- Altspec(alt.cte = alt.cte, n.sets = n.sets)
# Count the number of colums of the design matrix without constants
ncol.des.noconst <- sum(unlist(lapply(levels.list,ncol)))
#ncol.des.noconst <- sum(unlist(lapply(levels.list,function(x){
# if (is.matrix(x)) { ncol(x) }
# else{if (is.numeric(x)) { 1 }}
#})))
if (!identical(as.integer(ncol(par.draws)),
as.integer(n.cte + ncol.des.noconst))) {
stop("The sum of the number of columns in the components of 'par.draws' should equal the number of columns of design matrix (including alternative specific constants)")
}
### Random initial design.
if (is.null(start.des)) {
#create start designs
nr.starts <- n.start
start.des <- vector(mode = 'list', length = nr.starts)
okstart <- FALSE
while (okstart == FALSE) {
for (i in 1:nr.starts) {
# r is to know which levels to take in each attribute
r <- NULL
start <- NULL
for (j in 1:length(lvls)) {
r <- round(stats::runif((n.sets * n.alts), 1, lvls[j]))
start <- cbind(start, levels.list[[j]][r,])
}
colnames(start) <- c.nam
start.des[[i]] <- cbind(cte.des, start)
if (no.choice) {
start.des[[i]][ncsek, (ncol(cte.des) + 1):(ncol(cte.des) + ncol.des.noconst)] <- c(rep(0, ncol.des.noconst))
}
}
# Compute D-optimality for each design and each draw
d.start <- lapply(start.des, StartDB, par.draws, n.alts)
# If the DB-optimality of any starting desing is finite, continue
if (any(is.finite(unlist(lapply(d.start, mean, na.rm = TRUE))))) {
okstart <- TRUE
}
}
} else {
if (!is.list(start.des)) {
stop("'start.des' should be a list")
}
if (!(all(unlist(lapply(start.des, is.matrix))))) {
stop("'start.des' should contain matrices as components")
}
# Save the dimension of each starting design
dimstart <- as.matrix(lapply(start.des, dim))
# Save the number of random starts given
nr.starts <- length(dimstart)
if (nr.starts > 1.5) {
if (!isTRUE(all.equal(length(unique(unlist(dimstart))), 2))) {
stop("start designs have different dimensions")
}
}
if (!isTRUE(all.equal(n.alts * n.sets, unique(unlist(dimstart))[1]))) {
stop("number of rows of start design(s) does not match with 'n.alts' * 'n.sets'")
}
if (!isTRUE(all.equal(as.integer(n.cte + ncol.des.noconst),
unique(unlist(dimstart))[2]))) {
stop("number of columns of start design(s) does not match with the number of columns in the design matrix")
}
# Compute D-optimality for each design and each draw
d.start <- lapply(start.des, StartDB, par.draws, n.alts)
if (!any(is.finite(unlist(lapply(d.start, mean, na.rm = TRUE))))) {
stop("One or more of the provided start designs resulted in an unvalid db-error.")
}
}
### Improving the initial design
if (parallel) {
no_cores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(no_cores)
parallel::clusterExport(cl, c("n.sets", "par.draws", "n.alts", "n.cte",
"alt.cte", "no.choice", "max.iter","ncsek"),
envir = environment())
deslist <- parallel::parLapply(cl, start.des, CEAcore_ucpp, par.draws,
levels.list, n.alts, n.sets, n.cte, alt.cte,
no.choice, max.iter, ncsek)
parallel::stopCluster(cl)
} else {
deslist <- lapply(start.des, CEAcore_ucpp, par.draws, levels.list,
n.alts, n.sets, n.cte, alt.cte, no.choice,
max.iter = max.iter, ncsek)
}
bestdes <- deslist[[which.min(unlist(lapply(deslist, function(x) (x$error))))]]
ifelse(best, return(bestdes), return(deslist))
}
# Core of Coordinate Exchange
CEAcore_ucpp <- function(des, par.draws, levels.list, n.alts, n.sets, n.cte,
alt.cte, no.choice, max.iter, ncsek) {
converge <- FALSE # Boolean for convergence
change <- FALSE # Boolean for a change in the attribute
it <- 1 # Indicator of number of iterations
n.samples <- nrow(par.draws) # Number of samples from distribution of betas
n.par <- ncol(des) # Number of parameters
# Beginning of algorithm
while (!converge & it <= max.iter) {
# Compute DB-optimality of initial design
db.start <- mean(apply(par.draws, 1, Derr_ucpp, des = des,
n.alts = n.alts), na.rm = TRUE)
it <- it + 1 # Increase iteration
# Save design before iteration.
iter.des <- des
# For every row in the design.
sek <- 1:nrow(des)
if (no.choice) {
sek <- sek[-ncsek] # If no.choice is given, then it is not improved
}
# Loop for each row of the initial design
for (i in sek) {
# Loop for each attribute in the row
for (j in 1:length(levels.list)) {
#%%%%%%%
# This can be improved by removing mods object and just replace the
# initial design with the new attributes to compute the db.error and
# then if one of those is better than the initial db error. Replace
# finally the initial design (as it is in modfed)
# Initialize mods object to track changes of the design in each attribute
mods <- vector("list", nrow(levels.list[[j]]))
mods <- lapply(mods,function(x){des})
# Initialize db object for db error of each level
db <- numeric(nrow(levels.list[[j]]))
# Indicator of columns modified
ncol.ok <- sum(unlist(lapply(levels.list[(1:j - 1)],ncol)))
# Indicator of columns to change
ncol.mod <- ifelse( j == 1, (n.cte + j),
(n.cte + 1 + ncol.ok))
# ncol.left tracks which columns are left to be improved in the row
ncol.left <- sum(unlist(lapply(levels.list[-(1:j)],ncol)))
# Loop for each level in the attribute
for (k in 1:nrow(levels.list[[j]])) {
# Update design with new attribute
#mods[[k]][i, (n.cte + j):(n.par - ncol.left) ] <- levels.list[[j]][k,]
mods[[k]][i, ncol.mod:(n.par - ncol.left) ] <- levels.list[[j]][k,]
# Calculate D-optimality for each draw
d.errors <- apply(par.draws, 1, Derr_ucpp, des = mods[[k]],
n.alts = n.alts)
# Compute DB-optimality
db[k] <- mean(d.errors, na.rm = TRUE)
}
pr <- which.min(db)
db <- min(db)
# Change if lower db error.
if (!is.na(db) && !is.na(db.start)) {
if (db < db.start) {
des <- mods[[pr]]
db.start <- db
}
}
} # End loop for each attribute in the row
} # End loop for each row of the initial design
converge <- isTRUE(all.equal(des, iter.des)) # Convergence if no profile is swapped this iteration.
} # End while (after convergence)
#return(it, des)
# calculate percentage NA values.
na.percentage <- 0
d.errors <- apply(par.draws, 1, Derr_ucpp, des = des, n.alts = n.alts)
if (any(is.na(d.errors))) {
na.percentage <- scales::percent(sum(is.na(d.errors)) / n.samples)
}
# Utility balance.
# ub <- apply(par.draws, 1, Utbal, des = des, n.alts = n.alts)
# ub2 <- apply(par.draws, 1, InfoDes2, des = des, n.alts = n.alts, utbal = TRUE)
# ub_ucpp <- apply(par.draws, 1, InfoDes_cpp, des = des, n_alts = n.alts,
# utbal = TRUE)
# Utility balance using c++ function
ub <- apply(par.draws, 1, InfoDes_cpp, des = des, n_alts = n.alts,
utbal = TRUE)
pmat <- matrix(rowMeans(ub), ncol = n.alts, byrow = TRUE)
rownames(pmat) <- paste("set", 1:n.sets, sep = "")
colnames(pmat) <- paste(paste("Pr(", paste("alt", 1:n.alts, sep = ""),
sep = ""), ")", sep = "")
if (no.choice) {
colnames(pmat)[n.alts] <- "Pr(no choice)"
}
# Rownames design.
des.names <- Rcnames(n.sets = n.sets, n.alts = n.alts, alt.cte = alt.cte,
no.choice = no.choice)
rownames(des) <- des.names[[1]]
# Colnames alternative specific constants.
if (n.cte != 0 && !is.null(colnames(des))) {
colnames(des)[1:n.cte] <- des.names[[2]]
}
# Return design, D(B)error, percentage NA's, utility balance.
return(list("design" = des, "error" = db.start, "inf.error" = na.percentage,
"probs" = pmat))
}
#' Sequential Coordinate Exchange algorithm for MNL model.
#'
#' Selects the choice set that minimizes the DB-error when added to an initial
#' design, given (updated) parameter values.
#'
#' This algorithm is ideally used in an adaptive context. The algorithm will
#' select the next DB-efficient choice set given parameter values and possible
#' previously generated choice sets. In an adaptive context these parameter
#' values are updated after each observed response.
#'
#' Previously generated choice sets, which together form an initial design, can
#' be provided in \code{des}. When no design is provided, the algorithm will
#' select the most efficient choice set based on the fisher information of the
#' prior covariance matrix \code{prior.covar}.
#'
#' If \code{alt.cte = NULL}, \code{par.draws} should be a matrix in which each
#' row is a sample from the multivariate parameter distribution. In case that
#' \code{alt.cte} is not \code{NULL}, a list containing two matrices should be
#' provided to \code{par.draws}. The first matrix containing the parameter draws
#' for the alternative specific parameters. The second matrix containing the
#' draws for the rest of the parameters.
#'
#' The list of potential choice sets is created by selecting randomly a level for
#' each attribute in an alternative/profile. \code{n.cs} controls the number of
#' potential choice sets to consider. The default is \code{
#' NULL}, which means that the number of possible choice sets is the product of
#' attribute levels considered in the experiment. For instance, an experiment
#' with 3 attribute and 3 levels each will consider 3^3 = 27 possible choice sets.
#'
#' The \code{weights} argument can be used when the \code{par.draws} have
#' weights. This is for example the case when parameter values are updated using
#' \code{\link{ImpsampMNL}}.
#'
#' When \code{parallel} is \code{TRUE}, \code{\link[parallel]{detectCores}} will
#' be used to decide upon the number of available cores. That number minus 1
#' cores will be used to search for the optimal choice set. For small problems
#' (6 parameters), \code{parallel = TRUE} can be slower. For larger problems the
#' computation time will decrease significantly.
#'
#' *Note:* this function is faster than \code{\link[idefix]{SeqMOD}}, but
#' the output is not as stable. This happens because this function
#' makes a random search to get the choice set, whereas
#' \code{\link[idefix]{SeqMOD}} makes an exhaustive search.
#' @inheritParams CEA
#' @param par.draws A matrix or a list, depending on \code{alt.cte}.
#' @param des A design matrix in which each row is a profile. If alternative
#' specific constants are present, those should be included as the first
#' column(s) of the design. Can be generated with \code{\link{Modfed}} or
#' \code{\link{CEA}}
#' @param n.cs An integer indicating the number of possible random choice sets to
#' consider in the search for the next best choice set possible. The default is
#' \code{NULL}.
#' @param prior.covar Covariance matrix of the prior distribution.
#' @param weights A vector containing the weights of the draws. Default is
#' \code{NULL}. See also \code{\link{ImpsampMNL}}.
#' @param parallel Logical value indicating whether computations should be done
#' over multiple cores.
#' @param no.choice An integer indicating the no choice alternative. The default
#' is \code{NULL}.
#' @param reduce Logical value indicating whether the candidate set should be
#' reduced or not.
#' @return \item{set}{A matrix representing a DB efficient choice set.}
#' \item{error}{A numeric value indicating the DB-error of the whole
#' design.}
#' @importFrom Rdpack reprompt
#' @references \insertRef{idefix}{idefix}
#' @references \insertRef{ju}{idefix}
#' @references \insertRef{cea}{idefix}
#' @references \insertRef{cea_discrete}{idefix}
#' @examples
#' # DB efficient choice set, given a design and parameter draws.
#' # 3 attributes with 3 levels each
#' m <- c(0.3, 0.2, -0.3, -0.2, 1.1, 2.4) # mean (total = 6 parameters).
#' pc <- diag(length(m)) # covariance matrix
#' set.seed(123)
#' sample <- MASS::mvrnorm(n = 10, mu = m, Sigma = pc)
#' # Initial design.
#' des <- example_design
#' # Efficient choice set to add.
#' SeqCEA(des = des, lvls = c(3, 3, 3), coding = c("D", "D", "D"), n.alts = 2,
#' par.draws = sample, prior.covar = pc, parallel = FALSE)
#'
#' # DB efficient choice set, given parameter draws.
#' # with alternative specific constants
#' des <- example_design2
#' ac <- c(1, 1, 0) # Alternative specific constants.
#' m <- c(0.3, 0.2, -0.3, -0.2, 1.1, 2.4, 1.8, 1.2) # mean
#' pc <- diag(length(m)) # covariance matrix
#' pos <- MASS::mvrnorm(n = 10, mu = m, Sigma = pc)
#' sample <- list(pos[ , 1:2], pos[ , 3:8])
#' # Efficient choice set.
#' SeqCEA(des = des, lvls = c(3, 3, 3), coding = c("D", "D", "D"), n.alts = 3,
#' par.draws = sample, alt.cte = ac, prior.covar = pc, parallel = FALSE)
#' @export
SeqCEA <- function(des = NULL, lvls, coding, c.lvls = NULL, n.alts, par.draws,
prior.covar, alt.cte = NULL, no.choice = NULL,
weights = NULL, parallel = TRUE, reduce = TRUE, n.cs = NULL) {
# Error handling initial design
if (is.null(des)) {
n.sets <- 1L
} else {
if (!is.matrix(des)) {
stop("'des' should be a matrix or NULL")
}
if (!isTRUE(nrow(des) %% n.alts == 0)) {
stop("'n.alts' does not seem to match with the number of rows in 'des'")
}
n.sets <- nrow(des) / n.alts
}
### Error handling for design specifications
# No choice errors
if (!is.null(no.choice)) {
if (!is.wholenumber(no.choice)) {
stop("'no.choice' should be an integer or NULL")
}
if (any(isTRUE(no.choice > (n.alts + 0.2)), isTRUE(no.choice < 0.2))) {
stop("'no.choice' does not indicate one of the alternatives")
}
if (is.null(alt.cte)) {
stop("if there is a no choice alternative, 'alt.cte' should be specified")
}
if (!isTRUE(all.equal(alt.cte[no.choice], 1))) {
stop("the no choice alternative should correspond with a 1 in 'alt.cte'")
}
}
# Alternative constant errors
if (!is.null(alt.cte)) {
if (length(alt.cte) != n.alts) {
stop("'n.alts' does not match the 'alt.cte' vector")
}
if (!all(alt.cte %in% c(0, 1))) {
stop("'alt.cte' should only contain zero or ones.")
}
# alternative specific constants
n.cte <- length(which(alt.cte == 1L))
if (isTRUE(all.equal(n.cte, 0L))) {
alt.cte <- NULL
cte.des <- NULL
}
# Handling errors when there are alternative constants
if (!is.null(alt.cte)) {
if (n.cte > 0.2) {
if (!is.list(par.draws)) {
stop("'par.draws' should be a list when 'alt.cte' is not NULL")
}
if (!isTRUE(all.equal(length(par.draws), 2))) {
stop("'par.draws' should contain two components")
}
# If there is only one specific constant and is a vector, then it is
# transformed to a matrix
if (isTRUE(all.equal(n.cte, 1))) {
if (is.vector(par.draws[[1]])) {
par.draws[[1]] <- matrix(par.draws[[1]], ncol = 1)
}
}
if (!(all(unlist(lapply(par.draws, is.matrix))))) {
stop("'par.draws' should contain two matrices")
}
if (!isTRUE(all.equal(ncol(par.draws[[1]]), n.cte))) {
stop("the first component of 'par.draws' should contain the same number
of columns as there are non zero elements in 'alt.cte'")
}
dims <- as.data.frame(lapply(par.draws, dim))
if (!isTRUE(all.equal(dims[1, 1], dims[1, 2]))) {
stop("the number of rows in the components of 'par.draws' should be equal")
}
par.draws <- do.call("cbind", par.draws) # Transform par.draws to a matrix
}
# Create alternative specific design.
cte.des <- Altspec(alt.cte = alt.cte, n.sets = n.sets)
cte.set <- matrix(cte.des[1:n.alts, ], ncol = n.cte, byrow = FALSE)
}
} else {
cte.des <- NULL
n.cte <- 0
# if no alternative constants
if (!is.matrix(par.draws)) {
stop("'par.draws'should be a matrix when 'alt.cte' = NULL")
}
}
# Weights errors
n.par <- ncol(par.draws)
if (!is.null(weights)) {
if (!isTRUE(all.equal(length(weights), nrow(par.draws)))) {
stop("length of 'weights' does not match number total number of rows in 'par.draws'")
}
} else {
weights <- rep(1L, nrow(par.draws))
}
### Create all levels for each attribute
# Continuous attributes.
contins <- which(coding == "C")
n.contins <- length(contins)
# Change into correct coding.
coding <- dplyr::recode(coding, D = "contr.treatment", E = "contr.sum")
# Create all combinations of attribute levels.
levels.list <- lapply(X = as.list(lvls), function(x) (1:x))
# Replace continuous.
levels.list[contins] <- c.lvls
# Transform to matrix
levels.list[contins] <- lapply(X = levels.list[contins], as.matrix)
# Transform categorical attributes
categ <- which(coding %in% c("contr.treatment", "contr.sum"))
n.categ <- length(categ)
if (n.categ > 0) {
levels.list[categ] <- lapply(X = levels.list[categ], factor)
# Apply coding
if (n.contins > 0) {
for (i in 1:length(lvls)) {
if (!(i %in% contins)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
}else {
for (i in 1:length(lvls)) {
stats::contrasts(levels.list[[i]]) <- coding[i]
}
}
# Compute all possible values for each categorical attribute
levels.list[categ] <- lapply(X = levels.list[categ], stats::contrasts)
}
# Set colnames for the design matrix
c.nam = list()
for (i in 1:length(lvls)) {
if (coding[i] == "contr.treatment") {
c.nam[[i]] <- paste("Var", i, 2:lvls[i], sep = "")
} else {
if (coding[i] == "contr.sum") {
c.nam[[i]] <- paste("Var", i, 1:(lvls[i] - 1), sep = "")
} else {
c.nam[[i]] <- paste("Var", i, sep = "")
}
}
}
c.nam = unlist(c.nam)
# Count the number of colums of the design matrix without constants
ncol.des.noconst <- sum(unlist(lapply(levels.list,ncol)))
# if (!identical(as.integer(n.par), as.integer(n.cte + ncol.des.noconst))) {
# stop("The sum of the number of columns in the components of 'par.draws' should equal the number of columns of design matrix (including alternative specific constants)")
# }
if (!identical(as.integer(ncol(prior.covar)),
as.integer(n.cte + ncol.des.noconst))) {
stop("number of columns of 'prior.covar' does not equal the number of columns of design matrix (including alternative specific constants)")
}
## When a design is supplied
if (!is.null(des)) {
# Error par.draws
if (!isTRUE(all.equal(ncol(des), n.par))) {
stop("number of columns in 'par.draws' does not match the number of columns in 'des'")
}
# Error dimension of initial design and design matrix
if (!identical(as.integer(ncol(des)), as.integer(n.cte + ncol.des.noconst))) {
stop("number of columns in 'des' does not match the number of columns of design matrix (including alternative specific constants)")
}
# Starting and initializing values.
i.cov <- solve(prior.covar)
d.start <- apply(par.draws, 1, DerrC_ucpp, des = des, n.alts = n.alts,
i.cov = i.cov)
db.start <- mean(d.start, na.rm = TRUE)
full.comb <- Newsets_ucpp(levels.list = levels.list, n.alts = n.alts,
no.choice = no.choice, reduce = reduce,
c.names = c.nam, n.cs = n.cs)
# Adding alternative constants
if (!is.null(cte.des)) {
full.comb <- lapply(full.comb, function(x) cbind(cte.set, x))
}
# Adding these new choice sets to the initial design
full.des <- lapply(full.comb, function(x) rbind(des, x))
# Select the next best choice set
if (parallel) {
no_cores <- parallel::detectCores() - 1L
cl <- parallel::makeCluster(no_cores)
db.errors <- parallel::parLapply(cl, full.des, DBerrS.P_ucpp, par.draws,
n.alts, i.cov, weights)
parallel::stopCluster(cl)
} else {
# For each potential set, select best.
db.errors <- lapply(full.des, DBerrS.P_ucpp, par.draws, n.alts, i.cov,
weights)
}
dbs <- unlist(db.errors, use.names = FALSE)
set <- full.comb[[which.min(dbs)]]
row.names(set) <- NULL
colnames(set) <- colnames(des)
db <- min(dbs)
#return best set and db error design.
return(list("set" = set, "error" = db))
} # End if when a design is supplied
else {
# Starting and initializing values.
i.cov <- solve(prior.covar)
full.comb <- Newsets_ucpp(levels.list = levels.list, n.alts = n.alts,
no.choice = no.choice, reduce = reduce,
c.names = c.nam, n.cs = n.cs)
# Adding alternative constants
if (!is.null(cte.des)) {
full.comb <- lapply(full.comb, function(x) cbind(cte.set, x))
}
# Select the next best choice set
if (parallel) {
no_cores <- parallel::detectCores() - 1L
cl <- parallel::makeCluster(no_cores)
db.errors <- parallel::parLapply(cl, full.comb, DBerrS.P_ucpp, par.draws,
n.alts, i.cov, weights)
parallel::stopCluster(cl)
} else {
# For each potential set, select best.
db.errors <- lapply(full.comb, DBerrS.P_ucpp, par.draws, n.alts, i.cov,
weights)
}
dbs <- unlist(db.errors, use.names = FALSE)
set <- full.comb[[which.min(dbs)]]
row.names(set) <- NULL
if (n.cte > 0) {
colnames(set) <- c(paste("alt",1:n.cte,".cte", sep = ""), c.nam)
}
db <- min(dbs)
#return best set and db error design.
return(list("set" = set, "error" = db))
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.