This vignette demonstrates extracting parameter estimates and labels into a table that can be used for diagnostics, or to generate reports. Note, this vignette exclusively deals with NONMEM models.
If you are new to bbr, the “Getting
Started with bbr” vignette will take you through some basic
scenarios for modeling with NONMEM using bbr, introducing
you to its standard workflow and functionality.
We have a control stream and output directory, because the model has already been run.
MODEL_DIR <- system.file("test-refs", "param-labels", "example-model", package = "bbr")
mod1 <- read_model(file.path(MODEL_DIR, "510"))
class(mod1)
#> [1] "bbi_nonmem_model" "bbi_base_model" "bbi_model" "list"The model object you have just created can now be passed to the post-processing functions to create your tables.
Currently, the param_labels() function parses labels
from the comments in the control stream. Here is the relevant section of
our example control stream.
mod1 %>%
get_model_path() %>% # get control stream file path
read_lines(skip = 15, n_max = 18) %>% # read in only parameter section
paste(collapse = "\n") %>%
cat()
#> $THETA
#> (0, 13) ;[L/day] CL
#> (0, 75) ;[L] V
#> (0.001 );[L/day] CL_{CLCR}
#> $THETA
#> (0.001) ;[L/day] CL_{AGE}
#> (0.001) ;[L] V_{WT}
#> (0.001) ;[L] V_{AGE}
#>
#> $OMEGA BLOCK(2)
#> ; a comment
#> (0.04) ;[P] CL
#> (0.02) ;[R] CL-V
#> (0.04) ;[P] V
#>
#> $SIGMA
#> 0.04 ;[P] ResidualThis control stream is parsed into the tibble below, following the
syntax defined in the “Details” section of the
?param_labels documentation.
label_df <- mod1 %>%
param_labels() %>%
apply_indices(.omega = block(2))
label_df
#> # A tibble: 10 × 5
#> parameter_names label unit type param_type
#> <chr> <chr> <chr> <chr> <chr>
#> 1 THETA1 CL "L/day" "" THETA
#> 2 THETA2 V "L" "" THETA
#> 3 THETA3 CL_{CLCR} "L/day" "" THETA
#> 4 THETA4 CL_{AGE} "L/day" "" THETA
#> 5 THETA5 V_{WT} "L" "" THETA
#> 6 THETA6 V_{AGE} "L" "" THETA
#> 7 OMEGA(1,1) CL "" "[P]" OMEGA
#> 8 OMEGA(2,1) CL-V "" "[R]" OMEGA
#> 9 OMEGA(2,2) V "" "[P]" OMEGA
#> 10 SIGMA(1,1) Residual "" "[P]" SIGMANote, there are some subtleties to the apply_indices()
function that will be addressed in the next section.
The user can also extract parameter estimates using the
model_summary() and param_estimates()
functions.
param_df <- mod1 %>%
model_summary() %>%
param_estimates()
param_df
#> # A tibble: 10 × 8
#> parameter_names estimate stderr random_effect_sd random_effect_sdse fixed
#> <chr> <dbl> <dbl> <dbl> <dbl> <lgl>
#> 1 THETA1 12.3 0.961 NA NA FALSE
#> 2 THETA2 90.7 4.94 NA NA FALSE
#> 3 THETA3 0.163 0.335 NA NA FALSE
#> 4 THETA4 0.567 1.57 NA NA FALSE
#> 5 THETA5 0.403 0.342 NA NA FALSE
#> 6 THETA6 -0.395 2.32 NA NA FALSE
#> 7 OMEGA(1,1) 0.0517 0.0122 0.227 0.0269 FALSE
#> 8 OMEGA(2,1) 0.00345 0.00946 0.0654 0.178 FALSE
#> 9 OMEGA(2,2) 0.0538 0.0120 0.232 0.0259 FALSE
#> 10 SIGMA(1,1) 0.0450 0.00388 0.212 0.00914 FALSE
#> # ℹ 2 more variables: diag <lgl>, shrinkage <dbl>These two tibbles can be joined together to create a table for including in reports.
report_df <- inner_join(
label_df %>% select(-param_type),
param_df %>% select(parameter_names, estimate, stderr),
by = "parameter_names"
)
report_df
#> # A tibble: 10 × 6
#> parameter_names label unit type estimate stderr
#> <chr> <chr> <chr> <chr> <dbl> <dbl>
#> 1 THETA1 CL "L/day" "" 12.3 0.961
#> 2 THETA2 V "L" "" 90.7 4.94
#> 3 THETA3 CL_{CLCR} "L/day" "" 0.163 0.335
#> 4 THETA4 CL_{AGE} "L/day" "" 0.567 1.57
#> 5 THETA5 V_{WT} "L" "" 0.403 0.342
#> 6 THETA6 V_{AGE} "L" "" -0.395 2.32
#> 7 OMEGA(1,1) CL "" "[P]" 0.0517 0.0122
#> 8 OMEGA(2,1) CL-V "" "[R]" 0.00345 0.00946
#> 9 OMEGA(2,2) V "" "[P]" 0.0538 0.0120
#> 10 SIGMA(1,1) Residual "" "[P]" 0.0450 0.00388Because there are numerous ways of specifying the diagonal and
off-diagonal elements of an $OMEGA or $SIGMA
block in a control stream, automatically parsing the structure of these
blocks can be brittle and error prone. For this reason, indices are
not automatically added to the output of the
param_labels() function and are instead added with the
apply_indices() function.
By default apply_indices() assumes that all
$OMEGA and $SIGMA elements are diagonal. If
this is the case, you do not need to pass anything to
.omega or .sigma arguments described below.
However, be careful that you do not accidentally overlook this because
your indices will be incorrectly returned as simply (1,1),
(2,2), (3,3), etc.
Block structure is specified with two arguments, .omega
and .sigma which each take a logical vector. Each element
in this vector corresponds to a parameter specified in the control
stream and denotes whether that element is a diagonal in the relevant
matrix.
For example, an $OMEGA BLOCK(3) would be represented
with the vector
.omega = c(TRUE, FALSE, TRUE, FALSE, FALSE, TRUE) because
the first, third, and sixth elements are the diagonals and the others
represent the off-diagonals. However, the user doesn’t need to type this
explicitly because bbr has a helper function
block() that generates these vectors.
cat(paste("block(3): ", paste(block(3), collapse = ", "), "\n"))
#> block(3): TRUE, FALSE, TRUE, FALSE, FALSE, TRUE
cat(paste("block(4): ", paste(block(4), collapse = ", "), "\n"))
#> block(4): TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, TRUENotice the use of .omega = block(2) in the example from
the previous section.
More complicated block structures can be represented by concatenating these logical vectors together. For example, the following omega block:
cat(REF_OMEGA)
#>
#> $OMEGA BLOCK (3)
#> .1 ;[P] 5 P1NPF
#> .01 .1 ;[P] 6 CTFX
#> .01 .01 .1 ;[P] 7 LSF
#> $OMEGA BLOCK(2)
#> .1 ;[P] 8 FAKE1
#> .01 .1 ;[P] 9 FAKE2
#> $OMEGA BLOCK (1) 0.04 ; [P] IOV_{KA}
#> $OMEGA BLOCK(1) SAME
REF_OMEGA %>%
param_labels() %>%
apply_indices(
.omega = c(block(3), block(2), block(1), block(1))
)
#> # A tibble: 11 × 4
#> parameter_names label type param_type
#> <chr> <chr> <chr> <chr>
#> 1 OMEGA(1,1) "5 P1NPF" [P] OMEGA
#> 2 OMEGA(2,1) "" [A] OMEGA
#> 3 OMEGA(2,2) "6 CTFX" [P] OMEGA
#> 4 OMEGA(3,1) "" [A] OMEGA
#> 5 OMEGA(3,2) "" [A] OMEGA
#> 6 OMEGA(3,3) "7 LSF" [P] OMEGA
#> 7 OMEGA(4,4) "8 FAKE1" [P] OMEGA
#> 8 OMEGA(5,4) "" [A] OMEGA
#> 9 OMEGA(5,5) "9 FAKE2" [P] OMEGA
#> 10 OMEGA(6,6) "IOV_{KA}" [P] OMEGA
#> 11 OMEGA(7,7) "IOV_{KA}" [P] OMEGAIt is equally valid to replace any of these with logical vectors as
well. For instance, instead of c(..., block(1), block(1))
it may be easier to write c(..., rep(TRUE, 2)), especially
if there are a number of BLOCK(1) lines in a row in the
control stream.