Get initial parameter estimates and tweak initial parameter estimates

NAMESPACE

@@ -115,6 +115,7 @@ export(get_sigma)
 export(get_theta)
 export(get_yaml_path)
 export(inherit_param_estimates)
+export(initial_estimates)
 export(mod_ext)
 export(model_diff)
 export(model_summaries)
@@ -164,6 +165,7 @@ export(tags_diff)
 export(tail_lst)
 export(tail_output)
 export(test_threads)
+export(tweak_initial_estimates)
 export(update_model_id)
 export(use_bbi)
 export(wait_for_nonmem)

R/initial-estimates.R

@@ -0,0 +1,159 @@
+#' Retrieve and format initial parameter estimates
+#'
+#' @param .mod a `bbr` model object
+#' @param flag_fixed Logical (`TRUE`/`FALSE`). If `TRUE`, return which
+#'        parameters are fixed.
+#'
+#' @details
+#' `NA` values indicate that they were not specified in the control stream file.
+#' This is true for `THETA` bounds and whether a given parameter is `FIXED` or not.
+#'
+#' If you would like to format `OMEGA` or `SIGMA` records as *full matrices*,
+#' they are stored as attributes:
+#' ```
+#'  initial_est <- initial_estimates(.mod)
+#'  attr(initial_est, "omega_mat")
+#'  attr(initial_est, "sigma_mat")
+#' ```
+#'
+#' @examples
+#' mod1 <- read_model(
+#'   system.file("model", "nonmem", "basic", "1", package = "bbr")
+#' )
+#'
+#' initial_estimates(mod1)
+#'
+#' @export
+initial_estimates <- function(.mod, flag_fixed = FALSE){
+  initial_est <- get_initial_est(.mod, flag_fixed = flag_fixed)
+
+  # THETA labels
+  theta_inits <- initial_est$thetas
+  theta_inits$parameter_names <- sprintf("THETA(%d)", 1:nrow(theta_inits))
+  theta_inits <- theta_inits %>%
+    dplyr::rename(lower_bound = low, upper_bound = up)
+
+  # Format OMEGA
+  omega_inits <- matrix_to_df(initial_est$omegas, type = "omega")
+  if(isTRUE(flag_fixed)){
+    omega_inits <- omega_inits %>% dplyr::left_join(
+      matrix_to_df(attr(initial_est$omegas, "nmrec_flags")$fixed, type = "omega") %>%
+        dplyr::rename(fixed=init),
+      by = c("parameter_names")
+    )
+  }
+
+  # Format SIGMA
+  sigma_inits <- matrix_to_df(initial_est$sigmas, type = "sigma")
+  if(isTRUE(flag_fixed)){
+    sigma_inits <- sigma_inits %>% dplyr::left_join(
+      matrix_to_df(attr(initial_est$sigmas, "nmrec_flags")$fixed, type = "sigma") %>%
+        dplyr::rename(fixed=init),
+      by = c("parameter_names")
+    )
+  }
+
+  # Combine
+  initial_est_df <- dplyr::bind_rows(
+    theta_inits %>% dplyr::mutate(record_type = "theta"),
+    omega_inits %>% dplyr::mutate(record_type = "omega"),
+    sigma_inits %>% dplyr::mutate(record_type = "sigma")
+  ) %>% dplyr::relocate(parameter_names, record_type)
+
+  # Filter out NA values, as these were not specified in the control stream file
+  initial_est_df <- initial_est_df %>% dplyr::filter(!is.na(init))
+
+  # Add original matrices as attributes if needed
+  attr(initial_est_df, "omega_mat") <- get_initial_est(.mod) %>% purrr::pluck("omegas")
+  attr(initial_est_df, "sigma_mat") <- get_initial_est(.mod) %>% purrr::pluck("sigmas")
+
+  return(initial_est_df)
+}
+
+
+#' Retrieve initial parameter estimates
+#'
+#' @inheritParams initial_estimates
+#'
+#' @keywords internal
+get_initial_est <- function(.mod, flag_fixed = FALSE){
+
+  test_nmrec_version(.min_version = "0.4.0")
+  check_model_object(.mod, "bbi_nonmem_model")
+
+
+  ctl <- nmrec::read_ctl(get_model_path(.mod))
+
+  # Set flags
+  mark_flags <- if(isTRUE(flag_fixed)) "fix" else NULL
+
+  # Get initial estimates
+  theta_inits <- get_theta_inits(ctl, mark_flags = mark_flags)
+  omega_inits <- nmrec::get_omega(ctl, mark_flags = mark_flags)
+  sigma_inits <- nmrec::get_sigma(ctl, mark_flags = mark_flags)
+
+  # Format as list
+  inits <- list(
+    thetas = theta_inits,
+    omegas = omega_inits,
+    sigmas = sigma_inits
+  )
+
+  return(inits)
+}
+
+
+#' Helper function for getting the initial theta values, including bounds and
+#' the presence of `FIXED` parameters.
+#'
+#' @param ctl an `nmrec` ctl object.
+#' @param mark_flags TODO
+#'
+#' @keywords internal
+get_theta_inits <- function(ctl, mark_flags = NULL){
+
+  # Tabulate initial values and bounds
+  theta_inits <- purrr::map_dfc(c("init", "low", "up"), function(type){
+    theta_inits <- nmrec::get_theta(ctl, type = type)
+    tibble::tibble(!!rlang::sym(type) := theta_inits)
+  })
+
+  # Tabulate presence of FIXED parameters
+  theta_inits_fixed <- nmrec::get_theta(ctl, mark_flags = mark_flags)
+  theta_inits$fixed <- attr(theta_inits_fixed, "nmrec_flags")$fixed
+
+  return(theta_inits)
+}
+
+
+
+
+#' Format OMEGA or SIGMA matrix as data.frame
+#'
+#' @param mat a matrix
+#' @param type one of `c("omega","sigma")`
+#'
+#' @keywords internal
+matrix_to_df <- function(mat, type = c("omega","sigma")) {
+  type <- toupper(match.arg(type))
+
+  n <- nrow(mat)
+  indices <- lower.tri(mat, diag = TRUE)
+
+  row_indices <- row(indices)[indices]
+  col_indices <- col(indices)[indices]
+
+  parameter_names <- sprintf(paste0(type,"(%d,%d)"), row_indices, col_indices)
+
+  values <- mat[indices]
+
+  result_df <- tibble::tibble(
+    parameter_names = parameter_names,
+    init = values
+  )
+
+  return(result_df)
+}
+
+
+

R/tweak-initial-estimates.R

@@ -0,0 +1,162 @@
+
+#' Tweak the initial parameter estimates
+#'
+#' @param .mod model object to update.
+#' @param .p Percent to tweak the initial parameter estimates by. Represented as
+#'   a fraction.
+#' @param tweak type of estimates to tweak in the model. Defaults to
+#'   updating all of THETA, SIGMA, and OMEGA records.
+#' @param digits Number of significant digits to round estimates to.
+#'
+#' @details
+#'
+#' In the following cases, the initial estimate will *not* be updated:
+#'  - **Individual** `FIXED` `THETA` parameters
+#'     - e.g., `$THETA 1.2 FIX 1.5 0.2` --> would only skip the first value
+#'  - **Individual** `FIXED` `OMEGA` & `SIGMA` parameters for *diagonal* matrices
+#'     - e.g., `$OMEGA 0 FIX 1` --> would only skip the first value
+#'  - **Full** `FIXED` `OMEGA` & `SIGMA` *block* matrices
+#'     - e.g., `$OMEGA BLOCK(2) 0.1 0.001 0.1 FIX` --> would skip the full `OMEGA` record
+#'  - `THETA` parameters with no initial estimate
+#'     - e.g., `$THETA (0,,1)`
+#'
+#'  For bounded `THETA` estimates:
+#'   - If an initial `THETA` has bounds **and** an initial estimate
+#'  (e.g., `(0, 0.5, 1)`, `(0,1)`), the bounds will be respected when sampling
+#'  a percent to tweak by. If the tweaked value would fall below the lower bound,
+#'  the initial estimate will be set to the lower bound. The same is true for
+#'  upper bounds.
+#'     - e.g., `(0, 0.5, 1)` --> tweak initially, falls outside bound (`(0, 1.2, 1)`)
+#'     --> set to upper bound (`(0, 1, 1)`)
+#'
+#'
+#' @examples
+#' \dontrun{
+#' base_mod <- read_model(file.path(MODEL_DIR, "1"))
+#'
+#' mod2 <- copy_model_from(base_mod, "mod2") %>%
+#'   tweak_initial_estimates(.p = 0.2)
+#'
+#' # This function may be paired with `inherit_param_estimates()`:
+#' mod2 <- copy_model_from(base_mod, "mod2") %>%
+#'   inherit_param_estimates() %>% tweak_initial_estimates(.p = 0.2)
+#'
+#' # If you want to set the seed for reproducible results:
+#' mod2 <- withr::with_seed(1234, {
+#'    tweak_initial_estimates(mod2, .p = 0.2, digits = 3)
+#' })
+#'
+#' # The utilized `.Random.seed` is appended to the model object as an attribute:
+#' head(attr(mod2, "tweak_estimates.seed"))
+#' }
+#' @export
+tweak_initial_estimates <- function(
+    .mod,
+    .p = 0.1,
+    tweak = c("theta", "sigma", "omega"),
+    digits = 3
+){
+
+  # Initialize .Random.seed
+  # (.Random.seed does not exist until some randomness is used in R session)
+  if(!exists(".Random.seed")) set.seed(NULL)
+
+  # Assertions
+  test_nmrec_version(.min_version = "0.4.0")
+  check_model_object(.mod, "bbi_nonmem_model")
+  checkmate::assert_true(all(tweak %in% BBR_ESTIMATES_INHERIT))
+
+  # Get initial estimates
+  initial_est <- get_initial_est(.mod, flag_fixed = TRUE)
+
+  # Rounding
+  fmt_digits <- paste0("%.",digits,"G")
+
+  # Tweak THETA
+  init_thetas <- initial_est$thetas
+  init_thetas$new <- withr::with_preserve_seed(tweak_values(init_thetas$init, .p))
+
+  # Respect THETA bounds (if any) and FIXED estimates
+  init_thetas <- init_thetas %>% mutate(
+    new = dplyr::case_when(
+      fixed == TRUE ~ init,
+      !is.na(up) & new > up ~ up,
+      !is.na(low) & new < low ~ low,
+      TRUE ~ new
+    )
+  )
+  # Round
+  new_thetas <- init_thetas$new
+  new_thetas[!is.na(new_thetas)] <- new_thetas[!is.na(new_thetas)] %>%
+    signif(digits)
+
+
+  # Tweak OMEGA
+  init_omegas <- initial_est$omegas
+  fixed_omegas <- attr(init_omegas, "nmrec_flags")$fixed
+  fixed_omegas[is.na(fixed_omegas)] <- TRUE
+  new_omegas <- init_omegas[!fixed_omegas]
+  # Tweak values
+  new_omegas <- withr::with_preserve_seed(tweak_values(new_omegas, .p))
+  init_omegas[!fixed_omegas] <- new_omegas %>% signif(digits)
+
+
+  # Tweak SIGMA
+  init_sigmas <- initial_est$sigmas
+  fixed_sigmas <- attr(init_sigmas, "nmrec_flags")$fixed
+  fixed_sigmas[is.na(fixed_sigmas)] <- TRUE
+  new_sigmas <- init_sigmas[!fixed_sigmas]
+  # Tweak values
+  new_sigmas <- withr::with_preserve_seed(tweak_values(new_sigmas, .p))
+  init_sigmas[!fixed_sigmas] <- new_sigmas %>% signif(digits)
+
+
+  # nmrec model objects
+  mod_path <- get_model_path(.mod)
+  mod_lines <- nmrec::read_ctl(mod_path)
+
+  # Update THETA Block
+  if("theta" %in% tweak && !rlang::is_empty(new_thetas)){
+    nmrec::set_theta(
+      mod_lines, values = new_thetas, bounds = "keep",
+      fmt = fmt_digits
+    )
+  }
+
+  # Update OMEGA Block
+  if("omega" %in% tweak && !rlang::is_empty(init_omegas)){
+    nmrec::set_omega(
+      mod_lines, values = init_omegas, representation = "reset",
+      fmt = fmt_digits
+    )
+  }
+
+  # Update SIGMA Block
+  if("sigma" %in% tweak && !rlang::is_empty(init_sigmas)){
+    nmrec::set_sigma(
+      mod_lines, values = init_sigmas, representation = "reset",
+      fmt = fmt_digits
+    )
+  }
+
+  # Write out mod_lines to model
+  nmrec::write_ctl(mod_lines, mod_path)
+
+  # Capture seed for traceability
+  seed <- .Random.seed
+  attr(.mod, "tweak_estimates.seed") <- seed
+
+  return(.mod)
+}
+
+
+tweak_values <- function(values, .p){
+
+  # Sample percentages & Preserve seed if set
+  tweak_perc <- stats::runif(length(values), -.p, .p)
+
+  # Return new values
+  new_values <- values + (values * tweak_perc)
+
+  return(new_values)
+}

inst/model/nonmem/basic/mod_tweak.ctl

@@ -0,0 +1,47 @@
+$PROBLEM From bbr: see mod_tweak.yaml for details
+
+$INPUT ID TIME MDV EVID DV AMT  SEX WT ETN NUM
+$DATA ../../../../extdata/acop.csv IGNORE=@
+IGNORE(ID.EQ.2)
+$SUBROUTINES ADVAN2 TRANS2
+
+$PK
+ET=1
+IF(ETN.EQ.3) ET=1.3
+KA = THETA(1)
+CL = THETA(2)*((WT/70)**0.75)* EXP(ETA(1))
+V = THETA(3)*EXP(ETA(2))
+SC=V
+
+
+$THETA
+(0, 2)  ; KA
+(0, 3)  ; CL
+(0, 10) ; V2
+(0.02)  ; RUVp
+(1)     ; RUVa
+
+$OMEGA
+0.05    ; iiv CL
+0.2     ; iiv V2
+
+$SIGMA
+1 FIX
+
+$ERROR
+IPRED = F
+IRES = DV-IPRED
+W = IPRED*THETA(4) + THETA(5)
+IF (W.EQ.0) W = 1
+IWRES = IRES/W
+Y= IPRED+W*ERR(1)
+FAKE=1 ; for testing
+
+$EST METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT POSTHOC
+$COV
+$TABLE NUM IPRED NPDE CWRES NOPRINT ONEHEADER FILE=1.tab
+$TABLE NUM CL V KA ETAS(1:LAST) NOAPPEND NOPRINT ONEHEADER FILE=1par.tab
+$TABLE ID CL NOAPPEND NOPRINT ONEHEADER FIRSTONLY FILE=1first1.tab ; for testing firstonly
+$TABLE NUM CL NOAPPEND NOPRINT ONEHEADER FIRSTONLY FILE=1first2.tab ; for testing firstonly
+$TABLE NUM ID CL NOAPPEND NOPRINT ONEHEADER FIRSTONLY FILE=1first3.tab ; for testing firstonly
+$TABLE NUM CL KA FAKE NOAPPEND NOPRINT ONEHEADER FILE=1dups.tab ; for testing dropping duplicate columns

inst/model/nonmem/basic/mod_tweak.yaml

@@ -0,0 +1,7 @@
+model_type: nonmem
+description: Tweak estimates
+bbi_args:
+  overwrite: yes
+  threads: 4.0
+based_on:
+- '1'