Merge pull request #21 from edsandorf/v0.0.2-dev

V0.0.2 dev

DESCRIPTION

@@ -1,12 +1,12 @@
 Package: spdesign
 Type: Package
 Title: Designing Stated Preference Experiments
-Version: 0.0.1
+Version: 0.0.2
 Authors@R: c(
   person("Erlend Dancke", "Sandorf", email = "erlend.dancke.sandorf@nmbu.no", role = c("aut", "cre")),
   person("Danny", "Campbell", email = "danny.campbell@stir.ac.uk", role = c("aut")))
 Maintainer: Erlend Dancke Sandorf <erlend.dancke.sandorf@nmbu.no>
-Description: Contemporary software commonly used to design stated preference experiments are expensive and the code is closed source. This is a free software package with an easy to use interface to make flexible stated preference experimental designs using state-of-the-art methods. For an overview of stated choice experimental design theory, see e.g., Rose, J. M. & Bliemer, M. C. J. (2014) in Hess S. & Daly. A. <doi:10.4337/9781781003152>. The package website can be accessed at <https://spdesign.edsandorf.me>.
+Description: Contemporary software commonly used to design stated preference experiments are expensive and the code is closed source. This is a free software package with an easy to use interface to make flexible stated preference experimental designs using state-of-the-art methods. For an overview of stated choice experimental design theory, see e.g., Rose, J. M. & Bliemer, M. C. J. (2014) in Hess S. & Daly. A. <doi:10.4337/9781781003152>. The package website can be accessed at <https://spdesign.edsandorf.me>. We acknowledge funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant INSPiRE (Grant agreement ID: 793163).
 License: CC BY-SA 4.0
 Encoding: UTF-8
 LazyData: true

NAMESPACE

@@ -16,6 +16,7 @@ export(generate_design)
 export(make_draws)
 export(occurrences)
 export(priors)
+export(probabilities)
 export(rep_cols)
 export(rep_rows)
 export(set_default_options)

NEWS.md

@@ -1,2 +1,9 @@
-## spdesignR v0.0.1
+# spdesign v0.0.2
+* New function ´probabilities()´ will now return the choice probabilities by choice task. 
+* Suppress warnings when calculating the correlation between the blocking column and the attributes to avoid warning when calculating correlation with respect to a constant.
+* After a number of candidates without improvement try a new design candidate when using the RSC algorithm
+* Updated package load message
+* Fixed roxygen @docType issue
+
+# spdesign v0.0.1
 * This is the first working version of the `spdesign` package that is able to create simple efficient designs for the MNL model. 

R/block.R

@@ -74,7 +74,10 @@ block <- function(x, blocks, target = 0.0005, max_iter = 1000000) {
     # Calculate the correlation between the attributes and a random permutation
     # of the blocking variable.
     block <- sample(block)
-    correlation <- stats::cor(design, block)
+
+    # Suppress warnings when trying to calculate the correlation with respect to
+    # a constant.
+    correlation <- suppressWarnings(stats::cor(design, block))
     current <- mean(correlation ^ 2, na.rm = TRUE)
 
     if (current < blocked_design[["blocks_value"]]) {

R/design.R

@@ -57,15 +57,31 @@ generate_design <- function(utility,
                               cores = 1,
                               max_iter = 10000,
                               max_relabel = 10000,
+                              max_no_improve = 100000,
                               efficiency_threshold = 0.1,
                               sample_with_replacement = 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()
+  )
+
+  # Make sure that the best design candidate is always return if the loop is
+  # stopped prematurely Can on.exit have a function?
+  on.exit(
+    return(design_object),
+    add = TRUE
+  )
   ## Match arguments ----
-  model <- match.arg(model)
+  design_object[["model"]] <- model <- match.arg(model)
   efficiency_criteria <- match.arg(efficiency_criteria, several.ok = TRUE)
   algorithm <- match.arg(algorithm)
   draws <- match.arg(draws)
@@ -85,6 +101,7 @@ generate_design <- function(utility,
     max_iter = 10000,
     max_relabel = 10000,
     max_swap = 10000,
+    max_no_improve = 100000,
     efficiency_threshold = 0.1,
     sample_with_replacement = FALSE
   )
@@ -163,7 +180,7 @@ generate_design <- function(utility,
   # Prepare the list of priors ----
   cli_h2("Preparing the list of priors")
 
-  prior_values <- prepare_priors(utility, draws, R)
+  design_object[["prior_values"]] <- prior_values <- prepare_priors(utility, draws, R)
 
   cli_alert_success("Priors prepared successfully")
 
@@ -184,24 +201,6 @@ generate_design <- function(utility,
   # Evaluate designs ----
   cli_h1("Evaluating designs")
 
-  # Create the design object and make sure that the current status of the object
-  # is returned if the program is ended prematurely from clicking "stop"
-  design_object <- list()
-  class(design_object) <- "spdesign"
-
-  design_object[["utility"]] <- clean_utility(utility)
-  design_object[["prior_values"]] <- prior_values
-  design_object[["time"]] <- list(
-    time_start = Sys.time()
-  )
-
-  # Make sure that the best design candidate is always return if the loop is
-  # stopped prematurely Can on.exit have a function?
-  on.exit(
-    return(design_object),
-    add = TRUE
-  )
-
   # Optmization function!!!!!!!
   design_object <- switch(
     algorithm,

R/methods.R

@@ -81,7 +81,7 @@ print.spdesign <- function(x, ...) {
   cat("An 'spdesign' object\n\n")
   cat("Utility functions:\n")
   for (i in seq_along(x$utility)) {
-    cat(names(x$utility)[[i]], ":", x$utility[[i]], "\n")
+    cat(names(x$utility)[[i]], ":", clean_utility(x$utility)[[i]], "\n")
   }
 
   cat("\n\n")

R/probabilities.R

@@ -0,0 +1,56 @@
+#' Calculate the probabilities of the design
+#'
+#' Will take the design object and calculate the probabilities of each
+#' alternative and choice tasks.
+#'
+#' Using Bayesian priors the average across the prior distribution will be used.
+#'
+#' Using the specific type of model, either the MNL or RPL probs will be
+#' returned.
+#'
+#' @param x An 'spdesign' object.
+#'
+#' @return A matrix of probabilities for each alternative and choice task.
+#'
+#' @export
+probabilities <- function(x) {
+  pr_j <- switch(
+    x$model,
+    mnl = probabilities_mnl(x)
+  )
+
+  return(
+    pr_j
+  )
+}
+
+#' Calculate the MNL probabilities
+#'
+#' @inheritParams probabilities
+#'
+#' @return A matrix of probabilities for each alternative and choice task. With
+#' Bayesian priors the return is the average probabilites over the prior
+#' distribution
+probabilities_mnl <- function(x) {
+  db_env <- new.env()
+  list2env(as.list(as.data.frame(x$design)), envir = db_env)
+
+  pr_j <- lapply(x$prior_values, function(p) {
+    list2env(as.list(p), envir = db_env)
+
+    obs_utility <- lapply(update_utility(x$utility), function(v) eval(parse(text = v), envir = db_env))
+    exp_utility <- lapply(obs_utility, exp)
+    sum_utility <- Reduce("+", exp_utility)
+    pr_j <- lapply(exp_utility, function(v) {
+      v <- v / sum_utility
+      v[is.na(v)] <- 0
+      as.vector(v)
+    })
+
+    return(do.call(cbind, pr_j))
+  })
+
+  return(
+    Reduce("+", pr_j) / length(pr_j)
+  )
+}

R/rsc.R

@@ -16,14 +16,6 @@ rsc <- function(design_object,
                 candidate_set,
                 rows,
                 control) {
-  # Create a level balanced design candidate or near level balanced candidate
-  design_candidate <- generate_rsc_candidate(utility, rows)
-
-  # 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(design_candidate) %in% dummy_names(utility))) {
-    design_candidate[, i] <- as.factor(design_candidate[, i])
-  }
 
   # Set up the design environment
   design_env <- new.env()
@@ -42,12 +34,31 @@ rsc <- function(design_object,
   # Set iteration defaults
   iter <- 1
   iter_na <- 1
+  iter_no_improve <- 1
   alg <- "relabel"
   efficiency_current_best <- NA
 
   repeat {
+    # Create an initial design candidate OR a new candidate when a sufficient
+    # number of attempts have been made without any improvement.
+    if (iter == 1 || iter_no_improve > control$max_no_improve) {
+      # Create a level balanced design candidate or near level balanced candidate
+      design_candidate <- generate_rsc_candidate(utility, rows)
+
+      # 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(design_candidate) %in% dummy_names(utility))) {
+        design_candidate[, i] <- as.factor(design_candidate[, i])
+      }
+
+      iter_no_improve <- 1
+
+      if (iter > 1) {
+        cli_alert_info(paste0("We've tried ", control$max_no_improve, " candidates without improvement. Trying new base design candidate."))
+      }
+    }
 
-    # Swith algorithm every 10 000 iterations
+    # Swith algorithm every max_relabel iterations
     if (iter %% control$max_relabel == 0) {
       alg <- ifelse(alg == "relabel", "swap", "relabel")
     }
@@ -112,6 +123,9 @@ rsc <- function(design_object,
       design_object[["vcov"]] <- efficiency_outputs[["vcov"]]
       efficiency_current_best <- efficiency_current
 
+      # Reset iter_no_improve when we have an improvement.
+      iter_no_improve <- 1
+
     }
 
     # Check stopping conditions ----
@@ -132,6 +146,9 @@ rsc <- function(design_object,
     # Add to the iteration
     iter <- iter + 1
 
+    # Add to the no improvement iterator. It's reset upon improvement
+    iter_no_improve <- iter_no_improve + 1
+
   }
 
   # Return the design candidate

R/spdesign-package.R

@@ -0,0 +1,10 @@
+#' @keywords internal
+"_PACKAGE"
+
+## usethis namespace: start
+#' @import cli
+#' @import stringr
+#' @importFrom stats runif qnorm as.formula model.matrix
+#' @importFrom utils modifyList
+## usethis namespace: end
+NULL

R/zzz.R

@@ -58,12 +58,12 @@
         paste0(
           col_green(
             str_pad(
-              "Access development version: ",
+              "Package website with documentation and examples: ",
               pad_width,
               "right",
               " ")
           ),
-          col_white("devtools::install_github('edsandorf/spdesign')")
+          col_white("https://spdesign.edsandorf.me")
         ),
         paste0(
           col_green(

README.md

@@ -20,4 +20,4 @@ A development version of the package can be installed from Github.
 All software contains bugs and we would very much like to find these and root them out. If you find a bug or get an error message, please reach out so that we can try and improve the software. 
 
 # Acknowledgements
-We are grateful to Petr Mariel, Jürgen Meyerhoff and Ainhoa Vega for providing feedback and extensive testing of the package. The package comes with no warranty and the authors cannot be held liable for erros or mistakes resulting from use. 
+We are grateful to Petr Mariel, Jürgen Meyerhoff and Ainhoa Vega for providing feedback and extensive testing of the package. The package comes with no warranty and the authors cannot be held liable for erros or mistakes resulting from use. The authors acknowledge funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant INSPiRE (Grant agreement ID: 793163).

examples/mnl-design.R

@@ -18,7 +18,8 @@ design <- generate_design(utility, rows = 20,
                           model = "mnl", efficiency_criteria = "d-error",
                           algorithm = "rsc", draws = "scrambled-sobol",
                           control = list(
-                            max_iter = 21000
+                            max_iter = 21000,
+                            max_no_improve = 5000
                           ))
 
 # Add a blocking variable to the design with 4 blocks.