Use different weight strategies (like in Nova)

NAMESPACE

@@ -14,11 +14,13 @@ export(bootstrap_metrics)
 export(bootstrap_options)
 export(calculate_bayesian_impact)
 export(calculate_bootstrap_summ)
+export(calculate_fit_weights)
 export(calculate_shrinkage)
 export(calculate_stats)
 export(compare_psn_execute_results)
 export(compare_psn_proseval_results)
 export(fit_options)
+export(fit_weights)
 export(group_by_dose)
 export(group_by_time)
 export(install_default_literature_model)

R/calculate_fit_weights.R

@@ -0,0 +1,179 @@
+#' Construct a fit-weighting scheme specification
+#'
+#' Self-documenting helper that builds a `fit_weights` object describing how
+#' older observations should be downweighted relative to more recent ones
+#' during the iterative MAP Bayesian fitting step. The result can be passed as
+#' the `weights` argument to [run_eval()].
+#'
+#' Available schemes:
+#' - `"weight_all"`: all samples weighted equally (weight = 1).
+#' - `"weight_last_only"`: only the most recent sample is used (weight = 1),
+#'   all others are excluded (weight = 0).
+#' - `"weight_last_two_only"`: only the two most recent samples are used.
+#' - `"weight_gradient_linear"`: weights increase linearly from a minimum
+#'   (`w1`) for samples older than `t1` days to a maximum (`w2`) for samples
+#'   more recent than `t2` days. Accepts scheme parameter `gradient`, a list
+#'   with named elements `t1`, `w1`, `t2`, `w2`. Default:
+#'   `list(t1 = 7, w1 = 0, t2 = 2, w2 = 1)`.
+#' - `"weight_gradient_exponential"`: weights decay exponentially with the age
+#'   of the sample. Accepts scheme parameters `t12_decay` (half-life of decay
+#'   in hours, default 48) and `t_start` (delay in hours before decay starts,
+#'   default 0).
+#'
+#' @param scheme name of the weighting scheme (see Details).
+#' @param ... scheme-specific parameters, e.g. `t12_decay = 72` for
+#'   `"weight_gradient_exponential"`, or
+#'   `gradient = list(t1 = 5, w1 = 0.1, t2 = 1, w2 = 1)` for
+#'   `"weight_gradient_linear"`.
+#'
+#' @returns an object of class `fit_weights`.
+#' @examples
+#' fit_weights("weight_all")
+#' fit_weights("weight_gradient_exponential", t12_decay = 72)
+#' fit_weights("weight_gradient_linear", gradient = list(t1 = 5, w1 = 0.1, t2 = 1, w2 = 1))
+#' @export
+fit_weights <- function(
+  scheme = c(
+    "weight_all",
+    "weight_last_only",
+    "weight_last_two_only",
+    "weight_gradient_linear",
+    "weight_gradient_exponential"
+  ),
+  ...
+) {
+  scheme <- match.arg(scheme)
+  structure(
+    list(scheme = scheme, params = list(...)),
+    class = "fit_weights"
+  )
+}
+
+#' Calculate time-based sample weights for MAP Bayesian fitting
+#'
+#' Downweights older observations relative to more recent ones during the
+#' iterative MAP Bayesian fitting step. Can be passed as the `weights`
+#' argument to [run_eval()].
+#'
+#' `weights` may be a [fit_weights()] object, a string naming a scheme, or a
+#' named list with a `scheme` element plus optional scheme-specific parameters
+#' (e.g. `list(scheme = "weight_gradient_exponential", t12_decay = 72)`). See
+#' [fit_weights()] for the available schemes and their parameters.
+#'
+#' @param weights weighting scheme: a [fit_weights()] object, a string with the
+#'   scheme name, or a named list with a `scheme` element plus optional
+#'   scheme-specific parameters.
+#' @param t numeric vector of observation times (in hours)
+#'
+#' @returns numeric vector of weights the same length as `t`, or `NULL` if
+#'   `weights` is `NULL` or the scheme is not recognized.
+#' @export
+calculate_fit_weights <- function(weights = NULL, t = NULL) {
+  if (is.null(weights) || is.null(t)) return(NULL)
+
+  weights <- as_fit_weights(weights)
+  if (is.null(weights)) return(NULL)
+
+  weight_vec <- switch(
+    weights$scheme,
+    weight_gradient_linear      = .wt_gradient_linear(t, weights$params),
+    weight_gradient_exponential = .wt_gradient_exponential(t, weights$params),
+    weight_last_only            = .wt_last_only(t),
+    weight_last_two_only        = .wt_last_two_only(t),
+    weight_all                  = .wt_all(t)
+  )
+
+  if (!is.null(weight_vec)) {
+    weight_vec[t < 0] <- 0
+  }
+
+  weight_vec
+}
+
+# Normalize the various accepted `weights` inputs into a `fit_weights` object.
+# Returns NULL (with a warning) when the scheme cannot be recognized, so that
+# callers can cleanly ignore invalid input rather than error.
+as_fit_weights <- function(weights) {
+  if (inherits(weights, "fit_weights")) return(weights)
+
+  if (is.character(weights)) {
+    scheme <- weights
+    params <- list()
+  } else if (is.list(weights)) {
+    scheme <- weights$scheme
+    params <- weights[setdiff(names(weights), "scheme")]
+  } else {
+    warning("Weighting scheme not recognized, ignoring weights.")
+    return(NULL)
+  }
+
+  valid_schemes <- c(
+    "weight_gradient_linear",
+    "weight_gradient_exponential",
+    "weight_last_only",
+    "weight_last_two_only",
+    "weight_all"
+  )
+
+  if (length(scheme) != 1 || is.na(scheme) || !scheme %in% valid_schemes) {
+    warning("Weighting scheme not recognized, ignoring weights.")
+    return(NULL)
+  }
+
+  structure(list(scheme = scheme, params = params), class = "fit_weights")
+}
+
+.wt_gradient_linear <- function(t, params = list()) {
+  gradient <- list(t1 = 7, w1 = 0, t2 = 2, w2 = 1)
+  if (!is.null(params$gradient)) {
+    gradient[names(params$gradient)] <- params$gradient
+  }
+  if (gradient$t2 > gradient$t1) {
+    warning(
+      "weight_gradient_linear: t2 (", gradient$t2, ") > t1 (", gradient$t1,
+      "). t1 should be the older threshold and t2 the more recent one."
+    )
+  }
+  t_start <- max(c(0, max(t) - gradient$t1 * 24))
+  t_end   <- max(c(0, max(t) - gradient$t2 * 24))
+  if (t_end <= t_start) {
+    ifelse(t >= t_end, gradient$w2, gradient$w1)
+  } else {
+    ifelse(
+      t <= t_start, gradient$w1,
+      ifelse(
+        t >= t_end, gradient$w2,
+        gradient$w1 + (gradient$w2 - gradient$w1) * (t - t_start) / (t_end - t_start)
+      )
+    )
+  }
+}
+
+.wt_gradient_exponential <- function(t, params = list()) {
+  t12_decay <- if (!is.null(params$t12_decay)) params$t12_decay else 48
+  k_decay <- log(2) / t12_decay
+  t_diff <- max(t) - t
+  if (!is.null(params$t_start)) {
+    t_diff <- t_diff - params$t_start
+    t_diff <- ifelse(t_diff < 0, 0, t_diff)
+  }
+  exp(-k_decay * t_diff)
+}
+
+.wt_last_only <- function(t) {
+  weight_vec <- rep(0, length(t))
+  weight_vec[which.max(t)] <- 1
+  weight_vec
+}
+
+.wt_last_two_only <- function(t) {
+  weight_vec <- rep(0, length(t))
+  ranked <- order(t, decreasing = TRUE)
+  weight_vec[ranked[1]] <- 1
+  if (length(t) > 1) weight_vec[ranked[2]] <- 1
+  weight_vec
+}
+
+.wt_all <- function(t) {
+  rep(1, length(t))
+}

R/run_eval.R

@@ -14,6 +14,14 @@
 #'   this can be used to group peaks and troughs together, or to group
 #'   observations on the same day together. Grouping will be done prior to
 #'   running the analysis, so cannot be changed afterwards.
+#' @param weights optional sample downweighting scheme based on how long ago
+#'   observations were taken. Either a string naming the scheme (e.g.
+#'   `"weight_gradient_exponential"`), or a named list with a `scheme` element
+#'   plus any scheme-specific parameters (e.g.
+#'   `list(scheme = "weight_gradient_exponential", t12_decay = 72)`). See
+#'   [calculate_fit_weights()] for all available schemes and their parameters.
+#'   Default is `NULL` (no downweighting; all included samples are weighted
+#'   equally).
 #' @param censor_covariates with the `proseval` tool in PsN, there is “data
 #'   leakage” (of future covariates data): since the NONMEM dataset in each step
 #'   contains the covariates for the future, this is technically data leakage,
@@ -185,6 +193,7 @@ run_eval <- function(
     .x = data_parsed,
     .f = run_eval_core,
     mod_obj = mod_obj,
+    weights = weights,
     censor_covariates = censor_covariates,
     weight_prior = weight_prior,
     incremental = incremental,

R/run_eval_core.R

@@ -30,13 +30,14 @@ run_eval_core <- function(
   for(i in seq_along(iterations)) {
 
     ## Select which samples should be used in fit, for regular iterative
-    ## forecasting and incremental.
-    ## TODO: handle weighting down of earlier samples
-    weights <- handle_sample_weighting(
+    ## forecasting and incremental. Applies time-based downweighting if a
+    ## weighting scheme is provided via the `weights` argument.
+    sample_weights <- handle_sample_weighting(
       obs_data,
       iterations,
       incremental,
-      i
+      i,
+      weights = weights
     )
 
     ## Should covariate data be leaked? PsN::proseval does this,
@@ -69,7 +70,7 @@ run_eval_core <- function(
           covariates = cov_data,
           regimen = data$regimen,
           weight_prior = weight_prior,
-          weights = weights,
+          weights = sample_weights,
           iov_bins = mod_obj$bins,
           verbose = FALSE
         ),
@@ -102,6 +103,7 @@ run_eval_core <- function(
       par_dummy[, eta_names] <- NA_real_
       fit_pars <- dplyr::mutate(as.data.frame(par_dummy), id = obs_data$id[1])
     } else {
+      ## Data frame with predictive data
       ## Data frame with predictive data
       pred_data <- tibble::tibble(
         id = obs_data$id,
@@ -115,7 +117,7 @@ run_eval_core <- function(
         wres = fit$wres,
         cwres = fit$cwres,
         ofv = fit$fit$value,
-        ss_w = ss(fit$dv, fit$ipred, weights),
+        ss_w = ss(fit$dv, fit$ipred, sample_weights),
         `_iteration` = iterations[i],
         `_grouper` = obs_data$`_grouper`
       )
@@ -248,9 +250,9 @@ handle_covariate_censoring <- function(
 
 #' Handle weighting of samples
 #'
-#' This function is used to select the samples used in the fit (1 or 0),
-#' but also to select their weight, if a sample weighting strategy is
-#' selected.
+#' Binary selection of which samples are used in the fit (0 = excluded,
+#' 1 = included), combined with optional continuous downweighting of older
+#' samples via a time-based scheme (see [calculate_fit_weights()]).
 #'
 #' @inheritParams run_eval_core
 #' @param obs_data tibble or data.frame with observed data for individual
@@ -263,13 +265,24 @@ handle_sample_weighting <- function(
   obs_data,
   iterations,
   incremental,
-  i
+  i,
+  weights = NULL
 ) {
-  weights <- rep(0, nrow(obs_data))
+  binary_weights <- rep(0, nrow(obs_data))
   if(incremental) { # just fit current sample or group
-    weights[obs_data[["_grouper"]] %in% iterations[i]] <- 1
+    binary_weights[obs_data[["_grouper"]] %in% iterations[i]] <- 1
   } else { # fit all samples up until current sample
-    weights[obs_data[["_grouper"]] %in% iterations[1:i]] <- 1
+    binary_weights[obs_data[["_grouper"]] %in% iterations[1:i]] <- 1
+  }
+  if (!is.null(weights)) {
+    active_idx <- which(binary_weights == 1)
+    scheme_weights <- calculate_fit_weights(
+      weights = weights,
+      t = obs_data$t[active_idx]
+    )
+    if (!is.null(scheme_weights)) {
+      binary_weights[active_idx] <- scheme_weights
+    }
   }
-  weights
+  binary_weights
 }