R/design.R

In spdesign: Designing Stated Preference Experiments

#' Generate an efficient experimental design
#'
#' The function generates efficient experimental designs. The function takes
#' a set of indirect utility functions and generates efficient experimental
#' designs assuming that people are maximizing utility.
#'
#' No assumptions are made with respect to default values and it is up to the
#' user to specify optimization criteria, optmization routines, draws to use for
#' Bayesian priors and more.
#'
#' @param utility A named list of utility functions. See the examples and the
#' vignette for examples of how to define these correctly for different types
#' of experimental designs.
#' @param rows An integer giving the number of rows in the final design
#' @param model A character string indicating the model to optimize the design
#' for. Currently the only model programmed is the 'mnl' model and this is also
#' set as the default.
#' @param efficiency_criteria A character string giving the efficiency criteria
#' to optimize for. One of 'a-error', 'c-error', 'd-error' or 's-error'. No
#' default is set and argument must be specified. Optimizing for multiple
#' criteria is not yet implemented and will result in an error.
#' @param algorithm A character string giving the optimization algorithm to use.
#' No default is set and the argument must be specified to be one of 'rsc',
#' 'federov' or 'random'.
#' @param draws The type of draws to use with Bayesian priors. No default is set
#'  and must be specified even if you are not creating a Bayesian design. Can be
#' one of "pseudo-random", "mlhs", "standard-halton", "scrambled-halton",
#' "standard-sobol","scrambled-sobol".
#' @param R An integer giving the number of draws to use. The default is 100.
#' @param dudx A character string giving the name of the prior in the
#' denominator. Must be specified when optimizing for 'c-error'
#' @param candidate_set A matrix or data frame in the "wide" format containing
#' all permitted combinations of attributes. The default is NULL. If no
#' candidate set is provided, then it is built from the utility functions
#' subject to specified exclusions with \code{\link{build_candidate_set}}.
#' This is passed in as an object and not a character
#' string. The candidate set will be expanded to include zero columns to
#' consider alternative specific attributes.
#' @param exclusions A list of exclusions Often this list will be pulled
#' directly from the list of options or it is a modified list of exclusions
#' @param save_designs A boolean indicating whether to save up to 10 intermediate designs. The default value is FALSE.
#' @param control A list of control options. Set `allow_reversed_pairs = TRUE`
#' to include the profiles of exchangeable alternatives in every order when the
#' candidate set is built. See \code{\link{build_candidate_set}}.
#'
#' @return An object of class 'spdesign'. For the 'federov' algorithm, the list
#' element 'runs' holds the best design of each run of the search and its
#' efficiency criteria. The best design across all runs is the one returned as
#' the design.
#'
#' @export
generate_design <- function(
  utility,
  rows,
  model = "mnl",
  efficiency_criteria = c("a-error", "c-error", "d-error", "s-error"),
  algorithm = c("federov", "rsc", "random"),
  draws = c(
    "pseudo-random",
    "mlhs",
    "standard-halton",
    "scrambled-halton",
    "standard-sobol",
    "scrambled-sobol"
  ),
  R = 100,
  dudx = NULL,
  candidate_set = NULL,
  exclusions = NULL,
  save_designs = FALSE,
  control = list(
    cores = 1,
    max_iter = 10000,
    max_relabel = 10000,
    max_no_improve = 100000,
    efficiency_threshold = 0.000001,
    sample_with_replacement = FALSE,
    allow_reversed_pairs = FALSE
  )
) {
  # Match and check model arguments ----
  cli_h2("Checking function arguments")

  ## Create the design object ----
  design_object <- list()
  class(design_object) <- "spdesign"

  design_object[["utility"]] <- utility
  design_object[["time"]] <- list(
    time_start = Sys.time()
  )

  ## Match arguments ----
  design_object[["model"]] <- model <- match.arg(model)
  efficiency_criteria <- match.arg(efficiency_criteria, several.ok = TRUE)
  algorithm <- match.arg(algorithm)
  draws <- match.arg(draws)

  ## Check arguments ----
  stopifnot(is.list(utility) && !is.data.frame(utility))
  stopifnot(length(utility) > 1)
  stopifnot(all(do.call(c, lapply(utility, is_balanced, "[", "]"))))
  stopifnot(all(do.call(c, lapply(utility, is_balanced, "(", ")"))))
  # stopifnot(all_priors_and_levels_specified(utility))
  stopifnot(!any_duplicates(utility))
  stopifnot(!too_small(utility, rows))

  problem <- str_subset(attribute_names(utility), "_dummy$")

  if (length(problem) > 0) {
    stop(
      "Attribute names cannot end in '_dummy'. To dummy-code an attribute, ",
      "add '_dummy' to its parameter instead, e.g. ",
      "'b_x1_dummy[c(0.1, 0.2)] * x1[c(1, 2, 3)]'. Please rename: ",
      paste(problem, collapse = ", ")
    )
  }

  problem <- invalid_dummy_coding(utility)

  if (length(problem) > 0) {
    stop(
      "Dummy-coded attributes must have the levels 1, 2, ..., K, where 1 is ",
      "the base level, and K - 1 priors. Please check the levels and priors ",
      "of: ",
      paste(problem, collapse = ", ")
    )
  }

  # Set the default for control and replace the specified values in control
  default_control <- list(
    cores = 1,
    max_iter = 10000,
    max_relabel = 10000,
    max_swap = 10000,
    max_no_improve = 100000,
    efficiency_threshold = 0.000001,
    sample_with_replacement = FALSE,
    allow_reversed_pairs = FALSE
  )

  control <- modifyList(default_control, control)

  if (algorithm == "rsc") {
    cli_alert_info(
      "The cycling part of the algorithm is not used. It only applies to a
      small subset of designs. The algorithm swithes between relabeling of
      attribute levels and swapping of attributes."
    )
  }

  if (length(efficiency_criteria) > 1) {
    stop("Optimizing over multiple criteria is not implemented")
  }

  if (is.null(dudx) && efficiency_criteria == "c-error") {
    stop("The denominator index 'dudx' must be specified for c-error")
  }

  # Consider a core-check if relevant at a later point.
  if (control$cores > 1) {
    warning(
      "Multicore is not implemented yet. Design will be optimized using a single core."
    )
    control$cores <- 1
  }

  ## Candidate set ----
  # We are only creating the candidate set if we are using a random or modified federov algorithm
  if (!is.null(candidate_set) & algorithm == "rsc") {
    stop(
      "To use your supplied candidate set you must use either the 'random' or 'federov' algorithms."
    )
  }

  if (algorithm %in% c("random", "federov")) {
    cli_h2("Checking the candidate set and applying exclusions")

    # If no candidate set is supplied, build it with the exclusions applied.
    # Otherwise, check the supplied candidate set and apply the exclusions.
    if (is.null(candidate_set)) {
      cli_alert_info(
        "No candidate set supplied. The candidate set is built from the utility functions subject to supplied constraints."
      )

      candidate_set <- build_candidate_set(
        utility,
        exclusions,
        control$allow_reversed_pairs
      )

      cli_alert_success("Candidate set created")
    } else {
      stopifnot((is.matrix(candidate_set) || is.data.frame(candidate_set)))

      candidate_names_idx <- names(candidate_set) %in%
        names(expand_attribute_levels(utility))

      if (!all(candidate_names_idx)) {
        problem <- paste(
          names(candidate_set)[!candidate_names_idx],
          collapse = ", "
        )

        stop(
          paste0(
            "There are more attributes specified in the candidate set than are present in the utility functions. ",
            problem,
            " are not specified in the utility function. This could also be caused by a mismatch in the names. The names should be of the form <utility list element name>_<attribute name>. For example, in your case, they should correspond to: '",
            paste(names(expand_attribute_levels(utility)), collapse = ", "),
            "' The candidate set must be supplied in 'wide' format."
          )
        )
      }

      # Extract only the specified in the utility function to check
      regex <- as_whole_word(attribute_names(utility))
      utility_attributes <- vector(mode = "list", length = length(utility))
      for (i in seq_along(utility)) {
        idx <- str_detect(utility[[i]], regex)
        utility_attributes[[i]] <- paste(
          names(utility[i]),
          attribute_names(utility)[idx],
          sep = "_"
        )
      }

      utility_attributes <- do.call(c, utility_attributes)

      if (!all(utility_attributes %in% names(candidate_set))) {
        stop(
          paste0(
            "Not all attributes specified in the utility functions are specified in the candidate set. This could be caused by a mismatch in the names. The names should be of the form <utility list element name>_<attribute name>. For example, in your case, they should correspond to: '",
            paste(utility_attributes, collapse = ", "),
            "' The candidate set must be supplied in 'wide' format."
          )
        )
      }

      # Expand candidate set to be square, i.e., fill in zero columns, for non-specified. This in case of
      # Alternative specific attributes!
      expanded_names <- names(expand_attribute_levels(utility))

      # Skip expansion if no alternative specific attributes are present
      if (any(!(expanded_names %in% utility_attributes))) {
        expr <- paste(
          "cbind(candidate_set, ",
          paste(
            paste(
              expanded_names[!(expanded_names %in% utility_attributes)],
              0,
              sep = " = "
            ),
            collapse = ", "
          ),
          ")"
        )
        candidate_set <- eval(parse(text = expr))
      }

      candidate_set <- candidate_set[, expanded_names]

      # Levels not listed in the utility functions are allowed, except for
      # attributes with level occurrences specified
      if (level_occurrences_specified(utility)) {
        problem <- unlisted_levels(utility, candidate_set, rows)

        if (length(problem) > 0) {
          stop(
            "Level occurrences are specified for ",
            paste(problem, collapse = ", "),
            ", but the candidate set contains levels for these attributes ",
            "that are not listed in the utility functions."
          )
        }
      }

      candidate_set <- exclude(candidate_set, exclusions)
    }

    # Transform the candiate set such that attributes that are dummy coded
    # are turned into factors. This ensures that we can use the model.matrix()
    for (i in which(names(candidate_set) %in% dummy_names(utility))) {
      candidate_set[, i] <- as.factor(candidate_set[, i])
    }

    # candidate_set <- as.matrix(candidate_set)

    cli_alert_success("All exclusions successfully applied")
  }

  # Prepare the list of priors ----
  cli_h2("Preparing the list of priors")

  design_object[["prior_values"]] <- prior_values <- prepare_priors(
    utility,
    draws,
    R
  )

  cli_alert_success("Priors prepared successfully")

  # Set up parallel ----
  if (control$cores > 1) {
    cli_h2("Preparing multicore estimation")

    stop("Multicore not implmented")

    future::plan(
      future::multicore(workers = control$cores)
    )

    cli_alert_success("Multicore estimation prepared successfully")
  }

  # Evaluate designs ----
  cli_h1("Evaluating designs")

  # Optmization function!!!!!!!
  design_object <- switch(
    algorithm,
    random = random(
      design_object,
      model,
      efficiency_criteria,
      utility,
      prior_values,
      dudx,
      candidate_set,
      rows,
      save_designs,
      control
    ),
    federov = federov(
      design_object,
      model,
      efficiency_criteria,
      utility,
      prior_values,
      dudx,
      candidate_set,
      rows,
      save_designs,
      control
    ),
    rsc = rsc(
      design_object,
      model,
      efficiency_criteria,
      utility,
      prior_values,
      dudx,
      candidate_set,
      rows,
      save_designs,
      control
    )
  )

  if (is.null(design_object[["design"]])) {
    stop(
      "No design with a non-singular Fisher information matrix was found in ",
      control$max_iter,
      " iterations. If exclusions or level occurrences are tight, increasing ",
      "'max_iter' may help. Otherwise, check the utility functions for ",
      "perfect multicollinearity."
    )
  }

  design_object[["time"]][["time_end"]] <- Sys.time()

  # Turn the design object into a tibble to be tidyverse compatible
  design_object[["design"]] <- tibble::as_tibble(design_object[["design"]])

  # Print final closing messages
  cat("\n\n")
  cli_h1("Cleaning up design environment")
  cat(
    "Time spent searching for designs: ",
    Sys.time() - design_object$time$time_start,
    "\n"
  )

  return(
    design_object
  )
}