Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Get initial estimates and tweak initial theta estimates #646

Merged
merged 11 commits into from
Feb 7, 2024

Conversation

barrettk
Copy link
Collaborator

@barrettk barrettk commented Jan 31, 2024

See previous discussions here: #635

Getting initial parameter estimates

The goal was to return an object that looked similar to param_estimates(). The internal function get_initial_est(), returns all values (i.e. the full diagonally concatenated OMEGA matrix if multiple records are provided) in a list format, whereas the wrapper initial_estimates(), only displays values defined in the control stream file

`initial_estimates()` example
MOD1 <- read_model(
 system.file("model", "nonmem", "basic", "1", package = "bbr")
)

> initial_estimates(MOD1, flag_fixed = TRUE)
# A tibble: 8 × 7
  parameter_names  init lower_bound upper_bound fixed record_type record_number
  <chr>           <dbl>       <dbl>       <dbl> <lgl> <chr>               <int>
1 THETA(1)         2              0          NA FALSE theta                   1
2 THETA(2)         3              0          NA FALSE theta                   1
3 THETA(3)        10              0          NA FALSE theta                   1
4 THETA(4)         0.02          NA          NA FALSE theta                   1
5 THETA(5)         1             NA          NA FALSE theta                   1
6 OMEGA(1,1)       0.05          NA          NA FALSE omega                   1
7 OMEGA(2,2)       0.2           NA          NA FALSE omega                   1
8 SIGMA(1,1)       1             NA          NA TRUE  sigma                   1

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")

# Run
> attr(initial_est, "omega_mat")
     [,1] [,2]
[1,] 0.05   NA
[2,]   NA  0.2

Tweaking initial estimates

Usage

MOD1 <- read_model(
    system.file("model", "nonmem", "basic", "1", package = "bbr")
)

# Base usage
mod2 <- copy_model_from(MOD1, "mod2") %>%
  tweak_initial_estimates(.p = 0.2)

# This function may be paired with `inherit_param_estimates()`:
mod2 <- copy_model_from(MOD1, "mod2") %>%
  inherit_param_estimates() %>% tweak_initial_estimates(.p = 0.2)
Walkthrough

Starting Record

;Initial THETAs
$THETA
(0,,1) ;[LCLM2]
( 0.7 FIX) ;[LCLM]
(0.67, 0.7, 0.72) ;[LCLF]
( 2 )   ;[CLAM]
( 2.0);[CLAF]
( 0.7 ) ;[LV1M]
( 0.7 ) ;[LV1F]
( 2.0 )   ;[V1AM]
( 2.0 )   ;[V1AF]
( 0.7 ) ;[MU_3]
(  0.7 );[MU_4]
( 0.3 )     ;[SDSL]

Tweak values

tweak_initial_estimates(
    .mod
    .p = 0.1,
    digits = 2
)

After Tweaking

;Initial THETAs
$THETA
(0,,1) ;[LCLM2] # If no estimate is found, treat as fixed (skip)
( 0.7 FIX) ;[LCLM] # skip fixed parameters
(0.67, 0.72, 0.72) ;[LCLF] # If tweaked value would be outside of one of the bounds, set to the bound
( 1.9 )   ;[CLAM] # Otherwise tweak the value by the sampled percentage
( 2.1);[CLAF]
( 0.76 ) ;[LV1M]
( 0.67 ) ;[LV1F]
( 2.2 )   ;[V1AM]
( 1.8 )   ;[V1AF]
( 0.73 ) ;[MU_3]
(  0.68 );[MU_4]
( 0.27 )     ;[SDSL]

closes #632

@barrettk barrettk changed the title Feat/tweak initial estimates Get initial estimates and tweak initial theta estimates Jan 31, 2024
 - Returns the initial estimates as defined in the control stream file
 - additional metadata, including which record the estimate was defined in, are also appended
 - only theta records are currently supported, but omega and sigma will be added in a future PR
@barrettk barrettk force-pushed the feat/tweak_initial_estimates branch from 8211f79 to c6f3ef0 Compare February 1, 2024 16:13
@barrettk barrettk requested a review from seth127 February 1, 2024 16:16
seth127
seth127 previously requested changes Feb 2, 2024
Copy link
Collaborator

@seth127 seth127 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is looking great overall. I have a bunch of comments, but they're almost all very small things.

The one thing that I want to look into deeper next week is the adjust_theta_bounds part. See my comment on this part of the docs above

#' --> set to value within bounds (first iteration: `(1-0.5)*0.9)`: (`(0, 0.45, 1)`)

I think we're actually good, but I want to walk through this with you to make sure we're on the same page with the algorithm that's been implemented and the test cases for it.

initial_est_df <- initial_est_df %>% dplyr::filter(!is.na(.data$init))

# Add original matrices as attributes if needed
attr(initial_est_df, "omega_mat") <- suppressWarnings(get_initial_est(.mod)) %>%
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor thing here, but why do we call get_initial_est() again here (and two lines later)? Is there a reason we can't just pluck from the initial_est that we get on line 28 above?

I'm also a little concerned with the supressWarnings() wrapper. What's that about? Is there a specific warning we're suppressing?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I mostly did it as to only attach the matrix in this case (and not all the other attributes). I can rework this as to decrease the number of calls. The suppressWarnings was just to suppress prior warnings, that would have already been displayed during the first call.

#'
#' @inheritParams initial_estimates
#'
#' @keywords internal
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Small thing: let's start moving to using @noRd instead of @keywords internal. See message in fd3dedf for more context.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That sounds good for here. However, in some cases I do like generating helper docs for internal functions, in the event it makes it easier to use or recall how to specify things. Do you have any thoughts on still generating docs for some internal functions?

 - removed one redundant test (the second one covered more, so delete the first)
 - better commenting for added clarity
 - order initial estimates based on row, to reflect how `param_estimates()` works
 - add test for same blocks, and expecting a warning when old non-specific priors are used
 - informative --> specific
seth127
seth127 previously requested changes Feb 7, 2024
Copy link
Collaborator

@seth127 seth127 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is looking great. A couple tiny clean up comments, but I think we're good to merge very soon.


it("theta bounds - adjust_tweaked_theta()", {
# Note: the comments next to each theta record do not necessarily need to
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I follow (and appreciate) the rest of this comment, but the first sentence...

the comments next to each theta record do not necessarily need to be true.

I'm not really sure what you're getting at here. I think we should probably remove this sentence, or reword to make it clearer.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comments I had next to each record were originally true when I was setting the seed and used 3 digits. After making that test a loop those comments were no longer true. I thought they could add some clarity to how the function was intended to work, but I agree it sounds confusing/misleading. Opted to remove them.

 - Add additional checks when testing the parsing of non-specific priors
 - removed some comments and added clarity to tests
Copy link
Collaborator

@seth127 seth127 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great work here.

@barrettk barrettk merged commit 5e48bf5 into main Feb 7, 2024
4 checks passed
@barrettk barrettk deleted the feat/tweak_initial_estimates branch February 28, 2024 18:18
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

New Feature: jitter/tweak initial parameter estimates
2 participants