using Pumas, NCA, PharmaDatasets
using DataFramesMeta, CategoricalArrays
using CairoMakie, AlgebraOfGraphics
using SummaryTables
Introduction to Simulations in Pumas
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.
modelexpects anAbstractPumasModel, which (for now) refers to a model defined using the@modelmacro.populationaccepts aPopulationwhich was discussed in Module 4), or a singleSubjectwhich will be discussed later in the current module.paramshould be a single parameter set defined as aNamedTupleor 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.
randeffsis 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.ensemblealgis used to select the parallelization mode to be used for the simulation.diffeq_optionscan be used to pass additional options to the differential equation solver if the model does not have an analytical solution.rngcan be used to specify the random number generator to be used for the simulation.simulate_errorcan be used to disable (false) the inclusion of RUV in the value returned by the predictive distribution’s error model.
ImportantReproducibility
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.
NoteWarfarin Model
#* read dataset from PharmaDatasets
adeval = dataset("pumas/warfarin_pumas")
# Concentrations of exactly zero are below the limit of quantification
@rtransform! adeval @passmissing :conc = iszero(:conc) ? missing : :conc
# population
mypop = read_pumas(adeval; observations = [:conc, :pca], covariates = [:wtbl])
# warfarin model
mdl = @model begin
@metadata begin
desc = "Integrated Warfarin PK/PD model"
timeu = "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 ~ @. CombinedNormal(cp, σadd_pk, σprop_pk)
# 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.483322e+03 7.050748e+03 * time: 0.05883193016052246 1 1.715008e+03 1.254170e+03 * time: 46.81045699119568 2 1.403854e+03 6.276985e+02 * time: 47.98807191848755 3 1.288196e+03 2.125800e+02 * time: 48.80350399017334 4 1.268954e+03 1.661051e+02 * time: 49.562845945358276 5 1.254637e+03 1.317251e+02 * time: 50.31488084793091 6 1.242797e+03 1.279629e+02 * time: 50.9742169380188 7 1.236166e+03 9.206209e+01 * time: 51.66398787498474 8 1.232035e+03 5.421101e+01 * time: 52.44599986076355 9 1.229675e+03 2.814654e+01 * time: 53.17966985702515 10 1.228897e+03 2.765795e+01 * time: 53.85512900352478 11 1.228425e+03 2.753520e+01 * time: 54.47063899040222 12 1.227241e+03 3.763993e+01 * time: 55.14527487754822 13 1.224655e+03 6.404633e+01 * time: 55.81128191947937 14 1.219084e+03 9.969910e+01 * time: 56.355148792266846 15 1.211334e+03 1.138297e+02 * time: 56.85761380195618 16 1.204992e+03 7.571347e+01 * time: 57.39960694313049 17 1.200192e+03 7.939604e+01 * time: 57.891814947128296 18 1.197822e+03 7.391002e+01 * time: 58.43382000923157 19 1.193880e+03 6.120913e+01 * time: 58.93926191329956 20 1.188991e+03 6.241377e+01 * time: 59.399611949920654 21 1.184817e+03 4.196986e+01 * time: 59.9028639793396 22 1.182823e+03 4.796317e+01 * time: 60.3833429813385 23 1.182000e+03 3.801926e+01 * time: 60.8898138999939 24 1.181212e+03 2.029120e+01 * time: 61.36339092254639 25 1.180485e+03 1.840010e+01 * time: 61.86311197280884 26 1.179734e+03 2.924133e+01 * time: 62.35326600074768 27 1.178963e+03 4.387725e+01 * time: 62.82503890991211 28 1.177902e+03 4.949849e+01 * time: 63.29212784767151 29 1.176275e+03 4.717010e+01 * time: 63.795673847198486 30 1.174087e+03 3.144696e+01 * time: 64.3750069141388 31 1.172197e+03 4.057837e+01 * time: 64.94508981704712 32 1.171084e+03 4.237868e+01 * time: 65.6025857925415 33 1.170157e+03 3.848694e+01 * time: 66.13629698753357 34 1.168694e+03 3.461523e+01 * time: 66.7253029346466 35 1.165655e+03 4.449704e+01 * time: 67.32448983192444 36 1.160201e+03 4.525941e+01 * time: 67.83985590934753 37 1.154709e+03 2.787357e+01 * time: 68.3156168460846 38 1.150862e+03 1.899749e+01 * time: 68.79125094413757 39 1.149623e+03 1.237148e+01 * time: 69.26336693763733 40 1.149412e+03 1.129991e+01 * time: 69.74825692176819 41 1.149284e+03 1.143292e+01 * time: 70.21669888496399 42 1.149088e+03 1.149464e+01 * time: 70.67870688438416 43 1.148690e+03 1.254978e+01 * time: 71.13784694671631 44 1.147733e+03 1.848974e+01 * time: 71.58920001983643 45 1.145684e+03 2.548410e+01 * time: 72.08018779754639 46 1.142886e+03 2.834112e+01 * time: 72.58373594284058 47 1.141014e+03 3.293170e+01 * time: 73.05118680000305 48 1.139336e+03 3.534965e+01 * time: 73.58630895614624 49 1.138796e+03 3.556386e+01 * time: 74.08494091033936 50 1.138297e+03 3.456505e+01 * time: 74.56319499015808 51 1.137093e+03 3.204779e+01 * time: 75.03948187828064 52 1.134867e+03 2.076462e+01 * time: 75.49522280693054 53 1.132771e+03 1.781419e+01 * time: 75.95442581176758 54 1.131712e+03 1.000700e+01 * time: 76.4064838886261 55 1.131299e+03 7.295697e+00 * time: 76.90129494667053 56 1.131214e+03 7.207607e+00 * time: 77.38638687133789 57 1.131189e+03 7.170746e+00 * time: 77.84558987617493 58 1.131159e+03 7.132354e+00 * time: 78.30468201637268 59 1.131108e+03 7.121606e+00 * time: 78.79278182983398 60 1.130967e+03 7.467928e+00 * time: 79.23766779899597 61 1.130593e+03 1.459301e+01 * time: 79.71564388275146 62 1.129468e+03 2.906114e+01 * time: 80.17896580696106 63 1.125621e+03 6.675937e+01 * time: 80.71817183494568 64 1.121748e+03 1.469625e+02 * time: 81.38748598098755 65 1.117515e+03 8.094780e+01 * time: 81.89289784431458 66 1.090808e+03 5.887141e+01 * time: 82.80204200744629 67 1.090346e+03 5.763009e+01 * time: 83.81357979774475 68 1.082836e+03 2.114203e+01 * time: 84.40618896484375 69 1.081117e+03 1.419263e+01 * time: 85.01924896240234 70 1.080042e+03 7.506114e+00 * time: 85.62552499771118 71 1.079105e+03 6.966854e+00 * time: 86.28319501876831 72 1.078877e+03 6.199536e+00 * time: 86.98006200790405 73 1.078852e+03 5.824028e+00 * time: 87.67493200302124 74 1.078839e+03 5.522008e+00 * time: 88.42860293388367 75 1.078784e+03 4.536881e+00 * time: 89.13476181030273 76 1.078682e+03 3.592113e+00 * time: 89.85337901115417 77 1.078468e+03 5.852922e+00 * time: 90.57775688171387 78 1.078184e+03 6.285604e+00 * time: 91.28725981712341 79 1.077968e+03 3.961982e+00 * time: 91.96744799613953 80 1.077913e+03 3.197818e+00 * time: 92.61047983169556 81 1.077906e+03 2.606854e+00 * time: 93.23140501976013 82 1.077900e+03 2.222495e+00 * time: 93.88897800445557 83 1.077885e+03 2.247912e+00 * time: 94.50889682769775 84 1.077848e+03 2.279042e+00 * time: 95.17003393173218 85 1.077757e+03 3.260791e+00 * time: 95.8018147945404 86 1.077569e+03 4.666974e+00 * time: 96.46223783493042 87 1.077280e+03 5.069568e+00 * time: 97.12425088882446 88 1.077031e+03 3.446383e+00 * time: 97.77593398094177 89 1.076942e+03 1.206842e+00 * time: 98.46494889259338 90 1.076931e+03 1.219126e+00 * time: 99.14570879936218 91 1.076930e+03 1.159528e+00 * time: 99.83733582496643 92 1.076929e+03 1.095575e+00 * time: 100.5026068687439 93 1.076926e+03 1.096162e+00 * time: 101.17374801635742 94 1.076921e+03 1.080977e+00 * time: 101.82563781738281 95 1.076910e+03 1.299666e+00 * time: 102.51522588729858 96 1.076890e+03 1.952417e+00 * time: 103.238933801651 97 1.076858e+03 2.213834e+00 * time: 103.93375086784363 98 1.076832e+03 1.540694e+00 * time: 104.66967797279358 99 1.076821e+03 1.070195e+00 * time: 105.35756182670593 100 1.076819e+03 8.094975e-01 * time: 106.0415108203888 101 1.076819e+03 7.381349e-01 * time: 106.71964883804321 102 1.076818e+03 6.602328e-01 * time: 107.36917400360107 103 1.076817e+03 5.677929e-01 * time: 108.08423185348511 104 1.076813e+03 5.068602e-01 * time: 108.74154090881348 105 1.076806e+03 7.076927e-01 * time: 109.40183591842651 106 1.076790e+03 9.173744e-01 * time: 110.10334086418152 107 1.076769e+03 8.814664e-01 * time: 110.77104997634888 108 1.076753e+03 4.766855e-01 * time: 111.47770190238953 109 1.076749e+03 3.299655e-01 * time: 112.16301894187927 110 1.076748e+03 3.299949e-01 * time: 112.88556098937988 111 1.076748e+03 3.301208e-01 * time: 113.56006193161011 112 1.076748e+03 3.302106e-01 * time: 114.23961687088013 113 1.076748e+03 3.302321e-01 * time: 114.88364386558533 114 1.076748e+03 3.298995e-01 * time: 115.5705578327179 115 1.076746e+03 3.283989e-01 * time: 116.28547286987305 116 1.076743e+03 3.236283e-01 * time: 116.9319839477539 117 1.076736e+03 3.105718e-01 * time: 117.64321398735046 118 1.076719e+03 3.943687e-01 * time: 118.35608696937561 119 1.076682e+03 7.732129e-01 * time: 119.01122689247131 120 1.076620e+03 2.181546e+00 * time: 119.69675493240356 121 1.076543e+03 2.981177e+00 * time: 120.37044501304626 122 1.076500e+03 1.253603e+00 * time: 121.07064080238342 123 1.076492e+03 6.439285e-01 * time: 122.10538387298584 124 1.076491e+03 5.918611e-01 * time: 124.03505182266235 125 1.076490e+03 6.227171e-01 * time: 124.80739498138428 126 1.076489e+03 6.362415e-01 * time: 125.47609996795654 127 1.076489e+03 6.386330e-01 * time: 126.15632796287537 128 1.076489e+03 6.378581e-01 * time: 126.83844780921936 129 1.076488e+03 6.330405e-01 * time: 127.63648581504822 130 1.076487e+03 6.183830e-01 * time: 128.28257179260254 131 1.076485e+03 5.802552e-01 * time: 128.91200280189514 132 1.076478e+03 4.889058e-01 * time: 129.56687593460083 133 1.076467e+03 6.280202e-01 * time: 130.18684601783752 134 1.076452e+03 8.039889e-01 * time: 130.85084700584412 135 1.076440e+03 7.291919e-01 * time: 131.5578269958496 136 1.076436e+03 3.470876e-01 * time: 132.21841192245483 137 1.076436e+03 1.621274e-01 * time: 132.89499592781067 138 1.076435e+03 5.005089e-02 * time: 133.57436084747314 139 1.076434e+03 4.475264e-02 * time: 134.2335169315338 140 1.076433e+03 1.671055e-02 * time: 134.88869380950928 141 1.076433e+03 2.401701e-02 * time: 135.58629298210144 142 1.076433e+03 3.417469e-02 * time: 136.24390482902527 143 1.076433e+03 2.737273e-02 * time: 136.87294697761536 144 1.076433e+03 8.887351e-03 * time: 137.50707602500916 145 1.076433e+03 6.048890e-04 * time: 138.18009781837463
FittedPumasModel
Dynamical system type: Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)
Number of subjects: 32
Observation records: Active Missing
conc: 247 51
pca: 232 66
Total: 479 117
Number of parameters: Constant Optimized
0 16
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -1076.4328
----------------------
Estimate
----------------------
tvcl 0.1356
tvvc 7.9787
tvka 1.1818
tvalag 0.87767
tve0 96.615
tvimax 1.0
tvic50 1.1724
tvturn 18.835
Ωpk₁,₁ 0.069142
Ωpk₂,₂ 0.022038
Ωpk₃,₃ 0.84928
Ωpd₁,₁ 0.0028141
Ωpd₂,₂ 0.18092
σprop_pk 0.086547
σadd_pk 0.44227
σadd_pd 4.1414
----------------------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,
predictandsimobspredictnot technically simulation, but end result is the same, primary difference is the lack ofRUVcompared tosimobs. 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
TipApplication: 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)
# dosing rows must stay without a concentration
@rtransform :conc_new = :evid == 1 ? missing : :conc_new
# update route column for NCA
@transform :route = "ev"
# create nca data, specify obs col as conc_new
read_nca_pumas(_; concentration = :conc_new, dose = :amt, route = :route)
run_nca(_)
end32×48 DataFrame
7 rows omitted
| Row | id | profile | dose | route | duration | tau | assumed_tau | auc_extrap_percent_obs | auc_extrap_percent_pred | aucinf_dn_obs | aucinf_dn_pred | aucinf_obs | aucinf_pred | auclast | auclast_dn | aumc_extrap_percent_obs | aumc_extrap_percent_pred | aumcinf_obs | aumcinf_pred | aumclast | cl_f_obs | cl_f_pred | clast | clast_pred | cmax | cmax_dn | cmin | cmin_dn | half_life | lambdaz | lambdazadjr2 | lambdazintercept | lambdaznpoints | lambdazr2 | lambdaztimefirst | lambdaztimelast | mrtinf_obs | mrtinf_pred | mrtlast | n_samples | span | tfinal | tlag | tlast | tmax | tmin | vz_f_obs | vz_f_pred |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| String | Int64 | Float64 | String | Missing | Missing | Missing | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Int64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 1 | 1 | 100.0 | ev | missing | missing | missing | 1.10635 | 1.10635 | 3.51501 | 3.51501 | 351.501 | 351.501 | 347.612 | 3.47612 | 5.95913 | 5.95913 | 11406.7 | 11406.7 | 10726.9 | 0.284494 | 0.284494 | 0.126288 | 0.126288 | 10.8 | 0.108 | 0.0 | 0.0 | 21.3442 | 0.0324747 | 1.0 | 2.60717 | 3 | 1.0 | 96.0 | 144.0 | 32.4514 | 32.4514 | 30.8589 | 15 | 2.24885 | 144.0 | 0.5 | 144.0 | 9.0 | 0.5 | 8.76049 | 8.76049 |
| 2 | 2 | 2 | 100.0 | ev | missing | missing | missing | 8.18387 | 8.11264 | 7.5934 | 7.58752 | 759.34 | 758.752 | 697.197 | 6.97197 | 27.6183 | 27.4284 | 45091.0 | 44973.0 | 32637.6 | 0.131693 | 0.131795 | 1.1019 | 1.09146 | 11.2193 | 0.112193 | 0.0 | 0.0 | 39.0913 | 0.0177315 | 0.982873 | 2.64085 | 8 | 0.98532 | 12.0 | 144.0 | 59.3818 | 59.2724 | 46.8127 | 16 | 3.37671 | 144.0 | 0.5 | 144.0 | 9.0 | 0.0 | 7.42708 | 7.43284 |
| 3 | 3 | 3 | 120.0 | ev | missing | missing | missing | 15.9427 | 15.8916 | 8.96545 | 8.96 | 1075.85 | 1075.2 | 904.334 | 7.53612 | 44.3983 | 44.3041 | 85095.2 | 84951.2 | 47314.4 | 0.111539 | 0.111607 | 2.24883 | 2.24025 | 14.4 | 0.12 | 0.0 | 0.0 | 52.867 | 0.0131112 | 0.99836 | 2.69459 | 5 | 0.99877 | 48.0 | 144.0 | 79.0955 | 79.0096 | 52.3195 | 16 | 1.81588 | 144.0 | 0.5 | 144.0 | 9.0 | 0.0 | 8.5072 | 8.51238 |
| 4 | 4 | 4 | 60.0 | ev | missing | missing | missing | 10.8066 | 10.9356 | 13.6626 | 13.6823 | 819.754 | 820.941 | 731.166 | 12.1861 | 33.0929 | 33.3884 | 54571.1 | 54813.1 | 36511.9 | 0.0731927 | 0.0730869 | 1.47997 | 1.4998 | 11.9 | 0.198333 | 0.0 | 0.0 | 41.4901 | 0.0167063 | 0.993407 | 2.81104 | 3 | 0.996704 | 96.0 | 144.0 | 66.5701 | 66.7687 | 49.9366 | 16 | 1.1569 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 4.38114 | 4.37481 |
| 5 | 5 | 5 | 113.0 | ev | missing | missing | missing | 14.1411 | 14.1875 | 5.87913 | 5.88231 | 664.342 | 664.702 | 570.397 | 5.04776 | 42.7028 | 42.7962 | 48509.2 | 48588.4 | 27794.4 | 0.170093 | 0.170001 | 1.22807 | 1.23277 | 8.93906 | 0.0791068 | 0.0 | 0.0 | 53.0246 | 0.0130722 | 0.999112 | 2.09166 | 3 | 0.999556 | 96.0 | 144.0 | 73.0183 | 73.098 | 48.7281 | 16 | 0.905241 | 144.0 | 0.5 | 144.0 | 3.0 | 0.0 | 13.0118 | 13.0048 |
| 6 | 6 | 6 | 90.0 | ev | missing | missing | missing | 0.162159 | 0.26912 | 4.34934 | 4.354 | 391.44 | 391.86 | 390.806 | 4.34229 | 0.922718 | 1.5237 | 11557.9 | 11628.4 | 11451.2 | 0.22992 | 0.229674 | 0.0264344 | 0.0439179 | 13.4 | 0.148889 | 0.0 | 0.0 | 16.6441 | 0.0416451 | 0.939694 | 2.87146 | 10 | 0.946395 | 6.0 | 144.0 | 29.5265 | 29.6749 | 29.3016 | 15 | 8.2912 | 144.0 | 0.5 | 144.0 | 3.0 | 0.5 | 5.52094 | 5.51502 |
| 7 | 7 | 7 | 135.0 | ev | missing | missing | missing | 4.53414 | 5.56274 | 5.60647 | 5.66753 | 756.873 | 765.117 | 722.555 | 5.35226 | 17.7366 | 21.0983 | 37374.7 | 38967.1 | 30745.7 | 0.178365 | 0.176444 | 0.698002 | 0.865676 | 17.6 | 0.13037 | 0.0 | 0.0 | 34.079 | 0.0203394 | 0.969354 | 2.78463 | 11 | 0.972419 | 3.0 | 144.0 | 49.3805 | 50.9297 | 42.5514 | 16 | 4.13745 | 144.0 | 0.5 | 144.0 | 2.0 | 0.0 | 8.76945 | 8.67496 |
| 8 | 8 | 8 | 75.0 | ev | missing | missing | missing | 5.32736 | 4.75109 | 9.53126 | 9.47359 | 714.844 | 710.519 | 676.762 | 9.02349 | 20.793 | 18.8774 | 34938.6 | 34113.5 | 27673.8 | 0.104918 | 0.105557 | 0.814328 | 0.721846 | 12.9 | 0.172 | 0.0 | 0.0 | 32.4153 | 0.0213833 | 0.987479 | 2.75326 | 8 | 0.989268 | 12.0 | 144.0 | 48.8758 | 48.0121 | 40.8915 | 16 | 4.07215 | 144.0 | 0.5 | 144.0 | 9.0 | 0.0 | 4.90653 | 4.93639 |
| 9 | 9 | 9 | 105.0 | ev | missing | missing | missing | 14.4123 | 15.2131 | 8.45264 | 8.53247 | 887.527 | 895.909 | 759.614 | 7.23442 | 41.4362 | 42.9843 | 67599.5 | 69435.0 | 39588.9 | 0.118306 | 0.117199 | 1.70594 | 1.81773 | 11.5 | 0.109524 | 0.0 | 0.0 | 51.973 | 0.0133367 | 0.987148 | 2.51807 | 9 | 0.988754 | 9.0 | 144.0 | 76.1661 | 77.5022 | 52.1171 | 16 | 2.5975 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 8.87074 | 8.78774 |
| 10 | 11 | 10 | 123.0 | ev | missing | missing | missing | 9.08161 | 9.95779 | 7.58298 | 7.65677 | 932.706 | 941.782 | 848.002 | 6.89432 | 30.5925 | 32.7954 | 57074.0 | 58944.8 | 39613.6 | 0.131874 | 0.130603 | 1.3633 | 1.50938 | 17.5 | 0.142276 | 0.0 | 0.0 | 43.0665 | 0.0160948 | 0.976213 | 2.72935 | 9 | 0.979186 | 9.0 | 144.0 | 61.1918 | 62.5886 | 46.7141 | 16 | 3.13468 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 8.1936 | 8.11464 |
| 11 | 12 | 11 | 113.0 | ev | missing | missing | missing | 6.28718 | 6.68105 | 7.37587 | 7.407 | 833.473 | 836.991 | 781.071 | 6.91214 | 22.3203 | 23.467 | 46001.5 | 46690.8 | 35733.9 | 0.135577 | 0.135007 | 1.00889 | 1.07662 | 15.0 | 0.132743 | 0.0 | 0.0 | 36.0023 | 0.0192528 | 0.989375 | 2.84623 | 8 | 0.990893 | 12.0 | 144.0 | 55.1926 | 55.7841 | 45.7498 | 16 | 3.66643 | 144.0 | 0.5 | 144.0 | 9.0 | 0.0 | 7.04193 | 7.01234 |
| 12 | 13 | 12 | 113.0 | ev | missing | missing | missing | 7.13796 | 7.40025 | 5.70094 | 5.71709 | 644.206 | 646.031 | 598.223 | 5.29401 | 25.6683 | 26.4178 | 35656.3 | 36019.5 | 26504.0 | 0.17541 | 0.174914 | 0.835491 | 0.868646 | 12.9 | 0.114159 | 0.0 | 0.0 | 38.1489 | 0.0181695 | 0.984177 | 2.47559 | 9 | 0.986155 | 9.0 | 144.0 | 55.3492 | 55.7551 | 44.3045 | 16 | 3.53876 | 144.0 | 0.5 | 144.0 | 3.0 | 0.0 | 9.65407 | 9.6268 |
| 13 | 14 | 13 | 75.0 | ev | missing | missing | missing | 2.42909 | 2.40967 | 6.47357 | 6.47228 | 485.518 | 485.421 | 473.724 | 6.31632 | 11.2144 | 11.1327 | 19193.3 | 19175.7 | 17040.9 | 0.154474 | 0.154505 | 0.30629 | 0.30378 | 11.6 | 0.154667 | 0.0 | 0.0 | 26.6896 | 0.0259707 | 0.995124 | 2.54833 | 11 | 0.995612 | 3.0 | 144.0 | 39.5316 | 39.5031 | 35.9722 | 16 | 5.28296 | 144.0 | 0.5 | 144.0 | 2.0 | 0.0 | 5.94803 | 5.94921 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 21 | 22 | 21 | 87.0 | ev | missing | missing | missing | 8.37608 | 8.35081 | 10.0594 | 10.0567 | 875.171 | 874.929 | 801.866 | 9.21685 | 28.5563 | 28.4891 | 51762.6 | 51714.0 | 36981.1 | 0.0994092 | 0.0994366 | 1.27169 | 1.2675 | 13.6355 | 0.156729 | 0.0 | 0.0 | 39.9556 | 0.0173479 | 0.985457 | 2.73515 | 9 | 0.987275 | 9.0 | 144.0 | 59.1457 | 59.1064 | 46.1188 | 16 | 3.37875 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 5.73031 | 5.73189 |
| 22 | 23 | 22 | 110.0 | ev | missing | missing | missing | 12.144 | 12.1382 | 8.29193 | 8.29139 | 912.112 | 912.052 | 801.345 | 7.28496 | 37.3896 | 37.377 | 62673.7 | 62661.0 | 39240.2 | 0.120599 | 0.120607 | 1.63961 | 1.63872 | 12.5048 | 0.11368 | 0.0 | 0.0 | 46.8269 | 0.0148023 | 0.998688 | 2.62545 | 9 | 0.998852 | 9.0 | 144.0 | 68.7127 | 68.7033 | 48.9679 | 16 | 2.88296 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 8.14731 | 8.14784 |
| 23 | 24 | 23 | 115.0 | ev | missing | missing | missing | 18.3871 | 20.8861 | 9.05144 | 9.33735 | 1040.92 | 1073.79 | 849.521 | 7.38714 | 48.8603 | 52.8204 | 92279.7 | 1.00025e5 | 47191.5 | 0.11048 | 0.107097 | 2.08998 | 2.44902 | 11.6343 | 0.101168 | 0.0 | 0.0 | 63.4764 | 0.0109198 | 0.96369 | 2.46813 | 9 | 0.968229 | 9.0 | 144.0 | 88.6524 | 93.1512 | 55.5508 | 16 | 2.12678 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 10.1174 | 9.80761 |
| 24 | 25 | 24 | 112.0 | ev | missing | missing | missing | 11.9925 | 11.9673 | 9.58169 | 9.57894 | 1073.15 | 1072.84 | 944.452 | 8.4326 | 37.3318 | 37.2758 | 73580.0 | 73514.3 | 46111.3 | 0.104366 | 0.104396 | 1.85346 | 1.84903 | 14.6168 | 0.130507 | 0.0 | 0.0 | 48.1296 | 0.0144017 | 0.999712 | 2.6885 | 3 | 0.999856 | 96.0 | 144.0 | 68.5646 | 68.523 | 48.8233 | 16 | 0.997306 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 7.24678 | 7.24886 |
| 25 | 26 | 25 | 120.0 | ev | missing | missing | missing | 15.042 | 15.0328 | 8.40799 | 8.40708 | 1008.96 | 1008.85 | 857.191 | 7.14326 | 43.5163 | 43.4986 | 76970.1 | 76946.0 | 43475.5 | 0.118935 | 0.118947 | 1.9788 | 1.97738 | 12.5856 | 0.10488 | 0.0 | 0.0 | 53.1621 | 0.0130384 | 0.999968 | 2.5593 | 3 | 0.999984 | 96.0 | 144.0 | 76.2867 | 76.271 | 50.7186 | 16 | 0.902899 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 9.12188 | 9.12287 |
| 26 | 27 | 26 | 120.0 | ev | missing | missing | missing | 8.34752 | 9.18776 | 8.26742 | 8.34391 | 992.09 | 1001.27 | 909.275 | 7.5773 | 28.4085 | 30.594 | 59406.6 | 61277.3 | 42530.1 | 0.120957 | 0.119848 | 1.38519 | 1.53873 | 15.3615 | 0.128012 | 0.0 | 0.0 | 41.4405 | 0.0167263 | 0.991338 | 2.83955 | 9 | 0.992421 | 9.0 | 144.0 | 59.8803 | 61.1995 | 46.7736 | 16 | 3.25769 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 7.23151 | 7.16521 |
| 27 | 28 | 27 | 120.0 | ev | missing | missing | missing | 11.7711 | 11.4144 | 7.39845 | 7.36866 | 887.814 | 884.239 | 783.308 | 6.52757 | 36.5679 | 35.7643 | 59924.6 | 59175.0 | 38011.4 | 0.135163 | 0.13571 | 1.59104 | 1.53661 | 12.3499 | 0.102916 | 0.0 | 0.0 | 45.5286 | 0.0152244 | 0.991505 | 2.6219 | 9 | 0.992567 | 9.0 | 144.0 | 67.4968 | 66.9219 | 48.5268 | 16 | 2.96517 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 8.87805 | 8.91395 |
| 28 | 29 | 28 | 153.0 | ev | missing | missing | missing | 10.1961 | 11.2625 | 5.14084 | 5.20262 | 786.548 | 796.0 | 706.351 | 4.61667 | 32.7487 | 35.2481 | 51275.0 | 53254.2 | 34483.1 | 0.194521 | 0.192211 | 1.22659 | 1.37116 | 11.6258 | 0.0759856 | 0.0 | 0.0 | 45.3196 | 0.0152946 | 0.990591 | 2.51809 | 9 | 0.991767 | 9.0 | 144.0 | 65.1899 | 66.9022 | 48.8187 | 16 | 2.97884 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 12.7182 | 12.5672 |
| 29 | 30 | 29 | 105.0 | ev | missing | missing | missing | 10.1005 | 10.4871 | 8.0251 | 8.05977 | 842.636 | 846.275 | 757.526 | 7.21453 | 32.8292 | 33.7592 | 53713.6 | 54467.7 | 36079.9 | 0.124609 | 0.124073 | 1.34696 | 1.40456 | 12.4301 | 0.118382 | 0.0 | 0.0 | 43.798 | 0.015826 | 0.994098 | 2.61867 | 9 | 0.994835 | 9.0 | 144.0 | 63.7448 | 64.3617 | 47.6286 | 16 | 3.08233 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 7.87368 | 7.83982 |
| 30 | 31 | 30 | 125.0 | ev | missing | missing | missing | 12.4084 | 12.9299 | 7.10822 | 7.1508 | 888.528 | 893.85 | 778.276 | 6.22621 | 37.7782 | 38.8925 | 62344.8 | 63481.7 | 38792.1 | 0.140682 | 0.139845 | 1.58346 | 1.6599 | 12.0321 | 0.0962567 | 0.0 | 0.0 | 48.2617 | 0.0143622 | 0.997127 | 2.57492 | 9 | 0.997486 | 9.0 | 144.0 | 70.1664 | 71.0205 | 49.8436 | 16 | 2.79725 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 9.79527 | 9.73695 |
| 31 | 32 | 31 | 93.0 | ev | missing | missing | missing | 14.8495 | 14.8027 | 9.75247 | 9.74712 | 906.98 | 906.482 | 772.298 | 8.30428 | 42.6737 | 42.5831 | 68728.4 | 68620.0 | 39399.5 | 0.102538 | 0.102594 | 1.82584 | 1.81909 | 11.3276 | 0.121802 | 0.0 | 0.0 | 51.1296 | 0.0135567 | 0.999223 | 2.5505 | 3 | 0.999611 | 96.0 | 144.0 | 75.7772 | 75.6992 | 51.0159 | 16 | 0.938792 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 7.56366 | 7.56781 |
| 32 | 33 | 32 | 100.0 | ev | missing | missing | missing | 13.9197 | 13.7291 | 9.04179 | 9.02181 | 904.179 | 902.181 | 778.32 | 7.7832 | 40.9712 | 40.5847 | 66227.1 | 65796.4 | 39093.1 | 0.110598 | 0.110843 | 1.75804 | 1.73013 | 11.639 | 0.11639 | 0.0 | 0.0 | 49.6229 | 0.0139683 | 0.994638 | 2.55963 | 9 | 0.995308 | 9.0 | 144.0 | 73.2456 | 72.9304 | 50.2276 | 16 | 2.72052 | 144.0 | 0.5 | 144.0 | 6.0 | 0.0 | 7.91776 | 7.93529 |
- 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
SubjectandDosageRegimenconstructors without the additional complexity of creating a virtual population from scratch.
- Same population, different dosage regimen lets us examine the
NoteKey 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
Subjectsyntax by accessing data from the firstSubjectstored inmyfit.data[1].- Converting the mg/kg dose to mg requires
wtblwhich we extract from thecovariatesfield.
- Converting the mg/kg dose to mg requires
# 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 inmyfit.dataand modified them as we did above; this syntax below shows an equivalent approach that may be more intuitive to new users.
- Note, we could have iterated over all
# 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
layoutkwarg inmapping, and then separated using thepaginatefunction. x|yticklabelsizewas 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
Subjectfrom scratch along with thezero_randeffshelper function.
- We use this scenario to showcase creating a
NoteIntroduction 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,
)
endSimulated 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
Populationof 100Subjects 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
pcawithin the therapeutic range (20-35%) at any time during treatment. - The time needed to reach the first therapeutic
pcavalue. - The total time spent in the TR as a percentage of the dosing interval (i.e., 24 hours) at SS (Day 14).
- The probability of obtaining a
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
mapwith aRangebetween1:nand ado-block to create a vector of virtual subjects (i.e., aPopulation) - For simplicity, we will also use scalar literal values for the range of
wtblin theUniformcall instead of obtaining them programmatically withextremaor 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),
)
endPopulation
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:nas we did when creatingpop03. - The resulting vector of
SimulatedObservationobjects can be concatenated into a single object using thereduce(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
Subjectsyntax to add arep_idas a covariate for each subjectpop03inside themapbefore we callsimobs.
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)
end5-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, Nothing, Float64}, tvvc::RealDomain{Float64, Nothing, Float64}, tvka::RealDomain{Float64, Nothing, Float64}, tvalag::RealDomain{Float64, Nothing, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, Nothing, Float64}, tvturn::RealDomain{Float64, Nothing, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, Nothing, Float64}, σadd_pk::RealDomain{Float64, Nothing, Float64}, σadd_pd::RealDomain{Float64, Nothing, Int64}}}, ParamSet{@NamedTuple{ηpk::Pumas.RandomEffect{false, var"#16#17"}, ηpd::Pumas.RandomEffect{false, var"#18#19"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#20#21", var"#22#23"}}, Pumas.DCPObj{var"#28#29"}, var"#30#31", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, false), var"#32#33", var"#34#35"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), (:ic50, :vc, :imax), (:central, :e), (), var"#36#37"}, var"#38#39", ModelingToolkitBase.System, Pumas.PumasModelOptions{Pumas.AllSteadyStateStates}}, Subject{@NamedTuple{}, Pumas.ConstantInterpolationStructArray{Float64, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, 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::Vector{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, Val{true}}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.LUFactorization{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, Val{true}}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, @NamedTuple{depot::Vector{Float64}, central::Vector{Float64}, e::Vector{Float64}}, 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
goodsimsbecause this approach, while valid, is redundant because we are recreating thePopulationfrom scratch with each simulation. - We could simplify the code by using a nested
mapcall 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)
end5-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, Nothing, Float64}, tvvc::RealDomain{Float64, Nothing, Float64}, tvka::RealDomain{Float64, Nothing, Float64}, tvalag::RealDomain{Float64, Nothing, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, Nothing, Float64}, tvturn::RealDomain{Float64, Nothing, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, Nothing, Float64}, σadd_pk::RealDomain{Float64, Nothing, Float64}, σadd_pd::RealDomain{Float64, Nothing, Int64}}}, ParamSet{@NamedTuple{ηpk::Pumas.RandomEffect{false, var"#16#17"}, ηpd::Pumas.RandomEffect{false, var"#18#19"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#20#21", var"#22#23"}}, Pumas.DCPObj{var"#28#29"}, var"#30#31", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, false), var"#32#33", var"#34#35"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), (:ic50, :vc, :imax), (:central, :e), (), var"#36#37"}, var"#38#39", ModelingToolkitBase.System, Pumas.PumasModelOptions{Pumas.AllSteadyStateStates}}, Subject{@NamedTuple{}, Pumas.ConstantCovar{@NamedTuple{wtbl::Float64, rep_id::Int64}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, 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::Vector{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, Val{true}}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.LUFactorization{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, Val{true}}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, @NamedTuple{depot::Vector{Float64}, central::Vector{Float64}, e::Vector{Float64}}, 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, pcaThe
bettersimsyntax 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_warfarinthat accepts a single positional argument,dosethat we can use along withmap.dosewas also added as a covariate in theSubjectconstructor 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_warfarinto 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
endsimulate_warfarin (generic function with 1 method)Lastly, we perform the simulations using a
mapreducecall and save the result in a variable,sim03.mapreduceallows us to combine themapandreduce(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, Nothing, Float64}, tvvc::RealDomain{Float64, Nothing, Float64}, tvka::RealDomain{Float64, Nothing, Float64}, tvalag::RealDomain{Float64, Nothing, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, Nothing, Float64}, tvturn::RealDomain{Float64, Nothing, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, Nothing, Float64}, σadd_pk::RealDomain{Float64, Nothing, Float64}, σadd_pd::RealDomain{Float64, Nothing, Int64}}}, ParamSet{@NamedTuple{ηpk::Pumas.RandomEffect{false, var"#16#17"}, ηpd::Pumas.RandomEffect{false, var"#18#19"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#20#21", var"#22#23"}}, Pumas.DCPObj{var"#28#29"}, var"#30#31", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, false), var"#32#33", var"#34#35"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), (:ic50, :vc, :imax), (:central, :e), (), var"#36#37"}, var"#38#39", ModelingToolkitBase.System, Pumas.PumasModelOptions{Pumas.AllSteadyStateStates}}, Subject{@NamedTuple{}, Pumas.ConstantCovar{@NamedTuple{wtbl::Float64, rep_id::Int64, dose::Int64}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, 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::Vector{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, Val{true}}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.LUFactorization{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, Val{true}}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, @NamedTuple{depot::Vector{Float64}, central::Vector{Float64}, e::Vector{Float64}}, 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, pca4.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×34 DataFrame
21059975 rows omitted
| Row | id | time | conc | pca | evid | lags_depot | amt | cmt | rate | duration | ss | ii | wtbl | rep_id | dose | tad | dosenum | depot | central | e | ηpk₁ | ηpk₂ | ηpk₃ | ηpd₁ | ηpd₂ | cl | vc | ka | e0 | imax | ic50 | turn | kout | kin |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| String | Float64 | Float64? | Float64? | Int64? | Float64? | Float64? | Symbol? | Float64? | Float64? | Int8? | Float64? | 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.87767 | 5.0 | depot | 0.0 | 0.0 | 0 | 0.0 | 87.3753 | 1 | 5 | 0.0 | 1 | 0.0 | 0.0 | 95.9125 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 2 | 1 | 0.0 | 0.0 | 95.9125 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 0.0 | 1 | 0.0 | 0.0 | 95.9125 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 3 | 1 | 1.0 | 0.0155489 | 95.9074 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 1.0 | 1 | 4.84299 | 0.156899 | 95.9074 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 4 | 1 | 2.0 | 0.124914 | 95.5525 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 2.0 | 1 | 3.73117 | 1.26047 | 95.5525 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 5 | 1 | 3.0 | 0.207909 | 94.816 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 3.0 | 1 | 2.87459 | 2.09794 | 94.816 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 6 | 1 | 4.0 | 0.270601 | 93.8584 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 4.0 | 1 | 2.21466 | 2.73055 | 93.8584 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 7 | 1 | 5.0 | 0.317666 | 92.7754 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 5.0 | 1 | 1.70623 | 3.20547 | 92.7754 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 8 | 1 | 6.0 | 0.352704 | 91.6276 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 6.0 | 1 | 1.31452 | 3.55903 | 91.6276 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 9 | 1 | 7.0 | 0.378491 | 90.4544 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 7.0 | 1 | 1.01274 | 3.81924 | 90.4544 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 10 | 1 | 8.0 | 0.397164 | 89.2826 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 8.0 | 1 | 0.780244 | 4.00766 | 89.2826 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 11 | 1 | 9.0 | 0.41037 | 88.13 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 9.0 | 1 | 0.60112 | 4.14091 | 88.13 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 12 | 1 | 10.0 | 0.419376 | 87.0088 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 10.0 | 1 | 0.463119 | 4.23179 | 87.0088 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| 13 | 1 | 11.0 | 0.42516 | 85.9271 | 0 | missing | missing | missing | missing | missing | missing | missing | 87.3753 | 1 | 5 | 11.0 | 1 | 0.356799 | 4.29016 | 85.9271 | -0.344186 | 0.0131194 | -1.51102 | -0.00729419 | -0.224819 | 0.113499 | 10.0907 | 0.260812 | 95.9125 | 1.0 | 0.936322 | 18.8353 | 0.0530917 | 5.09216 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 21059989 | 100 | 325.0 | 3.99579 | 21.3603 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 13.0 | 14 | 1.01988e-5 | 36.7272 | 21.3603 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059990 | 100 | 326.0 | 3.92821 | 21.3722 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 14.0 | 14 | 3.16063e-6 | 36.106 | 21.3722 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059991 | 100 | 327.0 | 3.86178 | 21.3985 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 15.0 | 14 | 9.79486e-7 | 35.4954 | 21.3985 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059992 | 100 | 328.0 | 3.79646 | 21.4385 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 16.0 | 14 | 3.03545e-7 | 34.8951 | 21.4385 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059993 | 100 | 329.0 | 3.73225 | 21.4918 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 17.0 | 14 | 9.40693e-8 | 34.3049 | 21.4918 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059994 | 100 | 330.0 | 3.66913 | 21.5577 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 18.0 | 14 | 2.91523e-8 | 33.7247 | 21.5577 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059995 | 100 | 331.0 | 3.60708 | 21.6358 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 19.0 | 14 | 9.03435e-9 | 33.1543 | 21.6358 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059996 | 100 | 332.0 | 3.54607 | 21.7255 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 20.0 | 14 | 2.79977e-9 | 32.5936 | 21.7255 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059997 | 100 | 333.0 | 3.48609 | 21.8264 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 21.0 | 14 | 8.67654e-10 | 32.0423 | 21.8264 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059998 | 100 | 334.0 | 3.42713 | 21.9381 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 22.0 | 14 | 2.68887e-10 | 31.5004 | 21.9381 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21059999 | 100 | 335.0 | 3.36917 | 22.0602 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 23.0 | 14 | 8.33264e-11 | 30.9676 | 22.0602 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
| 21060000 | 100 | 336.0 | 3.31219 | 22.1922 | 0 | missing | missing | missing | missing | missing | missing | missing | 71.8961 | 200 | 15 | 24.0 | 14 | 2.58229e-11 | 30.4439 | 22.1922 | 0.125136 | 0.114773 | -0.00878139 | 0.0299923 | -0.066406 | 0.156784 | 9.19147 | 1.1715 | 99.5563 | 1.0 | 1.09704 | 18.8353 | 0.0530917 | 5.28561 |
NoteSanity 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_metricsfunction is a custom analysis fortable_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 | [21, 39] | [70, 85] | [89, 97] |
| Time to Target | |||
| Median | 122 | 83.6 | 63.1 |
| 95% CI | [104, 144] | [75.1, 95] | [57.6, 69.9] |
| Time in TR | |||
| Median | 95.6 | 96.9 | 97.2 |
| 95% CI | [87.8, 100] | [92.2, 99.6] | [91.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
pcapercentiles (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.