Simulate using Parameter Uncertainty

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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/RtmppduNoF/Rinst17e41f60bee3/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/RtmppduNoF/Rinst17e41f60bee3/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/RtmppduNoF/Rinst17e41f60bee3/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 addlKeepsCov=TRUE like NONMEM, specify directly to change
#> ℹ using addlDropSs=TRUE like NONMEM, specify directly to change
#> ℹ using ssAtDoseTime=TRUE like NONMEM, specify directly to change
#> ℹ using safeZero=FALSE since NONMEM does not use protection by default
#> ℹ using ss2cancelAllPending=FALSE since NONMEM does not cancel pending doses with SS=2
#> ℹ 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 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.32   4.20   1.33   3.89 0.201  0.132   0.136  -0.0170
#>  2      1 2       1.32   4.20   1.33   3.89 0.201 -0.213   0.0788 -0.297 
#>  3      1 3       1.32   4.20   1.33   3.89 0.201  0.0916  0.294   0.358 
#>  4      1 4       1.32   4.20   1.33   3.89 0.201 -0.270   0.333  -0.0127
#>  5      1 5       1.32   4.20   1.33   3.89 0.201  0.298  -0.0873 -0.234 
#>  6      1 6       1.32   4.20   1.33   3.89 0.201  0.196   0.826   0.228 
#>  7      1 7       1.32   4.20   1.33   3.89 0.201 -0.263   0.144   0.128 
#>  8      1 8       1.32   4.20   1.33   3.89 0.201  0.125  -0.610  -0.489 
#>  9      1 9       1.32   4.20   1.33   3.89 0.201  0.356  -0.483  -0.0748
#> 10      1 10      1.32   4.20   1.33   3.89 0.201  0.107  -0.275   0.484 
#> # ℹ 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.29  76.4  3.73  45.0  76.4   76.4 0.0830 0.0489 1417.
#> 2      1     1     2  4.29  76.4  3.73  45.0  76.4   76.4 0.0830 0.0489 1284.
#> 3      1     1     3  4.29  76.4  3.73  45.0  76.4   76.4 0.0830 0.0489 1168.
#> 4      1     1     4  4.29  76.4  3.73  45.0  76.4   76.4 0.0830 0.0489 1067.
#> 5      1     1     5  4.29  76.4  3.73  45.0  76.4   76.4 0.0830 0.0489  979.
#> 6      1     1     6  4.29  76.4  3.73  45.0  76.4   76.4 0.0830 0.0489  901.
#> # ℹ 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  87432.  94265.  98648. 2.5%      
#>  2 0.5        1 CENTRAL 105082. 106542. 108275. 50%       
#>  3 0.975      1 CENTRAL 111775. 113137. 114586. 97.5%     
#>  4 0.0250     2 CENTRAL  66092.  75242.  81946. 2.5%      
#>  5 0.5        2 CENTRAL  92729.  95053.  98169. 50%       
#>  6 0.975      2 CENTRAL 104349. 106948. 109566. 97.5%     
#>  7 0.0250     3 CENTRAL  51785.  60796.  68936. 2.5%      
#>  8 0.5        3 CENTRAL  82193.  85361.  89391. 50%       
#>  9 0.975      3 CENTRAL  97672. 101347. 104901. 97.5%     
#> 10 0.0250     4 CENTRAL  41524.  49825.  59256. 2.5%      
#> # ℹ 80 more rows

plot(sci)


plot(sci, log="y")

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.