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.04047703742980957
     1     1.762032e+03     1.253034e+03
 * time: 2.811324119567871
     2     1.449413e+03     6.274412e+02
 * time: 3.8159611225128174
     3     1.331405e+03     2.099880e+02
 * time: 4.813470125198364
     4     1.310843e+03     1.627399e+02
 * time: 5.8135600090026855
     5     1.295073e+03     1.376531e+02
 * time: 8.53337812423706
     6     1.282335e+03     1.324314e+02
 * time: 9.60521411895752
     7     1.275524e+03     9.521788e+01
 * time: 10.600411176681519
     8     1.271552e+03     5.354883e+01
 * time: 11.54372501373291
     9     1.269315e+03     2.407046e+01
 * time: 12.508581161499023
    10     1.268615e+03     2.368964e+01
 * time: 13.456582069396973
    11     1.268257e+03     2.366048e+01
 * time: 14.390335083007812
    12     1.267352e+03     3.190950e+01
 * time: 15.326550006866455
    13     1.265336e+03     5.399187e+01
 * time: 16.22954511642456
    14     1.260582e+03     8.775888e+01
 * time: 17.111937999725342
    15     1.252929e+03     1.110735e+02
 * time: 17.928004026412964
    16     1.246526e+03     8.328859e+01
 * time: 18.71481204032898
    17     1.240724e+03     8.711808e+01
 * time: 19.482946157455444
    18     1.236881e+03     7.439216e+01
 * time: 20.259448051452637
    19     1.232291e+03     6.878129e+01
 * time: 21.0018630027771
    20     1.227181e+03     5.604034e+01
 * time: 21.73366403579712
    21     1.223093e+03     3.856728e+01
 * time: 22.458101987838745
    22     1.221064e+03     5.186425e+01
 * time: 23.182790994644165
    23     1.220113e+03     4.306874e+01
 * time: 23.906623125076294
    24     1.218936e+03     2.178259e+01
 * time: 24.638844966888428
    25     1.217831e+03     2.144752e+01
 * time: 25.36596703529358
    26     1.216775e+03     2.201560e+01
 * time: 26.090232133865356
    27     1.215938e+03     2.969273e+01
 * time: 26.805341005325317
    28     1.214922e+03     3.438177e+01
 * time: 27.51904010772705
    29     1.212881e+03     4.013028e+01
 * time: 28.237519025802612
    30     1.207231e+03     5.214921e+01
 * time: 28.976107120513916
    31     1.197000e+03     6.567320e+01
 * time: 29.709606170654297
    32     1.189039e+03     6.178093e+01
 * time: 30.434771060943604
    33     1.173765e+03     7.561378e+01
 * time: 31.150236129760742
    34     1.171088e+03     2.385028e+01
 * time: 31.87293815612793
    35     1.170366e+03     2.359334e+01
 * time: 32.58821105957031
    36     1.170030e+03     2.394220e+01
 * time: 33.29846906661987
    37     1.169536e+03     2.287562e+01
 * time: 34.00157713890076
    38     1.168453e+03     1.921641e+01
 * time: 34.72014498710632
    39     1.166697e+03     1.661200e+01
 * time: 35.438637018203735
    40     1.164753e+03     2.146306e+01
 * time: 36.15893816947937
    41     1.163396e+03     1.589737e+01
 * time: 36.891162157058716
    42     1.162810e+03     1.879805e+01
 * time: 37.60688018798828
    43     1.162593e+03     1.922942e+01
 * time: 38.333057165145874
    44     1.162146e+03     1.865903e+01
 * time: 39.04715299606323
    45     1.161045e+03     2.060103e+01
 * time: 39.76336598396301
    46     1.158095e+03     4.011930e+01
 * time: 40.479026079177856
    47     1.150009e+03     6.883385e+01
 * time: 41.21449518203735
    48     1.147406e+03     7.331890e+01
 * time: 42.05571508407593
    49     1.144610e+03     7.504906e+01
 * time: 42.900644063949585
    50     1.138023e+03     6.294973e+01
 * time: 43.67860698699951
    51     1.132705e+03     4.283010e+01
 * time: 44.45889616012573
    52     1.129075e+03     2.327949e+01
 * time: 45.24944615364075
    53     1.128255e+03     2.383372e+01
 * time: 46.0436110496521
    54     1.127697e+03     2.546345e+01
 * time: 46.829432010650635
    55     1.126978e+03     2.642229e+01
 * time: 47.6283540725708
    56     1.126254e+03     2.679647e+01
 * time: 48.583171129226685
    57     1.124269e+03     2.481494e+01
 * time: 49.88638520240784
    58     1.121409e+03     3.925982e+01
 * time: 51.87938404083252
    59     1.118180e+03     3.554716e+01
 * time: 52.86155319213867
    60     1.116311e+03     1.691327e+01
 * time: 53.70243310928345
    61     1.115933e+03     1.725048e+01
 * time: 54.581916093826294
    62     1.115750e+03     1.702130e+01
 * time: 55.55463218688965
    63     1.115157e+03     1.748009e+01
 * time: 56.4709370136261
    64     1.113379e+03     2.744244e+01
 * time: 57.32971715927124
    65     1.107471e+03     4.999798e+01
 * time: 58.26607608795166
    66     1.099217e+03     7.740471e+01
 * time: 59.567793130874634
    67     1.098357e+03     6.707599e+01
 * time: 60.75169897079468
    68     1.093578e+03     3.715565e+01
 * time: 61.95774602890015
    69     1.090146e+03     1.668216e+01
 * time: 63.16390919685364
    70     1.088869e+03     1.082846e+01
 * time: 64.35664701461792
    71     1.087895e+03     9.101141e+00
 * time: 65.48611307144165
    72     1.086966e+03     1.225604e+01
 * time: 66.59634399414062
    73     1.086349e+03     1.200808e+01
 * time: 67.59993100166321
    74     1.086021e+03     9.321254e+00
 * time: 68.57169318199158
    75     1.085858e+03     8.113658e+00
 * time: 69.52956700325012
    76     1.085776e+03     8.127275e+00
 * time: 70.4882071018219
    77     1.085607e+03     7.969729e+00
 * time: 71.46570301055908
    78     1.085219e+03     7.491405e+00
 * time: 72.446613073349
    79     1.084323e+03     1.248692e+01
 * time: 73.44785404205322
    80     1.082506e+03     2.097437e+01
 * time: 74.46820712089539
    81     1.079393e+03     2.578588e+01
 * time: 75.50782918930054
    82     1.077539e+03     8.916753e+00
 * time: 76.53583407402039
    83     1.077366e+03     1.950844e+00
 * time: 77.53611207008362
    84     1.077341e+03     1.114665e+00
 * time: 78.54902911186218
    85     1.077339e+03     1.131449e+00
 * time: 79.54549908638
    86     1.077338e+03     1.121807e+00
 * time: 80.54905700683594
    87     1.077335e+03     1.078011e+00
 * time: 81.53591418266296
    88     1.077329e+03     1.045886e+00
 * time: 82.51691317558289
    89     1.077316e+03     1.385483e+00
 * time: 83.50935316085815
    90     1.077292e+03     1.538313e+00
 * time: 84.51403403282166
    91     1.077257e+03     1.233268e+00
 * time: 85.54174304008484
    92     1.077228e+03     5.728438e-01
 * time: 86.55996298789978
    93     1.077219e+03     5.059868e-01
 * time: 87.56942009925842
    94     1.077218e+03     5.033585e-01
 * time: 88.57844400405884
    95     1.077218e+03     5.028758e-01
 * time: 89.56705713272095
    96     1.077218e+03     5.022152e-01
 * time: 90.57410097122192
    97     1.077217e+03     5.008362e-01
 * time: 91.5552351474762
    98     1.077216e+03     4.977174e-01
 * time: 92.54172611236572
    99     1.077212e+03     4.905485e-01
 * time: 93.52854800224304
   100     1.077202e+03     7.608520e-01
 * time: 94.52535820007324
   101     1.077182e+03     1.030465e+00
 * time: 95.50503706932068
   102     1.077153e+03     1.039180e+00
 * time: 96.47164011001587
   103     1.077130e+03     6.005695e-01
 * time: 97.47609996795654
   104     1.077123e+03     5.073935e-01
 * time: 98.5042941570282
   105     1.077122e+03     4.817809e-01
 * time: 99.48327612876892
   106     1.077122e+03     4.675462e-01
 * time: 100.48115110397339
   107     1.077122e+03     4.482640e-01
 * time: 101.47480297088623
   108     1.077121e+03     4.167810e-01
 * time: 102.46560215950012
   109     1.077119e+03     4.160156e-01
 * time: 103.45274114608765
   110     1.077115e+03     7.459635e-01
 * time: 104.43729615211487
   111     1.077106e+03     1.228925e+00
 * time: 105.4296350479126
   112     1.077083e+03     1.876077e+00
 * time: 106.441721200943
   113     1.077043e+03     2.711759e+00
 * time: 107.44700908660889
   114     1.076990e+03     3.472557e+00
 * time: 108.44899201393127
   115     1.076942e+03     2.744137e+00
 * time: 109.46661615371704
   116     1.076925e+03     9.131430e-01
 * time: 110.48312616348267
   117     1.076924e+03     9.225454e-01
 * time: 111.48173999786377
   118     1.076924e+03     9.211516e-01
 * time: 112.46145606040955
   119     1.076923e+03     9.099628e-01
 * time: 113.47241806983948
   120     1.076920e+03     1.072665e+00
 * time: 114.48391008377075
   121     1.076913e+03     1.902228e+00
 * time: 115.49304103851318
   122     1.076897e+03     2.867368e+00
 * time: 116.49854111671448
   123     1.076863e+03     3.341982e+00
 * time: 117.51395416259766
   124     1.076807e+03     2.150147e+00
 * time: 118.53945899009705
   125     1.076761e+03     2.996847e-01
 * time: 119.55092716217041
   126     1.076749e+03     7.551790e-01
 * time: 120.55191111564636
   127     1.076747e+03     8.967744e-01
 * time: 121.5616180896759
   128     1.076746e+03     7.027976e-01
 * time: 122.56762409210205
   129     1.076744e+03     5.080790e-01
 * time: 123.57818007469177
   130     1.076743e+03     1.635753e-01
 * time: 124.59108901023865
   131     1.076741e+03     1.794543e-01
 * time: 125.61469507217407
   132     1.076740e+03     1.890049e-01
 * time: 126.60721898078918
   133     1.076740e+03     1.923774e-01
 * time: 127.62169599533081
   134     1.076740e+03     1.972350e-01
 * time: 128.6147871017456
   135     1.076740e+03     2.001884e-01
 * time: 129.59330797195435
   136     1.076740e+03     2.031750e-01
 * time: 130.58932304382324
   137     1.076740e+03     2.057419e-01
 * time: 131.57091116905212
   138     1.076740e+03     2.081866e-01
 * time: 132.56923508644104
   139     1.076740e+03     2.083821e-01
 * time: 133.56802415847778
   140     1.076740e+03     2.000800e-01
 * time: 134.56514716148376
   141     1.076739e+03     1.697404e-01
 * time: 135.57195401191711
   142     1.076738e+03     1.576578e-01
 * time: 136.69113898277283
   143     1.076737e+03     7.985707e-02
 * time: 138.14655900001526
   144     1.076737e+03     1.937434e-02
 * time: 139.71912217140198
   145     1.076737e+03     2.909237e-03
 * time: 140.88156914710999
   146     1.076737e+03     8.634669e-04
 * time: 141.97339797019958

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-09-24T08:39:52.332
     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	93.0704	1	5	0.0	1	0.0	0.0	98.8862	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
2	1	0.0	0.0	98.8862	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	0.0	1	0.0	0.0	98.8862	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
3	1	1.0	0.0459731	98.8716	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	1.0	1	4.53077	0.468697	98.8716	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
4	1	2.0	0.283193	98.1235	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	2.0	1	2.08009	2.88717	98.1235	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
5	1	3.0	0.38755	96.8915	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	3.0	1	0.95497	3.95109	96.8915	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
6	1	4.0	0.430987	95.5367	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	4.0	1	0.438428	4.39393	95.5367	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
7	1	5.0	0.446534	94.1834	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	5.0	1	0.201283	4.55243	94.1834	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
8	1	6.0	0.449353	92.8799	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	6.0	1	0.0924094	4.58117	92.8799	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
9	1	7.0	0.446403	91.6451	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	7.0	1	0.0424253	4.5511	91.6451	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
10	1	8.0	0.440879	90.4847	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	8.0	1	0.0194775	4.49478	90.4847	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
11	1	9.0	0.434245	89.3994	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	9.0	1	0.00894216	4.42714	89.3994	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
12	1	10.0	0.427173	88.3871	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	10.0	1	0.00410536	4.35504	88.3871	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
13	1	11.0	0.419969	87.445	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	93.0704	1	5	11.0	1	0.00188478	4.2816	87.445	0.178623	10.195	0.778484	98.8862	1.0	1.03506	18.8281	0.0531121	5.25206	0.0625663	-0.0405094	-0.411026	0.02323	-0.125685
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
21059989	100	325.0	4.16693	29.0372	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	13.0	14	7.98157e-12	32.3599	29.0372	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059990	100	326.0	4.08758	29.0631	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	14.0	14	7.7542e-13	31.7437	29.0631	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059991	100	327.0	4.00974	29.1081	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	15.0	14	7.3914e-14	31.1392	29.1081	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059992	100	328.0	3.93339	29.1715	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	16.0	14	6.95956e-15	30.5463	29.1715	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059993	100	329.0	3.85849	29.2524	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	17.0	14	3.76274e-16	29.9646	29.2524	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059994	100	330.0	3.78501	29.3501	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	18.0	14	-1.18625e-16	29.394	29.3501	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059995	100	331.0	3.71294	29.4638	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	19.0	14	1.55074e-16	28.8343	29.4638	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059996	100	332.0	3.64224	29.5929	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	20.0	14	-1.16109e-17	28.2852	29.5929	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059997	100	333.0	3.57288	29.7367	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	21.0	14	-2.40852e-15	27.7466	29.7367	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059998	100	334.0	3.50484	29.8946	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	22.0	14	-1.41987e-15	27.2182	29.8946	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21059999	100	335.0	3.4381	30.0659	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	23.0	14	1.82787e-16	26.6999	30.0659	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925
21060000	100	336.0	3.37263	30.2502	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	64.3371	200	15	24.0	14	2.53727e-16	26.1915	30.2502	0.149306	7.76589	2.33055	97.0326	1.0	1.80592	18.8281	0.0531121	5.15361	0.160206	0.0565492	0.685484	0.00430778	0.430925

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	30	78	94
95% CI	[23, 39]	[70, 86]	[89, 98]
Time to Target
Median	124	84.8	62.8
95% CI	[101, 142]	[75.5, 94.3]	[56.5, 70.9]
Time in TR
Median	95.5	96.7	97.3
95% CI	[88.2, 99.3]	[91.3, 99.6]	[92.1, 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