R/family-response.R

In sdmTMB: Spatial and Spatiotemporal SPDE-Based GLMMs with 'TMB'

.process_binomial_like_rows <- function(y_i, size, weights, rows,
  family_label, allow_counts, weighted_binary_counts = FALSE) {
  if (!any(rows)) return(list(y_i = y_i, size = size, weights = weights))
  values <- y_i[rows & !is.na(y_i)]
  if (!is.numeric(values) || any(values < 0)) {
    cli_abort("{family_label} rows must have non-negative numeric response values.")
  }

  classified <- .classify_binomial_like_numeric_rows(
    y_i, weights, allow_counts, weighted_binary_counts
  )
  counts <- rows & classified$count_rows
  bernoulli <- rows & classified$bernoulli_rows
  proportions <- rows & classified$proportion_rows

  # A numeric binomial response that contains proportions uses `weights` as
  # the number of trials. Binary rows with an available weight are then
  # boundary proportions (e.g. 3/3), which keeps the fit identical to an
  # equivalent cbind(success, failure) representation. Binary rows with a
  # missing weight stay Bernoulli, and the all-binary case keeps `weights`
  # as likelihood weights.
  if (!allow_counts && any(proportions) && !is.null(weights)) {
    weighted_binary <- bernoulli & !is.na(weights)
    proportions <- proportions | weighted_binary
    bernoulli <- bernoulli & !weighted_binary
  }

  if (!allow_counts && any(rows & !is.na(y_i) & y_i > 1)) {
    cli_abort("Binomial rows must have values between 0 and 1.")
  }
  if (allow_counts && any(proportions & !(y_i > 0 & y_i < 1))) {
    cli_abort("{family_label} rows must be integer counts or proportions in (0, 1).")
  }
  if (any(bernoulli) && !is.null(weights)) weights[bernoulli & is.na(weights)] <- 1

  converted <- counts | proportions
  if (any(converted)) {
    if (is.null(weights)) {
      kind <- if (allow_counts) "proportions or counts" else "proportions"
      cli_abort("{family_label} rows with {kind} require `weights` to supply binomial size.")
    }
    missing_size <- converted & is.na(weights)
    if (any(missing_size)) {
      # Match model-frame NA handling: a missing trial size makes the
      # observation unavailable to the likelihood, just like a missing
      # response. Use harmless placeholders for TMB's auxiliary vectors.
      y_i[missing_size] <- NA_real_
      size[missing_size] <- 1
      weights[missing_size] <- 1
      converted[missing_size] <- FALSE
    }
    if (any(weights[converted] <= 0)) {
      cli_abort("`weights` must be non-missing and > 0 for {tolower(family_label)} rows.")
    }
    if (any(counts & weights < y_i, na.rm = TRUE)) {
      cli_abort("{family_label} counts must be <= `weights` (binomial size).")
    }
    size[converted] <- weights[converted]
    y_i[proportions] <- y_i[proportions] * weights[proportions]
    weights[converted] <- 1
  }
  list(y_i = y_i, size = size, weights = weights)
}

.family_spec_process_response <- function(y_i, size, weights, family_spec) {
  row_family <- .family_spec_response_family_id(family_spec, y_i)
  component1 <- family_spec$components[family_spec$components$component == 1L, ]
  family_name <- component1$family_name[match(row_family, component1$family_id)]
  single <- family_spec$families$combine_kind[row_family] == "single"

  out <- .process_binomial_like_rows(
    y_i, size, weights, single & family_name == "binomial",
    "Binomial", allow_counts = FALSE
  )
  out <- .process_binomial_like_rows(
    out$y_i, out$size, out$weights, single & family_name == "censored_binomial",
    "Censored binomial", allow_counts = TRUE, weighted_binary_counts = TRUE
  )
  .process_binomial_like_rows(
    out$y_i, out$size, out$weights,
    single & family_name %in% c("betabinomial", "censored_betabinomial"),
    "Betabinomial", allow_counts = TRUE, weighted_binary_counts = TRUE
  )
}

# Censored binomial and beta-binomial counts, trials, and bounds on the count scale, as
# exact integers that both backends can use directly. Values within
# floating-point error of an integer (e.g., proportions times trials) are
# snapped to it. Otherwise counts and trials must be integers. A bound `U`
# can be any real number and means the integer count `floor(U)`; `NA`
# (right-censored) means the number of trials, the largest possible count.
# Returns `upr` as NULL if there are no such rows.
.normalize_censored_trials <- function(y_i, size, family_spec, upr = NULL) {
  row_family <- .family_spec_response_family_id(family_spec, y_i)
  component1 <- family_spec$components[family_spec$components$component == 1L, ]
  family_name <- component1$family_name[match(row_family, component1$family_id)]
  single <- family_spec$families$combine_kind[row_family] == "single"
  rows <- single & !is.na(y_i) &
    family_name %in% c("censored_binomial", "censored_betabinomial")
  if (is.null(upr) || !any(rows)) return(list(y_i = y_i, size = size, upr = NULL))
  snap <- function(x) ifelse(.is_whole_number(x), round(x), x)
  y_i[rows] <- snap(y_i[rows])
  size[rows] <- snap(size[rows])
  if (any(!is.finite(size[rows]) | size[rows] != round(size[rows]))) {
    cli_abort("Censored (beta-)binomial rows must have whole-number trial sizes.")
  }
  if (any(y_i[rows] != round(y_i[rows]))) {
    cli_abort(paste(
      "Censored (beta-)binomial rows must have whole-number counts",
      "(e.g., proportions times `weights` must be whole numbers)."
    ))
  }
  right <- rows & !is.finite(upr)
  upr[right] <- size[right]
  upr[rows] <- floor(snap(upr[rows]))
  if (any(upr[rows] < y_i[rows] | upr[rows] > size[rows])) {
    cli_abort(paste(
      "`censored_upper` must be between the observed count and the",
      "number of trials (from `weights`) for censored (beta-)binomial rows.",
      "Non-integer bounds are rounded down to the largest possible count."
    ))
  }
  list(y_i = y_i, size = size, upr = upr)
}

.family_spec_validate_response <- function(y_i, family_spec, upr = NULL) {
  row_family <- .family_spec_response_family_id(family_spec, y_i)
  component1 <- family_spec$components[family_spec$components$component == 1L, ]
  family_name <- component1$family_name[match(row_family, component1$family_id)]
  link_name <- component1$link_name[match(row_family, component1$family_id)]
  single <- family_spec$families$combine_kind[row_family] == "single"

  if (any(y_i[single & family_name %in% c("Gamma", "lognormal")] <= 0, na.rm = TRUE)) {
    cli_abort("Gamma and lognormal must have response values > 0.")
  }
  ordered_beta <- single & family_name == "ordbeta"
  if (any(y_i[ordered_beta] < 0 | y_i[ordered_beta] > 1, na.rm = TRUE)) {
    cli_abort("Ordered beta requires response values in [0, 1].")
  }
  if (any(y_i[single & link_name == "log"] < 0, na.rm = TRUE)) {
    cli_abort("`link = 'log'` but the response data include values < 0.")
  }
  censored <- single & family_name %in%
    c("censored_poisson", "censored_nbinom1", "censored_nbinom2")
  if (!is.null(upr) && any(y_i[censored] > upr[censored], na.rm = TRUE)) {
    cli_abort("Observed values must be <= `censored_upper` for censored Poisson and negative binomial rows.")
  }
  censored_nb <- single & family_name %in% c("censored_nbinom1", "censored_nbinom2")
  if (any(y_i[censored_nb] != round(y_i[censored_nb]), na.rm = TRUE)) {
    cli_abort("Censored negative binomial rows must have whole-number counts.")
  }
  invisible(NULL)
}

.family_spec_build_response <- function(y_i, family_spec) {
  row_family <- .family_spec_response_family_id(family_spec, y_i)
  response <- matrix(NA_real_, nrow = length(y_i), ncol = family_spec$n_m)
  single <- family_spec$families$combine_kind[row_family] == "single"
  response[single, 1L] <- y_i[single]
  response[!single, 1L] <- as.numeric(y_i[!single] > 0)
  if (family_spec$n_m == 2L) {
    positive <- !single & y_i > 0
    response[positive, 2L] <- y_i[positive]
  }
  response
}

.prepare_family_response <- function(y_i, weights, family_spec, upr = NULL) {
  component1 <- .family_spec_component_value(family_spec, 1L, 1L, "family_name")
  ordinary_binomial_like <- family_spec$n_f == 1L && family_spec$n_m == 1L &&
    component1 %in% c("binomial", "censored_binomial", "betabinomial",
      "censored_betabinomial")
  if (ordinary_binomial_like && (is.character(y_i) || is.factor(y_i))) {
    y_i <- factor(y_i)
    if (nlevels(y_i) > 2L) cli_abort("More than 2 levels detected for response")
    y_i <- pmin(as.numeric(y_i) - 1, 1)
  }
  if (ordinary_binomial_like && is.matrix(y_i)) {
    if (ncol(y_i) != 2L) cli_abort("Binomial matrix responses must have two columns.")
    processed <- list(y_i = y_i[, 1L], size = rowSums(y_i), weights = weights)
  } else {
    processed <- .family_spec_process_response(y_i, rep(1, NROW(y_i)),
      weights, family_spec)
  }
  .family_spec_validate_response(processed$y_i, family_spec, upr)
  censored <- .normalize_censored_trials(processed$y_i, processed$size,
    family_spec, upr)
  processed[names(censored)] <- censored
  processed$response <- .family_spec_build_response(processed$y_i, family_spec)
  processed
}

# Resolve `sdmTMB(censored_upper = )`: a vector or a column name in `data`.
.censored_upper_arg <- function(censored_upper, data) {
  if (is.null(censored_upper)) return(NULL)
  if (is.character(censored_upper)) {
    if (length(censored_upper) != 1L || !censored_upper %in% names(data)) {
      cli_abort("`censored_upper` must be a numeric vector or the name of a column in `data`.")
    }
    censored_upper <- data[[censored_upper]]
  }
  if (length(censored_upper) != nrow(data)) {
    cli_abort("`censored_upper` must have one value per row of `data`.")
  }
  if (anyNA(censored_upper)) {
    cli_warn(c(
      "Using `NA` in `censored_upper` for no upper bound is deprecated.",
      "i" = "Use `Inf` instead. `NA` is being treated as `Inf`."
    ), .frequency = "regularly", .frequency_id = "sdmTMB_censored_upper_na")
    censored_upper[is.na(censored_upper)] <- Inf
  }
  # Counts are integers, so bounds are rounded down; the likelihood sums the
  # counts below the bound and the gradients assume an integer bound. Values
  # within floating-point error of an integer are that integer.
  whole <- is.finite(censored_upper) & .is_whole_number(censored_upper)
  censored_upper[whole] <- round(censored_upper[whole])
  floor(censored_upper)
}