skpr v1.7.1: Add error bars to Monte Carlo eval output if detailedoutput = TRUE

Author: tylermorganwall

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: skpr
 Title: Design of Experiments Suite: Generate and Evaluate Optimal Designs
-Date: 2024-03-18
-Version: 1.7.0
+Date: 2024-03-25
+Version: 1.7.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>.

R/add_ci_bounds_mc_power.R

@@ -0,0 +1,18 @@
+#'@title Calculate CI bounds on Monte Carlo
+#'
+#'@description Calculates CI bounds for Monte Carlo power
+#'@return Power data.frame with conf intervals
+#'@keywords internal
+add_ci_bounds_mc_power = function(power_results, nsim, conf = 0.95) {
+  get_lcbs = function(power) {
+    return(unlist(lapply(power, \(x) (binom.test(x*nsim, n = nsim, conf.level = conf)[["conf.int"]][1]))))
+  }
+  get_ucbs = function(power) {
+    return(unlist(lapply(power, \(x) (binom.test(x*nsim, n = nsim, conf.level = conf)[["conf.int"]][2]))))
+  }
+  lcb = get_lcbs(power_results[["power"]])
+  ucb = get_ucbs(power_results[["power"]])
+  power_results$power_lcb = lcb
+  power_results$power_ucb = ucb
+  return(power_results)
+}

R/eval_design_custom_mc.R

@@ -25,10 +25,12 @@
 #'fit returned by `fitfunction` must have a method compatable with the car package.
 #'@param parameternames Vector of parameter names if the coefficients do not correspond simply to the columns in the model matrix
 #'(e.g. coefficients from an MLE fit).
+#'@param detailedoutput Default `FALSE`. If `TRUE`, return additional information about evaluation in results.
 #'@param advancedoptions Default `NULL`. Named list of advanced options. `advancedoptions$anovatype` specifies the Anova type in the car package (default type `III`),
 #'user can change to type `II`). `advancedoptions$anovatest` specifies the test statistic if the user does not want a `Wald` test--other options are likelyhood-ratio `LR` and F-test `F`.
 #'`advancedoptions$progressBarUpdater` is a function called in non-parallel simulations that can be used to update external progress bar.`advancedoptions$GUI` turns off some warning messages when in the GUI.
-#'If `advancedoptions$save_simulated_responses = TRUE`, the dataframe will have an attribute `simulated_responses` that contains the simulated responses from the power evaluation.
+#'If `advancedoptions$save_simulated_responses = TRUE`, the dataframe will have an attribute `simulated_responses` that contains the simulated responses from the power evaluation. `advancedoptions$ci_error_conf` will
+#'set the confidence level for power intervals, which are printed when `detailedoutput = TRUE`.
 #'@param anticoef The anticipated coefficients for calculating the power. If missing, coefficients will be
 #'automatically generated based on \code{effectsize}.
 #'@param effectsize The signal-to-noise ratio. Default 2. For a gaussian model, and for
@@ -102,7 +104,7 @@
 eval_design_custom_mc = function(design, model = NULL, alpha = 0.05,
                                  nsim, rfunction, fitfunction, pvalfunction,
                                  anticoef, effectsize = 2, contrasts = contr.sum,
-                                 coef_function = coef, calceffect = FALSE,
+                                 coef_function = coef, calceffect = FALSE, detailedoutput = FALSE,
                                  parameternames = NULL, advancedoptions = NULL, progress = TRUE,
                                  parallel = FALSE, parallel_pkgs = NULL, ...) {
   if (!is.null(advancedoptions)) {
@@ -135,6 +137,9 @@ eval_design_custom_mc = function(design, model = NULL, alpha = 0.05,
     progress = getOption("skpr_progress")
   }
 
+  if(is.null(advancedoptions$ci_error_conf)) {
+    advancedoptions$ci_error_conf = 0.95
+  }
   if (!is.null(advancedoptions$anovatype)) {
     anovatype = advancedoptions$anovatype
   } else {
@@ -351,6 +356,19 @@ eval_design_custom_mc = function(design, model = NULL, alpha = 0.05,
   attr(power_final, "runmatrix") = RunMatrixReduced
   attr(power_final, "anticoef") = anticoef
 
+  if (detailedoutput) {
+    if (nrow(power_final) != length(anticoef)){
+      power_final$anticoef = c(rep(NA, nrow(power_final) - length(anticoef)), anticoef)
+    } else {
+      power_final$anticoef = anticoef
+    }
+    power_final$alpha = alpha
+    power_final$trials = nrow(run_matrix_processed)
+    power_final$nsim = nsim
+    power_final = add_ci_bounds_mc_power(power_final, nsim = nsim, conf =  advancedoptions$ci_error_conf)
+    attr(power_final, "mc.conf.int") = advancedoptions$ci_error_conf
+  }
+
   if(!inherits(power_final,"skpr_eval_output")) {
     class(power_final) = c("skpr_eval_output", class(power_final))
   }

R/eval_design_mc.R

@@ -49,7 +49,8 @@
 #'@param advancedoptions Default `NULL`. Named list of advanced options. `advancedoptions$anovatype` specifies the Anova type in the car package (default type `III`),
 #'user can change to type `II`). `advancedoptions$anovatest` specifies the test statistic if the user does not want a `Wald` test--other options are likelyhood-ratio `LR` and F-test `F`.
 #'`advancedoptions$progressBarUpdater` is a function called in non-parallel simulations that can be used to update external progress bar.`advancedoptions$GUI` turns off some warning messages when in the GUI.
-#'If `advancedoptions$save_simulated_responses = TRUE`, the dataframe will have an attribute `simulated_responses` that contains the simulated responses from the power evaluation.
+#'If `advancedoptions$save_simulated_responses = TRUE`, the dataframe will have an attribute `simulated_responses` that contains the simulated responses from the power evaluation. `advancedoptions$ci_error_conf` will
+#'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 ... Additional arguments.
@@ -322,6 +323,9 @@ eval_design_mc = function(design, model = NULL, alpha = 0.05,
     }
     aliaspower = advancedoptions$aliaspower
   }
+  if(is.null(advancedoptions$ci_error_conf)) {
+    advancedoptions$ci_error_conf = 0.95
+  }
   alpha_adjust = FALSE
   if (adjust_alpha_inflation) {
     alpha_adjust = TRUE
@@ -966,6 +970,8 @@ eval_design_mc = function(design, model = NULL, alpha = 0.05,
     } else {
       retval$error_adjusted_alpha = alpha_parameter
     }
+    retval = add_ci_bounds_mc_power(retval, nsim = nsim, conf =  advancedoptions$ci_error_conf)
+    attr(retval, "mc.conf.int") = advancedoptions$ci_error_conf
   }
 
   colnames(estimates) = parameter_names

R/eval_design_survival_mc.R

@@ -34,7 +34,8 @@
 #' 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 detailedoutput Default `FALSE`. If `TRUE`, return additional information about evaluation in results.
 #'@param progress Default `TRUE`. Whether to include a progress bar.
-#'@param advancedoptions Default NULL. Named list of advanced options. Pass `progressBarUpdater` to include function called in non-parallel simulations that can be used to update external progress bar.
+#'@param advancedoptions Default `NULL`. Named list of advanced options. Pass `progressBarUpdater` to include function called in non-parallel simulations that can be used to update external progress bar.
+#'`advancedoptions$ci_error_conf` will set the confidence level for power intervals, which are printed when `detailedoutput = TRUE`.
 #'@param ... Any additional arguments to be passed into the \code{survreg} function during fitting.
 #'@return A data frame consisting of the parameters and their powers. The parameter estimates from the simulations are
 #'stored in the 'estimates' attribute. The 'modelmatrix' attribute contains the model matrix and the encoding used for
@@ -161,6 +162,9 @@ eval_design_survival_mc = function(design, model = NULL, alpha = 0.05,
     advancedoptions$GUI = FALSE
     progressBarUpdater = NULL
   }
+  if(is.null(advancedoptions$ci_error_conf)) {
+    advancedoptions$ci_error_conf = 0.95
+  }
   if (attr(terms.formula(model, data = design), "intercept") == 1) {
     nointercept = FALSE
   } else {
@@ -304,7 +308,7 @@ eval_design_survival_mc = function(design, model = NULL, alpha = 0.05,
       )
 
       #determine whether beta[i] is significant. If so, increment nsignificant
-      pvals = extractPvalues(fit)[1:ncol(ModelMatrix)]
+      pvals = extractPvalues(fit)[seq_len(ncol(ModelMatrix))]
       pvals = pvals[order(factor(names(pvals), levels = parameter_names))]
       pvals[is.na(pvals)] = 1
       stopifnot(all(names(pvals) == parameter_names))
@@ -360,7 +364,7 @@ eval_design_survival_mc = function(design, model = NULL, alpha = 0.05,
         )
 
         #determine whether beta[i] is significant. If so, increment nsignificant
-        pvals = extractPvalues(fit)[1:ncol(ModelMatrix)]
+        pvals = extractPvalues(fit)[seq_len(ncol(ModelMatrix))]
         pvals = pvals[order(factor(names(pvals), levels = parameter_names))]
         stopifnot(all(names(pvals) == parameter_names))
         pvals[is.na(pvals)] = 1
@@ -386,13 +390,17 @@ eval_design_survival_mc = function(design, model = NULL, alpha = 0.05,
   attr(retval, "alpha") = alpha
   attr(retval, "runmatrix") = RunMatrixReduced
 
-
   if (detailedoutput) {
-    retval$anticoef = anticoef
+    if (nrow(retval) != length(anticoef)){
+      retval$anticoef = c(rep(NA, nrow(retval) - length(anticoef)), anticoef)
+    } else {
+      retval$anticoef = anticoef
+    }
     retval$alpha = alpha
-    retval$distribution = distribution
     retval$trials = nrow(run_matrix_processed)
     retval$nsim = nsim
+    retval = add_ci_bounds_mc_power(retval, nsim = nsim, conf =  advancedoptions$ci_error_conf)
+    attr(retval, "mc.conf.int") = advancedoptions$ci_error_conf
   }
   if(!inherits(retval,"skpr_eval_output")) {
     class(retval) = c("skpr_eval_output", class(retval))

R/methods.R

@@ -80,13 +80,13 @@ print.skpr_eval_output = function(x, ...) {
   }
   if (!is.null(attr(x,"z.matrix.list")) && blocking) {
     number_blocks = unlist(lapply(attr(x,"z.matrix.list"),ncol))
-    block_str = paste(paste("Level ", 1:length(number_blocks), ": ", number_blocks, sep = ""),
+    block_str = paste(paste("Level ", seq_len(length(number_blocks)), ": ", number_blocks, sep = ""),
                       collapse=", ")
     cat(generate_text("Number of Blocks",block_str), sep = "\n")
   }
   if (!is.null(attr(x,"varianceratios")) && blocking) {
     vr = attr(x,"varianceratios")
-    vr_str =  paste(paste("Level ", 1:length(vr), ": ",as.character(vr), sep="", collapse=", "))
+    vr_str =  paste(paste("Level ", seq_len(length(vr)), ": ",as.character(vr), sep="", collapse=", "))
     cat(generate_text("Variance Ratios ",vr_str), sep = "\n")
   }
   if(!is.null(attr(x, "contrast_string"))) {
@@ -102,6 +102,9 @@ print.skpr_eval_output = function(x, ...) {
       cat(generate_text("Effect Analysis Method", attr(x, "effect_analysis_method_string")),sep="\n")
     }
   }
+  if(!is.null(attr(x, "mc.conf.int"))) {
+    cat(generate_text("MC Power CI Confidence", sprintf("%0.0f%%",100*attr(x, "mc.conf.int"))),sep="\n")
+  }
 }
 
 #'@title Print evaluation information

R/skprGUI.R

@@ -2428,7 +2428,7 @@ skprGUI = function(browser = FALSE, return_app = FALSE, multiuser = FALSE, progr
 
     filter_power_results = function(results) {
       col_results = colnames(results)
-      results[, !col_results %in% c("glmfamily",	"trials",	"nsim",	"blocking")]
+      results[, !col_results %in% c("glmfamily",	"trials",	"nsim",	"blocking" ,"power_ucb", "power_lcb")]
     }
 
     output$powerresults = gt::render_gt( {