Simulate using Parameter Uncertainty

This page shows a simple work-flow for directly simulating a different dosing paradigm than what was modeled taking into account the modeled uncertainty. This workflow is very similar to simply simulating without uncertainty in the parameters themselves.

Step 1: Import the model

library(nonmem2rx)
library(rxode2)
# its best practice to set the seed for the simulations
set.seed(42)
rxSetSeed(42)

# First we need the location of the nonmem control stream Since we are
# running an example, we will use one of the built-in examples in
# `nonmem2rx`
ctlFile <- system.file("mods/cpt/runODE032.ctl", package="nonmem2rx")
# You can use a control stream or other file. With the development
# version of `babelmixr2`, you can simply point to the listing file
mod <- nonmem2rx(ctlFile, lst=".res", save=FALSE, determineError=FALSE)
#> ℹ getting information from  '/tmp/RtmpTDyhzd/Rinst45f6fcdf3b4/nonmem2rx/mods/cpt/runODE032.ctl'
#> ℹ reading in xml file
#> ℹ done
#> ℹ reading in phi file
#> ℹ done
#> ℹ reading in lst file
#> ℹ abbreviated list parsing
#> ℹ done
#> ℹ done
#> ℹ splitting control stream by records
#> ℹ done
#> ℹ Processing record $INPUT
#> ℹ Processing record $MODEL
#> ℹ Processing record $THETA
#> ℹ Processing record $OMEGA
#> ℹ Processing record $SIGMA
#> ℹ Processing record $PROBLEM
#> ℹ Processing record $DATA
#> ℹ Processing record $SUBROUTINES
#> ℹ Processing record $PK
#> ℹ Processing record $DES
#> ℹ Processing record $ERROR
#> ℹ Processing record $ESTIMATION
#> ℹ Ignore record $ESTIMATION
#> ℹ Processing record $COVARIANCE
#> ℹ Ignore record $COVARIANCE
#> ℹ Processing record $TABLE
#> ℹ change initial estimate of `theta1` to `1.37034036528946`
#> ℹ change initial estimate of `theta2` to `4.19814911033061`
#> ℹ change initial estimate of `theta3` to `1.38003493562413`
#> ℹ change initial estimate of `theta4` to `3.87657341967489`
#> ℹ change initial estimate of `theta5` to `0.196446108190896`
#> ℹ change initial estimate of `eta1` to `0.101251418415006`
#> ℹ change initial estimate of `eta2` to `0.0993872449483344`
#> ℹ change initial estimate of `eta3` to `0.101302674763154`
#> ℹ change initial estimate of `eta4` to `0.0730497519364148`
#> ℹ read in nonmem input data (for model validation): /tmp/RtmpTDyhzd/Rinst45f6fcdf3b4/nonmem2rx/mods/cpt/Bolus_2CPT.csv
#> ℹ ignoring lines that begin with a letter (IGNORE=@)'
#> ℹ applying names specified by $INPUT
#> ℹ subsetting accept/ignore filters code: .data[-which((.data$SD == 0)),]
#> ℹ done
#> ℹ read in nonmem IPRED data (for model validation): /tmp/RtmpTDyhzd/Rinst45f6fcdf3b4/nonmem2rx/mods/cpt/runODE032.csv
#> ℹ done
#> ℹ changing most variables to lower case
#> ℹ done
#> ℹ replace theta names
#> ℹ done
#> ℹ replace eta names
#> ℹ done (no labels)
#> ℹ renaming compartments
#> ℹ done
#> ℹ solving ipred problem
#> ℹ done
#> ℹ solving pred problem
#> ℹ done

Step 2: Look at a different dosing paradigm

Lets say that in this case instead of a single dose, we want to see what the concentration profile is with a single day of BID dosing. In this case is done by creating a quick event table.

ev <- et(amt=120000, ii=12, until=24) %>%
  et(c(1:6, seq(8, 24, by=2))) %>%
  et(id=1:100)

Step 3: Solve using the uncertainty in the NONMEM model

To use the uncertainty in the model, it is a simple matter of telling how many times rxode2() should sample with nStud=X. In this case we will use 100.

s <- rxSolve(mod, ev, nStud=100)
#> ℹ using nocb interpolation like NONMEM, specify directly to change
#> ℹ using safeZero=FALSE since NONMEM does not use protection by default
#> ℹ using dfSub=120 from NONMEM
#> ℹ using dfObs=2280 from NONMEM
#> ℹ using thetaMat from NONMEM
#> ℹ using sigma from NONMEM
#> ℹ using NONMEM specified atol=1e-12
#> ℹ using NONMEM specified rtol=1e-06
#> ℹ using NONMEM specified ssRtol=1e-06
#> ℹ using NONMEM specified ssAtol=1e-12
#> ℹ thetaMat has too many items, ignored: 'omega.2.1', 'omega.3.1', 'omega.3.2', 'omega.4.1', 'omega.4.2', 'omega.4.3'
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> Warning: corrected 'thetaMat' to be a symmetric, positive definite matrix

s
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#> # A tibble: 10,000 × 11
#>    sim.id id    theta1 theta2 theta3 theta4   RSV    eta1    eta2    eta3
#>     <int> <fct>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>   <dbl>   <dbl>   <dbl>
#>  1      1 1       1.38   4.22   1.44   3.89 0.193  0.0233 -0.494  -0.505 
#>  2      1 2       1.38   4.22   1.44   3.89 0.193  0.0129  0.100  -0.240 
#>  3      1 3       1.38   4.22   1.44   3.89 0.193  0.180   0.493   0.0880
#>  4      1 4       1.38   4.22   1.44   3.89 0.193  0.748  -0.401   0.130 
#>  5      1 5       1.38   4.22   1.44   3.89 0.193 -0.164  -0.812  -0.237 
#>  6      1 6       1.38   4.22   1.44   3.89 0.193  0.0685  0.186   0.0503
#>  7      1 7       1.38   4.22   1.44   3.89 0.193  0.250   0.0895  0.490 
#>  8      1 8       1.38   4.22   1.44   3.89 0.193 -0.571   0.204   0.443 
#>  9      1 9       1.38   4.22   1.44   3.89 0.193  0.426  -0.394  -0.431 
#> 10      1 10      1.38   4.22   1.44   3.89 0.193  0.388  -0.236  -0.238 
#> # ℹ 9,990 more rows
#> # ℹ 1 more variable: eta4 <dbl>
#> ── Initial Conditions (x$inits): ──
#> CENTRAL    PERI 
#>       0       0 
#> 
#> Simulation with uncertainty in:
#> • parameters (x$thetaMat for changes)
#> • omega matrix (x$omegaList)
#> • sigma matrix (x$sigmaList)
#> 
#> ── First part of data (object): ──
#> # A tibble: 150,000 × 21
#>   sim.id    id  time    cl     v     q    v2    v1 scale1    k21    k12     f
#>    <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>
#> 1      1     1     1  4.05  41.5  2.54  60.5  41.5   41.5 0.0419 0.0611 2469.
#> 2      1     1     2  4.05  41.5  2.54  60.5  41.5   41.5 0.0419 0.0611 2115.
#> 3      1     1     3  4.05  41.5  2.54  60.5  41.5   41.5 0.0419 0.0611 1818.
#> 4      1     1     4  4.05  41.5  2.54  60.5  41.5   41.5 0.0419 0.0611 1569.
#> 5      1     1     5  4.05  41.5  2.54  60.5  41.5   41.5 0.0419 0.0611 1358.
#> 6      1     1     6  4.05  41.5  2.54  60.5  41.5   41.5 0.0419 0.0611 1181.
#> # ℹ 149,994 more rows
#> # ℹ 9 more variables: ipred <dbl>, rescv <dbl>, w <dbl>, ires <dbl>,
#> #   iwres <dbl>, y <dbl>, CENTRAL <dbl>, PERI <dbl>, DV <dbl>

Step 4: Summarize and plot

Since there is a bunch of data, a confidence band of the simulation with uncertainty would be helpful.

One way to do that is to select the interesting components, create a confidence interval and then plot the confidence bands:

sci <- confint(s, parm=c("CENTRAL", "PERI", "sim"))
#> summarizing data...done

sci
#> # A tibble: 90 × 7
#>        p1  time trt        p2.5     p50   p97.5 Percentile
#>     <dbl> <dbl> <fct>     <dbl>   <dbl>   <dbl> <fct>     
#>  1 0.0250     1 CENTRAL  87958.  93878.  98279. 2.5%      
#>  2 0.5        1 CENTRAL 104877. 106584. 108047. 50%       
#>  3 0.975      1 CENTRAL 111682. 113335. 114665. 97.5%     
#>  4 0.0250     2 CENTRAL  66598.  74518.  81202. 2.5%      
#>  5 0.5        2 CENTRAL  92436.  95159.  97666. 50%       
#>  6 0.975      2 CENTRAL 104230. 107250. 109672. 97.5%     
#>  7 0.0250     3 CENTRAL  51664.  60036.  67728. 2.5%      
#>  8 0.5        3 CENTRAL  82106.  85486.  88980. 50%       
#>  9 0.975      3 CENTRAL  97531. 101717. 104996. 97.5%     
#> 10 0.0250     4 CENTRAL  40963.  49433.  57207. 2.5%      
#> # ℹ 80 more rows

plot(sci)


plot(sci, log="y")