edsandorf · edsandorf · Jan 16, 2024 · Sep 8, 2023 · Sep 8, 2023 · Sep 8, 2023
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -11,3 +11,5 @@
 ^\.github$
 ^test.R$
 ^.DS_Store$
+^README\.Rmd$
+^test-function-call.R
diff --git a/.github/workflows/render-rmarkdown.yaml b/.github/workflows/render-rmarkdown.yaml
@@ -0,0 +1,37 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    paths: ['README.Rmd']
+
+name: render-rmarkdown
+
+jobs:
+  render-rmarkdown:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    steps:
+      - name: Checkout repo
+        uses: actions/checkout@v3
+        with:
+          fetch-depth: 0
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+
+      - name: Set repo
+        run: |
+          options(repos = structure(c(CRAN = "https://cloud.r-project.org")))
+          install.packages('rmarkdown', repos = 'https://cloud.r-project.org')
+        shell: Rscript {0}
+
+      - name: Render Rmarkdown files and Commit Results
+        run: |
+          RMD_PATH=($(git diff --name-only ${{ github.event.before }} ${{ github.sha }} | grep '[.]Rmd$'))
+          Rscript -e 'for (f in commandArgs(TRUE)) if (file.exists(f)) rmarkdown::render(f)' ${RMD_PATH[*]}
+          git config --local user.name "$GITHUB_ACTOR"
+          git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com"
+          git commit ${RMD_PATH[*]/.Rmd/.md} -m 'Re-build Rmarkdown files' || echo "No changes to commit"
+          git push origin || echo "No changes to commit"
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: spdesign
 Type: Package
 Title: Designing Stated Preference Experiments
-Version: 0.0.2
+Version: 0.0.3.9004
 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")))
@@ -16,7 +16,9 @@ Imports:
   cli,
   future,
   randtoolbox,
-  matrixStats
+  matrixStats,
+  dplyr,
+  tibble
 Suggests: 
   knitr,
   rmarkdown,

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,10 @@
+# spdesign v0.0.3
+* Fixed a bug related to optimizing for c-efficiency where it would sometimes fail to correctly identify the denominator. 
+* Fixed several bugs related to using a supplied candidate set with alternative specific constants and attributes. Checks have been updated. The code will now also add zero-columns for alternative specific constants and attributes in the utility functions where they are not present. This ensures that all matrices used when calculating the first- and second-order derivatives of the utility functions are square. 
+* Fixed an issue where it failed to catch a mismatch in naming between the supplied candidate set and the utilty functions which caused hard to debug situations. Error messages should now catch this and provide additional information to help find the cause. Syntax is updated to reflect this as well. 
+* Fixed an issue where the full factorial would be generated even when the "rsc" algorithm was used, which caused memory issues for large designs. It is now only generated for the "random" and "federov" algorithms. A small section is added to the syntax vignette to clarify this. 
+* Minor bug fixes
+
 # 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.

diff --git a/R/block.R b/R/block.R
@@ -56,7 +56,7 @@
   block <- rep(seq_len(blocks), nrow(design) / blocks)
 
   blocked_design[["blocks_value"]] <- 1
-  blocked_design[["design"]] <- cbind(design, block)
+  blocked_design[["design"]] <- dplyr::bind_cols(design, block = block)
   blocked_design[["blocks_correlation"]] <- stats::cor(design, block)
   blocked_design[["blocks_iter"]] <- 1
 
@@ -82,8 +82,8 @@
 
     if (current < blocked_design[["blocks_value"]]) {
       blocked_design[["blocks_value"]] <- current
-      blocked_design[["design"]] <- cbind(design, block)
-      blocked_design[["blocks_correlation"]] <- correlation
+      blocked_design[["design"]] <- dplyr::bind_cols(design, block = block)
+      blocked_design[["blocks_correlation"]] <- tibble::as_tibble(t(correlation))
       blocked_design[["blocks_iter"]] <- iter
 
     }

diff --git a/R/design.R b/R/design.R
@@ -27,12 +27,14 @@
 #' 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 parameter in the
+#' @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 the full factorial subject to specified
-#' exclusions will be used.
+#' exclusions will be used. 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 control A list of control options
@@ -80,6 +82,7 @@
     return(design_object),
     add = TRUE
   )
+
   ## Match arguments ----
   design_object[["model"]] <- model <- match.arg(model)
   efficiency_criteria <- match.arg(efficiency_criteria, several.ok = TRUE)
@@ -131,51 +134,85 @@
   }
 
   ## Candidate set ----
-  cli_h2("Checking the candidate set and applying exclusions")
+  # 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 no candidate set is supplied generate full factorial if not run simple
-  # checks
-  if (is.null(candidate_set)) {
-    cli_alert_info("No candidate set supplied. The design will use the full factorial subject to supplied constraints.")
+  if (algorithm %in% c("random", "federov")) {
+    cli_h2("Checking the candidate set and applying exclusions")
 
-    candidate_set <- full_factorial(expand_attribute_levels(utility))
+    # If no candidate set is supplied generate full factorial if not run simple
+    # checks
+    if (is.null(candidate_set)) {
+      cli_alert_info("No candidate set supplied. The design will use the full factorial subject to supplied constraints.")
 
-    cli_alert_success("Full factorial created")
+      candidate_set <- full_factorial(expand_attribute_levels(utility))
 
-  } else {
-    stopifnot((is.matrix(candidate_set) || is.data.frame(candidate_set)))
+      cli_alert_success("Full factorial created")
 
-    if (!all(names(candidate_set) %in% names(expand_attribute_levels(utility)))) {
-      stop(
-        "Not all attributes specified in the utility functions are specified in
-        the candidate set. Make sure that all attributes are specified and that
-        the names used in the utility functions correspond to the column names
-        of the supplied candidate set. The candidate set must be supplied
-        in 'wide' format."
-      )
-    }
+    } 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 <- paste0("\\b", 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.")
+        )
+      }
+
+      candidate_levels <- apply(candidate_set, 2, function(x) unique(sort(x)))
+      utility_levels <- lapply(expand_attribute_levels(utility), as.numeric)
+
+      # Subset utility levels to only correspond to the ones specified
+      utility_levels <- utility_levels[utility_attributes]
+
+      if (!identical(candidate_levels[sort(names(candidate_levels))], utility_levels[sort(names(utility_levels))])) {
+        problem <- paste(names(which(mapply(function(x, y) length(x) - length(y), candidate_levels, utility_levels) != 0)), collapse = ", ")
+
+        stop(
+          paste0("The attribute levels determined by the supplied candidate set differs from those supplied in the utility function. Please ensure that all specified levels are present in the candidate set. The error occurs because there are too few/many levels for: ", problem, " in the candidate set")
+        )
+      }
+
+      # Expand candidate set to be square, i.e., fill in zero columns, for non-specified
+      expanded_names <- names(expand_attribute_levels(utility))
+      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]
 
-    if (!identical(apply(candidate_set, 2, function(x) unique(sort(x))), lapply(expand_attribute_levels(utility), as.numeric))) {
-      stop(
-        "The attribute levels determined by the supplied candidate set differs
-        from those supplied in the utility function. Please ensure that all
-        specified levels are present in the candidate set. "
-      )
     }
-  }
 
-  # Apply the exclusions to the candidate set
-  candidate_set <- exclude(candidate_set, exclusions)
+    # Apply the exclusions to the candidate set
+    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])
-  }
+    # 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)
+    # candidate_set <- as.matrix(candidate_set)
 
-  cli_alert_success("All exclusions successfully applied")
+    cli_alert_success("All exclusions successfully applied")
+  }
 
   # Prepare the list of priors ----
   cli_h2("Preparing the list of priors")
@@ -235,6 +272,9 @@
 
   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")

diff --git a/R/efficiency-criteria.R b/R/efficiency-criteria.R
@@ -144,15 +144,16 @@ calculate_c_error <- function(design_vcov, p, dudx, return_all) {
 
   } else {
     # Local overwrite with respect to the actual position for correct subsetting
-    dudx <- which(names(p) == dudx)
+    # dudx <- which(names(p) == dudx)
+    dudx <- which(str_detect(names(p), dudx) == TRUE)
 
     c_eff <- p[-dudx]^-2 * (diag(design_vcov)[dudx] - 2 * p[dudx] * p[-dudx]^-1 * design_vcov[dudx, seq_len(nrow(design_vcov))[-dudx]] + (p[dudx] / p[-dudx])^2 * diag(design_vcov)[-dudx])
 
     # Check if all are to be returned
     if (return_all) {
       c_eff
     } else {
-      sum(c_eff)
+      sum(c_eff, na.rm = TRUE)
     }
   }
 }

diff --git a/R/methods.R b/R/methods.R
@@ -34,13 +34,14 @@
 summary.spdesign <- function(object, ...) {
   print(object, ...)
   cat("\n")
-  print(object[["design"]])
+  cat("Printing the first few rows of the design \n")
+  print(utils::head(object[["design"]]))
   cat("\n")
   cat("---------------------------------------------------------------------\n")
 
   if ("blocks_correlation" %in% names(object)) {
     cat("Correlation between the blocking vector and attributes: \n \n")
-    print(round(t(object$blocks_correlation), 3))
+    print(object$blocks_correlation)
     cat("\n")
     cat("---------------------------------------------------------------------\n")
   }
@@ -100,7 +101,7 @@
 #' @export
 coef.spdesign <- function(object, ...) {
   return(
-    object[["priors"]]
+    object[["prior_values"]]
   )
 }
 

diff --git a/R/random.R b/R/random.R
@@ -148,7 +148,7 @@
     fits <- fits_lvl_occurrences(utility, design_candidate, rows)
 
     if (show_warning && difftime(Sys.time(), time_start, units = "secs") > 60) {
-      cli_alert_info("No design candidate has been found that can achieve attribute level balance or that satisfies the constraints. You may want to let go of attribute level balance, reduce the number of constraints or increase the size of your design.")
+      cli_alert_info("No design candidate has been found. This could be because you have place too tight constraints on the design or that all design candidates result in a singular Fisher matrix. A singular Fisher matrix can happen if you have perfect multicollinearity in your utility functions.")
       show_warning <- FALSE
     }
   }

diff --git a/README.Rmd b/README.Rmd
@@ -0,0 +1,76 @@
+---
+output: github_document
+---
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.path = "man/figures/README-",
+  out.width = "100%"
+)
+```
+
+# spdesign: Designing Stated Preference Experiments
+# spdesign <a href="https://spdesign.edsandorf.me"><img src="man/figures/logo.svg" align="right" height="139" alt="spdesign website" /></a>
+
+<!-- badges: start -->
+[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version-last-release/spdesign)](https://cran.r-project.org/package=spdesign)
+[![CRAN RStudio mirror downloads](http://cranlogs.r-pkg.org/badges/spdesign)](http://www.r-pkg.org/pkg/spdesign)
+[![R build status](https://github.com/edsandorf/spdesign/workflows/R-CMD-check/badge.svg)](https://github.com/edsandorf/spdesign/actions?workflow=R-CMD-check)
+[![Coverage Status](https://codecov.io/github/edsandorf/spdesign/coverage.svg?branch=master)](https://codecov.io/github/edsandorf/spdesign?branch=master)
+<!-- badges: end -->
+
+Contemporary software commonly used to design stated preference experiments are expensive and the code is closed source. `spdesign` is a free software package with an easy to use interface to make flexible stated preference experimental designs using state-of-the-art methods. 
+
+## Installation
+
+The package can be installed from CRAN. 
+
+```{r install, eval=FALSE}
+install.packages("spdesign")
+```
+
+A development version of the package can be installed from Github. Remember to select the most recent development version (check the available branches).
+
+```{r install_github, eval=FALSE}
+devtools::install_github("edsandorf/spdesign")
+```
+
+## Example
+
+This is a basic example which shows you how to solve a common problem:
+
+```{r example, eval=FALSE}
+library(spdesign)
+
+#' Specifying a utility function with 3 attributes and a constant for the
+#' SQ alternative. The design has 20 rows.
+utility <- list(
+  alt1 = "b_x1[0.1]  * x1[1:5] + b_x2[0.4] * x2[c(0, 1)] + b_x3[-0.2] * x3[seq(0, 1, 0.25)]",
+  alt2 = "b_x1       * x1      + b_x2      * x2          + b_x3       * x3",
+  alt3 = "b_sq[0.15] * sq[1]"
+)
+
+# Generate designs ----
+design <- generate_design(utility,
+                          rows = 20,
+                          model = "mnl", 
+                          efficiency_criteria = "d-error",
+                          algorithm = "rsc")
+
+# Add a blocking variable to the design with 4 blocks.
+design <- block(design, 2)
+
+
+summary(design)
+
+```
+
+## Bugs and error reporting
+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. We also thank participants in the 2022 Summer School "Valuing options of adaption to climate change using choice experiments" at the University of Cape Town for valuable feedback on a beta version of the package. The package comes with no warranty and the authors cannot be held liable for errors 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).