skpr v1.9.1: Improve G-Efficiency calculation by explicitly sampling design space edges

Author: tylermorganwall

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: skpr
 Title: Design of Experiments Suite: Generate and Evaluate Optimal Designs
-Date: 2025-09-15
-Version: 1.9.0
+Date: 2025-09-19
+Version: 1.9.1
 Authors@R: c(person("Tyler", "Morgan-Wall", email = "tylermw@gmail.com", role = c("aut", "cre")),
              person("George", "Khoury", email = "george.m.khoury@gmail.com", 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>.

R/calculate_gefficiency.R

@@ -19,16 +19,16 @@ calculate_gefficiency = function(
   if (is.null(randsearches)) {
     randsearches = 10000
   }
-  variables = all.vars(get_attribute(design, "model"))
-  designmm = get_attribute(design, "model_matrix")
+  model = attr(design, "generating_model")
+  variables = all.vars(model)
+  designmm = attr(design, "model_matrix")
   modelentries = names(calculate_level_vector(
     run_matrix,
-    get_attribute(design, "model"),
+    model,
     FALSE
   ))
   model_entries_mul = gsub(":", "*", modelentries)
   model_entries_mul[1] = 1
-  fulllist = list()
   factorvars = attr(design, "contrastslist")
   variables = variables[!variables %in% names(factorvars)]
   rand_vector = function(factorlist, model_entries_mul) {
@@ -38,8 +38,16 @@ calculate_gefficiency = function(
     ))
     matrix(fulloutput, 1, length(fulloutput))
   }
-
+  invXtX = solve(t(designmm) %*% designmm)
   calculate_optimality_for_point = function(x, infval = FALSE) {
+    fulllist = list()
+    if (any(x > 1) || any(x < -1)) {
+      if (infval) {
+        return(Inf)
+      } else {
+        return(10000)
+      }
+    }
     for (i in seq_len(length(x))) {
       fulllist[[i]] = x[i]
     }
@@ -52,17 +60,9 @@ calculate_gefficiency = function(
     }
 
     mc_mm = rand_vector(fulllist, model_entries_mul)
-    if (any(x > 1) || any(x < -1)) {
-      if (infval) {
-        return(Inf)
-      } else {
-        return(10000)
-      }
-    }
-    100 *
-      (ncol(designmm)) /
-      (nrow(designmm) *
-        max(diag(mc_mm %*% solve(t(designmm) %*% designmm) %*% t(mc_mm))))
+    return(
+      as.numeric(mc_mm %*% invXtX %*% t(mc_mm))
+    )
   }
   modelentries = names(calculate_level_vector(
     run_matrix,
@@ -72,34 +72,43 @@ calculate_gefficiency = function(
   factorvars = attr(design, "contrastslist")
   variables = all.vars(get_attribute(design, "model"))
   variables = variables[!variables %in% names(factorvars)]
+
+  #First run edges
+  candset = attr(design, "candidate_set")
+  for (i in seq_len(ncol(candset))) {
+    if (is.numeric(candset[[i]])) {
+      num_col = candset[[i]]
+      min_max = range(candset)
+      is_exterior = num_col %in% min_max
+      candset = candset[is_exterior, ]
+    }
+  }
+  contrastslist = attr(design, "contrastslist")
+  candset_mm = model.matrix(model, data = candset, contrasts = contrastslist)
+  max_spv_verts = max(candset_mm %*% invXtX %*% t(candset_mm))
+
   if (calculation_type == "random") {
     vals = list()
-    lowest = 100
     for (i in seq_len(randsearches)) {
       vals[[i]] = calculate_optimality_for_point(
         2 * runif(length(variables)) - 1
       )
-      if (vals[[i]] < lowest) {
-        lowest = vals[[i]]
-      }
     }
-    return(min(unlist(vals)))
+    max_interior_val = max(unlist(vals))
+    max_val_overall = max(max_spv_verts, max_interior_val)
+    return(100 * ncol(designmm) / (nrow(designmm) * max_val_overall))
   }
 
   if (calculation_type == "optim") {
     vals = list()
-    lowest = 100
     for (i in seq_len(randsearches)) {
       vals[[i]] = optim(
         2 * runif(length(variables)) - 1,
         calculate_optimality_for_point,
         method = "SANN"
       )$value
-      if (vals[[i]] < lowest) {
-        lowest = vals[[i]]
-      }
     }
-    return(min(unlist(vals)))
+    return(100 * ncol(designmm) / (nrow(designmm) * max(unlist(vals))))
   }
   if (calculation_type == "custom" && !is.null(design_space_mm)) {
     return(GEfficiency(designmm, design_space_mm))

R/gen_design.R

@@ -1140,6 +1140,8 @@ gen_design = function(
     if (!is.null(custom_v)) {
       V = custom_v
     }
+  } else {
+    V = diag(trials) * varianceratio
   }
 
   initialreplace = FALSE
@@ -1975,6 +1977,8 @@ gen_design = function(
   attr(design, "generating.criterion") = optimality
   attr(design, "generating.contrast") = contrasts
   attr(design, "contrastslist") = contrastslist
+  attr(design, "variance.matrix") = V
+  attr(design, "candidate_set") = og_candidate_set
 
   if (!splitplot) {
     if (blocking) {
@@ -2097,7 +2101,6 @@ gen_design = function(
     error = function(e) {
     }
   )
-
   if (!splitplot && !blocking) {
     tryCatch(
       {
@@ -2197,6 +2200,9 @@ gen_design = function(
       }
     )
   }
+  if (is.null(attr(design, "I"))) {
+    attr(design, "I") = NA_real_
+  }
 
   #Re-order factors so levels with the lowest number of factors come first
   for (i in seq_len(ncol(design))) {
@@ -2327,7 +2333,6 @@ gen_design = function(
   } else {
     attr(design, "augmented") = FALSE
   }
-  attr(design, "candidate_set") = og_candidate_set
   modelmatrix_cor = model.matrix(
     model,
     design,

R/normalize_design.R

@@ -35,7 +35,7 @@ normalize_design = function(design, augmented = NULL) {
         "Design to be augmented and new design must have identical column classes"
       )
     }
-    for (column in 1:ncol(design)) {
+    for (column in seq_len(ncol(design))) {
       if (is.numeric(design[, column])) {
         midvalue = mean(c(
           max(c(design[, column], augmented[, column])),
@@ -46,7 +46,7 @@ normalize_design = function(design, augmented = NULL) {
       }
     }
   } else {
-    for (column in 1:ncol(design)) {
+    for (column in seq_len(ncol(design))) {
       if (is.numeric(design[, column])) {
         midvalue = mean(c(max(design[, column]), min(design[, column])))
         design[, column] = (design[, column] - midvalue) /

R/plot_fds.R

@@ -136,7 +136,7 @@ plot_fds = function(
       model = ~.
       model_matrix = model.matrix(
         model,
-        design,
+        data = design,
         contrasts.arg = temp_contrasts_list
       )
       factors = colnames(model_matrix)
@@ -145,13 +145,12 @@ plot_fds = function(
     new_model = TRUE
     model_matrix = model.matrix(
       model,
-      design,
+      data = design,
       contrasts.arg = temp_contrasts_list
     )
     factors = colnames(model_matrix)
   }
   #Need these inputs
-  model_matrix = attr(skpr_output, "model_matrix")
   if (!is.null(attr(skpr_output, "runmatrix"))) {
     design = attr(skpr_output, "runmatrix")
   } else {
@@ -187,6 +186,25 @@ plot_fds = function(
   }
   normalized_candidate_set = normalize_design(candidate_set)
 
+  candset_verts = candidate_set
+  for (i in seq_len(ncol(candset_verts))) {
+    if (is.numeric(candset_verts[[i]])) {
+      num_col = candset_verts[[i]]
+      min_max = range(candset_verts[[i]])
+      is_exterior = num_col %in% min_max
+      candset_verts = candset_verts[is_exterior, , drop = FALSE]
+    }
+  }
+  contrastslist = attr(design, "contrastslist")
+  candset_mm = model.matrix(
+    model,
+    data = candset_verts,
+    contrasts = contrastslist
+  )
+  invXtX = solve(t(model_matrix) %*% model_matrix)
+
+  spv_verts = diag(candset_mm %*% invXtX %*% t(candset_mm))
+
   #Need to generate this here
   Iopt = attr(skpr_output, "I")
   if (is.na(Iopt)) {
@@ -204,7 +222,7 @@ plot_fds = function(
     } else {
       model_matrix_cor = model.matrix(
         model,
-        design,
+        data = design,
         contrasts.arg = temp_contrasts_list
       )
     }
@@ -214,7 +232,7 @@ plot_fds = function(
       diag(nrow(model_matrix_cor))
     )
   }
-  if (is.null(Iopt)) {
+  if (is.na(Iopt) || is.null(Iopt)) {
     stop(
       "No I-optimality value found in design--was your design generated outside of skpr? If so, pass in a high resolution candidate set to `high_resolution_candidate_set` to ensure I-optimality is computed."
     )
@@ -238,7 +256,7 @@ plot_fds = function(
   }
   sample_list = list()
 
-  for (col in 1:ncol(skpr_output)) {
+  for (col in seq_len(ncol(skpr_output))) {
     if (inherits(skpr_output[, col], c("factor", "character"))) {
       vals = unique(skpr_output[, col])
     }
@@ -254,26 +272,15 @@ plot_fds = function(
   samples = as.data.frame(sample_list)
 
   #------Normalize/Center numeric columns ------#
-  for (column in 1:ncol(skpr_output)) {
-    if (is.numeric(skpr_output[, column])) {
-      midvalue = mean(c(max(skpr_output[, column]), min(skpr_output[, column])))
-      skpr_output[, column] = (skpr_output[, column] - midvalue) /
-        (max(skpr_output[, column]) - midvalue)
-    }
-  }
-  mm = model.matrix(model, skpr_output, contrasts.arg = contrastlist)
+  skpr_output_norm = normalize_design(skpr_output)
+  mm = model.matrix(model, skpr_output_norm, contrasts.arg = contrastlist)
   samplemm = model.matrix(model, samples, contrasts.arg = contrastlist)
 
   testcor = solve(t(mm) %*% solve(V) %*% mm)
 
-  v = list()
-
-  for (i in 1:nrow(samplemm)) {
-    xi = samplemm[i, ]
-    v[[i]] = t(xi) %*% testcor %*% xi
-  }
+  vals_interior = diag(samplemm %*% testcor %*% t(samplemm))
 
-  vars = do.call(rbind, v)
+  vars = c(spv_verts, vals_interior)
   varsordered = vars[order(vars)]
   meanindex = which(
     abs(mean(varsordered) - varsordered) ==