Created Augmented pred/ipred plots with `augPred()`

This is a simple process to create individual predictions augmented with more observations than was modeled. This allows smoother plots and a better examination of the observed concentrations for an individual and population.

Step 1: Convert the Monolix model to rxode2:

library(monolix2rx)

# First we need the location of the monolix mlxtran file. Since we are
# running an example, we will use one of the built-in examples in
# `monolix2rx`
pkgTheo <- system.file("theo/theophylline_project.mlxtran", package="monolix2rx")
# You can use a control stream or other file. With the development
# version of `babelmixr2`, you can simply point to the listing file

mod <- monolix2rx(pkgTheo)
#> ℹ integrated model file 'oral1_1cpt_kaVCl.txt' into mlxtran object
#> ℹ updating model values to final parameter estimates
#> ℹ done
#> ℹ reading run info (# obs, doses, Monolix Version, etc) from summary.txt
#> ℹ done
#> ℹ reading covariance from FisherInformation/covarianceEstimatesLin.txt
#> ℹ done
#> Warning in .dataRenameFromMlxtran(data, .mlxtran): NAs introduced by coercion
#> ℹ imported monolix and translated to rxode2 compatible data ($monolixData)
#> ℹ imported monolix ETAS (_SAEM) imported to rxode2 compatible data ($etaData)
#> ℹ imported monolix pred/ipred data to compare ($predIpredData)
#> ℹ solving ipred problem
#> ℹ done
#> ℹ solving pred problem
#> ℹ done

Step 2: convert the rxode2 model to nlmixr2

You can convert the model, mod, to a nlmixr2 fit object:

library(babelmixr2) # provides as.nlmixr2
#> Loading required package: nlmixr2
#> Loading required package: nlmixr2data
fit <- as.nlmixr2(mod)
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> rxode2 3.0.2.9000 using 2 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 7168
#> ℹ monolix parameter history integrated into fit object

fit
#> ── nlmixr² monolix2rx reading Monolix ver 5.1.1 ──
#> 
#>              OBJF     AIC      BIC Log-likelihood Condition#(Cov)
#> monolix  118.9368 355.482 377.7819       -169.741        21.26161
#>          Condition#(Cor)
#> monolix         1.383153
#> 
#> ── Time (sec fit$time): ──
#> 
#>            setup optimize covariance table compress as.nlmixr2
#> elapsed 0.031842    4e-06      5e-06 0.056    0.007      2.091
#> 
#> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
#> 
#>          Est.     SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> ka_pop  0.427  0.204 47.8       1.53 (1.03, 2.29)     75.4      1.05% 
#> V_pop  -0.786  0.045 5.72    0.456 (0.417, 0.497)     12.7      13.3% 
#> Cl_pop  -3.21 0.0837 2.61 0.0402 (0.0341, 0.0473)     27.6      2.65% 
#> a       0.433                               0.433                     
#> b      0.0543                              0.0543                     
#>  
#>   Covariance Type (fit$covMethod): monolix2rx
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink 
#>   Censoring (fit$censInformation): No censoring
#>   Minimization message (fit$message):  
#>     IPRED relative difference compared to Monolix IPRED: 0.04%; 95% percentile: (0%,0.52%); rtol=0.000379
#>     PRED relative difference compared to Monolix PRED: 0%; 95% percentile: (0%,0%); rtol=4.94e-07
#>     IPRED absolute difference compared to Monolix IPRED: atol=0.00253; 95% percentile: (0.000364, 0.00848)
#>     PRED absolute difference compared to Monolix PRED: atol=4.94e-07; 95% percentile: (1.13e-08, 0.000308) 
#> 
#> ── Fit Data (object fit is a modified tibble): ──
#> # A tibble: 120 × 20
#>   ID     TIME    DV  PRED    RES IPRED     IRES    IWRES omega_ka omega_V
#>   <fct> <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>    <dbl>   <dbl>
#> 1 1      0.25  2.84  2.78 0.0636  3.73 -0.887   -1.40       0.132  -0.183
#> 2 1      0.57  6.57  5.00 1.57    6.57  0.00239  0.00303    0.132  -0.183
#> 3 1      1.12 10.5   6.80 3.70    8.75  1.75     1.93       0.132  -0.183
#> # ℹ 117 more rows
#> # ℹ 10 more variables: omega_Cl <dbl>, CONC <dbl>, depot <dbl>, central <dbl>,
#> #   ka <dbl>, V <dbl>, Cl <dbl>, Cc <dbl>, tad <dbl>, dosenum <dbl>

Step 3: Create and plot an augmented prediction

ap <- augPred(fit)
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

head(ap)
#>     values        ind id  time Endpoint
#> 1 0.000000 Individual  1 0.000    depot
#> 2 3.726962 Individual  1 0.250    depot
#> 3 6.030526 Individual  1 0.493    depot
#> 4 6.567601 Individual  1 0.570    depot
#> 5 8.414706 Individual  1 0.986    depot
#> 6 8.746028 Individual  1 1.120    depot

# This augpred looks odd:
plot(ap)