R/design-utils.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

`%||%` <- function(x, y) if (is.null(x) || length(x) == 0L) y else x

.r4vn_design_require <- function(pkg, feature = NULL) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    msg <- sprintf("Package '%s' is required", pkg)
    if (!is.null(feature)) msg <- paste0(msg, " for ", feature)
    stop(paste0(msg, ". Install it with install.packages(\"", pkg, "\")."), call. = FALSE)
  }
  invisible(TRUE)
}

.r4vn_design_pkg_version <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) return(NA_character_)
  as.character(utils::packageVersion(pkg))
}

.r4vn_design_r4vn_version <- function() {
  tryCatch(as.character(utils::packageVersion("R4VN")), error = function(e) "development")
}

.r4vn_design_id <- function(prefix = "DS") {
  # Do not consume R's RNG merely to create an audit identifier. This keeps
  # design/sample-size bookkeeping from changing later randomization streams.
  now <- Sys.time()
  micro <- as.integer((as.numeric(now) * 1e6) %% 1e8)
  paste0(prefix, "-", format(now, "%Y%m%d-%H%M%S"), "-", sprintf("%08d", micro), "-", Sys.getpid())
}

.r4vn_design_hash <- function(x) {
  f <- tempfile(fileext = ".rds")
  on.exit(unlink(f), add = TRUE)
  saveRDS(x, f, version = 3)
  unname(tools::md5sum(f))
}

.r4vn_design_parse_numlist <- function(x) {
  if (is.numeric(x)) return(as.numeric(x))
  if (is.null(x) || !nzchar(trimws(x))) return(numeric())
  x <- gsub(";", ",", x, fixed = TRUE)
  out <- suppressWarnings(as.numeric(trimws(unlist(strsplit(x, ",", fixed = TRUE)))))
  out[is.finite(out)]
}

.r4vn_design_parse_named_counts <- function(x) {
  if (is.null(x) || !nzchar(trimws(x))) return(numeric())
  parts <- trimws(unlist(strsplit(x, ",", fixed = TRUE)))
  ans <- numeric()
  for (z in parts) {
    kv <- trimws(unlist(strsplit(z, "=", fixed = TRUE)))
    if (length(kv) == 2L) {
      v <- suppressWarnings(as.numeric(kv[2]))
      if (is.finite(v)) ans[kv[1]] <- v
    }
  }
  ans
}

.r4vn_design_fmt <- function(x, digits = 3L) {
  ifelse(is.na(x), "NA", formatC(x, digits = digits, format = "fg", flag = "#"))
}

.r4vn_design_pct <- function(x, digits = 1L) {
  paste0(formatC(100 * x, digits = digits, format = "f"), "%")
}

.r4vn_design_safe_filename <- function(x) {
  x <- gsub("[^A-Za-z0-9_-]+", "_", x)
  x <- gsub("_+", "_", x)
  x <- gsub("^_|_$", "", x)
  if (!nzchar(x)) x <- "R4VN_design"
  x
}

.r4vn_design_html_escape <- function(x) {
  x <- gsub("&", "&amp;", x, fixed = TRUE)
  x <- gsub("<", "&lt;", x, fixed = TRUE)
  x <- gsub(">", "&gt;", x, fixed = TRUE)
  x <- gsub('"', "&quot;", x, fixed = TRUE)
  x
}

.r4vn_design_project_new <- function(name = "Untitled Study") {
  structure(list(
    format = "R4VN design project",
    format_version = 1L,
    project_id = .r4vn_design_id("PROJECT"),
    name = name,
    created = as.character(Sys.time()),
    modified = as.character(Sys.time()),
    r_version = R.version.string,
    r4vn_version = .r4vn_design_r4vn_version(),
    study_design = NULL,
    sample_size = NULL,
    sampling = NULL,
    randomization = NULL
  ), class = "r4vn_design_project")
}

.r4vn_design_project_save <- function(project, file) {
  stopifnot(inherits(project, "r4vn_design_project"))
  project$modified <- as.character(Sys.time())
  saveRDS(project, file, version = 3)
  invisible(file)
}

.r4vn_design_project_load <- function(file) {
  x <- readRDS(file)
  if (!inherits(x, "r4vn_design_project")) {
    stop("The selected file is not an R4VN design project.", call. = FALSE)
  }
  x
}

.r4vn_ss_input <- function(id, label, default, help, min = -Inf, max = Inf,
                           step = NULL, type = "number", choices = NULL,
                           suffix = NULL, example = NULL) {
  list(id = id, label = label, default = default, help = help, min = min,
       max = max, step = step, type = type, choices = choices, suffix = suffix,
       example = example)
}


.r4vn_ss_parameter_guidance <- function(z) {
  stopifnot(is.list(z), !is.null(z$id))
  id <- z$id
  base <- z$help %||% "Planning parameter used by this sample-size method."

  friendly <- switch(id,
    p = "Your best estimate of the proportion of the target population with the outcome or characteristic of interest. Enter it as a number from 0 to 1.",
    p0 = "The proportion assumed under the null hypothesis or reference value.",
    p1 = "The expected proportion in group 1 (often the treatment/exposed group, depending on the method).",
    p2 = "The expected proportion in group 2 (often the control/unexposed group, depending on the method).",
    prevalence = "The expected proportion of recruited participants who have the target condition or outcome.",
    sensitivity = "The proportion of truly positive participants expected to test positive.",
    specificity = "The proportion of truly negative participants expected to test negative.",
    d = "The largest absolute margin of error you are willing to accept, expressed in the same scale as the parameter. For a proportion, 0.05 means \u00b15 percentage points.",
    relative_precision = "The margin of error expressed as a fraction of the expected value. For example, 20% relative precision around p = 0.20 gives an absolute precision of 0.04.",
    alpha = "The Type I error probability. For a conventional two-sided 95% confidence level, use 0.05.",
    power = "The probability of detecting the prespecified effect if that effect is truly present. Common planning values are 0.80 or 0.90.",
    ratio = "The planned sample-size ratio n2/n1. Use 1 for equal group sizes, 2 when group 2 will have twice as many participants as group 1.",
    sd = "The expected standard deviation of the outcome. Use a value from a previous study, pilot data, or a clinically justified planning assumption.",
    sd_diff = "The expected standard deviation of within-person or paired differences.",
    delta = "The smallest or expected difference that the study is designed to detect, in the original outcome units.",
    expected_diff = "The expected difference between groups or treatments under the alternative hypothesis.",
    margin = "The prespecified non-inferiority or equivalence margin. Its sign and interpretation depend on the selected method; justify it clinically before calculation.",
    conf_width = "The desired full width of the confidence interval. A width of 0.10 means an interval approximately \u00b10.05 around the estimate when symmetric.",
    r = "The expected correlation coefficient in the target population.",
    r0 = "The correlation assumed under the null hypothesis.",
    r1 = "The expected correlation under the alternative hypothesis or in group 1.",
    r2 = "The expected correlation in group 2.",
    auc = "The expected area under the ROC curve. 0.50 represents chance discrimination; larger values indicate better discrimination.",
    items = "The number of items/questions contributing to the scale or reliability coefficient.",
    calpha = "The Cronbach's alpha expected for the scale in the target population.",
    expected_alpha = "The Cronbach's alpha expected under the alternative hypothesis.",
    null_alpha = "The lowest or null Cronbach's alpha regarded as unacceptable for the hypothesis test.",
    alpha1 = "The expected Cronbach's alpha in group 1.",
    alpha2 = "The expected Cronbach's alpha in group 2.",
    icc = "The anticipated intraclass correlation coefficient. For cluster designs it describes within-cluster similarity; for reliability studies it represents expected reliability.",
    kappa = "The anticipated kappa agreement coefficient beyond chance.",
    raters = "The number of raters or repeated measurements contributing to the reliability estimate.",
    cluster_size = "The average number of participants expected within each cluster.",
    clusters_per_arm = "The number of clusters planned in each trial arm.",
    base_n = "The sample size that would be required under simple random sampling before applying a clustering/design-effect inflation.",
    hr = "The hazard ratio the study is designed to detect. A value below 1 indicates a lower hazard in the numerator group.",
    event_fraction = "The expected proportion of recruited participants who will experience the event during follow-up.",
    event_rate = "The expected event rate per unit of follow-up time.",
    rate1 = "The expected incidence/event rate in group 1, using the same time unit as group 2.",
    rate2 = "The expected incidence/event rate in group 2, using the same time unit as group 1.",
    rate_ratio = "The expected ratio of two incidence/event rates.",
    parameters = "The number of candidate predictor parameters, including dummy-variable degrees of freedom and nonlinear terms\u2014not merely the number of named variables.",
    rsquared = "The anticipated model R-squared on the scale required by the selected method. Use prior evidence or pilot modelling when possible.",
    shrinkage = "The target global shrinkage factor for prediction-model development. Values close to 1 indicate less desired overfitting.",
    mean_outcome = "The expected population mean of the continuous outcome.",
    timepoint = "The prediction or analysis time horizon, using the same time unit as the event rate/follow-up inputs.",
    mean_followup = "The expected mean follow-up duration in the same time units used elsewhere in the calculation.",
    measurements = "The number of repeated measurements planned per participant.",
    groups = "The number of independent study groups or treatment arms.",
    df = "The model degrees of freedom used by the selected goodness-of-fit or SEM/CFA calculation.",
    rmsea0 = "The RMSEA representing close/good fit under the null planning hypothesis.",
    rmsea1 = "The RMSEA representing the alternative level of lack of fit the study should be able to detect.",
    eta2 = "The expected eta-squared effect size: the proportion of outcome variance explained by the grouping factor.",
    f = "Cohen's f effect size for ANOVA-type designs.",
    f2 = "Cohen's f-squared effect size for multiple regression.",
    variance = "The expected outcome variance used by the selected model.",
    subjects_per_item = "A planning ratio of participants per scale item. This is a rule-of-thumb input, not a universal factor-analysis requirement.",
    minimum_n = "A minimum total sample imposed as a planning floor when using a rule-based EFA calculation.",
    base
  )

  example <- z$example
  if (is.null(example) || !nzchar(as.character(example))) {
    example <- switch(id,
      p = "If previous studies suggest a prevalence of 20%, enter 0.20.",
      prevalence = "If about 30% of recruited participants are expected to have the condition, enter 0.30.",
      d = "For \u00b15 percentage points of precision around a prevalence, enter 0.05.",
      alpha = "Use 0.05 for a two-sided 95% confidence level.",
      power = "Enter 0.80 for 80% power or 0.90 for 90% power.",
      ratio = "Enter 1 for a 1:1 allocation or 2 for n2:n1 = 2:1.",
      sd = "If a previous study reported SD = 12.4, enter 12.4.",
      delta = "If a 5-point mean difference is clinically important, enter 5.",
      conf_width = "For a desired CI width of 0.10, enter 0.10.",
      items = "For a 12-item questionnaire, enter 12.",
      calpha = "If Cronbach's alpha is expected to be about 0.80, enter 0.80.",
      icc = "If prior work suggests ICC = 0.03 for clustered outcomes, enter 0.03.",
      parameters = "Ten binary predictors plus one 4-level categorical predictor (3 parameters) means at least 13 candidate parameters.",
      paste0("Example planning value in R4VN: ", format(z$default, trim = TRUE, scientific = FALSE), ".")
    )
  }

  paste0(friendly, " ", example)
}

.r4vn_ss_method <- function(id, category, name, description, engine, inputs,
                            formula, reference, explore = character(),
                            method_type = "FORMULA-BASED", notes = NULL,
                            package = NULL, package_function = NULL) {
  list(
    id = id, category = category, name = name, description = description,
    engine = engine, inputs = inputs, formula = formula, reference = reference,
    explore = explore, method_type = method_type, notes = notes,
    package = package, package_function = package_function
  )
}

.r4vn_ss_z <- function(alpha = 0.05, two_sided = TRUE) {
  stats::qnorm(1 - if (two_sided) alpha / 2 else alpha)
}

.r4vn_ss_zpower <- function(power = 0.80) stats::qnorm(power)

.r4vn_ss_validate_prob <- function(x, name) {
  if (!is.finite(x) || x <= 0 || x >= 1) {
    stop(sprintf("%s must be strictly between 0 and 1.", name), call. = FALSE)
  }
  invisible(TRUE)
}

.r4vn_ss_validate_pos <- function(x, name, allow_zero = FALSE) {
  ok <- is.finite(x) && if (allow_zero) x >= 0 else x > 0
  if (!ok) stop(sprintf("%s must be %s.", name, if (allow_zero) "non-negative" else "positive"), call. = FALSE)
  invisible(TRUE)
}

.r4vn_ss_extract_scalar <- function(x, names = c("n", "N", "sample_size", "samplesize", "final", "final_n", "minimum_sample_size")) {
  if (is.numeric(x) && length(x) == 1L && is.finite(x)) return(as.numeric(x))
  if (is.data.frame(x) || is.matrix(x)) {
    nms <- colnames(x)
    if (!is.null(nms)) {
      hit <- intersect(names, nms)
      if (length(hit)) {
        z <- suppressWarnings(as.numeric(x[, hit[1]][1]))
        if (is.finite(z)) return(z)
      }
    }
  }
  if (is.list(x)) {
    nms <- names(x)
    if (!is.null(nms)) {
      for (nm in names) {
        if (nm %in% nms) {
          z <- suppressWarnings(as.numeric(x[[nm]][1]))
          if (is.finite(z)) return(z)
        }
      }
      # pmsampsize has changed internal names across versions; inspect likely final fields.
      idx <- grep("final|minimum|sample.*size|^n$", nms, ignore.case = TRUE)
      for (i in idx) {
        z <- suppressWarnings(as.numeric(x[[i]][1]))
        if (is.finite(z)) return(z)
      }
    }
    for (el in x) {
      z <- tryCatch(.r4vn_ss_extract_scalar(el, names), error = function(e) NA_real_)
      if (is.finite(z)) return(z)
    }
  }
  NA_real_
}

.r4vn_ss_adjust <- function(result, loss = 0, design_effect = 1, population = NA_real_) {
  if (is.null(result$group_n)) result$group_n <- result$base_n
  g0 <- as.numeric(result$group_n)
  if (!length(g0) || any(!is.finite(g0)) || any(g0 <= 0)) stop("Invalid base sample size.", call. = FALSE)
  g <- g0
  adj_steps <- character()
  adj_formula <- character()
  adj_substitution <- character()

  if (is.finite(population) && population > 0 && length(g) == 1L) {
    before <- g
    g <- g / (1 + (g - 1) / population)
    result$fpc_applied <- TRUE
    adj_formula <- c(adj_formula, "n_{FPC}=\\frac{n_0}{1+(n_0-1)/N}")
    adj_substitution <- c(adj_substitution, sprintf("n_{FPC}=\\frac{%.4f}{1+(%.4f-1)/%.0f}=%.4f",before,before,population,g))
    adj_steps <- c(adj_steps, sprintf("Finite-population adjusted n = %.4f", g))
  } else {
    result$fpc_applied <- FALSE
  }

  if (!is.finite(design_effect) || design_effect < 1) stop("Design effect must be at least 1.", call. = FALSE)
  if(design_effect != 1) {
    before <- g
    g <- g * design_effect
    adj_formula <- c(adj_formula, "n_{DEFF}=n\\times DEFF")
    adj_substitution <- c(adj_substitution, sprintf("n_{DEFF}=(%s)\\times %.4f=(%s)",paste(.r4vn_design_fmt(before,4),collapse=","),design_effect,paste(.r4vn_design_fmt(g,4),collapse=",")))
    adj_steps <- c(adj_steps, sprintf("After design effect %.4f: %s",design_effect,paste(.r4vn_design_fmt(g,4),collapse=" + ")))
  }

  if (!is.finite(loss) || loss < 0 || loss >= 1) stop("Loss/non-response must be between 0 and <1.", call. = FALSE)
  if(loss > 0) {
    before <- g
    g_unrounded <- g / (1 - loss)
    adj_formula <- c(adj_formula, "n_{final}=\\frac{n_{required}}{1-L}")
    adj_substitution <- c(adj_substitution, sprintf("n_{final}=\\frac{%s}{1-%.4f}=%s",paste(.r4vn_design_fmt(before,4),collapse=","),loss,paste(.r4vn_design_fmt(g_unrounded,4),collapse=",")))
    adj_steps <- c(adj_steps, sprintf("After %.1f%% loss/non-response: %s",100*loss,paste(.r4vn_design_fmt(g_unrounded,4),collapse=" + ")))
  } else g_unrounded <- g

  g_final <- ceiling(g_unrounded)
  adj_steps <- c(adj_steps, sprintf("Final rounded sample size%s = %s (total %d)",if(length(g_final)>1) " by group" else "",paste(g_final,collapse=" + "),sum(g_final)))

  result$base_group_n_unrounded <- g0
  result$adjusted_group_n_unrounded <- g_unrounded
  result$final_group_n <- g_final
  result$final_n <- sum(g_final)
  result$design_effect <- design_effect
  result$loss <- loss
  result$population <- population
  result$adjustment_formula <- adj_formula
  result$adjustment_substitution <- adj_substitution
  result$adjustment_steps <- adj_steps
  result
}

.r4vn_ss_result <- function(base_n, group_n = NULL, formula, substitution,
                            steps, details = NULL, raw = NULL) {
  if (is.null(group_n)) group_n <- base_n
  list(
    base_n = as.numeric(base_n), group_n = as.numeric(group_n), formula = formula,
    substitution = substitution, steps = steps, details = details, raw = raw
  )
}