R/measurement.R

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

#' Simulate item-response data
#'
#' Draws latent abilities from a standard multivariate normal (mean zero, unit
#' variance, correlation `ability_correlation`), projects them onto each item
#' through `dimensions`, and samples responses from a logistic item-response
#' model.
#'
#' @details
#' For the dichotomous models the probability of a correct response is
#' `guessing + (1 - guessing) * plogis(discrimination * (theta - difficulty))`,
#' so `difficulty` is on the ability scale (the point of inflection), not an
#' intercept. `"rasch"` fixes all discriminations at one and `"2pl"` sets
#' `guessing` to zero; only `"3pl"` uses `guessing`.
#'
#' For `model = "graded"` the cumulative category probabilities are
#' `plogis(discrimination * (theta - threshold))` for each threshold of an
#' item, and the category probabilities are their successive differences, so
#' the thresholds are also on the ability scale.
#'
#' @param n Number of respondents. A single whole number of at least 2.
#' @param discrimination Positive item discriminations, a scalar recycled
#'   across items (the default, `1`) or one value per item. Must equal one for
#'   `model = "rasch"`.
#' @param difficulty Item difficulties, on the same scale as the latent
#'   ability. A numeric vector, one difficulty per item, for the dichotomous
#'   models. `model = "graded"` instead requires an item-by-threshold matrix
#'   (or a tidy data frame with columns `item`, `dimension` and `difficulty`)
#'   whose thresholds increase within each item.
#' @param dimensions Item-by-dimension loading weights, or a tidy data frame
#'   with columns `item`, `dimension` and `loading`. `NULL` (the default) uses
#'   one dimension. Each row is rescaled to unit length before use, so only
#'   the relative weights within an item matter, not their scale. Row and
#'   column names, when present, name the items and the ability columns.
#' @param ability_correlation Latent ability correlation matrix, or a tidy
#'   data frame with columns `row`, `column` and `correlation`. `NULL` (the
#'   default) makes the abilities uncorrelated.
#' @param model Logistic model: `"2pl"` (the default), `"rasch"`, `"3pl"`, or
#'   `"graded"`.
#' @param guessing Lower-asymptote guessing parameters in `[0, 1)`, a scalar
#'   recycled across items (the default, `0.2`) or one value per item. It is
#'   silently reset to zero unless `model = "3pl"`.
#' @param seed Optional random seed.
#'
#' @return A `simulab_sim` base `data.frame` with one row per respondent and
#'   columns `id` and one integer column per item, coded `0`/`1` for the
#'   dichotomous models and `0` to the number of thresholds for
#'   `model = "graded"`. An `abilities` component holds one row per respondent
#'   with columns `id` and one true ability per dimension. A `parameters`
#'   component holds one row per item for the dichotomous models, with columns
#'   `item`, `parameter` (`"difficulty"`), `category` (`NA`), `value`,
#'   `discrimination` and `guessing`; for `model = "graded"` it holds one row
#'   per item-by-threshold combination, with `parameter` equal to
#'   `"threshold"`, `category` the threshold index, and no `guessing` column.
#' @inheritParams simulate_correlation
#' @export
#'
#' @examples
#' result <- simulate_irt(
#'   n = 200,
#'   discrimination = c(1, 1.2, 0.8),
#'   difficulty = c(-1, 0, 1),
#'   seed = 1
#' )
#' head(result)
#' as.data.frame(result, what = "parameters")
#'
#' # Three-parameter logistic with a guessing floor.
#' head(simulate_irt(
#'   n = 200,
#'   discrimination = c(1, 1.2, 0.8),
#'   difficulty = c(-1, 0, 1),
#'   model = "3pl",
#'   guessing = 0.2,
#'   seed = 1
#' ))
simulate_irt <- function(n, discrimination = 1, difficulty,
                         dimensions = NULL, ability_correlation = NULL,
                         model = c("2pl", "rasch", "3pl", "graded"),
                         guessing = 0.2, seed = NULL,
                         batch = NULL) {
  if (!is.null(batch)) return(.simulate_batch(batch, seed))
  if (is.data.frame(difficulty)) {
    difficulty <- .tidy_to_matrix(difficulty, "difficulty", "item", "dimension",
                                  "difficulty")
  }
  dimensions <- .tidy_to_matrix(dimensions, "dimensions", "item", "dimension",
                                "loading")
  ability_correlation <- .tidy_to_symmetric(
    ability_correlation, "ability_correlation", "row", "column", "correlation",
    diagonal = 1
  )
  # `difficulty` names the items; `dimensions`, `discrimination` and
  # `guessing` are matched to those names, and `ability_correlation` to the
  # dimension names, rather than by position.
  item_order <- if (is.matrix(difficulty)) rownames(difficulty) else names(difficulty)
  dimensions <- .align_matrix(dimensions, item_order, NULL, "dimensions")
  discrimination <- .align_vector(discrimination, item_order, "discrimination")
  guessing <- .align_vector(guessing, item_order, "guessing")
  if (is.matrix(dimensions)) {
    ability_correlation <- .align_matrix(ability_correlation, colnames(dimensions),
                                         colnames(dimensions), "ability_correlation")
  }
  model <- match.arg(model)
  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)),
    "`discrimination` must be a finite positive numeric vector" =
      is.numeric(discrimination) &&
        all(is.finite(discrimination)) &&
        all(discrimination > 0),
    "`difficulty` must be a numeric vector" =
      is.numeric(difficulty),
    "`guessing` must be a finite numeric vector" =
      is.numeric(guessing) &&
        all(is.finite(guessing)),
    "`dimensions` must be NULL or a matrix" =
      is.null(dimensions) || is.matrix(dimensions),
    "`ability_correlation` must be NULL or a matrix" =
      is.null(ability_correlation) || is.matrix(ability_correlation),
    "`seed` must be NULL or a single number" =
      is.null(seed) || (is.numeric(seed) && length(seed) == 1L)
  )
  items <- if (is.matrix(difficulty)) nrow(difficulty) else length(difficulty)
  if (length(discrimination) == 1L) discrimination <- rep(discrimination, items)
  if (length(guessing) == 1L) guessing <- rep(guessing, items)
  if (length(discrimination) != items) stop("discrimination must match the items.", call. = FALSE)
  if (length(guessing) != items || any(guessing < 0 | guessing >= 1)) {
    stop("guessing must contain values from zero up to, but not including, one.", call. = FALSE)
  }
  if (model != "3pl") guessing <- rep(0, items)
  if (model == "rasch" && any(abs(discrimination - 1) > 1e-12)) {
    stop("Rasch discrimination must equal one.", call. = FALSE)
  }
  if (is.null(dimensions)) dimensions <- matrix(1, items, 1L)
  if (!is.numeric(dimensions) || nrow(dimensions) != items ||
      any(rowSums(abs(dimensions)) == 0)) {
    stop("dimensions must contain one non-zero loading row per item.", call. = FALSE)
  }
  dimensions_count <- ncol(dimensions)
  normalized_dimensions <- dimensions / sqrt(rowSums(dimensions^2))
  if (is.null(ability_correlation)) ability_correlation <- diag(dimensions_count)
  ability_correlation <- .validate_correlation_matrix(ability_correlation)
  if (nrow(ability_correlation) != dimensions_count) {
    stop("ability_correlation must match dimensions.", call. = FALSE)
  }
  if (model == "graded" && !is.matrix(difficulty)) {
    stop("Graded items require an item-by-threshold difficulty matrix.", call. = FALSE)
  }
  ## Without this the matrix reaches `sweep()` as an over-long `STATS`, which
  ## warns rather than stops and yields a malformed parameter table.
  if (model != "graded" && is.matrix(difficulty)) {
    stop(sprintf(
      "An item-by-threshold difficulty matrix requires model = \"graded\"; model is \"%s\", which takes one difficulty per item.",
      model
    ), call. = FALSE)
  }
  if (is.matrix(difficulty) &&
      any(apply(difficulty, 1L, function(values) is.unsorted(values, strictly = TRUE)))) {
    stop("Graded thresholds must increase within item.", call. = FALSE)
  }
  item_names <- item_order %||% rownames(dimensions) %||% sprintf("item_%d", seq_len(items))
  ability_names <- if (!is.null(colnames(dimensions))) colnames(dimensions) else sprintf("theta_%d", seq_len(dimensions_count))
  generated <- .with_seed(seed, {
    abilities <- .draw_multivariate_normal(as.integer(n), rep(0, dimensions_count),
                                           rep(1, dimensions_count), ability_correlation)
    trait <- abilities %*% t(normalized_dimensions)
    responses <- if (model != "graded") {
      probabilities <- stats::plogis(sweep(trait, 2L, difficulty, `-`) *
                                        rep(discrimination, each = as.integer(n)))
      probabilities <- sweep(probabilities, 2L, 1 - guessing, `*`)
      probabilities <- sweep(probabilities, 2L, guessing, `+`)
      matrix(stats::rbinom(length(probabilities), 1L, probabilities),
             as.integer(n), items)
    } else {
      vapply(seq_len(items), function(item) {
        cumulative <- vapply(difficulty[item, ], function(threshold) {
          stats::plogis(discrimination[item] * (trait[, item] - threshold))
        }, numeric(as.integer(n)))
        probabilities <- cbind(1, cumulative) - cbind(cumulative, 0)
        vapply(seq_len(as.integer(n)), function(person) {
          sample.int(ncol(probabilities), 1L, prob = probabilities[person, ]) - 1L
        }, integer(1))
      }, integer(as.integer(n)))
    }
    list(abilities = abilities, responses = responses)
  })
  data <- data.frame(id = seq_len(as.integer(n)), generated$responses,
                     check.names = FALSE, row.names = NULL)
  names(data) <- c("id", item_names)
  scores <- data.frame(id = seq_len(as.integer(n)), generated$abilities,
                       check.names = FALSE, row.names = NULL)
  names(scores) <- c("id", ability_names)
  parameters <- if (model == "graded") {
    index <- expand.grid(item = seq_len(items), threshold = seq_len(ncol(difficulty)))
    data.frame(item = item_names[index$item], parameter = "threshold",
               category = index$threshold,
               value = difficulty[cbind(index$item, index$threshold)],
               discrimination = discrimination[index$item], row.names = NULL)
  } else {
    data.frame(item = item_names, parameter = "difficulty", category = NA_integer_,
               value = difficulty, discrimination = discrimination,
               guessing = guessing, row.names = NULL)
  }
  .new_simulab_sim(data, paste0("irt_", model), seed,
                   list(parameters = parameters, abilities = scores))
}

#' Simulate a hidden Markov model
#'
#' Draws an initial hidden state for each sequence from `initial`, evolves it
#' with `transition` for `chain_length` occasions, and emits one observed
#' category per occasion from `emission`.
#'
#' @param n Number of sequences. A single positive whole number.
#' @param transition Hidden-state transition matrix, or a tidy data frame with
#'   columns `from`, `to` and `probability`. Square, with rows indexing the
#'   current state and columns the next state; every row must sum to one.
#' @param chain_length Sequence length. A single positive whole number; every
#'   sequence has the same length.
#' @param emission Hidden-state-by-observed-category probability matrix, or a
#'   tidy data frame with columns `state`, `observation` and `probability`.
#'   Rows index the hidden states and must sum to one; columns index the
#'   observed categories.
#' @param initial Initial hidden-state probabilities, one per hidden state,
#'   summing to one. `NULL` (the default) starts every sequence in the first
#'   hidden state with probability one; supply a vector for a random start.
#' @param state_labels,observation_labels Optional labels, one unique value
#'   per hidden state and per observed category. They default to
#'   `"State 1"`, `"State 2"`, ... and `"Observation 1"`, `"Observation 2"`,
#'   ...
#' @param seed Optional random seed.
#'
#' @return A long-form `simulab_sim` base `data.frame` with one row per
#'   sequence-occasion (`n * chain_length` rows) and columns `id`, `occasion`,
#'   `state` (the true hidden state label) and `observation` (the emitted
#'   category label). Components `transitions` (columns `from`, `to`,
#'   `probability`), `emissions` (columns `state`, `observation`,
#'   `probability`) and `initial_probabilities` (columns `state`,
#'   `probability`) hold the generating parameters in tidy form.
#' @inheritParams simulate_correlation
#' @export
#'
#' @examples
#' result <- simulate_hmm(
#'   n = 50,
#'   transition = matrix(c(0.7, 0.3, 0.4, 0.6), nrow = 2, byrow = TRUE),
#'   chain_length = 20,
#'   emission = matrix(c(0.9, 0.1, 0.2, 0.8), nrow = 2, byrow = TRUE),
#'   seed = 1
#' )
#' head(result)
#' components(result)
simulate_hmm <- function(n, transition, chain_length, emission, initial = NULL,
                         state_labels = NULL, observation_labels = NULL,
                         seed = NULL,
                         batch = NULL) {
  if (!is.null(batch)) return(.simulate_batch(batch, seed))
  transition <- .tidy_to_square(transition, "transition", "from", "to",
                                "probability")
  emission <- .tidy_to_matrix(emission, "emission", "state", "observation",
                              "probability")
  # Emission rows and named initial probabilities are matched to the hidden
  # states by name, and a named transition's columns to its rows.
  if (is.matrix(transition)) {
    transition <- .align_matrix(transition, NULL, rownames(transition), "transition")
    emission <- .align_matrix(emission, rownames(transition), NULL, "emission")
    initial <- .align_vector(initial, rownames(transition), "initial")
  }
  stopifnot(
    "`n` must be a single positive whole number" =
      is.numeric(n) &&
        length(n) == 1L &&
        all(n >= 1) &&
        all(n == as.integer(n)),
    "`transition` must be a matrix" =
      is.matrix(transition) &&
        is.numeric(transition),
    "`chain_length` must be a single positive whole number" =
      is.numeric(chain_length) &&
        length(chain_length) == 1L &&
        all(chain_length >= 1) &&
        all(chain_length == as.integer(chain_length)),
    "`emission` must be a matrix" =
      is.matrix(emission) &&
        is.numeric(emission),
    "`initial` must be NULL or a numeric vector" =
      is.null(initial) || is.numeric(initial),
    "`seed` must be NULL or a single number" =
      is.null(seed) || (is.numeric(seed) && length(seed) == 1L)
  )
  states <- nrow(transition)
  if (ncol(transition) != states || nrow(emission) != states ||
      any(transition < 0) || any(emission < 0) ||
      any(abs(rowSums(transition) - 1) > 1e-8) ||
      any(abs(rowSums(emission) - 1) > 1e-8)) {
    stop("Transition and emission rows must be valid probability distributions.", call. = FALSE)
  }
  if (is.null(initial)) initial <- c(1, rep(0, states - 1L))
  if (length(initial) != states || any(initial < 0) || abs(sum(initial) - 1) > 1e-8) {
    stop("initial must contain one probability per hidden state.", call. = FALSE)
  }
  if (is.null(state_labels)) state_labels <- rownames(transition) %||% sprintf("State %d", seq_len(states))
  observed <- ncol(emission)
  if (is.null(observation_labels)) {
    observation_labels <- colnames(emission) %||% sprintf("Observation %d", seq_len(observed))
  }
  if (length(state_labels) != states || anyDuplicated(state_labels) ||
      length(observation_labels) != observed || anyDuplicated(observation_labels)) {
    stop("State and observation labels must match their matrices and be unique.", call. = FALSE)
  }
  chains <- .with_seed(seed, lapply(seq_len(as.integer(n)), function(id) {
    hidden <- .simulate_one_chain(transition, seq_len(states), as.integer(chain_length),
                                  sample.int(states, 1L, prob = initial))
    observations <- vapply(hidden, function(state) {
      sample.int(observed, 1L, prob = emission[state, ])
    }, integer(1))
    data.frame(id = id, occasion = seq_along(hidden),
               state = state_labels[hidden], observation = observation_labels[observations],
               stringsAsFactors = FALSE, row.names = NULL)
  }))
  data <- do.call(rbind, chains)
  transition_table <- .matrix_to_table(transition, "probability", state_labels, state_labels)
  names(transition_table)[1:2] <- c("from", "to")
  emission_table <- .matrix_to_table(emission, "probability", state_labels, observation_labels)
  names(emission_table)[1:2] <- c("state", "observation")
  initial_table <- data.frame(state = state_labels, probability = initial, row.names = NULL)
  .new_simulab_sim(data, "hmm", seed,
                   list(transitions = transition_table, emissions = emission_table,
                        initial_probabilities = initial_table))
}