R/statistical.R

In simulab: Simulates a Wide Variety of Data Shapes, Sizes and Distributions

#' Simulate a two-group design
#'
#' @param n_a,n_b Group sample sizes, each defaulting to `50`.
#' @param mean_a,mean_b Group means. A mean left `NULL` (the default) is drawn
#'   from the standard normal distribution and rounded to two decimals, inside
#'   the same seeded draw as the data, and is reported in the `parameters`
#'   table.
#' @param sd_a,sd_b Positive group standard deviations.
#' @param labels Two group labels.
#' @param outcome Name of the outcome variable.
#' @param seed Optional random seed.
#'
#' @return A `simulab_sim` base `data.frame` with one row per observation and
#'   columns `id`, `group` and `outcome` (renamed by `outcome`). Two tidy tables
#'   come from `as.data.frame()`: `what = "parameters"` has one row per group
#'   with columns `group`, `n`, `mean` and `sd`; `what = "effects"` has a single
#'   row with columns `contrast`, `mean_difference`, `pooled_sd` and `cohens_d`,
#'   the population Cohen's d formed from the sample-size-weighted pooled
#'   standard deviation.
#' @inheritParams simulate_correlation
#' @export
#'
#' @examples
#' result <- simulate_ttest(n_a = 40, n_b = 40, mean_a = 0, mean_b = 0.6, seed = 1)
#' head(result)
#' as.data.frame(result, what = "parameters")
#'
#' # With no arguments, the group means are drawn and reported.
#' drawn <- simulate_ttest(seed = 1)
#' as.data.frame(drawn, what = "parameters")
simulate_ttest <- function(n_a = 50L, n_b = 50L, mean_a = NULL, mean_b = NULL,
                           sd_a = 1, sd_b = 1,
                           labels = c("A", "B"), outcome = "outcome",
                           seed = NULL,
                           batch = NULL) {
  if (!is.null(batch)) return(.simulate_batch(batch, seed))
  stopifnot(
    "`n_a` must be a single whole number of at least 2" =
      is.numeric(n_a) &&
        length(n_a) == 1L &&
        all(n_a >= 2) &&
        all(n_a == as.integer(n_a)),
    "`n_b` must be a single whole number of at least 2" =
      is.numeric(n_b) &&
        length(n_b) == 1L &&
        all(n_b >= 2) &&
        all(n_b == as.integer(n_b)),
    "`mean_a` must be NULL or a single finite number" =
      is.null(mean_a) ||
        (is.numeric(mean_a) && length(mean_a) == 1L && all(is.finite(mean_a))),
    "`mean_b` must be NULL or a single finite number" =
      is.null(mean_b) ||
        (is.numeric(mean_b) && length(mean_b) == 1L && all(is.finite(mean_b))),
    "`sd_a` must be a single positive number" =
      is.numeric(sd_a) &&
        length(sd_a) == 1L &&
        all(sd_a > 0),
    "`sd_b` must be a single positive number" =
      is.numeric(sd_b) &&
        length(sd_b) == 1L &&
        all(sd_b > 0),
    "`labels` must be an unique atomic vector, of length 2" =
      is.atomic(labels) &&
        length(labels) == 2L &&
        !anyDuplicated(labels),
    "`outcome` must be a single non-empty string" =
      is.character(outcome) &&
        length(outcome) == 1L &&
        all(nzchar(outcome)),
    "`seed` must be NULL or a single number" =
      is.null(seed) || (is.numeric(seed) && length(seed) == 1L)
  )
  generated <- .with_seed(seed, {
    mean_a <- mean_a %||% .draw_means(1L)
    mean_b <- mean_b %||% .draw_means(1L)
    list(
      mean_a = mean_a,
      mean_b = mean_b,
      values = c(
        stats::rnorm(as.integer(n_a), mean_a, sd_a),
        stats::rnorm(as.integer(n_b), mean_b, sd_b)
      )
    )
  })
  mean_a <- generated$mean_a
  mean_b <- generated$mean_b
  values <- generated$values
  data <- data.frame(
    id = seq_along(values),
    group = rep(labels, times = c(as.integer(n_a), as.integer(n_b))),
    value = values,
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  names(data)[3L] <- outcome
  pooled_sd <- sqrt(
    ((n_a - 1) * sd_a^2 + (n_b - 1) * sd_b^2) /
      (n_a + n_b - 2)
  )
  parameters <- data.frame(
    group = labels,
    n = c(as.integer(n_a), as.integer(n_b)),
    mean = c(mean_a, mean_b),
    sd = c(sd_a, sd_b),
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  effects <- data.frame(
    contrast = sprintf("%s - %s", labels[2L], labels[1L]),
    mean_difference = mean_b - mean_a,
    pooled_sd = pooled_sd,
    cohens_d = (mean_b - mean_a) / pooled_sd,
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  .new_simulab_sim(data, type = "ttest", seed = seed,
                   tables = list(parameters = parameters, effects = effects))
}

#' Simulate a one-way group design
#'
#' @param n Group sample size or one size per group, defaulting to `30`.
#' @param means Group means. `NULL` (the default) draws one mean per group from
#'   the standard normal distribution, rounded to two decimals, inside the
#'   same seeded draw as the data. The number of groups is then the length of
#'   `n`, `sds` or `labels`, or `3` when all of them have length one.
#' @param sds Group standard deviations.
#' @param labels Optional group labels.
#' @param outcome Name of the outcome variable.
#' @param seed Optional random seed.
#'
#' @return A `simulab_sim` base `data.frame` with one row per observation and
#'   columns `id`, `group` and `outcome` (renamed by `outcome`). Two tidy tables
#'   come from `as.data.frame()`: `what = "parameters"` has one row per group
#'   with columns `group`, `n`, `mean` and `sd`; `what = "effects"` has a single
#'   row with columns `grand_mean`, `between_sum_squares`
#'   (`sum(n * (mean - grand_mean)^2)`), `within_sum_squares`
#'   (`sum((n - 1) * sd^2)`) and `eta_squared`, their ratio
#'   `between / (between + within)`. Because the within term uses the
#'   degrees of freedom `n - 1` rather than `n`, `eta_squared` is the
#'   design's sum-of-squares eta-squared and differs slightly from the
#'   variance-ratio population eta-squared.
#' @inheritParams simulate_correlation
#' @export
#'
#' @examples
#' result <- simulate_anova(n = 90, means = c(0, 0.4, 0.9), seed = 1)
#' head(result)
#' as.data.frame(result, what = "parameters")
simulate_anova <- function(n = 30L, means = NULL, sds = 1, labels = NULL,
                           outcome = "outcome", seed = NULL,
                           batch = NULL) {
  if (!is.null(batch)) return(.simulate_batch(batch, seed))
  stopifnot(
    "`n` must be a numeric vector, of whole numbers, with at least one element, each at least 2" =
      is.numeric(n) &&
        length(n) >= 1L &&
        all(n >= 2) &&
        all(n == as.integer(n)),
    "`means` must be NULL or a finite numeric vector with at least two elements" =
      is.null(means) ||
        (is.numeric(means) && length(means) >= 2L && all(is.finite(means))),
    "`sds` must be a finite positive numeric vector, with at least one element" =
      is.numeric(sds) &&
        length(sds) >= 1L &&
        all(is.finite(sds)) &&
        all(sds > 0),
    "`labels` must be NULL or an atomic vector" =
      is.null(labels) || is.atomic(labels),
    "`outcome` must be a single non-empty string" =
      is.character(outcome) &&
        length(outcome) == 1L &&
        all(nzchar(outcome)),
    "`seed` must be NULL or a single number" =
      is.null(seed) || (is.numeric(seed) && length(seed) == 1L)
  )
  groups <- if (is.null(means)) {
    sizes <- c(length(n), length(sds), length(labels))
    if (any(sizes > 1L)) max(sizes) else 3L
  } else {
    length(means)
  }
  if (groups < 2L) {
    stop("A one-way design needs at least two groups.", call. = FALSE)
  }
  if (length(n) == 1L) n <- rep(n, groups)
  if (length(sds) == 1L) sds <- rep(sds, groups)
  if (length(n) != groups || length(sds) != groups) {
    stop("n, means, and sds must describe the same number of groups.", call. = FALSE)
  }
  if (is.null(labels)) labels <- sprintf("Group %d", seq_len(groups))
  if (length(labels) != groups || anyDuplicated(labels)) {
    stop("labels must contain one unique value per group.", call. = FALSE)
  }
  generated <- .with_seed(seed, {
    means <- means %||% .draw_means(groups)
    list(means = means, values = unlist(Map(
      function(size, mean, sd) stats::rnorm(as.integer(size), mean, sd),
      n, means, sds
    ), use.names = FALSE))
  })
  means <- generated$means
  values <- generated$values
  data <- data.frame(
    id = seq_along(values),
    group = rep(labels, times = as.integer(n)),
    value = values,
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  names(data)[3L] <- outcome
  grand_mean <- sum(n * means) / sum(n)
  between <- sum(n * (means - grand_mean)^2)
  within <- sum((n - 1) * sds^2)
  parameters <- data.frame(
    group = labels,
    n = as.integer(n),
    mean = means,
    sd = sds,
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  effects <- data.frame(
    grand_mean = grand_mean,
    between_sum_squares = between,
    within_sum_squares = within,
    eta_squared = between / (between + within),
    row.names = NULL
  )
  .new_simulab_sim(data, type = "anova", seed = seed,
                   tables = list(parameters = parameters, effects = effects))
}

#' Simulate a linear regression design
#'
#' @param n Sample size, defaulting to `100`.
#' @param coefficients Named coefficients including optional `(Intercept)`.
#'   `NULL` (the default) draws an intercept from the standard normal
#'   distribution and `n_predictors` slopes, named `x1`, `x2`, ..., from the
#'   uniform distribution on \eqn{[-1, 1]}, each rounded to two decimals,
#'   inside the same seeded draw as the data. The drawn values are reported in
#'   the `coefficients` table.
#' @param n_predictors Number of predictors when `coefficients` is `NULL`, a
#'   single positive whole number defaulting to `2`. Ignored otherwise.
#' @param predictor_means,predictor_sds Predictor means and standard deviations.
#' @param correlation Optional predictor correlation matrix or tidy table.
#' @param error_sd Positive residual standard deviation.
#' @param outcome Name of the outcome variable.
#' @param seed Optional random seed.
#'
#' @return A `simulab_sim` base `data.frame` with one row per observation and
#'   columns `id`, one column per predictor, and `outcome` (renamed by
#'   `outcome`). Three tidy tables come from `as.data.frame()`:
#'   `what = "coefficients"` has one row per term with columns `term` and
#'   `coefficient`, the intercept first; `what = "effects"` has a single row
#'   with columns `signal_variance` (the quadratic form of the slopes in the
#'   predictor covariance matrix), `residual_variance` (`error_sd^2`) and the
#'   population `r_squared`
#'   `signal_variance / (signal_variance + residual_variance)`;
#'   `what = "predictor_correlation"` has one row per correlation-matrix cell
#'   with columns `row`, `column` and `correlation`.
#' @inheritParams simulate_correlation
#' @export
#'
#' @examples
#' # `coefficients` must be named; "(Intercept)" is optional.
#' result <- simulate_regression(
#'   n = 200,
#'   coefficients = c("(Intercept)" = 1, x1 = 0.5, x2 = -0.3),
#'   seed = 1
#' )
#' head(result)
#' as.data.frame(result, what = "coefficients")
#' as.data.frame(result, what = "effects")
#'
#' # With no coefficients, they are drawn and reported as the truth.
#' drawn <- simulate_regression(seed = 1)
#' as.data.frame(drawn, what = "coefficients")
simulate_regression <- function(n = 100L, coefficients = NULL, predictor_means = 0,
                                predictor_sds = 1, correlation = NULL,
                                error_sd = 1, outcome = "outcome",
                                seed = NULL,
                                batch = NULL, n_predictors = 2L) {
  if (!is.null(batch)) return(.simulate_batch(batch, seed))
  stopifnot(
    "`n` must be a single whole number of at least 2" =
      is.numeric(n) &&
        length(n) == 1L &&
        all(n >= 2) &&
        all(n == as.integer(n)),
    "`coefficients` must be NULL or a named numeric vector, with at least one element" =
      is.null(coefficients) ||
        (is.numeric(coefficients) &&
           length(coefficients) >= 1L &&
           !is.null(names(coefficients)) &&
           all(nzchar(names(coefficients)))),
    "`n_predictors` must be a single positive whole number" =
      is.numeric(n_predictors) &&
        length(n_predictors) == 1L &&
        all(n_predictors >= 1) &&
        all(n_predictors == as.integer(n_predictors)),
    "`predictor_means` must be a numeric vector" =
      is.numeric(predictor_means),
    "`predictor_sds` must be a numeric vector" =
      is.numeric(predictor_sds),
    "`correlation` must be NULL or a matrix or data frame" =
      is.null(correlation) || is.matrix(correlation) || is.data.frame(correlation),
    "`error_sd` must be a single positive number" =
      is.numeric(error_sd) &&
        length(error_sd) == 1L &&
        all(error_sd > 0),
    "`outcome` must be a single non-empty string" =
      is.character(outcome) &&
        length(outcome) == 1L &&
        all(nzchar(outcome)),
    "`seed` must be NULL or a single number" =
      is.null(seed) || (is.numeric(seed) && length(seed) == 1L)
  )
  draw_coefficients <- is.null(coefficients)
  if (draw_coefficients) {
    predictor_names <- sprintf("x%d", seq_len(as.integer(n_predictors)))
    # Placeholders of the right length for the validation below; the values
    # are drawn inside the seeded block together with the data.
    intercept <- 0
    slopes <- stats::setNames(numeric(length(predictor_names)), predictor_names)
  } else {
    intercept <- if ("(Intercept)" %in% names(coefficients)) coefficients["(Intercept)"] else 0
    slopes <- coefficients[names(coefficients) != "(Intercept)"]
    if (!length(slopes)) stop("At least one predictor coefficient is required.", call. = FALSE)
    predictor_names <- names(slopes)
  }
  if (length(predictor_means) == 1L) predictor_means <- rep(predictor_means, length(slopes))
  if (length(predictor_sds) == 1L) predictor_sds <- rep(predictor_sds, length(slopes))
  if (length(predictor_means) != length(slopes) ||
      length(predictor_sds) != length(slopes) || any(predictor_sds <= 0)) {
    stop("Predictor means and standard deviations must match coefficients.", call. = FALSE)
  }
  # Named predictor means and standard deviations, and a named or tidy
  # correlation, are matched to the coefficient names rather than by position.
  predictor_means <- .align_vector(predictor_means, predictor_names, "predictor_means")
  predictor_sds <- .align_vector(predictor_sds, predictor_names, "predictor_sds")
  if (is.null(correlation)) correlation <- diag(length(slopes))
  if (is.data.frame(correlation)) correlation <- .table_to_square_matrix(correlation, "correlation")
  correlation <- .validate_correlation_matrix(correlation)
  if (nrow(correlation) != length(slopes)) {
    stop("Predictor correlation dimensions must match coefficients.", call. = FALSE)
  }
  correlation <- .align_matrix(correlation, predictor_names, predictor_names, "correlation")
  dimnames(correlation) <- list(predictor_names, predictor_names)
  generated <- .with_seed(seed, {
    if (draw_coefficients) {
      intercept <- .draw_means(1L)
      slopes <- stats::setNames(.draw_slopes(length(predictor_names)), predictor_names)
    }
    predictors <- .draw_multivariate_normal(
      as.integer(n), predictor_means, predictor_sds, correlation
    )
    signal <- as.vector(intercept + predictors %*% slopes)
    response <- signal + stats::rnorm(as.integer(n), sd = error_sd)
    list(intercept = intercept, slopes = slopes, predictors = predictors,
         signal = signal, response = response)
  })
  intercept <- generated$intercept
  slopes <- generated$slopes
  data <- data.frame(
    id = seq_len(as.integer(n)),
    generated$predictors,
    outcome_value = generated$response,
    check.names = FALSE,
    row.names = NULL
  )
  names(data) <- c("id", predictor_names, outcome)
  scale_matrix <- diag(predictor_sds, nrow = length(predictor_sds))
  covariance <- scale_matrix %*% correlation %*% scale_matrix
  signal_variance <- as.numeric(t(slopes) %*% covariance %*% slopes)
  coefficient_table <- data.frame(
    term = c("(Intercept)", predictor_names),
    coefficient = c(intercept, slopes),
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  effects <- data.frame(
    signal_variance = signal_variance,
    residual_variance = error_sd^2,
    r_squared = signal_variance / (signal_variance + error_sd^2),
    row.names = NULL
  )
  .new_simulab_sim(
    data,
    type = "regression",
    seed = seed,
    tables = list(
      coefficients = coefficient_table,
      effects = effects,
      predictor_correlation = .matrix_to_table(correlation, "correlation")
    )
  )
}

#' Simulate Gaussian clusters
#'
#' @param n Sample size.
#' @param centers Matrix with one row per cluster and one column per variable,
#'   or a tidy data frame with columns `cluster`, `variable` and `center`.
#' @param sds Scalar, per-variable vector, cluster-by-variable matrix, or a
#'   tidy data frame with columns `cluster`, `variable` and `sd`.
#' @param proportions Optional cluster proportions.
#' @param labels Optional cluster labels.
#' @param seed Optional random seed.
#'
#' @return A `simulab_sim` base `data.frame` with one row per observation and
#'   columns `id`, `cluster` (the true membership label) and one column per
#'   variable. One tidy table comes from `as.data.frame()`:
#'   `what = "parameters"` has one row per cluster-variable combination with
#'   columns `cluster`, `variable`, `center`, `sd` and `proportion`.
#' @inheritParams simulate_correlation
#' @export
#'
#' @examples
#' # `centers` is one row per cluster and one column per variable.
#' result <- simulate_clusters(
#'   n = 150,
#'   centers = matrix(c(0, 0, 4, 4, 0, 4), nrow = 3, byrow = TRUE),
#'   sds = 1,
#'   seed = 1
#' )
#' head(result)
#' as.data.frame(result, what = "parameters")
simulate_clusters <- function(n, centers, sds = 1, proportions = NULL,
                              labels = NULL, seed = NULL,
                              batch = NULL) {
  if (!is.null(batch)) return(.simulate_batch(batch, seed))
  centers <- .tidy_to_matrix(centers, "centers", "cluster", "variable", "center")
  if (is.data.frame(sds)) {
    sds <- .tidy_to_matrix(sds, "sds", "cluster", "variable", "sd")
  }
  # `sds` is matched to `centers` by cluster and variable name, not position.
  sds <- .align_matrix(sds, rownames(centers), colnames(centers), "sds")
  if (is.matrix(centers)) sds <- .align_vector(sds, colnames(centers), "sds")
  stopifnot(
    "`n` must be a single whole number of at least 2" =
      is.numeric(n) &&
        length(n) == 1L &&
        all(n >= 2) &&
        all(n == as.integer(n)),
    "`centers` must be a matrix, with at least 2 rows" =
      is.matrix(centers) &&
        is.numeric(centers) &&
        nrow(centers) >= 2L,
    "`sds` must be a numeric vector, with at least one element" =
      is.numeric(sds) &&
        length(sds) >= 1L,
    "`proportions` must be NULL or a numeric vector" =
      is.null(proportions) || is.numeric(proportions),
    "`labels` must be NULL or an atomic vector" =
      is.null(labels) || is.atomic(labels),
    "`seed` must be NULL or a single number" =
      is.null(seed) || (is.numeric(seed) && length(seed) == 1L)
  )
  clusters <- nrow(centers)
  variables <- ncol(centers)
  if (is.null(proportions)) proportions <- rep(1 / clusters, clusters)
  if (length(proportions) != clusters || any(proportions <= 0)) {
    stop("proportions must contain one positive value per cluster.", call. = FALSE)
  }
  proportions <- proportions / sum(proportions)
  if (is.null(labels)) labels <- sprintf("Cluster %d", seq_len(clusters))
  if (length(labels) != clusters || anyDuplicated(labels)) {
    stop("labels must contain one unique value per cluster.", call. = FALSE)
  }
  sd_matrix <- if (length(sds) == 1L) {
    matrix(sds, nrow = clusters, ncol = variables)
  } else if (length(sds) == variables) {
    matrix(rep(sds, each = clusters), nrow = clusters)
  } else if (is.matrix(sds) && all(dim(sds) == dim(centers))) {
    sds
  } else {
    stop("sds must be scalar, per-variable, or cluster-by-variable.", call. = FALSE)
  }
  if (any(sd_matrix <= 0)) stop("Cluster standard deviations must be positive.", call. = FALSE)
  variable_names <- colnames(centers)
  if (is.null(variable_names)) variable_names <- sprintf("V%d", seq_len(variables))
  generated <- .with_seed(seed, {
    membership <- sample(seq_len(clusters), as.integer(n), replace = TRUE, prob = proportions)
    values <- vapply(seq_len(variables), function(variable) {
      stats::rnorm(
        as.integer(n),
        mean = centers[membership, variable],
        sd = sd_matrix[membership, variable]
      )
    }, numeric(as.integer(n)))
    list(membership = membership, values = values)
  })
  data <- data.frame(
    id = seq_len(as.integer(n)),
    cluster = labels[generated$membership],
    generated$values,
    check.names = FALSE,
    row.names = NULL
  )
  names(data) <- c("id", "cluster", variable_names)
  parameter_index <- expand.grid(
    cluster = seq_len(clusters),
    variable = seq_len(variables),
    KEEP.OUT.ATTRS = FALSE
  )
  parameters <- data.frame(
    cluster = labels[parameter_index$cluster],
    variable = variable_names[parameter_index$variable],
    center = centers[cbind(parameter_index$cluster, parameter_index$variable)],
    sd = sd_matrix[cbind(parameter_index$cluster, parameter_index$variable)],
    proportion = proportions[parameter_index$cluster],
    stringsAsFactors = FALSE,
    row.names = NULL
  )
  .new_simulab_sim(data, type = "clusters", seed = seed,
                   tables = list(parameters = parameters))
}