nlmixr

The broom and broom.mixed packages

broom and broom.mixed are packages that attempt to put standard model outputs into data frames. nlmixr supports the tidy and glance methods but does not support augment at this time.

Using a model with a covariance term, the Phenobarbital model, we can explore the different types of output that is used in the tidy functions.

To explore this, first we run the model:

library(nlmixr)
library(broom.mixed)

pheno <- function() {
  # Pheno with covariance
  ini({
    tcl <- log(0.008) # typical value of clearance
    tv <-  log(0.6)   # typical value of volume
    ## var(eta.cl)
    eta.cl + eta.v ~ c(1, 
                       0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
    # interindividual variability on clearance and volume
    add.err <- 0.1    # residual variability
  })
  model({
    cl <- exp(tcl + eta.cl) # individual value of clearance
    v <- exp(tv + eta.v)    # individual value of volume
    ke <- cl / v            # elimination rate constant
    d/dt(A1) = - ke * A1    # model differential equation
    cp = A1 / v             # concentration in plasma
    cp ~ add(add.err)       # define error model
  })
}

## We will run it two ways to allow comparisons
fit.s <- nlmixr(pheno, pheno_sd, "saem", control=list(logLik=TRUE, print=0),
                table=list(cwres=TRUE, npde=TRUE))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> Error in .fitFun(.ret) : 
#>   function 'rx_ae17520d67619efc8c6aa0c314d3ee5a__calc_lhs' not provided by package 'rx_ae17520d67619efc8c6aa0c314d3ee5a_'
#> Error in (function (data, inits, PKpars, model = NULL, pred = NULL, err = NULL,  : 
#>   Could not fit data.
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> 
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
fit.f <- nlmixr(pheno, pheno_sd, "focei",
                control=list(print=0), 
                table=list(cwres=TRUE, npde=TRUE))
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00

Glancing at the goodness of fit metrics

Often in fitting data, you would want to glance at the fit to see how well it fits. In broom, glance will give a summary of the fit metrics of goodness of fit:

glance(fit.s)
#> # A tibble: 1 x 5
#>    OBJF   AIC   BIC logLik conditionNumber
#>   <dbl> <dbl> <dbl>  <dbl>           <dbl>
#> 1  707. 1004. 1022.  -496.            6.32

Note in nlmixr it is possible to have more than one fit metric (based on different quadratures, FOCEi approximation etc). However, the glance only returns the fit metrics that are current.

If you wish you can set the objective function to the focei objective function (which was already calculated with CWRES).

setOfv(fit.s,"gauss3_1.6")

Now the glance gives the gauss3_1.6 values.

glance(fit.s)
#> # A tibble: 1 x 5
#>    OBJF   AIC   BIC logLik conditionNumber
#>   <dbl> <dbl> <dbl>  <dbl>           <dbl>
#> 1  707. 1004. 1022.  -496.            6.32

Of course you can always change the type of objective function that nlmixr uses:

setOfv(fit.s,"FOCEi") # Setting objective function to focei

By setting it back to the SAEM default objective function of FOCEi, the glance(fit.s) has the same values again:

glance(fit.s)
#> # A tibble: 1 x 5
#>    OBJF   AIC   BIC logLik conditionNumber
#>   <dbl> <dbl> <dbl>  <dbl>           <dbl>
#> 1  689.  986. 1004.  -487.            6.32

For convenience, you can do this while you glance at the objects:

glance(fit.s, type="FOCEi")
#> # A tibble: 1 x 5
#>    OBJF   AIC   BIC logLik conditionNumber
#>   <dbl> <dbl> <dbl>  <dbl>           <dbl>
#> 1  689.  986. 1004.  -487.            6.32

Tidying the model parameters

Tidying of overall fit parameters

You can also tidy the model estimates into a data frame with broom for processing. This can be useful when integrating into 3rd parting modeling packages. With a consistent parameter format, tasks for multiple types of models can be automated and applied.

The default function for this is tidy, which when applied to the fit object provides the overall parameter information in a tidy dataset:

tidy(fit.s)
#> # A tibble: 6 x 7
#>   effect   group         term              estimate std.error statistic  p.value
#>   <chr>    <chr>         <chr>                <dbl>     <dbl>     <dbl>    <dbl>
#> 1 fixed    <NA>          tcl                 -4.99     0.0729    -68.5   1   e+0
#> 2 fixed    <NA>          tv                   0.341    0.0550      6.21  2.52e-9
#> 3 ran_pars ID            sd__eta.cl           0.476   NA          NA    NA      
#> 4 ran_pars ID            sd__eta.v            0.403   NA          NA    NA      
#> 5 ran_pars ID            cor__eta.v.eta.cl    0.959   NA          NA    NA      
#> 6 ran_pars Residual(add) add.err              2.79    NA          NA    NA

Note by default these are the parameters that are actually estimated in nlmixr, not the back-transformed values in the table from the printout. Of course, with mu-referenced models, you may want to exponentiate some of the terms. The broom package allows you to apply exponentiation on all the parameters, that is:

## Transformation applied on every parameter
tidy(fit.s, exponentiate=TRUE) 
#> # A tibble: 6 x 7
#>   effect   group         term             estimate std.error statistic   p.value
#>   <chr>    <chr>         <chr>               <dbl>     <dbl>     <dbl>     <dbl>
#> 1 fixed    <NA>          tcl               0.00681  0.000496      13.7  1.60e-28
#> 2 fixed    <NA>          tv                1.41     0.0773        18.2  5.41e-40
#> 3 ran_pars ID            sd__eta.cl        0.476   NA             NA   NA       
#> 4 ran_pars ID            sd__eta.v         0.403   NA             NA   NA       
#> 5 ran_pars ID            cor__eta.v.eta.…  0.959   NA             NA   NA       
#> 6 ran_pars Residual(add) add.err           2.79    NA             NA   NA

Note:, in accordance with the rest of the broom package, when the parameters with the exponentiated, the standard errors are transformed to an approximate standard error by the formula: \(\textrm{se}(\exp(x)) \approx \exp(\textrm{model estimate}_x)\times \textrm{se}_x\). This can be confusing because the confidence intervals (described later) are using the actual standard error and back-transforming to the exponentiated scale. This is the reason why the default for nlmixr’s broom interface is exponentiate=FALSE, that is:

tidy(fit.s, exponentiate=FALSE) ## No transformation applied
#> # A tibble: 6 x 7
#>   effect   group         term              estimate std.error statistic  p.value
#>   <chr>    <chr>         <chr>                <dbl>     <dbl>     <dbl>    <dbl>
#> 1 fixed    <NA>          tcl                 -4.99     0.0729    -68.5   1   e+0
#> 2 fixed    <NA>          tv                   0.341    0.0550      6.21  2.52e-9
#> 3 ran_pars ID            sd__eta.cl           0.476   NA          NA    NA      
#> 4 ran_pars ID            sd__eta.v            0.403   NA          NA    NA      
#> 5 ran_pars ID            cor__eta.v.eta.cl    0.959   NA          NA    NA      
#> 6 ran_pars Residual(add) add.err              2.79    NA          NA    NA

If you want, you can also use the parsed back-transformation that is used in nlmixr tables (ie fit$parFixedDf). Please note that this uses the approximate back-transformation for standard errors on the log-scaled back-transformed values.

This is done by:

## Transformation applied to log-scaled population parameters
tidy(fit.s, exponentiate=NA) 
#> # A tibble: 6 x 7
#>   effect   group         term             estimate std.error statistic   p.value
#>   <chr>    <chr>         <chr>               <dbl>     <dbl>     <dbl>     <dbl>
#> 1 fixed    <NA>          tcl               0.00681  0.000496      13.7  1.60e-28
#> 2 fixed    <NA>          tv                1.41     0.0773        18.2  5.41e-40
#> 3 ran_pars ID            sd__eta.cl        0.476   NA             NA   NA       
#> 4 ran_pars ID            sd__eta.v         0.403   NA             NA   NA       
#> 5 ran_pars ID            cor__eta.v.eta.…  0.959   NA             NA   NA       
#> 6 ran_pars Residual(add) add.err           2.79    NA             NA   NA

Also note, at the time of this writing the default separator between variables is ., which doesn’t work well with this model giving cor__eta.v.eta.cl. You can easily change this by:

options(broom.mixed.sep2="..")
tidy(fit.s)
#> # A tibble: 6 x 7
#>   effect   group         term              estimate std.error statistic  p.value
#>   <chr>    <chr>         <chr>                <dbl>     <dbl>     <dbl>    <dbl>
#> 1 fixed    <NA>          tcl                 -4.99     0.0729    -68.5   1   e+0
#> 2 fixed    <NA>          tv                   0.341    0.0550      6.21  2.52e-9
#> 3 ran_pars ID            sd__eta.cl           0.476   NA          NA    NA      
#> 4 ran_pars ID            sd__eta.v            0.403   NA          NA    NA      
#> 5 ran_pars ID            cor__eta.v..eta.…    0.959   NA          NA    NA      
#> 6 ran_pars Residual(add) add.err              2.79    NA          NA    NA

This gives an easier way to parse value: cor__eta.v..eta.cl

Adding a confidence interval to the parameters

The default R method confint works with nlmixr fit objects:

confint(fit.s)
#>          model.est    estimate      2.5 %     97.5 %
#> tcl     -4.9889415 0.006812872 -5.1317729 -4.8461102
#> tv       0.3412029 1.406638665  0.2334961  0.4489098
#> add.err  2.7915933 2.791593287         NA         NA

This transforms the variables as described above. You can still use the exponentiate parameter to control the display of the confidence interval:

confint(fit.s, exponentiate=FALSE)
#>          model.est    estimate      2.5 %     97.5 %
#> tcl     -4.9889415 0.006812872 -5.1317729 -4.8461102
#> tv       0.3412029 1.406638665  0.2334961  0.4489098
#> add.err  2.7915933 2.791593287         NA         NA

However, broom has also implemented it own way to make these data a tidy dataset. The easiest way to get these values in a nlmixr dataset is to use:

tidy(fit.s, conf.level=0.9)
#> # A tibble: 6 x 9
#>   effect  group  term   estimate std.error statistic  p.value conf.low conf.high
#>   <chr>   <chr>  <chr>     <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
#> 1 fixed   <NA>   tcl      -4.99     0.0729    -68.5   1   e+0   -5.11     -4.87 
#> 2 fixed   <NA>   tv        0.341    0.0550      6.21  2.52e-9    0.251     0.432
#> 3 ran_pa… ID     sd__e…    0.476   NA          NA    NA         NA        NA    
#> 4 ran_pa… ID     sd__e…    0.403   NA          NA    NA         NA        NA    
#> 5 ran_pa… ID     cor__…    0.959   NA          NA    NA         NA        NA    
#> 6 ran_pa… Resid… add.e…    2.79    NA          NA    NA         NA        NA

The confidence interval is on the scale specified by exponentiate, by default the estimated scale.

If you want to have the confidence on the adaptive back-transformed scale, you would simply use the following:

tidy(fit.s, conf.level=0.9, exponentiate=NA)
#> # A tibble: 6 x 9
#>   effect  group term   estimate std.error statistic   p.value conf.low conf.high
#>   <chr>   <chr> <chr>     <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
#> 1 fixed   <NA>  tcl     0.00681  0.000496      13.7  1.60e-28  0.00604   0.00768
#> 2 fixed   <NA>  tv      1.41     0.0773        18.2  5.41e-40  1.29      1.54   
#> 3 ran_pa… ID    sd__e…  0.476   NA             NA   NA        NA        NA      
#> 4 ran_pa… ID    sd__e…  0.403   NA             NA   NA        NA        NA      
#> 5 ran_pa… ID    cor__…  0.959   NA             NA   NA        NA        NA      
#> 6 ran_pa… Resi… add.e…  2.79    NA             NA   NA        NA        NA

Extracting other model information with tidy

The type of information that is extracted can be controlled by the effects argument.

Extracting only fixed effect parameters

The fixed effect parameters can be extracted by effects="fixed"

tidy(fit.s, effects="fixed")
#> # A tibble: 2 x 6
#>   effect term  estimate std.error statistic       p.value
#>   <chr>  <chr>    <dbl>     <dbl>     <dbl>         <dbl>
#> 1 fixed  tcl     -4.99     0.0729    -68.5  1            
#> 2 fixed  tv       0.341    0.0550      6.21 0.00000000252

Extracting only random parameters

The random standard deviations can be extracted by effects="ran_pars":

tidy(fit.s, effects="ran_pars")
#> # A tibble: 4 x 4
#>   effect   group         term               estimate
#>   <chr>    <chr>         <chr>                 <dbl>
#> 1 ran_pars ID            sd__eta.cl            0.476
#> 2 ran_pars ID            sd__eta.v             0.403
#> 3 ran_pars ID            cor__eta.v..eta.cl    0.959
#> 4 ran_pars Residual(add) add.err               2.79

Extracting random values (also called ETAs)

The random values, or in NONMEM the ETAs, can be extracted by effects="ran_vals" or effects="random"

head(tidy(fit.s, effects="ran_vals"))
#> # A tibble: 6 x 5
#>   effect   group level term   estimate
#>   <chr>    <chr> <int> <fct>     <dbl>
#> 1 ran_vals ID        1 eta.cl  -0.0991
#> 2 ran_vals ID        2 eta.cl  -0.229 
#> 3 ran_vals ID        3 eta.cl   0.239 
#> 4 ran_vals ID        4 eta.cl  -0.510 
#> 5 ran_vals ID        5 eta.cl   0.318 
#> 6 ran_vals ID        6 eta.cl  -0.150

This duplicate method of running effects is because the broom package supports effects="random" while the broom.mixed package supports effects="ran_vals".

Extracting random coefficients

Random coefficients are the population fixed effect parameter + the random effect parameter, possibly transformed to the correct scale.

In this case we can extract this information from a nlmixr fit object by:

head(tidy(fit.s, effects="ran_coef"))
#> # A tibble: 6 x 5
#>   effect   group level term  estimate
#>   <chr>    <chr> <int> <fct>    <dbl>
#> 1 ran_coef ID        1 tcl      -5.09
#> 2 ran_coef ID        2 tcl      -5.22
#> 3 ran_coef ID        3 tcl      -4.75
#> 4 ran_coef ID        4 tcl      -5.50
#> 5 ran_coef ID        5 tcl      -4.67
#> 6 ran_coef ID        6 tcl      -5.14

This can also be changed by the exponentiate argument:

head(tidy(fit.s, effects="ran_coef", exponentiate=NA))
#> # A tibble: 6 x 5
#>   effect   group level term  estimate
#>   <chr>    <chr> <int> <fct>    <dbl>
#> 1 ran_coef ID        1 tcl    0.00617
#> 2 ran_coef ID        2 tcl    0.00542
#> 3 ran_coef ID        3 tcl    0.00865
#> 4 ran_coef ID        4 tcl    0.00409
#> 5 ran_coef ID        5 tcl    0.00936
#> 6 ran_coef ID        6 tcl    0.00586
head(tidy(fit.s, effects="ran_coef", exponentiate=TRUE))
#> # A tibble: 6 x 5
#>   effect   group level term  estimate
#>   <chr>    <chr> <int> <fct>    <dbl>
#> 1 ran_coef ID        1 tcl    0.00617
#> 2 ran_coef ID        2 tcl    0.00542
#> 3 ran_coef ID        3 tcl    0.00865
#> 4 ran_coef ID        4 tcl    0.00409
#> 5 ran_coef ID        5 tcl    0.00936
#> 6 ran_coef ID        6 tcl    0.00586

Example of using a tidy model estimates for other packages

Dotwhisker

As explained above, this standard format makes it easier for tidyverse packages to interact with model information. An example of this is piping the tidy information to dplyr to filter the effects and then to the dotwhisker package to plot the model parameter confidence intervals.

options(broom.mixed.sep2=": ", broom.mixed.sep2=", ")
library(ggplot2)
library(dotwhisker)
library(dplyr)
fit.s %>%
    tidy(exponentiate=NA) %>%
    filter(effect=="fixed") %>%
    dwplot()

You may also compare models easily by using the dotwhisker package

dwplot(list("SAEM"=fit.s, "FOCEi"=fit.f), exponentiate=NA)

Huxtable

This allows easy creation of report ready tables in many formats including word.

Huxtable relies on the broom implementation

library(huxtable)
tbl <- huxreg('Phenobarbitol'=fit.s)

tbl
Phenobarbitol
tcl -4.989    
(0.073)   
tv 0.341 ***
(0.055)   
sd__eta.cl 0.476    
(NA)        
sd__eta.v 0.403    
(NA)        
cor__eta.v, eta.cl 0.959    
(NA)        
add.err 2.792    
(NA)        
N 155        
logLik -486.921    
AIC 985.842    
*** p < 0.001; ** p < 0.01; * p < 0.05.

Like dwplot, you can also use huxtable to compare runs:

huxreg('SAEM'=fit.s, 'FOCEi'=fit.f)
SAEM FOCEi
tcl -4.989     -5.000    
(0.073)    (0.086)   
tv 0.341 *** 0.335 ***
(0.055)    (0.062)   
sd__eta.cl 0.476     0.497    
(NA)         (NA)        
sd__eta.v 0.403     0.396    
(NA)         (NA)        
cor__eta.v, eta.cl 0.959     0.980    
(NA)         (NA)        
add.err 2.792     2.826    
(NA)         (NA)        
N 155         155        
logLik -486.921     -486.744    
AIC 985.842     985.487    
*** p < 0.001; ** p < 0.01; * p < 0.05.

A word-based table can also be easily created with the tool:

library(officer)
library(flextable)

ft  <- huxtable::as_flextable(tbl);
    
read_docx() %>%
    flextable::body_add_flextable(ft)  %>%
    print(target="pheno.docx")

Which produces the following word document.

Happy tidying!