skpr v1.8.1: Update moment matrix API for evaluation functions to require candidate set if user wants I optimality calculated. Use non-constrained by default.

-Fix skprGUI table color bug

Author: tylermorganwall

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: skpr
 Title: Design of Experiments Suite: Generate and Evaluate Optimal Designs
 Date: 2025-02-24
-Version: 1.8.0
+Version: 1.8.1
 Authors@R: c(person("Tyler", "Morgan-Wall", email = "[email protected]", role = c("aut", "cre")),
              person("George", "Khoury", email = "[email protected]", role = c("aut")))
 Description: Generates and evaluates D, I, A, Alias, E, T, and G optimal designs. Supports generation and evaluation of blocked and split/split-split/.../N-split plot designs. Includes parametric and Monte Carlo power evaluation functions, and supports calculating power for censored responses. Provides a framework to evaluate power using functions provided in other packages or written by the user. Includes a Shiny graphical user interface that displays the underlying code used to create and evaluate the design to improve ease-of-use and make analyses more reproducible. For details, see Morgan-Wall et al. (2021) <doi:10.18637/jss.v099.i01>.
@@ -26,11 +26,11 @@ Imports:
     lmerTest,
     methods,
     progress,
-    scales,
     doRNG,
     doFuture,
     progressr,
-    geometry
+    geometry,
+    digest
 LinkingTo: Rcpp, RcppEigen
 Suggests: testthat, mbest, ggplot2, lmtest, cli, gridExtra, rintrojs, shinythemes, shiny, shinyjs, gt, shinytest2
 Encoding: UTF-8

R/eval_design.R

@@ -37,6 +37,12 @@
 #'The reordering will be presenting in the output when `detailedoutput = TRUE`.
 #'@param advancedoptions Default `NULL`. A named list with parameters to specify additional attributes to calculate. Options: `aliaspower`
 #'gives the degree at which the Alias matrix should be calculated.
+#'@param candidate_set Default `NA`. If you generated your design externally from skpr and there are disallowed combinations in your design,
+#'the calculated I-optimality values will not be correct, as it assumes an unconstrained unit hypercube. To calculate the restricted regions,
+#'skpr will find the convex hull of the point set and generate a higher density of points in that region. Note that this only supports
+#'convex constraints.
+#'@param mm_sample_density Default `20`. The density of points to sample when calculating the moment matrix to compute I-optimality. Only
+#'required if the design was generated outside of skpr and there are disallowed combinations.
 #'@param ... Additional arguments.
 #'@return A data frame with the parameters of the model, the type of power analysis, and the power. Several
 #'design diagnostics are stored as attributes of the data frame. In particular,
@@ -173,6 +179,8 @@ eval_design = function(
   reorder_factors = FALSE,
   detailedoutput = FALSE,
   advancedoptions = NULL,
+  candidate_set = NA,
+  mm_sample_density = 20,
   ...
 ) {
   if (missing(design)) {
@@ -501,24 +509,30 @@ eval_design = function(
   if(!is.null(attr(input_design, "generating.model"))) {
     og_design_factors = attr(terms.formula(attr(input_design, "generating.model")),"factors")
     new_model_factors = attr(terms.formula(model),"factors")
-    identical_main_effects = all(rownames(og_design_factors) == rownames(new_model_factors))
-    identical_interactions = all(colnames(og_design_factors) == colnames(new_model_factors))
-    if(identical_interactions && identical_main_effects) {
-      mm = attr(input_design, "moments.matrix")
-      imported_mm = TRUE
+    if(all(dim(og_design_factors) == dim(new_model_factors))) {
+      identical_main_effects = all(rownames(og_design_factors) == rownames(new_model_factors))
+      identical_interactions = all(colnames(og_design_factors) == colnames(new_model_factors))
+      if(identical_interactions && identical_main_effects) {
+        mm = attr(input_design, "moments.matrix")
+        imported_mm = TRUE
+      }
     }
   }
   if(!imported_mm) {
-    if(all(classvector)) {
-      mm = gen_momentsmatrix(colnames(attr(run_matrix_processed, "modelmatrix")), levelvector, classvector)
-    } else {
-      if(!is.null(attr(design, "candidate_set"))) {
-        mm = gen_momentsmatrix_continuous(formula = model,
-                                          candidate_set = attr(input_design, "candidate_set"),
-                                          n_samples_per_dimension = 20)
+    mm = gen_momentsmatrix(colnames(attr(run_matrix_processed, "modelmatrix") ), levelvector, classvector)
+    if(!is.na(candidate_set)) {
+      if(!all(classvector)) {
+        hash_mm = digest::digest(list(model, candidate_set, mm_sample_density))
+        if(!exists(hash_mm, envir = skpr_moment_matrix_cache)) {
+          mm = gen_momentsmatrix_continuous(formula = model,
+                                            candidate_set = candidate_set,
+                                            n_samples_per_dimension = mm_sample_density)
+          assign(hash_mm, mm, envir = skpr_moment_matrix_cache)
+        } else {
+          mm = get(hash_mm,envir = skpr_moment_matrix_cache)
+        }
       } else {
-        warning("Candidate set not included with design: assuming no disallowed combinations when calculating moment matrix.")
-        mm = gen_momentsmatrix(colnames(attr(run_matrix_processed, "modelmatrix")), levelvector, classvector)
+        mm = gen_momentsmatrix(colnames(attr(run_matrix_processed, "modelmatrix") ), levelvector, classvector)
       }
     }
   }

R/eval_design_mc.R

@@ -53,6 +53,12 @@
 #'set the confidence level for power intervals, which are printed when `detailedoutput = TRUE`.
 #'@param parallel Default `FALSE`. If `TRUE`, the Monte Carlo power calculation will use all but one of the available cores. If the user wants to set the number of cores manually, they can do this by setting `options("cores")` to the desired number (e.g. `options("cores" = parallel::detectCores())`).
 #' NOTE: If you have installed BLAS libraries that include multicore support (e.g. Intel MKL that comes with Microsoft R Open), turning on parallel could result in reduced performance.
+#'@param candidate_set Default `NA`. If you generated your design externally from skpr and there are disallowed combinations in your design,
+#'the calculated I-optimality values will not be correct, as it assumes an unconstrained unit hypercube. To calculate the restricted regions,
+#'skpr will find the convex hull of the point set and generate a higher density of points in that region. Note that this only supports
+#'convex constraints.
+#'@param mm_sample_density Default `20`. The density of points to sample when calculating the moment matrix to compute I-optimality. Only
+#'required if the design was generated outside of skpr and there are disallowed combinations.
 #'@param ... Additional arguments.
 #'@return A data frame consisting of the parameters and their powers, with supplementary information
 #'stored in the data frame's attributes. The parameter estimates from the simulations are stored in the "estimates"
@@ -254,6 +260,8 @@ eval_design_mc = function(
   detailedoutput = FALSE,
   progress = TRUE,
   advancedoptions = NULL,
+  candidate_set = NA,
+  mm_sample_density = 20,
   ...
 ) {
   if (!firth || glmfamily != "binomial") {
@@ -1294,25 +1302,31 @@ eval_design_mc = function(
   imported_mm = FALSE
   if(!is.null(attr(design, "generating.model"))) {
     og_design_factors = attr(terms.formula(attr(design, "generating.model")),"factors")
-    new_model_factors = attr(terms.formula(model),"factors")
-    identical_main_effects = all(rownames(og_design_factors) == rownames(new_model_factors))
-    identical_interactions = all(colnames(og_design_factors) == colnames(new_model_factors))
-    if(identical_interactions && identical_main_effects) {
-      mm = attr(design, "moments.matrix")
-      imported_mm = TRUE
+    new_model_factors = attr(terms.formula(model_formula),"factors")
+    if(all(dim(og_design_factors) == dim(new_model_factors))) {
+      identical_main_effects = all(rownames(og_design_factors) == rownames(new_model_factors))
+      identical_interactions = all(colnames(og_design_factors) == colnames(new_model_factors))
+      if(identical_interactions && identical_main_effects) {
+        mm = attr(design, "moments.matrix")
+        imported_mm = TRUE
+      }
     }
   }
   if(!imported_mm) {
-    if(all(classvector)) {
-      mm = gen_momentsmatrix(colnames(attr(run_matrix_processed, "modelmatrix")), levelvector, classvector)
-    } else {
-      if(!is.null(attr(design, "candidate_set"))) {
-        mm = gen_momentsmatrix_continuous(formula = model,
-                                          candidate_set = attr(design, "candidate_set"),
-                                          n_samples_per_dimension = 20)
+    mm = gen_momentsmatrix(colnames(ModelMatrix), levelvector, classvector)
+    if(!is.na(candidate_set)) {
+      if(!all(classvector)) {
+        hash_mm = digest::digest(list(model_formula, candidate_set, mm_sample_density))
+        if(!exists(hash_mm, envir = skpr_moment_matrix_cache)) {
+          mm = gen_momentsmatrix_continuous(formula = model_formula,
+                                            candidate_set = candidate_set,
+                                            n_samples_per_dimension = mm_sample_density)
+          assign(hash_mm, mm, envir = skpr_moment_matrix_cache)
+        } else {
+          mm = get(hash_mm,envir = skpr_moment_matrix_cache)
+        }
       } else {
-        warning("Candidate set not included with design: assuming no disallowed combinations when calculating moment matrix.")
-        mm = gen_momentsmatrix(colnames(attr(run_matrix_processed, "modelmatrix")), levelvector, classvector)
+        mm = gen_momentsmatrix(colnames(ModelMatrix), levelvector, classvector)
       }
     }
   }

R/skprGUI.R

@@ -2727,10 +2727,7 @@ skprGUI = function(
       display_table = display_table %>%
         gt::data_color(
           columns = "power",
-          palette = scales::col_numeric(
-            palette = colorRampPalette(c("white", "darkgreen"))(100),
-            domain = c(0, 1)
-          ),
+          palette = colorRampPalette(c("white", "darkgreen"))(100),
           domain = c(0, 1),
           alpha = 0.3,
           autocolor_text = FALSE
@@ -2766,7 +2763,7 @@ skprGUI = function(
           ) %>%
           gt::tab_source_note(
             source_note = sprintf(
-              "Note: Power values marked in %s are within the simulation uncertainty for user-specified Type-I error (increase the number of simulations)",
+              "Note: Power values marked in %s are within the simulation uncertainty for user-specified minimum power: increase the number of simulations.",
               color_maybe
             )
           )
@@ -2784,9 +2781,9 @@ skprGUI = function(
           ) %>%
           gt::tab_source_note(
             source_note = sprintf(
-              "Note: Power values marked in %s fall below the user-specified Type-I error (%0.2f)",
+              "Note: Power values marked in %s fall below the user-specified minimum power (%0.2f)",
               color_bad,
-              alpha
+              desired_power_level
             )
           )
       }

R/zzz.R

@@ -46,3 +46,6 @@ assign(
   detect_multicore_support(),
   envir = skpr_system_setup_env
 )
+
+skpr_moment_matrix_cache = new.env(parent = emptyenv())
+

man/eval_design.Rd

@@ -4,7 +4,8 @@ library(shiny)
 test_that("{shinytest2} recording: apps", {
   options("in_skpr_test_environment" = TRUE)
   on.exit(options("in_skpr_test_environment" = NULL), add = TRUE)
-  app <- AppDriver$new(app_dir = testthat::test_path("apps"),
+  app <- AppDriver$new(app_dir = testthat::test_path("apps"), screenshot_args = FALSE,
+                       expect_values_screenshot_args  = FALSE,
     variant = platform_variant(), name = "apps", height = 923,
       width = 1619, timeout = 10*1000, wait = TRUE)
   app$set_inputs(setseed = TRUE)