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 parameter estimates and tweak initial parameter estimates #635

Closed
wants to merge 19 commits into from

Conversation

barrettk
Copy link
Collaborator

@barrettk barrettk commented Jan 10, 2024

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: 9 × 6
  parameter_names record_type  init lower_bound upper_bound fixed
  <chr>           <chr>       <dbl>       <dbl>       <dbl> <lgl>
1 THETA(1)        theta        2              0          NA FALSE
2 THETA(2)        theta        3              0          NA FALSE
3 THETA(3)        theta       10              0          NA FALSE
4 THETA(4)        theta        0.02          NA          NA FALSE
5 THETA(5)        theta        1             NA          NA FALSE
6 OMEGA(1,1)      omega        0.05          NA          NA FALSE
8 OMEGA(2,2)      omega        0.2           NA          NA FALSE
9 SIGMA(1,1)      sigma        1             NA          NA TRUE 

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

The details below are also found in the function docs:

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

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)


# If you want to set the seed for reproducible results:
mod2 <- withr::with_seed(1234, {
  copy_model_from(MOD1, "mod2") %>%
    tweak_initial_estimates(.p = 0.2, digits = 3)
})

# The utilized `.Random.seed` is appended to the model object as an attribute:
head(attr(mod2, "tweak_estimates.seed"))
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,
    tweak = c("theta", "sigma", "omega"),
    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

…ith other relevant information

 - The information returned from this function will be needed for the new tweak/jitter function. It currently relies on internal `nmrec` functions, and likely needs to be reworked to be more readable.
 - I think a function like this belongs in `nmrec`, but wanted to develop a prototype in `bbr` to facilitate the conversation in terms of what information I need
 - In the end, we want to export a function *similar* to what I have so far, but with improved formatting; perhaps as a table, similar to `param_estimates()`. The `bbr` function could be a wrapper of the `nmrec` function we develop, but I think the core work should be done in `nmrec`
@barrettk
Copy link
Collaborator Author

barrettk commented Jan 10, 2024

Examples so far:

Setup

source("/data/Projects/package_dev/bbr/tests/testthat/setup-workflow-ref.R", echo=FALSE)
MODEL_DIR_C <- system.file("model", "nonmem", "complex",   package = "bbr")

Example models

These are the examples that have been tested thus far

.mod <- MOD1
.mod <- read_model(file.path(MODEL_DIR_C, "1001"))
.mod <- read_model(file.path(MODEL_DIR_C, "example2_saemimp"))

Example output

Base Model. Bounds must be tabulated so that we can respect them when tweaking the initial estimates (see conversation here)

MOD1
> .mod <- MOD1
> get_param_inits(.mod)
$theta
# A tibble: 5 × 3
  low   init  index
  <chr> <chr> <int>
1 0     2         1
2 0     3         2
3 0     10        3
4 NA    0.02      4
5 NA    1         5

$omega
$omega$matrices
$omega$matrices[[1]]
     [,1] [,2]
[1,] 0.05  0.0
[2,] 0.00  0.2


$omega$fixed
$omega$fixed[[1]]
[1] FALSE



$sigma
$sigma$matrices
$sigma$matrices[[1]]
     [,1]
[1,]    1


$sigma$fixed
$sigma$fixed[[1]]
[1] TRUE

Model with multiple records: here we create a combined matrix, which will facilitate passing jittered values to set_omega(). If one of the matrices was FIXED, we would need additional logic for not overwriting those values

1001
> .mod <- read_model(file.path(MODEL_DIR_C, "1001"))
> get_param_inits(.mod)
$theta
# A tibble: 6 × 2
  init               index
  <chr>              <int>
1 -0.574219127214817     1
2 5.67154052944669       2
3 1.03568637508119       3
4 0.241769673098703      4
5 -11.3632385211334      5
6 5.11272713931817       6

$omega
$omega$matrices
$omega$matrices[[1]]
     [,1]
[1,] 0.09

$omega$matrices[[2]]
     [,1] [,2] [,3]
[1,] 0.09 0.00 0.00
[2,] 0.01 0.09 0.00
[3,] 0.01 0.01 0.09


$omega$fixed
$omega$fixed[[1]]
[1] FALSE

$omega$fixed[[2]]
[1] FALSE


$omega$full_matrix
     [,1] [,2] [,3] [,4]
[1,] 0.09 0.00 0.00 0.00
[2,] 0.00 0.09 0.00 0.00
[3,] 0.00 0.01 0.09 0.00
[4,] 0.00 0.01 0.01 0.09


$sigma
$sigma$matrices
$sigma$matrices[[1]]
     [,1] [,2]
[1,] 0.04 0.00
[2,] 0.00 0.04


$sigma$fixed
$sigma$fixed[[1]]
[1] FALSE

Note: this model as prior blocks, so the output here is technically incorrect, or at least lacking in terms of differentiating between actual records and priors. Additional logic would be needed to filter these out or error if the expected length wasnt found

example2_saemimp
> .mod <- read_model(file.path(MODEL_DIR_C, "example2_saemimp"))
> get_param_inits(.mod)
$theta
# A tibble: 12 × 3
   init  fixed index
   <chr> <chr> <int>
 1 0.7   NA        1
 2 0.7   NA        2
 3 2     NA        3
 4 2.0   NA        4
 5 0.7   NA        5
 6 0.7   NA        6
 7 2.0   NA        7
 8 2.0   NA        8
 9 0.7   NA        9
10 0.7   NA       10
11 0.3   NA       11
12 4     TRUE     12

$omega
$omega$matrices
$omega$matrices[[1]]
      [,1]  [,2]  [,3] [,4]
[1,] 0.500 0.000 0.000  0.0
[2,] 0.001 0.500 0.000  0.0
[3,] 0.001 0.001 0.500  0.0
[4,] 0.001 0.001 0.001  0.5

$omega$matrices[[2]]
     [,1] [,2] [,3] [,4]
[1,] 0.01 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00
[3,] 0.00 0.00 0.01 0.00
[4,] 0.00 0.00 0.00 0.01


$omega$fixed
$omega$fixed[[1]]
[1] FALSE

$omega$fixed[[2]]
[1] TRUE


$omega$full_matrix
      [,1]  [,2]  [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.500 0.000 0.000  0.0 0.00 0.00 0.00 0.00
[2,] 0.001 0.500 0.000  0.0 0.00 0.00 0.00 0.00
[3,] 0.001 0.001 0.500  0.0 0.00 0.00 0.00 0.00
[4,] 0.001 0.001 0.001  0.5 0.00 0.00 0.00 0.00
[5,] 0.000 0.000 0.000  0.0 0.01 0.00 0.00 0.00
[6,] 0.000 0.000 0.000  0.0 0.00 0.01 0.00 0.00
[7,] 0.000 0.000 0.000  0.0 0.00 0.00 0.01 0.00
[8,] 0.000 0.000 0.000  0.0 0.00 0.00 0.00 0.01


$sigma
$sigma$matrices
$sigma$matrices[[1]]
     [,1]
[1,]    1


$sigma$fixed
$sigma$fixed[[1]]
[1] TRUE

 - includes a refactored get_param_inits(), which now utilizes new WIP `nmrec` `get_{theta,omega,sigma}` functions.
 - Prototype supports all of the cases we have mentioned so far. Have manually tested other cases brought up in passing.
 - Code could be improved a bit to be more readable; perhaps a few helper functions could be used.
 - docs and tests are currently lacking
@barrettk
Copy link
Collaborator Author

barrettk commented Jan 11, 2024

Some examples with the most recent refactor:

.mod <- read_model(file.path(MODEL_DIR_C, "example2_saemimp"))

Getting initial estimates

Internal function examples

Base case

> get_initial_est(.mod)
$thetas
# A tibble: 13 × 3
    init   low    up
   <dbl> <dbl> <dbl>
 1  NA    0     1   
 2   0.7 NA    NA   
 3   0.7  0.67  0.72
 4   2   NA    NA   
 5   2   NA    NA   
 6   0.7 NA    NA   
 7   0.7 NA    NA   
 8   2   NA    NA   
 9   2   NA    NA   
10   0.7 NA    NA   
11   0.7 NA    NA   
12   0.3 NA    NA   
13   4   NA    NA   

$omegas
      [,1]  [,2]  [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.500    NA    NA   NA   NA   NA   NA   NA
[2,] 0.001 0.500    NA   NA   NA   NA   NA   NA
[3,] 0.001 0.001 0.500   NA   NA   NA   NA   NA
[4,] 0.001 0.001 0.001  0.5   NA   NA   NA   NA
[5,]    NA    NA    NA   NA 0.01   NA   NA   NA
[6,]    NA    NA    NA   NA 0.00 0.01   NA   NA
[7,]    NA    NA    NA   NA 0.00 0.00 0.01   NA
[8,]    NA    NA    NA   NA 0.00 0.00 0.00 0.01

$sigmas
     [,1]
[1,]    1

Flag fixed parameters

> get_initial_est(.mod, flag_fixed = TRUE)
$thetas
# A tibble: 13 × 4
    init   low    up fixed
   <dbl> <dbl> <dbl> <lgl>
 1  NA    0     1    FALSE
 2   0.7 NA    NA    TRUE 
 3   0.7  0.67  0.72 FALSE
 4   2   NA    NA    FALSE
 5   2   NA    NA    FALSE
 6   0.7 NA    NA    FALSE
 7   0.7 NA    NA    FALSE
 8   2   NA    NA    FALSE
 9   2   NA    NA    FALSE
10   0.7 NA    NA    FALSE
11   0.7 NA    NA    FALSE
12   0.3 NA    NA    FALSE
13   4   NA    NA    TRUE 

$omegas
      [,1]  [,2]  [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.500    NA    NA   NA   NA   NA   NA   NA
[2,] 0.001 0.500    NA   NA   NA   NA   NA   NA
[3,] 0.001 0.001 0.500   NA   NA   NA   NA   NA
[4,] 0.001 0.001 0.001  0.5   NA   NA   NA   NA
[5,]    NA    NA    NA   NA 0.01   NA   NA   NA
[6,]    NA    NA    NA   NA 0.00 0.01   NA   NA
[7,]    NA    NA    NA   NA 0.00 0.00 0.01   NA
[8,]    NA    NA    NA   NA 0.00 0.00 0.00 0.01
attr(,"nmrec_flags")
attr(,"nmrec_flags")$fixed
      [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7] [,8]
[1,] FALSE    NA    NA    NA   NA   NA   NA   NA
[2,] FALSE FALSE    NA    NA   NA   NA   NA   NA
[3,] FALSE FALSE FALSE    NA   NA   NA   NA   NA
[4,] FALSE FALSE FALSE FALSE   NA   NA   NA   NA
[5,]    NA    NA    NA    NA TRUE   NA   NA   NA
[6,]    NA    NA    NA    NA TRUE TRUE   NA   NA
[7,]    NA    NA    NA    NA TRUE TRUE TRUE   NA
[8,]    NA    NA    NA    NA TRUE TRUE TRUE TRUE


$sigmas
     [,1]
[1,]    1
attr(,"nmrec_flags")
attr(,"nmrec_flags")$fixed
     [,1]
[1,] TRUE

Exported function: initial_estimates()

Here I was thinking about filtering out initial values with NA, as this indicates the value was not specified in the control stream file, and would likely reflect the user's expectation more appropriately. Edit: decided to do this for now

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

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

Example

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,
    tweak = c("theta", "sigma", "omega"),
    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]

 - add tests for getting initial estimates
 - update some docs
 - need to put more thought into settting up template models. Couldnt use the approach I initially had in mind
 - did a lot of testing with setting seeds both within the function (mostly via `withr`), and set outside. I think it's pretty robust now, but would like a second pair of eyes to ensure nothing else should be done

 - Pulled out random sampling into a separate helper function. Originally with the intention of inheriting the seed of the parent environment, though this did not work as desired. Decided to leave it as a separate function for now for code readability, and in the event we pass a `seed` argument within these functions
@barrettk barrettk requested a review from seth127 January 16, 2024 21:39
@barrettk
Copy link
Collaborator Author

barrettk commented Jan 16, 2024

Hey @seth127, I just marked you for an intermediate review. Some things that still need to be done or talked about:

Updates

  • Update the required nmrec version in a few functions/tests, along with the drone yaml, once a new version of nmrec is released
    • This may require a slight refactor if anything changes on the nmrec side
    • Wanted to talk about whether it made sense to update drone vs wait till we migrated to github actions.

Discussions

  • I didnt right up any tests for the initial_estimates wrapper given the nature of it, but let me know if you think I should
  • Let me know how you feel about the current tests for getting the initial estimates (pretty minimal given that nmrec will handle most operations)
  • Discussion regarding the test framework for jittering the initial estimates
  • Any thoughts on setting a seed? (see my commit messages d169643 & a916af3)
  • Argument and function names

Im going to spend a little more time on documentation, but after that I was planning on waiting for a first review/some other parts to move along

Comment on lines 11 to 31
#' @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)`)
#'
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@kyleam when you have a minute, im curious if you have any pointers/opinions on these docs

Copy link
Collaborator

Choose a reason for hiding this comment

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

I know you asked Kyle, but these make sense to me.

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 good, just few minor comments from this review.

I also want to note that we still have two outstanding questions about the omega (and also sigma?) matrices:

  1. Do we tweak the off-diagonals (if they're defined in the control stream) or not?
  2. How should we handle checking for (and possibly forcing) positive definiteness?

Either way, as discussed with @barrettk offline, our current plan for that is:

  • Continue with this PR and merge it to main once review comments are addressed and the relevant nmrec changes have landed (add functions for getting inits nmrec#30)
  • @barrettk to create a new issue to specifically address the outstanding matrix questions mentioned above.
  • Address that new issue in a separate PR.
  • Plan to not make a bbr release until all of that ^ has landed on main.

Comment on lines 11 to 31
#' @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)`)
#'
Copy link
Collaborator

Choose a reason for hiding this comment

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

I know you asked Kyle, but these make sense to me.

@barrettk barrettk marked this pull request as ready for review January 19, 2024 18:28
@barrettk barrettk requested a review from seth127 January 19, 2024 18:41
@barrettk
Copy link
Collaborator Author

As @seth127 noted, we will need a separate PR to address these issues. The new tweak_initial_estimates() function works for THETA parameters, but is flawed in its handling of matrices. Still needed:

  • evaluate positive-definiteness of the tweaked matrix, and adjust if needed
  • In the event of a correlation matrix, the tweaked matrix must first be converted to a covariance matrix before evaluating positive-definiteness/making necessary adjustments, and then converting it back to the correct format (variance on diagonal, correlation on off-diagonal)

 - fix test in R 4.0
@barrettk
Copy link
Collaborator Author

@seth127 something I just thought about. This is related to a conversation @kyleam and I had regarding the handling of prior blocks that follow the old naming convention (see final thoughts).

Basically, prior blocks currently get "grouped" when creating the combined matrix. This would also occur for THETA records:

Example
> MODEL_DIR_X
[1] "inst/model/nonmem/complex"
> .mod <- read_model(file.path(MODEL_DIR_X, "example2_saemimp"))
ctl <- nmrec::read_ctl(get_model_path(.mod))
> nmrec::extract_omega(ctl)
      [,1]  [,2]  [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.500    NA    NA   NA   NA   NA   NA   NA
[2,] 0.001 0.500    NA   NA   NA   NA   NA   NA
[3,] 0.001 0.001 0.500   NA   NA   NA   NA   NA
[4,] 0.001 0.001 0.001  0.5   NA   NA   NA   NA
[5,]    NA    NA    NA   NA 0.01   NA   NA   NA
[6,]    NA    NA    NA   NA 0.00 0.01   NA   NA
[7,]    NA    NA    NA   NA 0.00 0.00 0.01   NA
[8,]    NA    NA    NA   NA 0.00 0.00 0.00 0.01

In the $OMEGA records below, only the first one is an actual record:

;Initial OMEGAs
$OMEGA BLOCK(4)
0.5  ;[p]
0.001  ;[f]
0.5  ;[p]
0.001 ;[f]
0.001 ;[f]
0.5  ;[p]
0.001 ;[f]
0.001 ;[f]
0.001 ;[f]
0.5 ;[p]

; Degrees of freedom to OMEGA prior matrix:
$THETA 4 FIX
; Prior OMEGA matrix
$OMEGA BLOCK(4)
0.01 FIX
0.0  0.01
0.0 0.0 0.01
0.0 0.0 0.0 0.01

I was reminded about this conversation when doing some preliminary testing for the positive-definiteness, but now im not sure if we should address this in this PR, given that it also implicates THETA records.

Im not really sure how to handle this. In the case of inherit_param_estimates(), Kyle compares the length of the passed-in values (from model_summary output) to the detected number of record values. We cant do that here though, as the length of the passed-in values is determined by the length of the fetched initial estimates. I dont see a way around parsing the PRIOR block or allowing priors to be jittered.

@kyleam
Copy link
Collaborator

kyleam commented Jan 23, 2024

Basically, prior blocks currently get "grouped" when creating the combined matrix.

My opinion for how to proceed is what I said in the issue you linked: we should document that using non-informative record names for priors is not supported with this feature.

@seth127
Copy link
Collaborator

seth127 commented Jan 24, 2024

@kyleam said

we should document that using non-informative record names for priors is not supported with this feature.

and I think I agree with that, in principle. However, in practice, we never know that they're using non-informative priors in this case.

Above, @barrettk said

[nmrec] compares the length of the passed-in values (from model_summary output) to the detected number of record values. We cant do that here though, as the length of the passed-in values is determined by the length of the fetched initial estimates.

I've looked through the code a bit and I'm pretty sure that's right. My evidence:

  • @kyleam I think you're saying we should just let this case fall through to this error
  • But that error only gets triggered when length(values) != len_expected
  • In the inherit_param_estimates() call, values is ultimately coming from a model_summary() on a previously executed model (e.g. here). At that point, we know they're all real parameters (and not a non-informative priors), and so the mismatch/error in set_param() gets triggered.
  • However, in the tweak_initial_estimates() call, the thing that is eventually passed to values is coming from the nmrec::extract_* functions (here, which leads to here). And since the extract_* don't differentiate between real parameters and non-informative priors... we never get the mismatch.
I also tried this empirically and saw the same thing:

I tried with example2_saemimp model, but first I had to remove the FIX from the final THETA (which is actually a prior) to verify that it would indeed get tweaked. And it does:

.mod <- read_model(file.path(MODEL_DIR_X, "example2_saemimp"))
tweak_initial_estimates(.mod)

results in this

Screenshot 2024-01-23 at 10 33 28 PM

Soooooo... that means we can't just leave it as is and expect that error to bubble up.

@seth127
Copy link
Collaborator

seth127 commented Jan 24, 2024

In effect, this means that the current state actually tweaks non-informative priors too. The big question from me: is it ok for us to tweak the non-informative priors the same way that we tweak the initial parameter estimates?

If so, then we can just document that behavior and move on.

If that's not ok... then I'm not sure what the next step is. I guess I'd just say "let's cross that bridge when we come to it" (and hope we don't come to it).

@curtisKJ @timwaterhouse could either of you weigh in on my "big question" at the top of this comment?

One final note (in favor of leaving as is): anecdotally, in most cases I've seen, the non-informative priors use FIX. In that case, they wouldn't be tweaked anyway. I'm not sure if that's "standard practice" or not, but I thought it was worth noting.

@kyleam
Copy link
Collaborator

kyleam commented Jan 24, 2024

Thanks @seth127. I haven't reviewed this PR or been involved in many of the design discussions, so it's helpful to see you expand on those details.

@kyleam I think you're saying we should just let this case fall through to this error

In the linked thread, I had assumed that would trigger, but for the latest comment here I took for granted @barrettk's claim that that doesn't get triggered. (Thank you both for pointing that out.)

So, in the comment above, I was just casting my vote for documenting that this feature is not meant to be used with uninformative prior names, accepting that it will jitter them if there (i.e. your "document that behavior and move on").

A few more thoughts:

  • IIUC the common workflow here is "inherit_param_estimates then tweak", in which case the first step would flag control streams that are using non-informative names. So in practice it seems to me that this problem would often be flagged.

  • If bbr wanted to give an error in most cases, could it abort if there is a $PRIOR record and not any informative record names? That'd hopefully not have false positives, and I think it might catch most cases, though it would miss some (e.g. subroutine used instead of prior record).

  • Or it's perhaps worth considering reworking things so that jittering is coupled to estimates from a previous model. You could still use nmrec::extract_* to get information to intersect, but the core init values could come from a previous model, and you could reliably error on a shape mismatch. However, I guess that's probably being avoided in anticipation for cases where "tweak estimates" might be used to blow out a series of models for assessing sensitivity of the results to the initial values, in which case requiring a finished model as the starting point would be a bit odd (but maybe not too bad in practice).

  • It's perhaps worth emphasizing in the docs that these are tools meant to help modelers modify the control stream, but they should review the result very carefully, as if they were editing the model by hand or reviewing changes a person was proposing. (I haven't taken a look at the documentation that exists, so apologies if something along those lines is already there.)

@seth127
Copy link
Collaborator

seth127 commented Jan 24, 2024

Thanks @kyleam I appreciate those thoughts. You said

IIUC the common workflow here is "inherit_param_estimates then tweak", in which case the first step would flag control streams that are using non-informative names. So in practice it seems to me that this problem would often be flagged.

I actually had this same thought this morning. That's reassuring, and I think will likely cover a lot of cases.

To be clear (for anyone else reading this), we're talking about a pattern like this, which several people have mentioned wanting to be able to do:

mod2 <- mod1 %>%
    copy_model_from() %>%
    inherit_param_estimates() %>%
    tweak_initial_estimates()

In that case, if mod1 was using non-informative priors, the error I mentioned above would be triggered by the inherit_param_estimates() call, before we even get to the tweak_initial_estimates() function that we're discussing here.


We're going to discuss with a few scientific folks offline. Unless they object, I'm leaning towards this

documenting that this feature is not meant to be used with uninformative prior names, accepting that it will jitter them if there

@seth127
Copy link
Collaborator

seth127 commented Jan 25, 2024

I had said

We're going to discuss with a few scientific folks offline. Unless they object, I'm leaning towards

documenting that this feature is not meant to be used with uninformative prior names, accepting that it will jitter them if there

Unfortunately, @timwaterhouse and @curtisKJ do object. Essentially, they said we definitely don't want to tweak the priors. Two things we're going to look into:

  1. Look at the NONMEM docs to verify whether you need to specify FIX for a parameter record to function as non-informative prior. If the answer is "yes", then this solves our problem, because we don't tweak those values anyway.
  2. Look into ways that we can reliably tell whether a given parameter record is an real parameter or a non-informative prior. Basically, restart this discussion about whether the $PRIOR block is sufficient and, if not, what else we would need to check.

@barrettk I believe you're taking this ^ on, right?

@kyleam
Copy link
Collaborator

kyleam commented Jan 25, 2024

Unfortunately, @timwaterhouse and @curtisKJ do object. Essentially, they said we definitely don't want to tweak the priors.

The proposal isn't that anyone wants to tweak them. It's just accepting that the caller can violate the documented use in a way that leads to them being tweaked. Was the conclusion "we're going to insist that the outcome is severe enough that we really want to prevent that mistake"?

Two things we're going to look into:

Perhaps it's worth considering rethinking the design to always use an executed model as the starting point (my "reworking things so that jittering is coupled to estimates from a previous model" bullet point above). That'd give you a reliable answer to what the shape should be.

@kyleam
Copy link
Collaborator

kyleam commented Jan 25, 2024

Just a few comments, in case it's helpful...

Look at the NONMEM docs to verify whether you need to specify FIX for a parameter record to function as non-informative prior. If the answer is "yes", then this solves our problem, because we don't tweak those values anyway.

At least based on the nwpri docs, it looks like that could be the case. I see a lot of "values on these records must be fixed", but perhaps a closer look would show cases where that's not said or perhaps there are other more relevant docs.

[...] about whether the $PRIOR block is sufficient

My conclusion that the prior block isn't sufficient is based on the nwpri docs saying "a $PRIOR record may be used instead of a user-written PRIOR subroutine, in which case the input arguments listed above may be specified as options of $PRIOR".

@barrettk
Copy link
Collaborator Author

barrettk commented Jan 26, 2024

In response to @seth127's comment:

Look at the NONMEM docs to verify whether you need to specify FIX for a parameter record to function as non-informative prior. If the answer is "yes", then this solves our problem, because we don't tweak those values anyway.

Seth and I talked offline about this. It is still something worth looking into, but we decided it isnt really a complete solution, as we want to be able to identify prior records in order to remove them from initial_estimates() output. This may not be totally feasible, but the motivation for this is strong, and it was determined that we should see what can be done here.

In response to @kyleam's comments:

The proposal isn't that anyone wants to tweak them. It's just accepting that the caller can violate the documented use in a way that leads to them being tweaked. Was the conclusion "we're going to insist that the outcome is severe enough that we really want to prevent that mistake"?

The conclusion is that we needed to guarantee priors wouldnt be altered, even using the old method. It's apparently a large deal so we need some method of determining if priors are being used at the minimum. Identifying the specific records that are priors is ideal, but we need to at least know if they're being used.

Perhaps it's worth considering rethinking the design to always use an executed model as the starting point (my "reworking things so that jittering is coupled to estimates from a previous model" bullet point above). That'd give you a reliable answer to what the shape should be.

The jittering feature has to be separate from the inherit_param_estimates() function unfortunately, as we want tweak_initial_estimates() to be usable on its own. Ideally we dont want to require it to be a previously run model. There have been discussions on creating a wrapper around tweak_initial_estimates() to quickly create a bunch of new models with different initial estimates.

  • These reasons dont technically prohibit what you're suggesting though. I think that is worth considering if we get stuck, though I think there is a strong preference not to require a previously run model. I'd like to hear @seth127's opinion on it, but I know we want to look into how feasible parsing the various ways of specifying priors would be. @curtisKJ suggested there was one really old method that few people were even aware of, but that one or two methods cover the vast majority of cases (one being the $PRIOR block, need to look more into the other). I think if we could look for the presence of those two blocks, and make as many inferences as we can (with good documentation on its limitations/warnings), that would be the most preferable option by far.
  • I think some of this logic could be worth revisiting, though I know it's far from complete (see how it's used here).

@kyleam
Copy link
Collaborator

kyleam commented Jan 26, 2024

I think if we could look for the presence of those two blocks, and make as many inferences as we can (with good documentation on its limitations/warnings), that would be the most preferable option by far.

It seems to me that that line of thinking is in direct contradiction to claiming you must "guarantee priors wouldnt be altered".

The jittering feature has to be separate from the inherit_param_estimates() function unfortunately, as we want tweak_initial_estimates() to be usable on its own

I understand that (I think, at least). My suggestion is that you reconsider that design decision due to the desire to reliably know if there is a shape mismatch due to prior use. If we decide not to go that direction, that's of course fine, but here it feels like you're just dismissing the suggestion by reasserting what's already a given.

@barrettk
Copy link
Collaborator Author

barrettk commented Jan 26, 2024

It seems to me that that line of thinking is in direct contradiction to claiming you must "guarantee priors wouldnt be altered".

Perhaps I didnt provide enough detail here. We need to know if there are priors with certainty (so worst case we could error out and not support it). The line you referenced would be the ideal state, and would handle almost all cases of being able to filter out priors. It's fine if we cant guarantee the identification of specific prior records.

I understand that (I think, at least). My suggestion is that you reconsider that design decision due to the desire to reliably know if there is a shape mismatch due to prior use. If we decide not to go that direction, that's of course fine, but here it feels like you're just dismissing the suggestion by reasserting what's already a given.

It's more that this is not the direction we want to go first. That is something we would consider if parsing priors seems too finicky/unreliable, but we want to attempt that direction first. We could only go that direction if we required a previously executed model run to pull out previous estimates, and there is at least a notable preference to not do that based on my discussions with Seth, Curtis and Tim (easier to have these discussions over zoom). That's something I can discuss with @seth127 next week when he gets back, but for now I have directions to look into the parsing of prior records.

@seth127
Copy link
Collaborator

seth127 commented Jan 30, 2024

Proposed path forward on priors

I talked with @barrettk for bit yesterday and I have a proposed path forward that I'd love some feedback on. I'll say upfront that my main motivation here is not wanting to try to identify and filter out the non-informative priors. Given discussion above, and some further research, this seems like a huge pain. I also don't like the idea of making this feature rely on a previously run model. (I'm happy to go into detail on either of those points, if anyone is interested.)

So then, my proposal:

  • We are pretty confident from the NONMEM docs that non-informative priors must be fixed. This is great news because it means that our function would not tweak them, regardless of whether it "knows" they are priors and not regular parameters. (@barrettk can maybe find more, but one example is in the nwpri.utl page in NM75 manual, which says "these records must be fixed" in three different places.)
  • Given this ^, I have less concern about not being able know which are non-informative priors, because our primary constraint was that we didn't want to tweak them erroneously. Now we know that we won't.
  • I propose that we have a fairly simple heuristic (details below) to check if we think they are using non-informative priors, and then give a warning that these are not supported (notably, initial_estimates() will return them as if they're regular parameters) and suggest they switch to informative priors.
  • Then we move on with the tweak because, as noted above, we won't actually tweak the priors (because they'll be fixed).

Thoughts on that? Anything that I'm missing?

Details

Here is some pseudo-code I threw together with the heuristic I have in mind:

bbr:::using_old_priors
#' This function uses a simple heuristic to decide
#' whether the control stream appears to be using the old
#' "non-informative" method of specifying priors
bbr:::using_old_priors(.ctl) {
  # check if PRIORS are being used
  prior_recs <- nmrec::select_records(.ctl, "prior")
  prior_subs <- nmrec::select_records(.ctl, "sub")
  prior_subs <- any(stringr::str_detect(prior_subs, "PRIOR")) # probably want nmrec to parse to see if PRIOR shows up in here...

  # check if INFORMATIVE style is present
  informative_priors <- purrr::walk(
    c("thetap", "thetapv", "omegap", "omegapd", "sigmap", "sigmapd"), # ARE THERE OTHERS???
    ~ {nmrec::select_records(.ctl, .x)}
  )    
  
  # if priors are specified AND _not_ using informative style, then assume using old style
  if (
      ((length(prior_recs) || length(prior_subs)) > 0) && 
      length(informative_priors) == 0   
    ) {
    return(TRUE)
  }
  return(FALSE)
}

Then we call that in initial_estimates() (maybe here) with a message something like this:

if (isTRUE(using_old_priors(ctl))) {
  rlang::warn(paste(
    "This model appears to be using the old 'non-informative' method of specifying priors.",
    "That will cause this function to extract and display priors as if they are initial parameter estimates.",
    "Consider changing the model to use informative prior record names (such as THETAP and THETAPV)."
  ))
}

and then we go on with our day...

@barrettk
Copy link
Collaborator Author

barrettk commented Jan 30, 2024

Still working on some matrix stuff before I commit the latest changes, but here are some examples of getting the initial estimates, based on some of the above discussion (including @seth127's most recent comment):

Example with no priors (`MOD1`)
> initial_estimates(.mod)
# A tibble: 8 × 6
  parameter_names record_type  init lower_bound upper_bound record_number
  <chr>           <chr>       <dbl>       <dbl>       <dbl>         <int>
1 THETA(1)        theta        2              0          NA             1
2 THETA(2)        theta        3              0          NA             1
3 THETA(3)        theta       10              0          NA             1
4 THETA(4)        theta        0.02          NA          NA             1
5 THETA(5)        theta        1             NA          NA             1
6 OMEGA(1,1)      omega        0.05          NA          NA             1
7 OMEGA(2,2)      omega        0.2           NA          NA             1
8 SIGMA(1,1)      sigma        1             NA          NA             1
Example when using informative priors
> initial_estimates(.mod)
# A tibble: 15 × 6
   parameter_names record_type    init lower_bound upper_bound record_number
   <chr>           <chr>         <dbl>       <dbl>       <dbl>         <int>
 1 THETA(1)        theta        -0.574          NA          NA             1
 2 THETA(2)        theta         5.67           NA          NA             1
 3 THETA(3)        theta         1.04           NA          NA             1
 4 THETA(4)        theta         0.242          NA          NA             1
 5 THETA(5)        theta       -11.4            NA          NA             1
 6 THETA(6)        theta         5.11           NA          NA             1
 7 OMEGA(1,1)      omega         0.09           NA          NA             1
 8 OMEGA(2,2)      omega         0.09           NA          NA             2
 9 OMEGA(3,2)      omega         0.01           NA          NA             2
10 OMEGA(4,2)      omega         0.01           NA          NA             2
11 OMEGA(3,3)      omega         0.09           NA          NA             2
12 OMEGA(4,3)      omega         0.01           NA          NA             2
13 OMEGA(4,4)      omega         0.09           NA          NA             2
14 SIGMA(1,1)      sigma         0.04           NA          NA             1
15 SIGMA(2,2)      sigma         0.04           NA          NA             1
Example when ***not*** using informative priors
> initial_estimates(.mod)
# A tibble: 33 × 6
   parameter_names record_type  init lower_bound upper_bound record_number
   <chr>           <chr>       <dbl>       <dbl>       <dbl>         <int>
 1 THETA(1)        theta         0.7          NA          NA             1
 2 THETA(2)        theta         0.7          NA          NA             1
 3 THETA(3)        theta         2            NA          NA             1
 4 THETA(4)        theta         2            NA          NA             1
 5 THETA(5)        theta         0.7          NA          NA             1
 6 THETA(6)        theta         0.7          NA          NA             1
 7 THETA(7)        theta         2            NA          NA             1
 8 THETA(8)        theta         2            NA          NA             1
 9 THETA(9)        theta         0.7          NA          NA             1
10 THETA(10)       theta         0.7          NA          NA             1
# ℹ 23 more rows
# ℹ Use `print(n = ...)` to see more rows
Warning message:
! This model appears to be using the old 'non-informative' method of specifying priors.- That will cause this function to extract and display priors *as if* they are initial parameter estimates.- Consider changing the model to use informative prior record names (such as THETAP and THETAPV). 

Function for checking for non-informative priors:

cc @kyleam to check out this modified function as well as Seth's comment above

Function
#' Check if PRIORS are being used
using_old_priors <- function(ctl) {
  # pull prior related records
  prior_recs <- nmrec::select_records(ctl, "prior")
  subroutine_recs <- nmrec::select_records(ctl, "sub")
  if(length(subroutine_recs) > 0){
    prior_subs <- purrr::map_lgl(subroutine_recs, function(prior_sub){
      prior_sub_str <- prior_sub$get_lines()
      # filter out comments
      prior_sub_str <- gsub(";.*", "", prior_sub_str)
      # look for priors
      any(stringr::str_detect(prior_sub_str, "(?i)PRIOR"))
    })
    prior_subs <- any(prior_subs)
  }else{
    prior_subs <- FALSE
  }

  # check if INFORMATIVE style is present
  prior_names <- c("thetap", "thetapv", "omegap", "omegapd", "sigmap", "sigmapd")
  informative_priors <- purrr::map(prior_names, function(.x){
    recs <- nmrec::select_records(ctl, .x)
    if(length(recs) == 0) return(NULL) else return(recs)
  }) %>% purrr::list_c()

  # if priors are specified AND _not_ using informative style, then assume using old style
  if ((length(prior_recs) > 0 || isTRUE(prior_subs)) && length(informative_priors) == 0){
    return(TRUE)
  }
  return(FALSE)
}

@kyleam
Copy link
Collaborator

kyleam commented Jan 30, 2024

@seth127

So then, my proposal:
[...]
Thoughts on that? Anything that I'm missing?

Sounds like a reasonable path forward to me.

 - `initial_estimates()` now tabuluates the record number
 - We check if non-informative priors are being used
 - tweaking each parameter type is now pulled out into separate helper functions. The `tweak_matrix` function will be altered in a later PR to address positive-definiteness
@timwaterhouse
Copy link
Collaborator

@seth127 Your proposal looks reasonable to me as well. My only suggestion would be to change the language for "informative" and "non-informative" record names, since that has specific meaning for Bayesian priors. I know the NONMEM documentation calls them that, but I think something like this might be more clear:

if (isTRUE(using_old_priors(ctl))) {
  rlang::warn(paste(
    "This model appears to be using $THETA, $OMEGA, and/or $SIGMA to specify priors.",
    "That will cause this function to extract and display priors as if they are initial parameter estimates.",
    "Consider changing the model to use the more specific prior record names (such as $THETAP and $THETAPV)."
  ))
}

@barrettk
Copy link
Collaborator Author

barrettk commented Jan 30, 2024

@timwaterhouse The only issue with changing the naming convention, is that nmrec also uses those names. I think we'd want those to be consistent. So if this is important I think we'd want to make the change in both places.

@kyleam
Copy link
Collaborator

kyleam commented Jan 30, 2024

I agree with the suggestion to avoid "non-informative" and will update nmrec.

kyleam added a commit to metrumresearchgroup/nmrec that referenced this pull request Jan 30, 2024
The set-param and extract-param docs call records like THETAP and
THETAPV as "informative prior record names" because NONMEM refers to
records like THETAP and THETAPV as "informative record names".
However, that invites confusion with the unrelated statistical concept
of informative/non-informative priors (even more so because the nmrec
docs unnecessarily insert "prior" between "informative" and "record
names").

Drop "informative" entirely to avoid confusion.  While at it, reduce
"prior record names" to "record names" because "prior" makes it harder
to parse and is unnecessary given the context.

Thanks to @timwaterhouse for the suggestion.

Re: metrumresearchgroup/bbr#635 (comment)
@barrettk barrettk closed this Jan 31, 2024
@barrettk barrettk deleted the feature/tweak-inits branch February 28, 2024 18:17
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
4 participants