R/lr_age.R

In mispitools: Missing Person Identification Tools

#' Likelihood Ratio for Age Variable
#'
#' @description
#' Simulates age observations and optionally computes likelihood ratios (LRs)
#' under either H1 (unidentified person is the missing person) or H2
#' (unidentified person is not the missing person).
#'
#' Ages are categorized into two groups based on whether they fall within
#' the missing person's estimated age range:
#' \itemize{
#'   \item T1: Age within range (MPa - MPr) to (MPa + MPr)
#'   \item T0: Age outside range
#' }
#'
#' @param MPa Numeric. Missing person's estimated age in years. Default: 40.
#' @param MPr Numeric. Age range tolerance (plus/minus years). The matching
#'   interval is (MPa - MPr) to (MPa + MPr). Default: 6.
#' @param UHRr Numeric. Additional uncertainty range for the unidentified
#'   person's age estimation. Default: 1.
#' @param gam Numeric. Gamma parameter for age uncertainty scaling.
#'   The uncertainty interval is age +/- (gam * age + UHRr). Default: 0.07.
#' @param numsims Integer. Number of simulations to perform. Default: 1000.
#' @param nsims Deprecated. Use \code{numsims} instead.
#' @param epa Numeric (0-1). Error rate for age categorization. Default: 0.05.
#' @param erRa Numeric (0-1). Error rate in the reference/database.
#'   Defaults to \code{epa}.
#' @param H Integer (1 or 2). Hypothesis to simulate under:
#'   \itemize{
#'     \item 1: H1 (Related) - Unidentified person IS the missing person
#'     \item 2: H2 (Unrelated) - Unidentified person is NOT the missing person
#'   }
#'   Default: 1.
#' @param modelA Character. Reference age distribution model:
#'   \itemize{
#'     \item "uniform": Assumes uniform age distribution (default)
#'     \item "custom": Uses empirical frequencies from simulations
#'   }
#' @param LR Logical. If TRUE, compute and return LR values. Default: FALSE.
#' @param seed Integer. Random seed for reproducibility. Default: 1234.
#'
#' @return A data.frame with columns:
#'   \itemize{
#'     \item \code{group}: Age group classification ("T1" or "T0")
#'     \item \code{age}: Simulated age value
#'     \item \code{UHRmin}: Lower bound of uncertainty interval
#'     \item \code{UHRmax}: Upper bound of uncertainty interval
#'     \item \code{LRa}: Likelihood ratio (only if LR = TRUE)
#'   }
#'
#' @details
#' \strong{Under H1 (Related):}
#' Age is sampled to fall within the MP's range with probability (1 - erRa),
#' outside with probability erRa.
#'
#' \strong{Under H2 (Unrelated):}
#' Age is sampled uniformly from 1-80, then categorized.
#'
#' \strong{LR Calculation (uniform model):}
#' \itemize{
#'   \item LR(T1) = (1 - epa) / P(T1), where P(T1) = 2*MPr/80
#'   \item LR(T0) = epa / P(T0), where P(T0) = 1 - P(T1)
#' }
#'
#' @section Deprecation:
#' Soft-deprecated in mispitools 2.0. This feature is generalised by
#' \code{\link{nongenetic_feature}} (a continuous feature); collapsing
#' the population grid to two cells (within tolerance vs outside)
#' reproduces the T1 / T0 LR exactly. The legacy function still works
#' for the 2.0 release-candidate cycle and will be removed afterwards.
#'
#' @seealso
#' \code{\link{nongenetic_feature}} for the unified replacement,
#' \code{\link{sim_lr_prelim}} for unified preliminary LR simulations,
#' \code{\link{lr_sex}}, \code{\link{lr_hair_color}} for other variables.
#'
#' @references
#' Marsico FL, et al. (2023). "Likelihood ratios for non-genetic evidence
#' in missing person cases." \emph{Forensic Science International: Genetics},
#' 66, 102891. \doi{10.1016/j.fsigen.2023.102891}
#'
#' @export
#' @import dplyr
#' @examples
#' # Simulate under H1 (related)
#' sim_h1 <- lr_age(MPa = 40, MPr = 6, H = 1, numsims = 100)
#' table(sim_h1$group)
#'
#' # Simulate under H2 with LR values
#' sim_h2 <- lr_age(MPa = 40, MPr = 6, H = 2, numsims = 100, LR = TRUE)
#' head(sim_h2)
#'
#' # Narrower age range (more discriminating)
#' sim_narrow <- lr_age(MPa = 35, MPr = 3, numsims = 500, LR = TRUE)
#' summary(sim_narrow$LRa)

lr_age <- function(MPa = 40,
                   MPr = 6,
                   UHRr = 1,
                   gam = 0.07,
                   numsims = 1000,
                   epa = 0.05,
                   erRa = epa,
                   H = 1,
                   modelA = c("uniform", "custom")[1],
                   LR = FALSE,
                   seed = 1234,
                   nsims = NULL) {

  ng_soft_deprecate("lr_age",
    "Use nongenetic_feature(type = \"age\", ...) with the per-feature engine.")

  # Handle deprecated nsims parameter
  if (!is.null(nsims)) {
    warning("Parameter 'nsims' is deprecated. Use 'numsims' instead.",
            call. = FALSE)
    numsims <- nsims
  }

  # Input validation
  if (!is.numeric(MPa) || MPa < 1 || MPa > 80) {
    stop("MPa must be between 1 and 80")
  }
  if (!is.numeric(MPr) || MPr < 0) {
    stop("MPr must be non-negative")
  }
  if (MPa - MPr < 1 || MPa + MPr > 80) {
    warning("Age range (MPa +/- MPr) extends beyond [1, 80]; results may be affected")
  }
  if (!is.numeric(gam) || gam < 0) {
    stop("gam must be non-negative")
  }
  if (!is.numeric(numsims) || numsims < 1) {
    stop("numsims must be a positive integer")
  }
  if (!is.numeric(epa) || epa < 0 || epa > 1) {
    stop("epa must be between 0 and 1")
  }
  if (!is.numeric(erRa) || erRa < 0 || erRa > 1) {
    stop("erRa must be between 0 and 1")
  }
  if (!H %in% c(1, 2)) {
    stop("H must be 1 or 2")
  }
  if (!modelA %in% c("uniform", "custom")) {
    stop("modelA must be 'uniform' or 'custom'")
  }

  set.seed(seed)
  sims <- list()
  Age <- seq(1, 80)
  MPmin <- max(1, MPa - MPr)  # Ensure valid range
  MPmax <- min(80, MPa + MPr)

  # Calculate probabilities under uniform model
  if (modelA == "uniform") {
    T1p <- (MPmax - MPmin) / length(Age)
    T0p <- 1 - T1p
    LR1 <- (1 - epa) / T1p
    LR0 <- epa / T0p
  }

  # Define age groups (inclusive of boundaries)
  T1a <- Age[Age <= MPmax & Age >= MPmin]
  T0a <- Age[Age < MPmin | Age > MPmax]

  if (H == 1) {
    # H1: Sample with high probability of being in range
    group <- unlist(sample(c("T1", "T0"), size = numsims, prob = c(1 - erRa, erRa), replace = TRUE))
    ages <- unlist(lapply(group, function(x) ifelse(x == "T1", sample(T1a, 1), sample(T0a, 1))))

    sims <- as.data.frame(cbind(group, ages))
    names(sims) <- c("group", "age")

    sims <- dplyr::mutate(sims, UHRmin = as.numeric(ages) - gam * as.numeric(ages) - UHRr)
    sims <- dplyr::mutate(sims, UHRmax = as.numeric(ages) + gam * as.numeric(ages) + UHRr)
  }
  else if (H == 2) {
    # H2: Sample uniformly from population
    ages <- unlist(sample(Age, numsims, replace = TRUE))
    group <- unlist(lapply(ages, function(x) ifelse(x >= MPmin & x <= MPmax, "T1", "T0")))

    sims <- as.data.frame(cbind(group, ages))
    names(sims) <- c("group", "age")
    sims <- dplyr::mutate(sims, UHRmin = as.numeric(ages) - gam * as.numeric(ages) - UHRr)
    sims <- dplyr::mutate(sims, UHRmax = as.numeric(ages) + gam * as.numeric(ages) + UHRr)
  }

  # Custom model: calculate probabilities from empirical data
  if (modelA == "custom") {
    T1p <- length(subset(sims$group, sims$group == "T1")) / length(sims$group)
    T0p <- 1 - T1p
    LR1 <- (1 - epa) / T1p
    LR0 <- epa / T0p
  }

  if (LR == TRUE) {
    sims <- dplyr::mutate(sims, LRa = ifelse(group == "T1", LR1, LR0))
    names(sims) <- c("group", "Age", "UHRmin", "UHRmax", "LRa")
    return(sims)
  }

  return(as.data.frame(sims))
}