Get initial parameter estimates and tweak initial parameter estimates

[barrettk]

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

```r 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 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: ```r 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 ```r 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 https://github.com/metrumresearchgroup/bbr/issues/632

[barrettk]

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

```r > .mod <- MOD1 > get_param_inits(.mod) $theta # A tibble: 5 × 3 low init index 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

```r > .mod <- read_model(file.path(MODEL_DIR_C, "1001")) > get_param_inits(.mod) $theta # A tibble: 6 × 2 init index 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

```r > .mod <- read_model(file.path(MODEL_DIR_C, "example2_saemimp")) > get_param_inits(.mod) $theta # A tibble: 12 × 3 init fixed index 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 ```

[barrettk]

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 ```r > get_initial_est(.mod) $thetas # A tibble: 13 × 3 init low up 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 ```r > get_initial_est(.mod, flag_fixed = TRUE) $thetas # A tibble: 13 × 4 init low up fixed 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

```r > initial_estimates(MOD1, flag_fixed = TRUE) # A tibble: 9 × 6 parameter_names record_type init lower_bound upper_bound fixed 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: ```r 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 ```r 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] ```

[barrettk]

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

[barrettk]

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)

[barrettk]

@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

```r > MODEL_DIR_X [1] "inst/model/nonmem/complex" > .mod <- read_model(file.path(MODEL_DIR_X, "example2_saemimp")) ``` ```r 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]

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]

@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: ```r .mod <- read_model(file.path(MODEL_DIR_X, "example2_saemimp")) tweak_initial_estimates(.mod) ``` results in this 

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

[seth127]

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]

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]

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]

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]

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]

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]

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]

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]

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]

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

```r #' 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]

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

```r > initial_estimates(.mod) # A tibble: 8 × 6 parameter_names record_type init lower_bound upper_bound record_number 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

```r > initial_estimates(.mod) # A tibble: 15 × 6 parameter_names record_type init lower_bound upper_bound record_number 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

```r > initial_estimates(.mod) # A tibble: 33 × 6 parameter_names record_type init lower_bound upper_bound record_number 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

```r #' 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]

@seth127

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

Sounds like a reasonable path forward to me.

[timwaterhouse]

@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]

@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]

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