Introduction to Simulations in Pumas

Author

Haden Bunn

using Pumas, PharmaDatasets
using DataFramesMeta, CategoricalArrays
using CairoMakie, AlgebraOfGraphics
using SummaryTables

1 Introduction

Simulation is a fundamental part of most pharmacometrics analyses. Fortunately, Pumas provides a powerful interface for simulating a variety of endpoints without the need for additional software. This tutorial will serve as a basic introduction to simulation in Pumas and will focus on simulating additional observations or alternative dosing scenarios in an existing population. Advanced simulation (e.g., clinical trial simulation) will be discussed in later tutorials.

2 Getting Started

Simulations are performed using the simobs function which has multiple methods, all of which are detailed in the Pumas documentation. Checking ?simobs in the REPL provides the following function signature:

simobs(
    model::AbstractPumasModel,
    population::Union{Subject,Population},
    param,
    randeffs = sample_randeffs(model, param, population);
    obstimes = nothing,
    ensemblealg = EnsembleSerial(),
    diffeq_options = NamedTuple(),
    rng = Random.default_rng(),
    simulate_error = Val(true),
)

The first three positional arguments (model, population, param) are all required.

model expects an AbstractPumasModel, which (for now) refers to a model defined using the @model macro.
population accepts a Population which was discussed in Module 4), or a single Subject which will be discussed later in the current module.
param should be a single parameter set defined as a NamedTuple or a vector of such parameter sets.

The remaining arguments have default values and need not be defined explicitly; however, it is worth knowing how the defaults affect each simulation.

randeffs is used to specify the random effect (“eta”) values for each subject. If left to the default, these values are generated by sampling the associated prior distributions defined in the model.
ensemblealg is used to select the parallelization mode to be used for the simulation.
diffeq_options can be used to pass additional options to the differential equation solver if the model does not have an analytical solution.
rng can be used to specify the random number generator to be used for the simulation.
simulate_error can be used to disable (false) the inclusion of RUV in the value returned by the predictive distribution’s error model.

Reproducibility

Many users are likely familiar with the concept of a random number generator (RNG) and the role they play in computational exercises where values are randomly sampled from a distribution. Using an RNG will make it (nearly) impossible to reproduce the results of a simulation unless steps are taken at the start to ensure reproducibility. In short, your results will differ slightly from those in the tutorial if you are executing the code locally, and that is to be expected. We will discuss this topic in greater detail later; for now, just focus on understanding the simulation workflow.

3 Setup

The examples below were created using the final integrated PK/PD model for warfarin. If you have downloaded this tutorial and are working through it locally, make sure you execute the code in the setup block before continuing. The example code assumes that the warfarin dataset (adeval), model (mdl), and fitted pumas model (“fpm”, myfit) all exist in your current session.

Warfarin Model

#* read dataset from PharmaDatasets
adeval = dataset("pumas/warfarin_pumas")

# population
mypop = read_pumas(adeval; observations = [:conc, :pca], covariates = [:wtbl])

# warfarin model
mdl = @model begin

    @metadata begin
        desc = "Integrated Warfarin PK/PD model"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        # Clearance, L/hr
        tvcl ∈ RealDomain(lower = 0.0, init = 0.134)
        # Volume of distribution, central, L
        tvvc ∈ RealDomain(lower = 0.0, init = 8.11)
        # absorption rate constant, hr^-1
        tvka ∈ RealDomain(lower = 0.0, init = 1.32)
        # absorption lag, hr
        tvalag ∈ RealDomain(lower = 0.0, init = 0.1)

        # PD parameters
        # Baseline, %
        tve0 ∈ RealDomain(lower = 0.0, init = 95, upper = 100)
        # Imax, %
        tvimax ∈ RealDomain(lower = 0.0, init = 0.8, upper = 1)
        # IC50, mg/L
        tvic50 ∈ RealDomain(lower = 0.0, init = 1.0)
        # Turnover
        tvturn ∈ RealDomain(lower = 0.0, init = 14.0)
        # Inter-individual variability
        """
          - Ωcl
          - Ωvc
          - Ωka
        """
        Ωpk ∈ PDiagDomain([0.01, 0.01, 0.01])
        """
          - Ωe0
          - Ωic50
        """
        Ωpd ∈ PDiagDomain([0.01, 0.01])
        # Residual variability
        # proportional error, pk
        σprop_pk ∈ RealDomain(lower = 0.0, init = 0.2)
        # additive error, pk, mg/L
        σadd_pk ∈ RealDomain(lower = 0.0, init = 0.2)
        # additive error, pca, %
        σadd_pd ∈ RealDomain(lower = 0.0, init = 1)
    end

    @random begin
        # mean = 0, covariance = Ωpk
        ηpk ~ MvNormal(Ωpk)
        # mean = 0, covariance = Ωpd
        ηpd ~ MvNormal(Ωpd)
    end

    @covariates wtbl

    @pre begin
        # PK
        cl = tvcl * (wtbl / 70)^0.75 * exp(ηpk[1])
        vc = tvvc * (wtbl / 70) * exp(ηpk[2])
        ka = tvka * exp(ηpk[3])
        # PD
        e0 = tve0 * exp(ηpd[1])
        imax = tvimax
        ic50 = tvic50 * exp(ηpd[2])
        turn = tvturn
        #kout = log(2) / turn
        kout = 1 / turn
        kin = e0 * kout
    end

    @dosecontrol begin
        lags = (depot = tvalag,)
    end

    @init begin
        e = e0
    end

    @vars begin
        # inhibitory model
        imdl := 1 - (imax * (central / vc)) / (ic50 + (central / vc))
    end

    @dynamics begin
        depot' = -ka * depot
        central' = ka * depot - (cl / vc) * central
        e' = kin * imdl - kout * e
    end

    @derived begin
        cp := @. central / vc
        # warfarin concentration, mg/L
        conc ~ @. Normal(cp, sqrt((σprop_pk * cp)^2 + σadd_pk^2))
        # prothrombin complex activity, % of normal
        pca ~ @. Normal(e, σadd_pd)
    end
end

# fitted pumas model, fpm
myfit = fit(mdl, mypop, init_params(mdl), FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     5.536943e+03     7.052675e+03
 * time: 0.034646034240722656
     1     1.762032e+03     1.253034e+03
 * time: 3.983558177947998
     2     1.449413e+03     6.274412e+02
 * time: 4.9340980052948
     3     1.331405e+03     2.099880e+02
 * time: 5.822710037231445
     4     1.310843e+03     1.627399e+02
 * time: 6.697432994842529
     5     1.295073e+03     1.376531e+02
 * time: 7.574280023574829
     6     1.282335e+03     1.324314e+02
 * time: 8.445069074630737
     7     1.275524e+03     9.521788e+01
 * time: 9.304376125335693
     8     1.271552e+03     5.354883e+01
 * time: 10.171005010604858
     9     1.269315e+03     2.407046e+01
 * time: 11.024396181106567
    10     1.268615e+03     2.368964e+01
 * time: 11.88334608078003
    11     1.268257e+03     2.366048e+01
 * time: 12.730183124542236
    12     1.267352e+03     3.190950e+01
 * time: 13.55699610710144
    13     1.265336e+03     5.399187e+01
 * time: 14.37921404838562
    14     1.260582e+03     8.775888e+01
 * time: 15.138981103897095
    15     1.252929e+03     1.110735e+02
 * time: 15.883120059967041
    16     1.246526e+03     8.328859e+01
 * time: 16.585114002227783
    17     1.240724e+03     8.711808e+01
 * time: 17.249306201934814
    18     1.236881e+03     7.439216e+01
 * time: 17.92052912712097
    19     1.232291e+03     6.878129e+01
 * time: 18.57830309867859
    20     1.227181e+03     5.604034e+01
 * time: 19.224945068359375
    21     1.223093e+03     3.856728e+01
 * time: 19.863075017929077
    22     1.221064e+03     5.186425e+01
 * time: 20.502418041229248
    23     1.220113e+03     4.306874e+01
 * time: 21.142597198486328
    24     1.218936e+03     2.178259e+01
 * time: 21.78066921234131
    25     1.217831e+03     2.144752e+01
 * time: 22.408284187316895
    26     1.216775e+03     2.201560e+01
 * time: 23.041134119033813
    27     1.215938e+03     2.969273e+01
 * time: 23.667880058288574
    28     1.214922e+03     3.438177e+01
 * time: 24.295493125915527
    29     1.212881e+03     4.013028e+01
 * time: 24.924723148345947
    30     1.207231e+03     5.214921e+01
 * time: 25.557337045669556
    31     1.197000e+03     6.567320e+01
 * time: 26.183232069015503
    32     1.189039e+03     6.178093e+01
 * time: 26.8185031414032
    33     1.173765e+03     7.561378e+01
 * time: 27.46991801261902
    34     1.171088e+03     2.385028e+01
 * time: 28.108517169952393
    35     1.170366e+03     2.359334e+01
 * time: 28.738765001296997
    36     1.170030e+03     2.394220e+01
 * time: 29.35625910758972
    37     1.169536e+03     2.287562e+01
 * time: 29.98232412338257
    38     1.168453e+03     1.921641e+01
 * time: 30.617645978927612
    39     1.166697e+03     1.661200e+01
 * time: 31.26036500930786
    40     1.164753e+03     2.146306e+01
 * time: 31.93475317955017
    41     1.163396e+03     1.589737e+01
 * time: 32.59160399436951
    42     1.162810e+03     1.879805e+01
 * time: 33.27702212333679
    43     1.162593e+03     1.922942e+01
 * time: 33.927324056625366
    44     1.162146e+03     1.865903e+01
 * time: 34.5688362121582
    45     1.161045e+03     2.060103e+01
 * time: 35.2326021194458
    46     1.158095e+03     4.011930e+01
 * time: 35.91108202934265
    47     1.150009e+03     6.883385e+01
 * time: 36.57724118232727
    48     1.147406e+03     7.331890e+01
 * time: 37.32923102378845
    49     1.144610e+03     7.504906e+01
 * time: 38.07316207885742
    50     1.138023e+03     6.294973e+01
 * time: 38.818812131881714
    51     1.132705e+03     4.283010e+01
 * time: 39.51548409461975
    52     1.129075e+03     2.327949e+01
 * time: 40.283414125442505
    53     1.128255e+03     2.383372e+01
 * time: 41.070449113845825
    54     1.127697e+03     2.546345e+01
 * time: 41.86984419822693
    55     1.126978e+03     2.642229e+01
 * time: 42.59227204322815
    56     1.126254e+03     2.679647e+01
 * time: 43.32767105102539
    57     1.124269e+03     2.481494e+01
 * time: 44.069902181625366
    58     1.121409e+03     3.925982e+01
 * time: 44.86302709579468
    59     1.118180e+03     3.554716e+01
 * time: 45.65375804901123
    60     1.116311e+03     1.691327e+01
 * time: 46.45603609085083
    61     1.115933e+03     1.725048e+01
 * time: 47.295788049697876
    62     1.115750e+03     1.702130e+01
 * time: 48.12573218345642
    63     1.115157e+03     1.748009e+01
 * time: 48.97201108932495
    64     1.113379e+03     2.744244e+01
 * time: 49.854093074798584
    65     1.107471e+03     4.999798e+01
 * time: 50.819382190704346
    66     1.099217e+03     7.740471e+01
 * time: 52.14011001586914
    67     1.098357e+03     6.707599e+01
 * time: 53.208301067352295
    68     1.093578e+03     3.715565e+01
 * time: 54.23335313796997
    69     1.090146e+03     1.668216e+01
 * time: 55.2799551486969
    70     1.088869e+03     1.082846e+01
 * time: 56.35427403450012
    71     1.087895e+03     9.101141e+00
 * time: 57.46889519691467
    72     1.086966e+03     1.225604e+01
 * time: 58.58069920539856
    73     1.086349e+03     1.200808e+01
 * time: 59.6416130065918
    74     1.086021e+03     9.321254e+00
 * time: 60.59870505332947
    75     1.085858e+03     8.113658e+00
 * time: 61.512824058532715
    76     1.085776e+03     8.127275e+00
 * time: 62.383790016174316
    77     1.085607e+03     7.969729e+00
 * time: 63.29163098335266
    78     1.085219e+03     7.491405e+00
 * time: 64.25534105300903
    79     1.084323e+03     1.248692e+01
 * time: 65.23602604866028
    80     1.082506e+03     2.097437e+01
 * time: 66.2549991607666
    81     1.079393e+03     2.578588e+01
 * time: 67.16749620437622
    82     1.077539e+03     8.916753e+00
 * time: 68.16218209266663
    83     1.077366e+03     1.950844e+00
 * time: 69.0757200717926
    84     1.077341e+03     1.114665e+00
 * time: 69.95699620246887
    85     1.077339e+03     1.131449e+00
 * time: 70.84955501556396
    86     1.077338e+03     1.121807e+00
 * time: 71.71126103401184
    87     1.077335e+03     1.078011e+00
 * time: 72.59285306930542
    88     1.077329e+03     1.045886e+00
 * time: 73.46510100364685
    89     1.077316e+03     1.385483e+00
 * time: 74.3443820476532
    90     1.077292e+03     1.538313e+00
 * time: 75.22292399406433
    91     1.077257e+03     1.233268e+00
 * time: 76.1412000656128
    92     1.077228e+03     5.728438e-01
 * time: 77.05811715126038
    93     1.077219e+03     5.059868e-01
 * time: 77.95202112197876
    94     1.077218e+03     5.033585e-01
 * time: 78.82808303833008
    95     1.077218e+03     5.028758e-01
 * time: 79.69551205635071
    96     1.077218e+03     5.022152e-01
 * time: 80.58601498603821
    97     1.077217e+03     5.008362e-01
 * time: 81.48142218589783
    98     1.077216e+03     4.977174e-01
 * time: 82.36715316772461
    99     1.077212e+03     4.905485e-01
 * time: 83.24948811531067
   100     1.077202e+03     7.608520e-01
 * time: 84.1238341331482
   101     1.077182e+03     1.030465e+00
 * time: 85.00410509109497
   102     1.077153e+03     1.039180e+00
 * time: 85.87290215492249
   103     1.077130e+03     6.005695e-01
 * time: 86.75855708122253
   104     1.077123e+03     5.073935e-01
 * time: 87.64319515228271
   105     1.077122e+03     4.817809e-01
 * time: 88.53810620307922
   106     1.077122e+03     4.675462e-01
 * time: 89.49876117706299
   107     1.077122e+03     4.482640e-01
 * time: 90.5514931678772
   108     1.077121e+03     4.167810e-01
 * time: 91.54678106307983
   109     1.077119e+03     4.160156e-01
 * time: 92.55010318756104
   110     1.077115e+03     7.459635e-01
 * time: 93.52128314971924
   111     1.077106e+03     1.228925e+00
 * time: 94.46232104301453
   112     1.077083e+03     1.876077e+00
 * time: 95.41552400588989
   113     1.077043e+03     2.711759e+00
 * time: 96.31253600120544
   114     1.076990e+03     3.472557e+00
 * time: 97.2081789970398
   115     1.076942e+03     2.744137e+00
 * time: 98.1005961894989
   116     1.076925e+03     9.131430e-01
 * time: 98.99382305145264
   117     1.076924e+03     9.225454e-01
 * time: 99.86799812316895
   118     1.076924e+03     9.211516e-01
 * time: 100.74650716781616
   119     1.076923e+03     9.099628e-01
 * time: 101.62293410301208
   120     1.076920e+03     1.072665e+00
 * time: 102.50113821029663
   121     1.076913e+03     1.902228e+00
 * time: 103.38058805465698
   122     1.076897e+03     2.867368e+00
 * time: 104.27910113334656
   123     1.076863e+03     3.341982e+00
 * time: 105.1949691772461
   124     1.076807e+03     2.150147e+00
 * time: 106.10362720489502
   125     1.076761e+03     2.996847e-01
 * time: 107.02980613708496
   126     1.076749e+03     7.551790e-01
 * time: 107.95058298110962
   127     1.076747e+03     8.967744e-01
 * time: 108.86090803146362
   128     1.076746e+03     7.027976e-01
 * time: 109.80578017234802
   129     1.076744e+03     5.080790e-01
 * time: 110.75959801673889
   130     1.076743e+03     1.635753e-01
 * time: 111.67248106002808
   131     1.076741e+03     1.794543e-01
 * time: 112.63213300704956
   132     1.076740e+03     1.890049e-01
 * time: 113.5797131061554
   133     1.076740e+03     1.923774e-01
 * time: 114.51833605766296
   134     1.076740e+03     1.972350e-01
 * time: 115.47011017799377
   135     1.076740e+03     2.001884e-01
 * time: 116.39256715774536
   136     1.076740e+03     2.031750e-01
 * time: 117.29740619659424
   137     1.076740e+03     2.057419e-01
 * time: 118.2215940952301
   138     1.076740e+03     2.081866e-01
 * time: 119.1143901348114
   139     1.076740e+03     2.083821e-01
 * time: 120.02702498435974
   140     1.076740e+03     2.000800e-01
 * time: 120.93974304199219
   141     1.076739e+03     1.697404e-01
 * time: 121.85793614387512
   142     1.076738e+03     1.576578e-01
 * time: 122.77054715156555
   143     1.076737e+03     7.985707e-02
 * time: 123.65238809585571
   144     1.076737e+03     1.937434e-02
 * time: 124.55515003204346
   145     1.076737e+03     2.909237e-03
 * time: 125.4507851600647
   146     1.076737e+03     8.634669e-04
 * time: 126.33820700645447

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             16

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1076.7369

----------------------
           Estimate
----------------------
tvcl        0.13551
tvvc        7.9849
tvka        1.1742
tvalag      0.87341
tve0       96.616
tvimax      1.0
tvic50      1.1737
tvturn     18.828
Ωpk₁,₁      0.069236
Ωpk₂,₂      0.021885
Ωpk₃,₃      0.83994
Ωpd₁,₁      0.0028192
Ωpd₂,₂      0.18162
σprop_pk    0.08854
σadd_pk     0.4189
σadd_pd     4.1536
----------------------

4 Simulation Basics

Simulation complexity increases as the conditions in a given scenario diverge from the set of conditions (e.g., population, dosage regimen) used to develop the underlying model.
We examine four scenarios with increasing complexity to introduce the user to simulations in Pumas.

4.1 Scenario 1

We begin with the simple goal of generating complete profiles for each subject in the original dataset.
- Simple because underlying population and dosage regimen are unchanged.
Two approaches, predict and simobs
- predict not technically simulation, but end result is the same, primary difference is the lack of RUV compared to simobs. The random effects will additionally be set to the empirical bayes estimates.

# sampling times from studies
stimes = [0.5, 1, 1.5, 2, 3, 6, 9, 12, 24, 36, 48, 72, 96, 120, 144]

# predictions based on individual EBEs and obstimes
# if passed, obstimes are merged with existing observation times vector
mypred = predict(myfit, myfit.data; obstimes = stimes)

Prediction
  Subjects: 32
  Predictions: conc, pca
  Covariates: wtbl

Application: NCA

With a few simple modifications, the mypred object can used for NCA. Refer to the documentation for more information on performing NCA in Pumas.

mynca = @chain DataFrame(mypred) begin
    # merge conc and conc_ipred into new column, conc_new
    transform([:conc, :conc_ipred] => ByRow((c, i) -> coalesce(c, i)) => :conc_new)
    # update route column for NCA
    @transform :route = "ev"
    # create nca population, specify obs col as conc_new
    read_nca(_; observations = :conc_new)
    run_nca(_)
end

NCA Report
     Timestamp: 2025-07-23T15:36:45.129
     Version number: 0.1.0

Output Parameters DataFrame
32×38 DataFrame
 Row │ id      dose     tlag     tmax     cmax      tlast    clast      clast_ ⋯
     │ String  Float64  Float64  Float64  Float64   Float64  Float64    Float6 ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ 1         100.0      0.5      9.0  10.8        144.0  0.124115    0.124 ⋯
   2 │ 2         100.0      0.5      9.0  11.2304     144.0  1.10587     1.093
   3 │ 3         120.0      0.5      9.0  14.4        144.0  2.25022     2.241
   4 │ 4          60.0      0.5      6.0  11.9        144.0  1.48714     1.505
   5 │ 5         113.0      0.5      3.0   8.93322    144.0  1.22051     1.200 ⋯
   6 │ 6          90.0      0.5      3.0  13.4        144.0  0.0272853   0.044
   7 │ 7         135.0      0.5      2.0  17.6        144.0  0.71204     0.873
   8 │ 8          75.0      0.5      9.0  12.9        144.0  0.808874    0.719
  ⋮  │   ⋮        ⋮        ⋮        ⋮        ⋮         ⋮         ⋮          ⋮  ⋱
  26 │ 27        120.0      0.5      6.0  15.3014     144.0  1.40167     1.548 ⋯
  27 │ 28        120.0      0.5      6.0  12.3473     144.0  1.58808     1.535
  28 │ 29        153.0      0.5      6.0  11.5899     144.0  1.24201     1.379
  29 │ 30        105.0      0.5      6.0  12.4077     144.0  1.35182     1.407
  30 │ 31        125.0      0.5      6.0  12.0114     144.0  1.59054     1.663 ⋯
  31 │ 32         93.0      0.5      6.0  11.3265     144.0  1.82322     1.816
  32 │ 33        100.0      0.5      6.0  11.638      144.0  1.75404     1.728
                                                  31 columns and 17 rows omitted

An example of a similar analysis using simobs.

simobs(myfit.model, myfit.data, coef(myfit), empirical_bayes(myfit); obstimes = stimes)

Simulated population (Vector{<:Subject})
  Simulated subjects: 32
  Simulated variables: conc, pca

Example post-processing for mypred

@chain DataFrame(mypred) begin
    # only observations
    filter(df -> df.evid == 0, _)
    # if conc missing, 1, else 0
    transform(:conc => ByRow(c -> ismissing(c) ? 1 : 0) => :isnew)
    # merge conc and conc_ipred into new column, conc_new
    transform([:conc, :conc_ipred] => ByRow((c, i) -> coalesce(c, i)) => :conc_new)
    # CT scatter plot
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc_new => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
        color = :isnew => renamer(0 => "No", 1 => "Yes") => "Predicted",
    ) *
    visual(Scatter)
    draw(
        _;
        axis = (;
            title = "Individual Concentration-Time Profiles",
            subtitle = "Warfarin ~ TAD",
        ),
    )
end

4.2 Scenario 2

In this scenario we generate observation (conc, pca) time profiles for the trial population following a one-time LD of 0.75 mg/kg PO.
- Same population, different dosage regimen lets us examine the Subject and DosageRegimen constructors without the additional complexity of creating a virtual population from scratch.

Key Concept: Subject Constructor

The Subject constructor is a fundamental part of most simulation workflows in Pumas. If you have not reviewed the corresponding tutorial, we recommend doing so before proceeding here.

Often, the best approach to building a simulation in Pumas is to focus on a single subject workflow, then, once everything is working, use a repeated-evaluation construct to complete the analysis.
In the initial setup, we showcase the mutating Subject syntax by accessing data from the first Subject stored in myfit.data[1].
- Converting the mg/kg dose to mg requires wtbl which we extract from the covariates field.

# first subject in population used in model fit
sub01 =
    Subject(myfit.data[1]; events = DosageRegimen(0.75 * myfit.data[1].covariates(0).wtbl))

sim01 = simobs(
    mdl,                        # model
    sub01,                      # subject or population of subjects                      
    coef(myfit),                # parameter estimates
    empirical_bayes(myfit)[1];  # random effects (i.e., EBEs)
    obstimes = stimes,          # obstimes for full study profile
    simulate_error = false,      # set RUV=0
)

SimulatedObservations
  Simulated variables: conc, pca
  Time: [0.5, 1.0, 1.5, 2.0, 3.0, 6.0, 9.0, 12.0, 24.0, 36.0, 48.0, 72.0, 96.0, 120.0, 144.0]

With a working single subject simulation, we can move on to simulating observations for a population.
Here, we take a slightly different approach by creating a dataframe of subject-level covariates and iterating over each row to create our population.
- Note, we could have iterated over all Subjects store in myfit.data and modified them as we did above; this syntax below shows an equivalent approach that may be more intuitive to new users.

# df with one row for each unique patient in original dataset
_patients = combine(groupby(adeval, :id), first)

# iterate over _patients creating 1 subject per row
pop02 = map(eachrow(_patients)) do r
    Subject(
        id = r.id[1],
        events = DosageRegimen(0.75 * r.wtbl[1]),
        covariates = (; wtbl = r.wtbl[1]),
    )
end

sim02 = simobs(
    mdl,
    pop02,
    coef(myfit),
    empirical_bayes(myfit);
    obstimes = stimes,
    simulate_error = false,
)

Simulated population (Vector{<:Subject})
  Simulated subjects: 32
  Simulated variables: conc, pca

A bit of additional post-processing
The figure below shows individual warfarin CT profiles for all subjects.

@chain DataFrame(sim02) begin
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
    ) *
    visual(Scatter)
    draw(
        _;
        axis = (;
            title = "Individual Concentration-Time Profiles",
            subtitle = "0.75 mg/kg x1",
        ),
    )
end

Individual CT profiles can be created by stratifying the data using the layout kwarg in mapping, and then separated using the paginate function.
x|yticklabelsize was adjusted to improve readability along those axes.

# plot layers
_plt = @chain DataFrame(sim02) begin
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
        layout = :id => nonnumeric,
    ) *
    visual(Scatter)
end

# draw(paginate(...)) returns a vector of `FigureGrid` objects
_pgrid = draw(
    paginate(_plt, layout = 16);
    figure = (;
        size = (6.5, 6.5) .* 96,
        title = "Individual Concentration-Time Profiles",
        subtitle = "0.75 mg/kg x1",
    ),
    axis = (; xticks = 0:24:144, xticklabelsize = 12, yticklabelsize = 12),
)

2-element Vector{AlgebraOfGraphics.FigureGrid}:
 FigureGrid()
 FigureGrid()

The result is a Vector{FigureGrid} with figures that can be accessed via indexing. The first 16 subjects are shown in the panel below.

_pgrid[1]

4.3 Scenario 3

In this scenario we assess the impact of augmented clearance on target attainment after a one-time LD of 1.5 mg/kg PO.
- We use this scenario to showcase creating a Subject from scratch along with the zero_randeffs helper function.

Introduction to Julia Callback

These last two scenarios should reinforce why Julia fundamentals are so important and why they were chosen for Module 1. We encourage the reader to revisit that tutorial if any of the code that follows is unclear.

# final parameter estimates
_params = coef(myfit)

sim03 = map([0.8, 1, 1.2]) do i
    simobs(
        mdl,
        Subject(
            id = "CL: $i",
            events = DosageRegimen(1.5 * 70),
            covariates = (; wtbl = 70),
        ),
        merge(_params, (; tvcl = _params.tvcl * i)),
        zero_randeffs(mdl, _params);
        obstimes = 0.5:0.5:144,
        simulate_error = false,
    )
end

Simulated population (Vector{<:Subject})
  Simulated subjects: 3
  Simulated variables: conc, pca

Visualize PK profile

@chain DataFrame(sim03) begin
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
        color = :id => "Scenario",
    ) *
    visual(Lines)
    draw(
        _;
        axis = (;
            title = "Population Concentration-Time Profiles with Augmented CL",
            subtitle = "1.5 mg/kg x1",
        ),
    )
end

Visual PD profile
More complex figures in AoG can be easier to manage if their respective layers are stored in separate variables.

# band for therapeutic range
tr_layer = mapping(0:144, 20, 35) * visual(Band; color = (:blue, 0.2))

# profiles
profiles =
    data(DataFrame(sim03)) *
    mapping(
        :time => "Time After Dose, hours",
        :pca => "PCA, % of Normal",
        group = :id => nonnumeric,
        color = :id => "Scenario",
    ) *
    visual(Lines)

draw(
    tr_layer + profiles;
    axis = (;
        title = "Population Concentration-Time Profiles with Augmented CL",
        subtitle = "1.5 mg/kg x1",
    ),
)

4.4 Scenario 4

In this scenario we combine the concepts discussed above to evaluate three alternative dosage regimens: 5, 10, or 15 mg PO daily for 14 days.
- The estimated population half-life for warfarin per our model is ~41 hours which means it should take roughly 9 days (on average) to achieve steady-state; we extend this to 14 days to ensure each of our virtual subjects is at SS prior to evaluation.
We will generate a Population of 100 Subjects and use it simulate a total of 600 trials (200 per dosage regimen).
- We will not include RUV, since most variability comes from BSV and RUV can make results difficult to interpret.
We are interested in three metrics.
- The probability of obtaining a pca within the therapeutic range (20-35%) at any time during treatment.
- The time needed to reach the first therapeutic pca value.
- The total time spent in the TR as a percentage of the dosing interval (i.e., 24 hours) at SS (Day 14).

4.4.1 Setup

We begin, as before, by developing the code for a single simulation that we can then reuse for the remaining dosage regimens. While working through the setup, we will keep our code simple by limiting our “population” and replicates to 5. We will also focus on one dosage regimen, 5 mg PO daily for 14 days. This will allow us to spot check our code and the results to ensure the output it what we expect instead of trying to troubleshoot for the full population and profile. In order, we must:

Create a population of subjects that has a single covariate (wtbl) that is sampled from a uniform distribution of observed values (40-102 kg).
Simulate an appropriate number of observations (hourly observations will be sufficient).
Repeat the simulation in #2 for a total of 5 simulations.
Store the output in a format that will make post-processing and evaluation as easy as possible.
Process and evaluate the result then present our findings in a meaningful way.

4.4.1.1 Population

We can combine map with a Range between 1:n and a do-block to create a vector of virtual subjects (i.e., a Population)
For simplicity, we will also use scalar literal values for the range of wtbl in the Uniform call instead of obtaining them programmatically with extrema or some other function.

pop03 = map(1:5) do i
    _wtbl = rand(Uniform(40, 120))
    Subject(
        id = i,
        events = DosageRegimen(5; ii = 24, addl = 13),
        covariates = (; wtbl = _wtbl),
    )
end

Population
  Subjects: 5
  Covariates: wtbl
  Observations:

4.4.1.2 Simulation

Simulating observations for 2 subjects is comparatively simple, and we can repeat that simulation by mapping over a range 1:n as we did when creating pop03.
The resulting vector of SimulatedObservation objects can be concatenated into a single object using the reduce(vcat, myvectorofsims), then converted into a data frame for post-processing.
Since we need to summarize values from each simulation, we will need to include a variable to track the iteration number for each simulation. We can do this by leveraging the mutating Subject syntax to add a rep_id as a covariate for each subject pop03 inside the map before we call simobs.

goodsim = map(1:5) do i
    # rebuild pop03 using mutating Subject to add rep_id
    _pop = map(pop03) do s
        Subject(s; covariates = (; rep_id = i))
    end
    # simulate values for _pop 
    simobs(mdl, _pop, coef(myfit); obstimes = 0:1:(24*14), simulate_error = false)
end

5-element Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvvc = 1, tvka = 1, tvalag = 1, tve0 = 1, tvimax = 1, tvic50 = 1, tvturn = 1, Ωpk = 3, Ωpd = 2, σprop_pk = 1, σadd_pk = 1, σadd_pd = 1), 5, (:depot, :central, :e), (:cl, :vc, :ka, :e0, :imax, :ic50, :turn, :kout, :kin), ParamSet{@NamedTuple{tvcl::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvvc::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvka::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvalag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvturn::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pd::RealDomain{Float64, TransformVariables.Infinity{true}, Int64}}}, Pumas.RandomObj{(), var"#1#12"}, Pumas.TimeDispatcher{var"#2#13", var"#3#14"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#5#16", var"#6#17"}}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, StepRange{Int64, Int64}, @NamedTuple{tvcl::Float64, tvvc::Float64, tvka::Float64, tvalag::Float64, tve0::Float64, tvimax::Float64, tvic50::Float64, tvturn::Float64, Ωpk::PDMats.PDiagMat{Float64, Vector{Float64}}, Ωpd::PDMats.PDiagMat{Float64, Vector{Float64}}, σprop_pk::Float64, σadd_pk::Float64, σadd_pd::Float64}, @NamedTuple{ηpk::Vector{Float64}, ηpd::Vector{Float64}}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}, @NamedTuple{saveat::StepRange{Int64, Int64}, ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, SciMLBase.ODESolution{Float64, 2, Vector{StaticArraysCore.SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{StaticArraysCore.SVector{3, Float64}}}, Nothing, Any, Any, Any, SciMLBase.DEStats, Vector{Int64}, Nothing, Nothing, Nothing}, @NamedTuple{lags::Vector{@NamedTuple{depot::Float64}}}, @NamedTuple{cl::Vector{Float64}, vc::Vector{Float64}, ka::Vector{Float64}, e0::Vector{Float64}, imax::Vector{Float64}, ic50::Vector{Float64}, turn::Vector{Float64}, kout::Vector{Float64}, kin::Vector{Float64}}, @NamedTuple{conc::Vector{Float64}, pca::Vector{Float64}}, @NamedTuple{}}}}:
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca

We used the label goodsims because this approach, while valid, is redundant because we are recreating the Population from scratch with each simulation.
We could simplify the code by using a nested map call to create the population during each iteration.

bettersim = map(1:5) do i
    _pop = map(1:5) do s
        Subject(
            id = s,
            events = DosageRegimen(5, ii = 24, addl = 13),
            covariates = (; wtbl = rand(Uniform(40, 102)), rep_id = i),
        )
    end
    simobs(mdl, _pop, coef(myfit); obstimes = 0:1:(24*14), simulate_error = false)
end

5-element Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvvc = 1, tvka = 1, tvalag = 1, tve0 = 1, tvimax = 1, tvic50 = 1, tvturn = 1, Ωpk = 3, Ωpd = 2, σprop_pk = 1, σadd_pk = 1, σadd_pd = 1), 5, (:depot, :central, :e), (:cl, :vc, :ka, :e0, :imax, :ic50, :turn, :kout, :kin), ParamSet{@NamedTuple{tvcl::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvvc::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvka::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvalag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvturn::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pd::RealDomain{Float64, TransformVariables.Infinity{true}, Int64}}}, Pumas.RandomObj{(), var"#1#12"}, Pumas.TimeDispatcher{var"#2#13", var"#3#14"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#5#16", var"#6#17"}}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{}, Pumas.ConstantCovar{@NamedTuple{wtbl::Float64, rep_id::Int64}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, StepRange{Int64, Int64}, @NamedTuple{tvcl::Float64, tvvc::Float64, tvka::Float64, tvalag::Float64, tve0::Float64, tvimax::Float64, tvic50::Float64, tvturn::Float64, Ωpk::PDMats.PDiagMat{Float64, Vector{Float64}}, Ωpd::PDMats.PDiagMat{Float64, Vector{Float64}}, σprop_pk::Float64, σadd_pk::Float64, σadd_pd::Float64}, @NamedTuple{ηpk::Vector{Float64}, ηpd::Vector{Float64}}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}, @NamedTuple{saveat::StepRange{Int64, Int64}, ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, SciMLBase.ODESolution{Float64, 2, Vector{StaticArraysCore.SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{StaticArraysCore.SVector{3, Float64}}}, Nothing, Any, Any, Any, SciMLBase.DEStats, Vector{Int64}, Nothing, Nothing, Nothing}, @NamedTuple{lags::Vector{@NamedTuple{depot::Float64}}}, @NamedTuple{cl::Vector{Float64}, vc::Vector{Float64}, ka::Vector{Float64}, e0::Vector{Float64}, imax::Vector{Float64}, ic50::Vector{Float64}, turn::Vector{Float64}, kout::Vector{Float64}, kin::Vector{Float64}}, @NamedTuple{conc::Vector{Float64}, pca::Vector{Float64}}, @NamedTuple{}}}}:
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca

The bettersim syntax provides a reasonable solution for simulating a single dose, now we just need to abstract that code out into a function so that we can apply it to our two remaining doses.
We create a function, simulate_warfarin that accepts a single positional argument, dose that we can use along with map.
- dose was also added as a covariate in the Subject constructor so that we can use it for stratification during post-processing.
- The values for number of subjects (100), and number of samples (200) were hard-coded for simplicity. In a real-world application it would be better to pass those parameters as arguments to simulate_warfarin to improve its overall utility.

function simulate_warfarin(dose)
    #! using literal values for n samples
    _sim = map(1:200) do s
        # create a population
        #! using literal for n subjects
        _pop = map(1:100) do p
            Subject(
                id = p,
                #! using literals for dosing frequency and duration
                events = DosageRegimen(dose; ii = 24, addl = 13),
                #! using literals for wtbl range, adding rep_id and dose as covariates
                covariates = (; wtbl = rand(Uniform(40, 102)), rep_id = s, dose),
            )
        end
        # simulation
        simobs(
            mdl,
            _pop,
            coef(myfit);
            #! using literal for n days in obstimes range
            obstimes = 0:1:(24*14),
            #! no RUV
            simulate_error = false,
        )
    end
end

simulate_warfarin (generic function with 1 method)

Lastly, we perform the simulations using a mapreduce call and save the result in a variable, sim03.
- mapreduce allows us to combine the map and reduce(vcat) steps into a single function call.

sim03 = mapreduce(simulate_warfarin, vcat, [5, 10, 15])

600-element Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvvc = 1, tvka = 1, tvalag = 1, tve0 = 1, tvimax = 1, tvic50 = 1, tvturn = 1, Ωpk = 3, Ωpd = 2, σprop_pk = 1, σadd_pk = 1, σadd_pd = 1), 5, (:depot, :central, :e), (:cl, :vc, :ka, :e0, :imax, :ic50, :turn, :kout, :kin), ParamSet{@NamedTuple{tvcl::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvvc::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvka::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvalag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvturn::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pd::RealDomain{Float64, TransformVariables.Infinity{true}, Int64}}}, Pumas.RandomObj{(), var"#1#12"}, Pumas.TimeDispatcher{var"#2#13", var"#3#14"}, Pumas.DCPChecker{Pumas.TimeDispatcher{var"#5#16", var"#6#17"}}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{}, Pumas.ConstantCovar{@NamedTuple{wtbl::Float64, rep_id::Int64, dose::Int64}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, StepRange{Int64, Int64}, @NamedTuple{tvcl::Float64, tvvc::Float64, tvka::Float64, tvalag::Float64, tve0::Float64, tvimax::Float64, tvic50::Float64, tvturn::Float64, Ωpk::PDMats.PDiagMat{Float64, Vector{Float64}}, Ωpd::PDMats.PDiagMat{Float64, Vector{Float64}}, σprop_pk::Float64, σadd_pk::Float64, σadd_pd::Float64}, @NamedTuple{ηpk::Vector{Float64}, ηpd::Vector{Float64}}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}, dose::Vector{Int64}}, @NamedTuple{saveat::StepRange{Int64, Int64}, ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, SciMLBase.ODESolution{Float64, 2, Vector{StaticArraysCore.SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{StaticArraysCore.SVector{3, Float64}}}, Nothing, Any, Any, Any, SciMLBase.DEStats, Vector{Int64}, Nothing, Nothing, Nothing}, @NamedTuple{lags::Vector{@NamedTuple{depot::Float64}}}, @NamedTuple{cl::Vector{Float64}, vc::Vector{Float64}, ka::Vector{Float64}, e0::Vector{Float64}, imax::Vector{Float64}, ic50::Vector{Float64}, turn::Vector{Float64}, kout::Vector{Float64}, kin::Vector{Float64}}, @NamedTuple{conc::Vector{Float64}, pca::Vector{Float64}}, @NamedTuple{}}}}:
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 ⋮
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca

4.4.1.3 Post-processing

The starting point for post-processing will depend on the output needed to answer the question of interest; in this case a simple tabular summary of metrics with 95% CIs and a graphical summary of 90% PIs for each regimen will suffice.
We will start from a data frame (sim03df)

#! takes ~3-5min on 16vCPU
sim03df = DataFrame(reduce(vcat, (sim03)))

21060000×35 DataFrame

21059975 rows omitted

Row	id	time	conc	pca	evid	lags_depot	amt	cmt	rate	duration	ss	ii	route	wtbl	rep_id	dose	tad	dosenum	depot	central	e	cl	vc	ka	e0	imax	ic50	turn	kout	kin	ηpk_1	ηpk_2	ηpk_3	ηpd_1	ηpd_2
	String	Float64	Float64?	Float64?	Int64	Float64?	Float64?	Symbol?	Float64?	Float64?	Int8?	Float64?	NCA.Route?	Float64?	Int64?	Int64?	Float64	Int64	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64	Float64	Float64	Float64	Float64
1	1	0.0	missing	missing	1	0.873414	5.0	depot	0.0	0.0	0	0.0	NullRoute	59.489	1	5	0.0	1	0.0	0.0	96.5153	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
2	1	0.0	0.0	96.5153	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	0.0	1	0.0	0.0	96.5153	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
3	1	1.0	0.0638109	96.5077	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	1.0	1	4.52001	0.479522	96.5077	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
4	1	2.0	0.390486	96.0931	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	2.0	1	2.0365	2.9344	96.0931	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
5	1	3.0	0.532183	95.3803	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	3.0	1	0.917548	3.99922	95.3803	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
6	1	4.0	0.590621	94.5822	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	4.0	1	0.413403	4.43837	94.5822	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
7	1	5.0	0.611629	93.7794	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	5.0	1	0.186259	4.59624	93.7794	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
8	1	6.0	0.615853	93.004	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	6.0	1	0.0839196	4.62798	93.004	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
9	1	7.0	0.612595	92.2689	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	7.0	1	0.0378101	4.6035	92.2689	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
10	1	8.0	0.606044	91.578	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	8.0	1	0.0170354	4.55427	91.578	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
11	1	9.0	0.598087	90.9319	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	9.0	1	0.00767534	4.49447	90.9319	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
12	1	10.0	0.589572	90.3294	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	10.0	1	0.00345814	4.43048	90.3294	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
13	1	11.0	0.58088	89.7687	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	59.489	1	5	11.0	1	0.00155807	4.36516	89.7687	0.114856	7.51474	0.797282	96.5153	1.0	2.73226	18.8281	0.0531121	5.12613	-0.0433602	0.10202	-0.387166	-0.0010379	0.844984
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
21059989	100	325.0	2.8743	25.2125	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	13.0	14	0.178912	24.8079	25.2125	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059990	100	326.0	2.80564	25.1936	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	14.0	14	0.124171	24.2152	25.1936	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059991	100	327.0	2.73685	25.1992	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	15.0	14	0.0861784	23.6216	25.1992	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059992	100	328.0	2.66853	25.2288	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	16.0	14	0.0598106	23.0319	25.2288	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059993	100	329.0	2.60107	25.2817	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	17.0	14	0.0415105	22.4496	25.2817	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059994	100	330.0	2.53473	25.3573	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	18.0	14	0.0288097	21.877	25.3573	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059995	100	331.0	2.46967	25.455	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	19.0	14	0.0199948	21.3155	25.455	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059996	100	332.0	2.406	25.574	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	20.0	14	0.0138771	20.766	25.574	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059997	100	333.0	2.34377	25.7136	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	21.0	14	0.00963115	20.2289	25.7136	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059998	100	334.0	2.28302	25.8731	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	22.0	14	0.00668433	19.7046	25.8731	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21059999	100	335.0	2.22375	26.0518	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	23.0	14	0.00463914	19.193	26.0518	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908
21060000	100	336.0	2.16595	26.2491	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	98.9213	200	15	24.0	14	0.00321972	18.6941	26.2491	0.227941	8.63091	0.365236	95.3328	1.0	1.00107	18.8281	0.0531121	5.06333	0.260652	-0.268029	-1.16783	-0.0133653	-0.15908

Sanity Check

We can quickly check the observation-time profile(s) for a single subject to limit the risk of down-stream errors as we continue our analysis. The profile in the figure below appears reasonable.

@chain sim03df begin
    filter(df -> df.id == "1" && df.rep_id == 1, _)
    filter(df -> df.evid == 0, _)
    select(:time, :conc, :pca)
    # default col names for stack are variable and value
    stack(Not([:time]))
    data(_) * mapping(:time, :value, row = :variable) * visual(Lines)
    draw(_; facet = (; linkyaxes = false))
end

Our tabular summary will include three metrics of interest (TA, TTA, TTR) which will be evaluated for each subject, in each simulation.
We will take the average for each metric per simulation and then report the relevant percentiles (2.5, 50, 97.5).
Since there is effectively one assessment per subject, we can make use of the split-apply-combine design for data frames.
The war_metrics function is a custom analysis for table_one; see the SummaryTables.jl documentation for details.

# custom analysis function for table_one
function war_metrics(col)
    all(ismissing, col) && return ("-" => "Median", "-" => "95% CI")
    (
        median(col) => "Median",
        Concat("[", quantile(col, 0.025), ", ", quantile(col, 0.975), "]") => "95% CI",
    )
end

@chain sim03df begin
    # drop records where pca is missing
    dropmissing(_, :pca)
    # first combine step evaluates metrics for individual subjects
    combine(groupby(_, [:dose, :rep_id, :id])) do gdf
        #! metrics 1 and 2
        # find index of first pca value in TR; returns index or nothing 
        i = findfirst(x -> 20 ≤ x < 35, gdf.pca)
        # if `i` was found, return 1 (true), else 0 (false)
        ta_i = Int(!isnothing(i))
        # if no index was found, return missing, else return corresponding time
        tta_i = isnothing(i) ? missing : gdf.time[i]

        #! metric 3
        # temporary df of SS obs from start of Day 14 (312 hours) that are in TR
        _ssdf = filter(df -> df.time >= 312 && 20 ≤ df.pca < 35, gdf)
        # if no obs found, return missing, else return TTR as percentage of ii
        ttr_i =
            iszero(nrow(_ssdf)) ? missing :
            ((last(_ssdf.time) - first(_ssdf.time)) / 24) * 100

        # return a named tuple of the 3 metrics for each subject
        return (; ta_i, tta_i, ttr_i)
    end
    # second combine summarizes each metric per simulation (rep_id)
    combine(
        groupby(_, [:dose, :rep_id]),
        # mean(0|1) * 100 = TA percentage
        :ta_i => (x -> mean(x) * 100) => :ta,
        # applies anonymous function to tta_i and ttr_i cols
        # possible all values could be missing, else could have just used `mean`
        [:tta_i, :ttr_i] .=> function (c)
            all(ismissing, c) && return missing
            mean(skipmissing(c))
        end .=> [:tta, :ttr],
    )
    # from SummaryTables.jl
    table_one(
        _,
        [
            :ta => war_metrics => "Probability of TA",
            :tta => war_metrics => "Time to Target",
            :ttr => war_metrics => "Time in TR",
        ],
        sort = false,
        groupby = :dose => "Dose, mg",
        show_total = false,
    )
end


	Dose, mg
	5	10	15

Probability of TA
Median	31	78	94
95% CI	[23, 39]	[70, 86]	[88, 98]
Time to Target
Median	124	85.4	62.9
95% CI	[105, 145]	[74.5, 95.8]	[56.9, 69.9]
Time in TR
Median	94.9	96.7	97.2
95% CI	[87.6, 100]	[91.7, 99.6]	[92, 100]

The tabular summary focuses on average response, the graphical summary provide a better understanding of the range of predicted values that we might expect.
We will summarize the relevant pca percentiles (5, 50, 90%) at each time point following the Day 14 dose.

_tbl = @chain sim03df begin
    dropmissing(_, :pca)
    # Day 14 (SS) observations only
    filter(df -> df.time >= 312, _)
    # Summarize by dose and tad; l, m, h are 5th,50th,95th percentile
    combine(
        groupby(_, [:dose, :tad]),
        :pca => (x -> quantile(x, [0.05, 0.5, 0.9])') => [:l, :m, :h],
    )
end

# plot layers
median_layer = mapping(:tad, :m) * visual(Lines; linewidth = 2)
pi_layer = mapping(:tad, :l, :h) * visual(Band, alpha = 0.2)
tr_layer =
    mapping([20, 35], color = "TR" => AlgebraOfGraphics.scale(:secondary)) *
    visual(HLines; linestyle = :dash, alpha = 0.5)

# color and facet map
cf_map = mapping(
    color = :dose => renamer(5 => "5mg", 10 => "10mg", 15 => "15mg") => "",
    col = :dose => renamer(5 => "5mg", 10 => "10mg", 15 => "15mg"),
)

# combine layers and draw
(data(_tbl) * (pi_layer + median_layer) * cf_map) + tr_layer |> draw(
    scales(;
        Y = (; label = "PCA, % of normal"),
        X = (; label = "Time after previous dose, hours"),
        secondary = (; palette = [:gray30]),
    );
    figure = (;
        size = (6, 4) .* 96,
        title = "Predicted PCA-Time Profiles at Steady-state (Day 14) by Dose",
        subtitle = "Median (line), 90%PI (band), TR (dash-line)",
        titlealign = :left,
    ),
    axis = (;
        limits = (0, 24, 0, 80),
        xticks = [0, 12, 24],
        xlabelpadding = 10,
        yticks = 0:20:80,
    ),
    legend = (; orientation = :horizontal, framevisible = false, position = :bottom),
)

4.5 Evaluation

10 mg PO daily dosage regimen offers a reasonable balance of TA, TTA, TTR.

5 Conclusion

Presented the basics of simulation in Pumas using several examples that utilize built-in functionality and user-defined functions.

Reuse

CC BY-SA 4.0