Pumas Docstrings

DataFrame(
  obs::SimulatedPopulation;
  include_events::Bool = true,
  include_covariates::Bool = true,
  include_ebes::Bool = true,
  include_icoefs::Bool = true,
)

Create a DataFrame with the information in a simobs result obs. The data follows Pumas Data Format ready to be imported with read_pumas.

DataFrame(
  pop::Population;
  include_covariates = true,
  include_observations = true,
  include_events = true,
)

Create a DataFrame with the Population information. The data follows Pumas Data Format ready to be imported with read_pumas.

Examples:

julia> events = DosageRegimen(1.0; time = 1.0)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     1.0      1      1.0     1      0.0      0      0.0       0.0     0

julia> pop = [Subject(; id = 1, events), Subject(; id = 2, events)]
Population
  Subjects: 2
  Observations:

julia> DataFrame(pop)
2×11 DataFrame
 Row │ id      time     amt       evid  cmt    rate      duration  ss     ii        tad       dosenum
     │ String  Float64  Float64?  Int8  Int64  Float64?  Float64?  Int8?  Float64?  Float64?  Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 1           1.0       1.0     1      1       0.0       0.0      0       0.0       0.0        1
   2 │ 2           1.0       1.0     1      1       0.0       0.0      0       0.0       0.0        1

DataFrame(dr::DosageRegimen)

Create a DataFrame with the information in the dosage regimen. The data follows Pumas Data Format ready to be imported with read_pumas.

Examples:

julia> DataFrame(DosageRegimen(100))
1×9 DataFrame
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1      0.0      0      0.0       0.0     0

julia> DataFrame(DosageRegimen(100; ii = 24, addl = 2))
1×9 DataFrame
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      2      0.0       0.0     0

DataFrame(subject::AbstractSubject)

Create a DataFrame with the Subject information. The data follows Pumas Data Format ready to be imported with read_pumas.

Examples:

julia> DataFrame(Subject(; events = DosageRegimen(1.0; time = 1.0)))
1×11 DataFrame
 Row │ id      time     amt       evid  cmt    rate      duration  ss     ii        tad       dosenum
     │ String  Float64  Float64?  Int8  Int64  Float64?  Float64?  Int8?  Float64?  Float64?  Int64
─────┼────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 1           1.0       1.0     1      1       0.0       0.0      0       0.0       0.0        1

julia> DataFrame(
         Subject(;
           id = "AKJ491",
           events = DosageRegimen(1.0; time = 1.0),
           observations = (; glucose = [8.31, 7.709, 5.19]),
           time = [0.3, 0.9, 3.0],
           covariates = (; bloodpressure = [(143, 95.01), (141.3, 94.8), (130.4, 85.1)]),
           covariates_time = [0.0, 0.5, 5.0],
         ),
       )
7×13 DataFrame
 Row │ id      time     evid  glucose      amt        cmt      rate       duration   ss       ii         bloodpressure ⋯
     │ String  Float64  Int8  Float64?     Float64?   Int64?   Float64?   Float64?   Int8?    Float64?   Tuple…?       ⋯
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ AKJ491      0.0     2  missing      missing    missing  missing    missing    missing  missing    (143, 95.01)  ⋯
   2 │ AKJ491      0.3     0        8.31   missing    missing  missing    missing    missing  missing    (143, 95.01)
   3 │ AKJ491      0.5     2  missing      missing    missing  missing    missing    missing  missing    (141.3, 94.8)
   4 │ AKJ491      0.9     0        7.709  missing    missing  missing    missing    missing  missing    (141.3, 94.8)
   5 │ AKJ491      1.0     1  missing            1.0        1        0.0        0.0        0        0.0  (141.3, 94.8) ⋯
   6 │ AKJ491      3.0     0        5.19   missing    missing  missing    missing    missing  missing    (141.3, 94.8)
   7 │ AKJ491      5.0     2  missing      missing    missing  missing    missing    missing  missing    (130.4, 85.1)
                                                                                                       2 columns omitted

DataFrame(
  fpmi::FittedPumasModelInspection;
  include_covariates = true,
  include_observations = true,
  include_events = true,
  include_icoef = true,
  include_dosecontrol = true,
  include_ebes = true,
  include_dynamics = true,
)

Create a DataFrame with the information in an inspect result fpmi. The data follows Pumas Data Format ready to be imported with read_pumas.

DataFrame(
  obs::SimulatedObservations;
  include_events::Bool = true,
  include_covariates::Bool = true,
  include_ebes::Bool = true,
  include_icoefs::Bool = true,
)

Create a DataFrame with the information in a simobs result obs. The data follows Pumas Data Format ready to be imported with read_pumas.

Vector-valued simulated variables must provide exactly one value for each observation time point. Other variables, e.g. scalars, are repeated for every time point.

DataFrame(pred::Tuple{Vararg{SubjectPrediction}})

Create a DataFrame with the information in the subject's per-component predictions pred of a mixture model, with the component in a mixture_component column. The data follows Pumas Data Format ready to be imported with read_pumas.

DataFrame(vpred::Vector{<:SubjectPrediction})

Create a DataFrame with the information in a population prediction vpred. The data follows Pumas Data Format ready to be imported with read_pumas.

DataFrame(vresid::Vector{<:SubjectResidual})

Create a DataFrame with the information in a population weighted residuals vresid. The data follows Pumas Data Format ready to be imported with read_pumas.

DataFrame(vpred::Vector{<:Tuple{Vararg{SubjectPrediction}}})

Create a DataFrame with the information in a population's per-component predictions vpred of a mixture model, with the component in a mixture_component column. The data follows Pumas Data Format ready to be imported with read_pumas.

DataFrame(pred::SubjectPrediction)

Create a DataFrame with the information in the subject's prediction pred. The data follows Pumas Data Format ready to be imported with read_pumas.

MvNormal(pmi::FittedPumasModelInference)

Returns an MvNormal with a mean equal to the optimal parameter values in the fitted model and a covariance equal to the variance-covariance estimate.

MvNormal(model::PumasModel, means::NamedTuple, vcov::NamedTuple)

Returns an MvNormal with mean equal to the vectorized means argument and a covariance equal to the block diagonal matrix of the covariance matrices in the vcov argument. The mean argument should be a named tuple of values in their natural domain. vcov should be a named tuple of numbers (for real parameters) or matrices for vector and matrix parameters. For parameters that are themselves symmetric matrices, e.g:

Ω = [
  Ω₁₁ Ω₁₂ Ω₁₃
  Ω₂₁ Ω₂₂ Ω₂₃
  Ω₃₁ Ω₃₂ Ω₃₃
]

the covariance matrix should be the covariance of the vector: [Ω₁₁, Ω₁₂, Ω₂₂, Ω₁₃, Ω₂₃, Ω₃₃]. Note that because the matrix is symmetric, only the upper triangle is used iterating over the columns first in an outer loop then the over the rows of each column in an inner loop.

Example:

means = (θ1 = 1.0, θ2 = [1.0, 2.0], θ3 = [1.0 0.1; 0.1 1.0])
vcov =
  (θ1 = 0.1, θ2 = [0.5 -0.2; -0.2 0.6], θ3 = [0.5 -0.2 0.1; -0.2 0.6 -0.3; 0.1 -0.3 0.45])
dist = MvNormal(model, means, vcov)

Chains(b::BayesMCMCResults; subject::Union{Nothing, Int} = nothing)

Constructs an MCMCChains.Chains object for the samples of the population or subject-specific parameters. If subject is nothing (the default value) a Chains object of the population parameters' samples is returned. If subject is an integer, the samples of the subject-specific parameters corresponding to the subject with the given index will be returned instead.

AnalyticalPKProblem(pkprob, prob2, sys2)

A problem that is partially an analytical problem that can be evaluated independently of the rest of the system.

Fields:

• pkprob: the analytical part of the problem
• prob2: a problem that represents the rest of the system
• sys2 (optional): the system of prob2 - used for latexification.

BayesMCMC

An instance of the Hamiltonian Monte Carlo (HMC) based No-U-Turn Sampler (NUTS) algorithm. All the options of the algorithm can be passed to the constructor of BayesMCMC as keyword arguments. Example:

alg = BayesMCMC(; nsamples = 2000, nadapts = 1000)

The main options that can be set are:

• target_accept (default: 0.8): Target acceptance ratio for the NUTS algorithm.
• nsamples (default: 2000): Number of Markov Chain Monte Carlo (MCMC) samples to generate, including adaptation samples.
• nadapts (default: nsamples ÷ 2): Number of adaptation steps in the NUTS algorithm.
• nchains (default: 4): Number of MCMC chains to sample.
• ess_per_chain (default: ∞): Target effective sample size (ESS) per chain, sampling terminates if the target is reached.
• check_every (default: max(2, ess_per_chain ÷ 5)): Number of samples after which the ESS per chain is checked; only relevant when ess_per_chain is finite.
• time_limit (default: ∞): A time limit for sampling in seconds, sampling terminates if the time limit is reached.
• check_initial (default: true): A Bool to check if the model's initial parameter values have a valid log-likelihood before begin fitting.

To customize the differential equation solver, the parallelism, the random number generator, or to fix some parameters during sampling, use the diffeq_options, ensemblealg, parallel_chains, parallel_subjects, rng, and constantcoef keyword arguments of fit.

Additionally, there are options to fine-tune the NUTS sampler behavior:

• adapt_init_buffer (default: 0): Width of initial fast adaptation interval for the NUTS mass matrix adaptation.
• adapt_term_buffer (default: max(50, nadapts ÷ 20)): Width of final fast adaptation interval for the NUTS mass matrix adaptation.
• adapt_window_size (default: 5): Initial width of slow adaptation interval for the NUTS mass matrix adaptation.

For the NUTS mass matrix adaptation, if you want the same behavior as Stan, set the following values:

BayesMCMC(; ..., adapt_init_buffer = 75, adapt_term_buffer = 50, adapt_window_size = 25)

Bootstrap(;
    rng::Random.AbstractRNG = Random.default_rng(),
    samples::Integer = 200,
    stratify_by::Union{Nothing,Vector{Symbol}} = nothing,
    ensemblealg::SciMLBase.EnsembleAlgorithm = EnsembleThreads(),
)

Construct a Bootstrap object that can be supplied to infer to perform bootstrap-based inference with the provided settings.

See ?infer and the Inference documentation for more information about obtaining inference results.

Arguments

• rng: Pseudo-random number generator used for resampling the population.
• samples: Number of resampled populations.
• stratify_by: Vector of names of covariate(s) used for stratification in the resampling scheme. By default (stratify_by = nothing), no stratification is applied. The vector must contain at least one covariate; pass nothing to disable stratification.
• ensemblealg: Mode of parallelization across fit calls of the resampled populations. Possible options are: EnsembleSerial(), EnsembleThreads() and EnsembleDistributed(). Note that each resampled population is fit using the same mode of parallelization as the fit of the original population.

ByObservation(; allsubjects::Bool = true)

Using this method, each observation or collection of observations is treated as a single data point. When computing predictive loglikelihoods using this method, the predictive loglikelihood computed is the conditional loglikelihood of one or more observations for one or more subjects.

This method has one option, allsubjects. If allsubjects is set to true (the default value), the ith observation of each subject are all grouped together into a single data point. If allsubjects is set to false, then each observation per subject is its individual data point.

Examples:

Assume there are 2 subjects and 3 observations per subject. When using LeaveK(K = 1) as the splitting method together with ByObservation(allsubjects = false), the training and validation splittings are:

On the other hand if allsubjects is set to true, the training and validation splittings are:

BySubject(; marginal = LaplaceI())

Using this method, each subject is treated as a single data point. The predictive loglikelihood computed for each subject can be either the marginal loglikelihood or conditional loglikelihood.

This method has one option, marginal. If marginal is set to nothing, the predictive loglikelihood computed for each subject is the conditional loglikelihood using the posterior samples for the subject-specific parameters. Otherwise, the predictive loglikelihood computed for each subject is the marginal loglikelihood using marginal as the marginalization algorithm.

The default value of marginal is LaplaceI() which uses the Laplace method to integrate out the subject-specific parameters. Other options include: FOCE() and FO().

Central1

A one compartment model with dosing into Central. Equivalent to

Central' = -CL/Vc*Central

where clearance, CL, and volume, Vc, are required to be defined in the model through @param, @covariates, or @pre.

Central1Meta1

A model with a central compartment, Central, and a metabolite, Metabolite. Equivalent to

Central'    = -(CL+CLfm)/Vc*Central
Metabolite' =       CLfm/Vc*Central  - CLm/Vm*Metabolite

where clearances (CL and CLm) and volumes (Vc and Vm) and formation clearance of metabolite CLfm are required to be defined in the model through @param, @covariates, or @pre.

Central1Periph1

A two-compartment model with a central compartment, Central and a peripheral compartment, Peripheral. Equivalent to

Central'    = -(CL+Q)/Vc*Central + Q/Vp*Peripheral
Peripheral' =       Q/Vc*Central - Q/Vp*Peripheral

where clearance, CL, and volumes, Vc and Vp, and distribution clearance, Q, are required to be defined in the model through @param, @covariates, or @pre.

Central1Periph1Meta1

A two compartment model with a central compartment, Central, with a peripheral compartment, Peripheral, and a metabolite compartment, Metabolite. Equivalent to

Central'     = -(CL+Q+CLfm)/Vc*Central + Q/Vp*CPeripheral
CPeripheral' =            Q/Vc*Central - Q/Vp*CPeripheral
Metabolite'  = -CLm/Vm*Metabolite + CLfm/Vc*Central

where clearances (CL and CLm) and volumes (Vc, Vp and Vm), distribution clearance (Q), and formation clearance of metabolite CLfm are required to be defined in the model through @param, @covariates, or @pre.

Central1Periph1Meta1Periph1

A two compartment model with a central compartment, Central, with a peripheral compartment, Peripheral, and a metabolite compartment, Metabolite, with a peripheral compartment, MPeripheral. Equivalent to

Central'     = -(CL+Q+CLfm)/Vc*Central + Q/Vp*CPeripheral
CPeripheral' =            Q/Vc*Central - Q/Vp*CPeripheral
Metabolite'  = -(CLm+Qm)/Vm*Metabolite + Qm/Vmp*MPeripheral + CLfm/Vc*Central
MPeripheral' =        Qm/Vm*Metabolite - Qm/Vmp*MPeripheral

where clearances (CL and CLm) and volumes (Vc, Vp, Vm and Vmp), distribution clearances (Q and Qm) and formation clearance of metabolite CLfm are required to be defined in the model through @param, @covariates, or @pre.

Central1Periph2

A three-compartment model with a central compartment, Central and two peripheral compartments, Peripheral1 and Peripheral2. Equivalent to

Central'    = -(CL+Q1+Q2)/Vc*Central + Q1/Vp1*Peripheral1 + Q2/Vp2*Peripheral2
Peripheral1' =     Q1    /Vc*Central - Q1/Vp1*Peripheral1
Peripheral2' =        Q2 /Vc*Central                      - Q2/Vp2*Peripheral2

where clearance, CL, and volumes, Vc, Vp1 and Vp2, and distribution clearances, Q1 and Q2, are required to be defined in the model through @param, @covariates, or @pre.

CombinedCorrelatedNormal(μ::Real, σadd::Real, σprop::Real, ρ::Real)

Construct a normal distribution with mean μ and variance σadd^2 + σprop^2 * μ^2 + 2 * ρ * σadd * σprop * abs(μ) defined in terms of additive standard deviation σadd ≥ 0, proportional standard deviation σprop ≥ 0, and correlation ρ ∈ [-1, 1].

See also ProportionalNormal, CombinedNormal.

CombinedNormal(μ::Real, σadd::Real, σprop::Real)

Construct a normal distribution with mean μ and variance σadd^2 + σprop^2 * μ^2 defined in terms of additive standard deviation σadd ≥ 0 and proportional standard deviation σprop ≥ 0.

See also ProportionalNormal, CombinedCorrelatedNormal.

Constrained

Constrain a Distribution within a RealDomain or a VectorDomain. The most common case is an MvNormal constrained within a VectorDomain. The only supported method for Constrained is logpdf. Notice that the result does not represent a probability distribution since the remaining probability mass is not scaled by the mass excluded by the constraints. For that reason, it is only safe to use Constrained for priors in the @param block.

Example

julia> d = Constrained(MvNormal(fill(1.0, 1, 1)), lower = -1, upper = 1)
Constrained(
  ZeroMeanFullNormal(
  dim: 1
  μ: Zeros(1)
  Σ: [1.0;;]
  );
  lower = [-1],
  upper = [1],
)

julia> logpdf(d, [0])
-0.9189385332046728

julia> logpdf(d, [-2])
-Inf

@param x ∈ CorrDomain(n::Int)
@param x ∈ CorrDomain(; init::AbstractMatrix{<:Real})

Specify a parameter as an n-by-n correlation matrix by providing either the number n of rows and columns or the initial correlation matrix init.

If the number n of rows and columns in the matrix is specified, the initial correlation matrix is set to the identity matrix.

Depots1Central1

A one compartment model with a central compartment, Central, and a depot, Depot. Equivalent to

Depot'   = -Ka*Depot
Central' =  Ka*Depot - CL/Vc*Central

where absoption rate, Ka, clearance, CL, and volume, Vc, are required to be defined in the model through @param, @covariates, or @pre.

Depots1Central1Periph1

A two-compartment model with a central compartment, Central, a peripheral compartment, Peripheral, and a depot Depot. Equivalent to

Depot'      = -Ka*Depot
Central'    =  Ka*Depot - (CL+Q)/Vc*Central + Q/Vp*Peripheral
Peripheral' =                 Q/Vc*Central - Q/Vp*Peripheral

where absorption rate, Ka, clearance, CL, and volumes, Vc and Vp, and distribution clearance, Q, are required to be defined in the model through @param, @covariates, or @pre.

Depots2Central1

A one compartment model with a central compartment, Central, and two depots, Depot1 and Depot2. Equivalent to

Depot1'  = -Ka1*Depot1
Depot2'  = -Ka2*Depot2
Central' =  Ka1*Depot1 + Ka2*Depot2 - CL/Vc*Central

where absorption rates, Ka1 and Ka2, clearance, CL, and volume, Vc, are required to be defined in the model through @param, @covariates, or @pre.

When using this model during simulation or estimation, it is preferred to have 2 dosing rows for each subject in the dataset, where the first dose goes into cmt = :Depot1 (or cmt = 1) and the second dose goes into cmt = :Depot2 (or cmt = 2). Central compartment is cmt = :Central or (cmt = 3), e.g.,

ev = DosageRegimen([100, 100], cmt = [:Depot1, :Depot2])
s1 = Subject(id = 1, events = ev)

Depots2Central1Periph1

A two-compartment model with a central compartment, Central, a peripheral compartment, Peripheral, and two depots Depot1 and Depot2. Equivalent to

Depot1'      = -Ka1*Depot1
Depot2'      = -Ka2*Depot2
Central'    =  Ka1*Depot1 + Ka2*Depot2 - (CL+Q)/Vc*Central + Q/Vp*Peripheral
Peripheral' =                 Q/Vc*Central - Q/Vp*Peripheral

where absorption rates, Ka1 and Ka2, clearance, CL, and volumes, Vc and Vp, and distribution clearance, Q, are required to be defined in the model through @param, @covariates, or @pre.

DosageRegimen(
  amt::Numeric;
  time::Numeric = 0,
  cmt::Union{Numeric,Symbol} = 1,
  evid::Numeric = 1,
  ii::Numeric = zero.(time),
  addl::Numeric = 0,
  rate::Numeric = zero.(amt)./oneunit.(time),
  duration::Numeric = zero(amt)./oneunit.(time),
  ss::Numeric = 0,
)

Lazy representation of a series of events.

Examples

julia> DosageRegimen(100; ii = 24, addl = 6)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0

julia> DosageRegimen(50; ii = 12, addl = 13)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0

julia> DosageRegimen(200; ii = 24, addl = 2)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    200.0     1     24.0      2      0.0       0.0     0

You can also create a new DosageRegimen from various existing DosageRegimens:

evs = DosageRegimen(
  regimen1::DosageRegimen,
  regimen2::DosageRegimen;
  offset = nothing
)

offset specifies if regimen2 should start after an offset following the end of the last event in regimen1.

Examples

julia> e1 = DosageRegimen(100; ii = 24, addl = 6)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0

julia> e2 = DosageRegimen(50; ii = 12, addl = 13)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0

julia> evs = DosageRegimen(e1, e2)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0
   2 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0

julia> DosageRegimen(e1, e2; offset = 10)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8
─────┼────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0
   2 │   178.0      1     50.0     1     12.0     13      0.0       0.0     0

ExactCrossvalidation(; split_method, split_by)

Defines an instance of the ExactCrossvalidation algorithm for crossvalidation. In this algorithm, the fitting is re-run while leaving out a subset of the data each time.

The parallelism used to fit the splits is controlled by the ensemblealg keyword argument of crossvalidate.

The way by which the data is split between training and validation sets is determined using the keyword arguments split_method and split_by. The split_method argument can be of any of the following types:

• LeaveK for leave-K-out cross-validation
• KFold for K-fold cross-validation
• LeaveFutureK for leaving K future points at a time

while split_by can be of any of the following types:

• BySubject for leaving out subjects
• ByObservation for leaving out individual observations per subject

Please refer to the documentation of each of the above types for the parameters of each split method, e.g. typing ?LeaveK or ?BySubject.

Examples:

Assume there are 5 subjects and 10 observations per subject and that ft is the result of the fit function. The following are some of the combinations in which the above inputs can be used:

• Leave-1-observation-out cross-validation, leaving 1 observation for all the subjects at a time. allsubjects = true means that the same observation index is removed for all the subjects, e.g. the 10th observation for all the subjects is used for validation in the first run, then the 9th observation is used for validation in the second run, etc.

      split_method = LeaveK(K = 1)
      split_by = ByObservation(allsubjects = true)
      cv_method = ExactCrossvalidation(; split_method = split_method, split_by = split_by)
      r = crossvalidate(ft, cv_method; ensemblealg = EnsembleThreads())

• Leave-2-future-observations-out cross-validation, leaving 2 future observations per subject such that no less than 4 observations are used in training. allsubjects = false means that only the data of one subject at a time will be split. So in the first run, only observations 1 to 4 of subject 1 are used for training and observations 5 and 6 of subject 1 are used for validation. In the second run, observations 1 to 6 of subject 1 are used for training and observations 7 and 8 are used for validation. In the third run, observations 1 to 8 of subject 1 get used for training and observations 9 and 10 are used for validation. For all 3 runs, all the observations of subjects 2 to 5 are used in training. Then in the fourth to sixth runs, subject 2's data gets split. In the seventh to ninth runs, subject 3's data gets split, etc.

      split_method = LeaveFutureK(K = 2, minimum = 4)
      split_by = ByObservation(allsubjects = false)
      cv_method = ExactCrossvalidation(; split_method = split_method, split_by = split_by)
      r = crossvalidate(ft, cv_method; ensemblealg = EnsembleThreads())

• Leave-1-subject-out cross-validation, marginalizing the subject-specific parameters using FOCE when computing predictive likelihoods. Using this method, subjects {1, 2, 3, 4} are used for training in the first run and subject {5} is used for validation. In the second run, subjects {1, 2, 3, 5} are used for training and subject {4} is used for validation, etc. The predictive likelihoods computed are the marginal likelihood of the subject approximated using FOCE. FO and LaplaceI can also be used instead.

      split_method = LeaveK(K = 1)
      split_by = BySubject(marginal = FOCE())
      cv_method = ExactCrossvalidation(; split_method = split_method, split_by = split_by)
      r = crossvalidate(ft, cv_method; ensemblealg = EnsembleThreads())

• Leave-1-future-subject-out cross-validation, where the predictive likelihoods are conditional likelihoods computed using 0s for the subject-specific parameters instead of marginalizing. The minimum number of subjects is 3 so there are only 2 runs. In the first run, the training subjects are {1, 2, 3} and the validation subject is {4}. In the second run, the training subjects are {1, 2, 3, 4} and the validation subject is {5}.

      split_method = LeaveFutureK(K = 1, minimum = 3)
      split_by = BySubject(marginal = nothing)
      cv_method = ExactCrossvalidation(; split_method = split_method, split_by = split_by)
      r = crossvalidate(ft, cv_method; ensemblealg = EnsembleThreads())

FO()

A first-order approximation method of the marginal likelihood.

The method linearizes the model around the central value of the random effects and uses this linearization to approximate the marginal likelihood integral analytically.

Since FO is not based on and hence does not compute empirical Bayes estimates, it is faster but typically less accurate than FOCE or LaplaceI.

See also FOCE, LaplaceI, SAEM and NaivePooled.

FOCE(; kwargs...)

First-order conditional estimation of the marginal likelihood.

The method approximates the marginal likelihood based on a linearization around the conditional mode (empirical Bayes estimate) of each subject's random effects.

Keyword Arguments

• optim_alg (default: Optim.NewtonTrustRegion(use_fg = false)): optimization algorithm used for the optimization of the empirical Bayes estimates.
• optim_options::NamedTuple (default: (;)): optimization options used for the optimization of the empirical Bayes estimates, with names and values that match the Optim.Options type of the optimization algorithm optim_alg. By default, options show_trace = false (optimization trace is not printed), allow_f_increases = false (optimization is aborted if the objective function value increases) and g_tol = 1e-5 (optimization is stopped if the norm of the gradient is smaller than 1e-5) are used, unless specified otherwise in optim_options.
• randeffs_initialization::RandeffsInitialization.T (default: RandeffsInitialization.PreviousOrth): strategy for initializing random effects at each iteration when fitting a model. See RandeffsInitialization for available options.

See also FO, LaplaceI, SAEM and NaivePooled.

JointMAP

An instance of the maximum a-posteriori (MAP) algorithm where the mode of the joint posterior of the population and subject-specific parameters is found using the L-BFGS optimization algorithm.

alg = JointMAP()

The differential equation solver options (diffeq_options), the parallelism (ensemblealg), the optimization algorithm and its settings (optim_alg, optim_options, optimize_fn), and the fixed parameters (constantcoef) can be customized by passing them as keyword arguments to fit.

KFold(; K = 5, shuffle = false, rng = nothing)

K-fold splitting method. In this method, the data is split K times into 2 disjoint groups, each time starting from the full data set. The 2 groups are typically called training and validation subsets, where the validation subset has floor(N / K) data points, N being the total number of data points. In the next iteration, the whole data set is then re-split using another disjoint validation subset of floor(N / K) different points, disjoint from the previous validation subsets. This process is done iteratively until almost each data point has shown up in 1 and only 1 validation subset. If N is divisible by K, each point will show up in 1 and only 1 validation subset. Otherwise, the remainder points will be part of the training subset for all the splittings and will not show up in any validation subset. This method can be used for computing predictive loglikelihoods, in-sample loglikelihoods and doing cross-validation using subsets of out-of-sample data.

The data is typically a vector of some sort, e.g. observations or subjects, and the splittings are order-dependent. Before performing the splittings, you can randomly shuffle the data vector by setting the shuffle keyword argument to true (default is false) getting rid of the sensitivity to the original order of the data. You can additionally pass an optional pseudorandom number generator rng to control the pseudo-randomness for maximum reproducibility.

Example:

Assume the original data is ["A", "B", "C", "D"]. 4-fold splitting without shuffling results in the following splittings of the data:

where each data point showed once and only once in a validation subset.

2-fold splitting without shuffling results in the following splittings:

LaplaceI(; kwargs...)

Laplace approximation of the marginal likelihood.

The method approximates the marginal likelihood based on a second-order approximation around the conditional mode (empirical Bayes estimate) of each subject's random effects.

Keyword Arguments

• optim_alg (default: Optim.NewtonTrustRegion(use_fg = false)): optimization algorithm used for the optimization of the empirical Bayes estimates.
• optim_options::NamedTuple (default: (;)): optimization options used for the optimization of the empirical Bayes estimates, with names and values that match the Optim.Options type of the optimization algorithm optim_alg. By default, options show_trace = false (optimization trace is not printed), allow_f_increases = false (optimization is aborted if the objective function value increases) and g_tol = 1e-5 (optimization is stopped if the norm of the gradient is smaller than 1e-5) are used, unless specified otherwise in optim_options.
• randeffs_initialization::RandeffsInitialization.T (default: RandeffsInitialization.PreviousOrth): strategy for initializing random effects at each iteration when fitting a model. See RandeffsInitialization for available options.

See also FO, FOCE, SAEM and NaivePooled.

LeaveFutureK(; K = 1, minimum = 2)

Leave-future-K-out splitting method. In this method, the data is assumed to be a timeseries. The goal is to split the data into "past" and "future". Using this method, the data is split multiple times into 3 disjoint groups where the third group is discarded, each time starting from the full data set. The first 2 groups are typically called the past/training subset and the future/validation subset, where the future validation subset has K future data points. In the next iteration, the whole data set is then re-split using another disjoint subset of K data points as the future validation subset. This process is done iteratively until the training set has no less than minimum number of points remaining. Using this method, each data point can show up in at most 1 future validation subset. This method can be used for computing predictive loglikelihoods, in-sample loglikelihoods and doing cross-validation using subsets of future out-of-sample data. The default values of K and minimum are 1 and 2 respectively.

Example:

Assume the original data is ["A", "B", "C", "D", "E", "F"]. Leave-1-future-out splitting with minimum = 2 results in the following splittings of the data where the third column is discarded:

Leave-2-future-out splitting with minimum = 2 results in the following splittings:

LeaveK(; K = 5, shuffle = false, rng = nothing)

Leave-K-out splitting method. In this method, the data is split multiple times into 2 disjoint groups, each time starting from the full data set. The 2 groups are typically called training and validation subsets, where the validation subset has K data points. In the next iteration, the whole data set is then re-split using another disjoint subset of K data points as the validation subset. This process is done iteratively until almost each data point has shown up in 1 and only 1 validation subset. This method can be used for computing predictive loglikelihoods, in-sample loglikelihoods and doing cross-validation using subsets of out-of-sample data.

The data is typically a vector of some sort, e.g. observations or subjects, and the splittings are order-dependent. Before performing the splittings, you can randomly shuffle the data vector by setting the shuffle keyword argument to true (default is false) getting rid of the sensitivity to the original order of the data. You can additionally pass an optional pseudorandom number generator rng to control the pseudo-randomness for maximum reproducibility.

Example:

Assume the original data is ["A", "B", "C", "D"]. Leave-1-out splitting without shuffling results in the following splittings of the data:

where each data point showed once and only once in a validation subset.

Leave-2-out splitting without shuffling results in the following splittings:

MarginalMCMC

An instance of the marginal Hamiltonian Monte Carlo (HMC) based No-U-Turn Sampler (NUTS) algorithm where the subject-specific parameters are marginalized and only the population parameters are sampled. All the options of the algorithm can be passed to the constructor of MarginalMCMC as keyword arguments. Example:

alg = MarginalMCMC(; approx = LaplaceI(), nsamples = 2000, nadapts = 1000)

The main options that can be set are:

• approx: The likelihood-approximation algorithm used to marginalize out the subject-specific parameters. It can be set to FO(), FOCE(), or LaplaceI(), for example.
• target_accept (default: 0.8): Target acceptance ratio for the NUTS algorithm.
• nsamples (default: 2000): Number of Markov Chain Monte Carlo (MCMC) samples to generate, including adaptation samples.
• nadapts (default: nsamples ÷ 2): Number of adaptation steps in the NUTS algorithm.
• nchains (default: 4): Number of MCMC chains to sample.
• ess_per_chain (default: ∞): Target effective sample size (ESS) per chain. Sampling terminates if the target is reached.
• check_every (default: max(2, ess_per_chain ÷ 5)): Number of samples after which the ESS per chain is checked; only relevant when ess_per_chain is finite.
• time_limit (default: ∞): Time limit for sampling in seconds. Sampling terminates if the time limit is reached.

To customize the differential equation solver, the parallelism, the random number generator, or to fix some parameters during sampling, use the diffeq_options, ensemblealg, parallel_chains, parallel_subjects, rng, and constantcoef keyword arguments of fit.

NaivePooled()

A likelihood approximation method that ignores any contribution from random effects to the likelihood.

NaivePooled() is particularly useful in models that do not have any random effects, and for which therefore FO, FOCE, LaplaceI, or SAEM are inappropriate as there are no unobserved variables to marginalize out.

If there are random effects present in the model that should be ignored, it is necessary to set the dispersion of the random effects to zero by specifying parameters from @params appropriately. When fitting a model, the relevant parameter values have to be fixed with the constantcoef keyword argument.

See also FO, FOCE, LaplaceI, and SAEM.

NaturalScale()

The parameters in their original domain, i.e., the scale that coef and variance-covariance matrices are expressed on.

OrderedLogit(μ::AbstractVector{<:Real})

An ordered logit distribution is parameterized by a vector μ of ordered logits (of length K-1 implying a support of 1:K).

The cumulative probabilities of a random variable X with ordered logit distribution are given by

\[P(X \leq k) = \begin{cases} \frac{1}{1 + \exp{-\mu_k}} & \text{for } k \in \{1, \ldots, K-1\},\\ 1 & \text{for } k = K. \end{cases}\]

OrderedLogit(μ)  # OrderedLogit distribution with ordered logits `μ`
ncategories(d)     # Get the number of categories, i.e. K

Here, μ must be a real vector in monotonically increasing order.

External links:

• Ordered Logit on Wikipedia

OrderedVectorDomain(n::Int; lower=-∞)
OrderedVectorDomain(; init::AbstractVector{<:Real}, lower=-∞)

Return an ordered vector domain of dimension n, optionally lower-bounded by lower.

Any n-dimensional vector with first element ≥ lower and increasing elements is in this domain. Instead of passing n, an initial value init can be passed as a keyword argument.

@param x ∈ PDiagDomain(n::Int)
@param x ∈ PDiagDomain(; init::AbstractVector{<:Real})

Specify a parameter as an n-by-n positive semi-definite diagonal matrix by providing either the number n of rows and columns or the initial diagonal (variance) elements init.

If the number n of rows and columns in the matrix is specified, the initial diagonal elements are all set to 1.

@param x ∈ PSDDomain(n::Int)
@param x ∈ PSDDomain(; init=Matrix{Float64}(I, n, n))

Specifies a parameter as a symmetric n-by-n positive semi-definite matrix. The user must either specify n as the number of rows and columns in the matrix or init as the initial value of the parameter. init should be a positive semi- definite Matrix of Float64s.

PSISCrossvalidation(; split_method, split_by, filter_nans = false, pareto_shape_threshold = Inf)

Define an instance of the PSISCrossvalidation algorithm for crossvalidation.

In this algorithm, the fitting is not re-run but the MCMC samples are re-weighted using importance sampling (IS) while leaving out a subset of the data each time. The weights are then smoothed using a generalized Pareto distribution. The smoothed weights are then used to estimate the predictive likelihood of the out-of-sample data. This algorithm is known as the Pareto-smoothed importance sampling, leave-one-out (PSIS-LOO) crossvalidation method. In Pumas, we generalize the approach to handle any training-validation splitting of the data.

The shape parameter k of the generalized Pareto distribution fitted to each set of importance weights is a diagnostic for the reliability of the corresponding predictive-likelihood estimate: values k > 0.7 indicate that the PSIS approximation is unreliable for that split. For any such split, a warning is emitted automatically.

By default, no splits are removed and elpd returns the sum of pointwise log predictive densities over all held-out splits. The optional keyword arguments filter_nans and pareto_shape_threshold can be used for filtering splits instead: setting filter_nans = true drops splits whose PSIS smoothing fails (produces NaN), and setting pareto_shape_threshold to a finite value additionally drops splits with pareto shape parameter greater than that threshold. Splits removed in this way are recorded in the problematic_indices field of the returned PSISCVResult, and elpd is then computed from the remaining "good" splits only.

The way by which the data is split between training and validation sets is determined using the keyword arguments split_method and split_by. The split_method argument can be of any of the following types:

• LeaveK for leave-K-out cross-validation
• KFold for K-fold cross-validation
• LeaveFutureK for leaving K future points at a time

while split_by can be of any of the following types:

• BySubject for leaving out subjects
• ByObservation for leaving out individual observations per subject

Please refer to the documentation of each of the above types for the parameters of each split method, e.g. typing ?LeaveK or ?BySubject.

Examples:

Assume there are 5 subjects and 10 observations per subject and that ft is the result of the fit function. The following are some of the combinations in which the above inputs can be used:

• Leave-1-observation-out cross-validation, leaving 1 observation for all the subjects at a time. allsubjects = true means that the same observation index is removed for all the subjects, e.g. the 10th observation for all the subjects is used for validation in the first run, then the 9th observation is used for validation in the second run, etc.

      split_method = LeaveK(K = 1)
      split_by = ByObservation(allsubjects = true)
      cv_method = PSISCrossvalidation(; split_method, split_by)
      r = crossvalidate(ft, cv_method)

• Leave-2-future-observations-out cross-validation, leaving 2 future observations per subject such that no less than 4 observations are used in training. allsubjects = false means that only the data of one subject at a time will be split. So in the first run, only observations 1 to 4 of subject 1 are used for training and observations 5 and 6 of subject 1 are used for validation. In the second run, observations 1 to 6 of subject 1 are used for training and observations 7 and 8 are used for validation. In the third run, observations 1 to 8 of subject 1 get used for training and observations 9 and 10 are used for validation. For all 3 runs, all the observations of subjects 2 to 5 are used in training. Then in the fourth to sixth runs, subject 2's data gets split. In the seventh to ninth runs, subject 3's data gets split, etc.

      split_method = LeaveFutureK(K = 2, minimum = 4)
      split_by = ByObservation(allsubjects = false)
      cv_method = PSISCrossvalidation(; split_method, split_by)
      r = crossvalidate(ft, cv_method)

• Leave-1-subject-out cross-validation, marginalizing the subject-specific parameters using FOCE when computing predictive likelihoods. Using this method, subjects {1, 2, 3, 4} are used for training in the first run and subject {5} is used for validation. In the second run, subjects {1, 2, 3, 5} are used for training and subject {4} is used for validation, etc. The predictive likelihoods computed are the marginal likelihood of the subject approximated using FOCE. FO and LaplaceI can also be used instead.

      split_method = LeaveK(K = 1)
      split_by = BySubject(marginal = FOCE())
      cv_method = PSISCrossvalidation(; split_method, split_by)
      r = crossvalidate(ft, cv_method)

• Leave-1-future-subject-out cross-validation, where the predictive likelihoods are conditional likelihoods computed using 0s for the subject-specific parameters instead of marginalizing. The minimum number of subjects is 3 so there are only 2 runs. In the first run, the training subjects are {1, 2, 3} and the validation subject is {4}. In the second run, the training subjects are {1, 2, 3, 4} and the validation subject is {5}.

      split_method = LeaveFutureK(K = 1, minimum = 3)
      split_by = BySubject(marginal = nothing)
      cv_method = PSISCrossvalidation(; split_method, split_by)
      r = crossvalidate(ft, cv_method)

ParamSet(params::NamedTuple)

Construct a Pumas.jl parameter set.

params should be a NamedTuple that maps a parameter name to a Domain or a Distribution.

Example

ParamSet((;
  tvKa = RealDomain(lower = 0.0),
  tvCL = RealDomain(lower = 0.0),
  tvω = PDiagDomain(2),
))

A Population is an AbstractVector of Subjects.

ProportionalGamma(μ::Real, σ::Real)

Construct a gamma distribution with mean μ > 0 and variance (σ * μ)^2, defined in terms of the constant coefficient of variation σ ≥ 0.

Equivalently a Gamma with shape 1 / σ^2 and scale μ * σ^2.

See also ProportionalNormal.

ProportionalNormal(μ::Real, σprop::Real)

Construct a normal distribution with mean μ and standard deviation σprop * abs(μ) defined in terms of proportional standard deviation σprop ≥ 0.

See also CombinedNormal, CombinedCorrelatedNormal.

@param x ∈ RealDomain(;lower=-∞,upper=∞,init=0)

Specifies a parameter as a real value. lower and upper are the respective bounds, init sets the initial value of the parameter and is used for init_params.

SAEM(; kwargs...)

Stochastic approximation expectation-maximization (SAEM) algorithm for maximum-likelihood estimation of population parameters.

Keyword Arguments

• iters::Tuple{Int,Int,Int} (default: (200, 100, 100)): a tuple indicating the number of iterations to run in the exploratory, middle, and smoothing phases.
• α::Real (default: 0.7): the learning rate for the middle phase. On the kth update, α^k weight is placed on the kth iteration, and 1 - α^k on the previous.
• λ::Real (default: 0.975): maximum change rate of standard deviation parameters of random effects and error models in consecutive iterations. A value of 0.95 means that the value x of a standard deviation parameter is restricted to the interval [0.95x, x/0.95] in the following iteration. Slowing the decline in standard deviation parameters increases exploration, while limiting the growth prevents excessive flattening of the log-likelihood.
• repeat::Int (default: 1): number of times the iterations (exploratory, middle and smoothing phases) are repeated.
• verbosity::Int (default: 0): verbosity level. Available options are 0 (silent), 1 (printing once per repetition), and 2 (printing once per iteration).

Log-Likelihood

Since SAEM is an expectation-maximization algorithm, the (log-)likelihood is not calculated during fitting. To approximate the log-likelihood of a Pumas model fpm_saem fitted with SAEM, you can use likelihood approximation methods such as FOCE or LaplaceI:

loglikelihood(fpm_saem, FOCE())
loglikelihood(fpm_saem, LaplaceI())

See also FO, FOCE, LaplaceI, and NaivePooled.

SIR(;
    rng::Random.AbstractRNG = Random.default_rng(),
    samples::Union{Int,Vector{Int}},
    resamples::Union{Int,Vector{Int}},
    proposal_scale::Union{NaturalScale,UnconstrainedScale} = UnconstrainedScale(),
    proposal_shape::Union{Nothing,YeoJohnson} = nothing,
)

Construct an SIR object that can be used in infer function calls to perform sampling importance re-sampling (SIR) with the supplied settings. See ?infer and the Inference documentation for more information about obtaining inference results.

Confidence intervals based on SIR respect bounds on the parameters (e.g., the confidence interval for a strictly positive parameter cannot contain negative values).

In each iteration, parameter samples are drawn from a multivariate normal proposal distribution whose mean is set to the parameter estimates of the fitted Pumas model. In the initial iteration, the covariance matrix of the proposal distribution is either supplied by the user or derived from the inference results that the SIR procedure is based on: It is the variance-covariance matrix if the SIR procedure is based on the results of an infer call, and the covariance matrix of the parameter samples if it is based on the results of a Bootstrap, SIR, or MarginalMCMC inference (see infer). In subsequent iterations, the covariance matrix of the proposal distribution is set to the empirical covariance matrix of the parameter samples resampled in the previous iteration. In each iteration, parameter samples are resampled by weighted sampling without replacement.

The proposal_scale and proposal_shape keyword arguments select the space in which a proposal distribution is fitted to parameter samples. They apply to the proposal distributions of the second and all subsequent iterations, and to the one of the initial iteration if the SIR procedure is based on the parameter samples of a Bootstrap, SIR, or MarginalMCMC inference. A proposal distribution that is built from a covariance matrix is used in the space of that matrix, and neither keyword argument affects it.

Keyword Arguments

• rng: the pseudo-random number generator used within the SIR procedure.

• samples: the number(s) of parameter samples simulated from the proposal distribution in each iteration.

• resamples: the number(s) of parameter samples to resample in each iteration. In each iteration, the number of resampled parameter samples must not be greater than the number of originally sampled parameters. In every iteration but the final one, it must be greater than the number of estimated parameters as well since otherwise the covariance matrix of the proposal distribution of the next iteration is singular.

• proposal_scale: the scale of the parameters that a proposal distribution is fitted on. Supported values:

  • UnconstrainedScale() (default): the unconstrained parameterization that Pumas uses internally. Every sample corresponds to a valid parameter value, and hence no sample is rejected.
  • NaturalScale(): the parameters in their original domain, which is also the scale that covariance matrices are expressed on. Samples that violate the domain constraints of the parameters are rejected.

• proposal_shape: how the marginals of the parameter samples are reshaped before a proposal distribution is fitted to them. Supported values:

  • nothing (default): not at all.
  • YeoJohnson(): with a Yeo-Johnson transformation, one parameter per marginal, estimated from the samples. Samples outside of the image of the transformation are rejected.

Both samples and resamples must be supplied since there are no default values for these keyword arguments.

SimplexDomain(n::Int)

Specifies a parameter as a probability vector of length n constrained to the unit simplex (all elements non-negative, summing to 1). The initial value is a uniform distribution fill(1/n, n).

Pumas.SimulatedObservations

Simulation results mainly obtained through a simobs call. Possible uses:

- Check ODE solution with [`isvalid`](@ref).
- [`Subject`](@ref) construction through `Subject(sim::SimulatedObservations)`. Also works with [`Population`](@ref) through `Subject.(sim::SimulatedObservations)`.
- Conversion to `DataFrame` with `DataFrame(sim::SimulatedObservations)`.

- Post-processed by [`postprocess`](@ref) to obtain queries and statistics.

Pumas.SimulatedPopulation

Vector of Pumas.SimulatedObservations returned by simobs. Possible uses:

- Post-processing by [`postprocess`](@ref).
- Generate a visual predictive check with [`vpc`](@ref).
- Organized in a `DataFrame` with `DataFrame(sims::SimulatedPopulation)`.
- Getting statistical summaries through `Statistics.stat(func, sims::SimulatedPopulation)`, where `stat` can be `cov`, `cor`, `mean`, `var`, `std`. For the last 3 options, the `func` argument can be omitted.

Subject

The data corresponding to a single subject:

Fields:

• id: identifier
• observations: a named tuple of the dependent variables
• covariates: a named tuple containing the covariates, or nothing.
• events: a DosageRegimen, or nothing.
• time: a vector of time stamps for the observations

When there are time varying covariates, each covariate is interpolated with a common covariate time support. The interpolated values are then used to build a multi-valued interpolant for the complete time support.

From the multi-valued interpolant, certain discontinuities are flagged in order to use that information for the differential equation solvers and to correctly apply the analytical solution per region as applicable.

Constructor

Subject(;id = "1",
         observations = NamedTuple(),
         events = nothing,
         time = observations isa AbstractDataFrame ? observations.time : nothing,
         covariates::Union{Nothing, NamedTuple} = nothing,
         covariates_time = observations isa AbstractDataFrame ? observations.time : nothing,
         covariates_direction = :left)

Subject may be constructed from an <:AbstractDataFrame with the appropriate schema or by providing the arguments directly through separate DataFrames / structures.

Examples:

julia> Subject()
Subject
  ID: 1

julia> Subject(
         id = 20,
         events = DosageRegimen(200, ii = 24, addl = 2),
         covariates = (WT = 14.2, HT = 5.2),
       )
Subject
  ID: 20
  Events: 3
  Covariates: WT, HT

julia> Subject(covariates = (WT = [14.2, 14.7], HT = fill(5.2, 2)), covariates_time = [0, 3])
Subject
  ID: 1
  Covariates: WT, HT

Subject

Constructor

Subject(simsubject::SimulatedObservations)

Roundtrip the result of simobs, i.e. SimulatedObservations to a Subject/Population

Example:

sims = simobs(model, pop, params)

To convert sims to a Population, broadcast as below

Subject.(sims)

Subject(subject::Subject; id::AbstractString, events)

Construct a Subject from an existing subject while replacing information according to the keyword arguments. The possible keyword arguments are

• id, sets the subject identifier to the string id. Defaults to the identifier already in subject.
• events, sets the subject events to the input DosageRegimen. Defaults to the events already present in subject.
• covariates and covariates_times, used as input to add covariate information

Examples:

julia> s1 = Subject(; id = "AKJ491", events = DosageRegimen(1.0; time = 1.0))
Subject
  ID: AKJ491
  Events: 1

julia> Subject(s1; id = "AKJ492")
Subject
  ID: AKJ492
  Events: 1

TimeToEvent{T}

Distribution-like struct to store the hazard and the cumulative hazard in a time-to-event model. The dependent variable in a model that uses TimeToEvent should be a censoring variable that is one when the variable is observed (not censored) zero when the variable is right censored. Currently, no other censoring types are supported.

Example

...
@pre begin
  θeff = θ*DOSE
  λ = λ₀*exp(θeff)
end

@dynamics begin
  Λ' = λ
end

@derived begin
  DV ~ @. TimeToEvent(λ, Λ)
end
...

UnconstrainedScale()

The unconstrained parameterization of the parameters that Pumas uses internally, in which every vector of real numbers corresponds to a valid parameter value.

@param x ∈ VectorDomain(n::Int; lower=-∞,upper=∞,init=0)

Specifies a parameter as a real vector of length n. lower and upper are the respective bounds, init sets the initial value of the parameter and is used for init_params. The keyword arguments can all take either numbers or vectors. If numbers then the same value will be applied to all n vector elements. If you specify a vector (of length n) then you can adjust the parameters of the VectorDomain elements individually.

YeoJohnson()

Yeo-Johnson transformation of the marginals of the parameter samples, with one parameter $\lambda$ per marginal that is estimated from the samples. See SIR.

The transformation is

\[\psi(x, \lambda) = \begin{cases} \bigl((1 + x)^{\lambda} - 1\bigr) / \lambda, & x \ge 0, \lambda \ne 0, \\ \log(1 + x), & x \ge 0, \lambda = 0, \\ -\bigl((1 - x)^{2 - \lambda} - 1\bigr) / (2 - \lambda), & x < 0, \lambda \ne 2, \\ -\log(1 - x), & x < 0, \lambda = 2. \end{cases}\]

On the non-negative branch it is the Box-Cox transformation of $1 + x$, and on the negative branch the negative of the Box-Cox transformation of $1 - x$ with $2 - \lambda$ in place of $\lambda$. Reflecting the second branch in this way is what extends Box-Cox, which is defined for positive values only, to the whole real line: the result is an increasing continuously differentiable bijection of the real line onto an open interval for every $\lambda$, and the identity for $\lambda = 1$. Values $\lambda < 1$ compress the right tail of the marginal and $\lambda > 1$ its left tail.

$\lambda$ is estimated separately for every marginal by maximizing the profile likelihood of a normal distribution over the transformed marginal, i.e., by minimizing

\[f(\lambda) = \log \widehat{\sigma}^2(\lambda) - \frac{2 (\lambda - 1)}{n} \sum_{i=1}^n \operatorname{sign}(x_i) \log(1 + |x_i|),\]

where $\widehat{\sigma}^2(\lambda)$ is the variance of the transformed marginal and the sum is its mean log-Jacobian, which is what makes the objective comparable across $\lambda$.

References

Yeo, I.-K. and Johnson, R. A. (2000). A new family of power transformations to improve normality or symmetry. Biometrika 87(4), 954-959. https://doi.org/10.1093/biomet/87.4.954

ZeroSumDomain(n; init = zeros(n))
ZeroSumDomain(init)
ZeroSumDomain(; init)

Return a zero-sum domain of dimension n. Any n-dimensional vector with a sum of 0 is in this domain. Instead of passing n, an initial value init can be passed instead either as a positional or keyword argument.

isvalid(s::SimulatedObservations) -> Bool

Check if ODE solution of a simulation was successful.

Examples

mdl = @model begin
  @options begin
    checklinear = false # Force numerical integration of ODE
  end
  @param begin
    θkₑ ∈ RealDomain(lower = 0.0)
  end
  @pre begin
    kₑ = θkₑ
    cᵥ = 0.1
    V = 1.0
  end
  @dynamics begin
    Central' = -kₑ * Central
  end
  @derived begin
    μ = Central ./ V
    y ~ @. Gamma(1 / cᵥ^2, abs(μ) * cᵥ^2)
  end
end

sim = simobs(
  mdl,
  Subject(events = DosageRegimen(100), time = range(1.0, 40.0, length = 12)),
  (; θkₑ = 1e20),
)

julia> isvalid(sim)
false

sol = solve(
  model::AbstractPumasModel,
  population::Union{Population, Subject},
  param::NamedTuple,
  randeffs::NamedTuple=sample_randeffs(rng, model, population, param);
  ensemblealg=EnsembleSerial(),
  diffeq_options=NamedTuple(),
  rng=default_rng()
)

Solve the model applied to the Subject(s) within population using parameters param and random effects randeffs. By default, the times at which to save the solution are taken to be the observation times for each Subject within population. This can be overridden by supplying a vector or range as saveat to diffeq_options - e.g. diffeq_options=(; saveat=[0., 0.5, 1.]).

Arguments

• model may either be a PumasModel or a PumasEMModel.

• population may either be a Population of Subjects or a single Subject.

• param is parameter set in the form of a NamedTuple, e.g. (; tvCL=1., tvKa=2., tvV=1.).

• randeffs is an optional argument that, if used, should specify the random effects for each subject in population. Such random effects are specified by NamedTuples for PumasModels (e.g. (; tvCL=1., tvKa=2.)) and by Tuples for PumasEMModels (e.g. (1., 2.)). If population is a single Subject (without being enclosed in a vector) then a single such random effect specifier should be passed. If, however, population is a Population of multiple Subjects then randeffs should be a vector containing one such specifier for each Subject. The functions init_randeffs, zero_randeffs, and sample_randeffs are sometimes convenient for generating randeffs:

randeffs = zero_randeffs(model, population, param)
solve(model, population, param, randeffs)

If no randeffs is provided, then random ones are generated according to the distribution in the model.

• ensemblealg is a keyword argument that allows you to choose between different modes of parallelization. Options include EnsembleSerial(), EnsembleThreads() and EnsembleDistributed().

• diffeq_options is a keyword argument that takes a NamedTuple of options to pass on to the differential equation solver. The absolute tolerance, relative tolerance, and the maximum number of iterations of the numerical solvers for computing steady states can be specified with the ss_abstol, ss_reltol, and ss_maxiters options. By default, ss_abstol and ss_reltol are set to the absolute and relative tolerances of the differential equation solver, and ss_maxiters is set to 1000.

• rng is a keyword argument where you can specify which random number generator to use.

logpdf(d::Pumas.Constrained, x)

Evaluate the logarithm of the probability density of the constrained distribution, d, at x.

• If d is a constrained univariate distribution then x should be a scalar.
• If d is a constrained multivariate distribution then x should be a vector.

Evaluations of x outside of d's constraints returns -Inf.

Note that d itself is not a true distribution since its probability mass is not rescaled to 1.

cond(pmi::FittedPumasModelInference; correlation = true)

Return the condition number of the correlation matrix (or variance-covariance matrix with correlation = false) of the free parameters stored in pmi. Throws an error if pmi is the result of a call to infer with Pumas.Bootstrap or if the variance-covariance calculation failed.

gewekediag(b::BayesMCMCResults; subject::Union{Nothing,Int} = nothing, first::Real = 0.1, last::Real = 0.5)

Computes the Geweke diagnostic ^[Geweke1991] for each chain outputting a p-value per parameter. If subject is nothing (default), the chains diagnosed are those of the population parameters. If subject is set to an integer index, the chains diagnosed are those of the subject-specific parameters corresponding to the subject with the input index.

The Geweke diagnostic compares the sample means of two disjoint sub-chains X₁ and X₂ of the entire chain using a normal difference of means hypothesis test where the null and alternative hypotheses are defined as:

• H₀: μ₁ = μ₂
• H₁: μ₁ ≂̸ μ₂ where μ₁ and μ₂ are the population means. The first sub-chain X₁ is taken as the first (first * 100)% of the samples in the chain, where first is a keyword argument defaulting to 0.1. The second sub-chain X₂ is taken as the last (last * 100)% of the samples in the chain, where last is a keyword argument defaulting to 0.5.

The test statistic used is: z₀ = x̅₁ - x̅₂ / √(s₁² + s₂²) where x̅₁ and x̅₂ are the sample means of X₁ and X₂ respectively, and s₁ and s₂ are the Markov Chain standard error (MCSE) estimates of X₁ and X₂ respectively. Auto-correlation is assumed within the samples of each individual sub-chain, but the samples in X₁ are assumed to be independent of the samples in X₂. The p-value output is an estimate of P(|z| > |z₀|), where z is a standard normally distributed random variable.

Low p-values indicate one of the following:

• The first and last parts of the chain are sampled from distributions with different means, i.e. non-convergence,
• The need to discard some initial samples as burn-in, or
• The need to run the sampling for longer due to lack of samples or high auto-correlation.

High p-values indicate the inability to conclude that the means of the first and last parts of the chain are different with statistical significance. However, this alone does not guarantee convergence to a fixed posterior distribution because:

• Either the standard deviations or higher moments of X₁ and X₂ may be different, or
• The independence assumption between X₁ and X₂ may not be satisfied when high auto-correlation exists.

heideldiag(b::BayesMCMCResults; subject::Union{Nothing,Int} = nothing, alpha::Real = 0.05, eps::Real = 0.1, start::Integer = 1)

Compute the Heidelberger and Welch diagnostic ^{[Heidelberger1983]} for each chain. If subject is nothing (default), the chains diagnosed are those of the population parameters. If subject is set to an integer index, the chains diagnosed are those of the subject-specific parameters corresponding to the subject with the input index.

The Heidelberger diagnostic attempts to:

• Identify a cutoff point for the initial transient phase for each parameter, after which the samples can be assumed to come from a steady-state posterior distribution. The initial transient phase can be removed as burn-in. The cutoff point for each parameter is given in the burnin column of the output dataframe.
• Estimate the relative confidence interval (confidence interval divided by the mean) for the mean of the steady-state posterior distribution of each parameter, assuming such steady-state distribution exists in the samples. A large confidence interval implies either the lack of convergence to a stationary distribution or lack of samples. Half the relative confidence interval is given in the halfwidth column of the output dataframe. The test column will be true (1) if the halfwidth is less than the input target eps (default is 0.1) and false (0) otherwise. Note that parameters with mean value close to 0 can have erroneous relative confidence intervals because of the division by the mean. The test value can therefore be expected to be false (0) for those parameters without concluding lack of convergence.
• Quantify the extent to which the distribution of the samples is stationary using statistical testing. The returned p-value, shown in the pvalue column of the output dataframe, can be considered a measure of mean stationarity. A p-value lower than the input threshold alpha (default is 0.05) implies lack of stationarity of the mean, i.e. the posterior samples did not converge to a steady-state distribution with a fixed mean.

The Heidelberger diagnostic only tests for the mean of the distribution. Therefore, it can only be used to detect lack of convergence and not to prove convergence. In other words, even if all the numbers seem normal, one cannot conclude that the chain converged to a stationary distribution or that it converged to the true posterior.

rafterydiag(b::BayesMCMCResults; subject::Union{Nothing,Int} = nothing, q = 0.025, r = 0.005, s = 0.95, eps = 0.001)

Compute the Raftery and Lewis diagnostic ^{[Raftery1992]}. This diagnostic is used to determine the number of iterations required to estimate a specified quantile q within a desired degree of accuracy. The diagnostic is designed to determine the number of autocorrelated samples required to estimate a specified quantile $\theta_q$, such that $\Pr(\theta \le \theta_q) = q$, within a desired degree of accuracy. In particular, if $\hat{\theta}_q$ is the estimand and $\Pr(\theta \le \hat{\theta}_q) = \hat{P}_q$ the estimated cumulative probability, then accuracy is specified in terms of r and s, where $\Pr(q - r < \hat{P}_q < q + r) = s$. Thinning may be employed in the calculation of the diagnostic to satisfy its underlying assumptions. However, users may not want to apply the same (or any) thinning when estimating posterior summary statistics because doing so results in a loss of information. Accordingly, sample sizes estimated by the diagnostic tend to be conservative (too large).

Furthermore, the argument r specifies the margin of error for estimated cumulative probabilities and s the probability for the margin of error. eps specifies the tolerance within which the probabilities of transitioning from initial to retained iterations are within the equilibrium probabilities for the chain. This argument determines the number of samples to discard as a burn-in sequence and is typically left at its default value.

Censored(distribution::UnivariateDistribution, lower::Real, upper::Real)

Construct a censored distribution based on distribution and the censoring limits lower and upper.

apply_stat(stat, v::AbstractVector)

Apply the summarizing function stat to the contents of the vector in a way that is equivalent to calling stat on a vector of all the elements that are in the ith position of each element of v. For example, this can be used to create a NamedTuple of the mean across a vector of NamedTuples such that the key k in the output is the mean of all the values belonging to the key k in each element of the vector v.

Example

julia> v = [(a = 3, b = 4), (a = 1, b = 9), (a = 6, b = 0)]
3-element Vector{NamedTuple{(:a, :b), Tuple{Int64, Int64}}}:
 (a = 3, b = 4)
 (a = 1, b = 9)
 (a = 6, b = 0)

julia> apply_stat(mean, v)
(a = 3.3333333333333335,
 b = 4.333333333333333,)

julia> apply_stat(var, v)
(a = 6.333333333333334,
 b = 20.333333333333336,)

julia> apply_stat(std, v)
(a = 2.5166114784235836,
 b = 4.509249752822894,)

bootstrap_fits(fpm::FittedPumasModelInference{<:Any, <:Bootstraps})
bootstrap_fits(fpm::Bootstraps)

Return a vector of the fits of the resampled populations in a Bootstrap inference.

censored_latent(
    distribution::UnivariateDistribution,
    lower::Union{Real,Nothing},
    upper::Union{Real,Nothing},
)
censored_latent(
    distribution::UnivariateDistribution;
    lower::Union{Real,Nothing} = nothing,
    upper::Union{Real,Nothing} = nothing,
)

Construct a censored distribution based on distribution and the censoring limits lower and upper. When fitting a model with a censored_latent dependent variable, the log probability of the censored distribution is used in the likelihood calculation. However, when calculating the mean or other summary statistics, the underlying uncensored distribution is used.

Examples

When censoring a distribution it is necessary to provide lower and upper bounds. Keyword arguments can be used to specify the bounds. The default is to use nothing for both lower and upper bounds to indicate no censoring at either end. The keyword for can be used whenever the bound is fixed for all subjects and observations. If the cp is the modelled concentration and σ is the standard deviation of the error, then the following code will create a censored normal distribution with a lower bound of 1e-3 similar to how the functions would be used in a @derived block.

julia> cp = [10.0, 4.0, 1.0]
3-element Vector{Float64}:
 10.0
  4.0
  1.0

julia> σ = 0.35
0.35

julia> @. censored_latent(Normal(cp, σ); lower = 1e-3);

This can be used to set a lower censoring point, and upper censoring point (with the upper keyword) or both. However, if the censoring point depends on a covariate we cannot broadcast (@.) if the keyword is a vector. Then we can use the first version mentioned at the top where limits are given af positional arguments.

julia> cp = [10.0, 4.0, 1.0]
3-element Vector{Float64}:
 10.0
  4.0
  1.0

julia> lower = [1e-3, 1e-3, 1e-3]
3-element Vector{Float64}:
 0.001
 0.001
 0.001

julia> σ = 0.35
0.35

julia> @. censored_latent(Normal(cp, σ), lower, nothing);

Similarly, the second positional argument can be nothing and the third argument can be a vector or all three arguments can be vectors if both lower and upper limits are relevant.

center_randeffs(
  model::Union{PumasModel,PumasEMModel},
  [subject_or_population::Union{AbstractSubject,AbstractPopulation},]
  param::NamedTuple,
)

Return a set of central random effects for model model using parameters param.

By default, random effects are set to the mean of their distributions. For some distributions, however, the median or the mode are used instead.

The optional argument subject_or_population takes a subject or population for which the random effects are returned. If a population of subjects is provided, a vector of random effects is returned with one element for each subject within the population.

julia> Pumas.center_stat(Normal)
mean (generic function with 130 methods)

julia> Pumas.center_stat(LogNormal)
median (generic function with 52 methods)

conditional_nll(
    model::AbstractPumasModel,
    subject::Subject,
    param::NamedTuple,
    randeffs::NamedTuple;
    diffeq_options::NamedTuple = NamedTuple(),
)

Compute the conditional negative log-likelihood of model m for subject with parameters param and random effects randeffs. diffeq_options is a NamedTuple of options passed to the ODE solver. Requires that the derived produces distributions.

correlation_diagnostic(
  pmi::FittedPumasModelInference;
  high_cor_threshold::Real = 0.7,
  medium_cor_threshold::Real = high_cor_threshold / 2,
)

Extract lists of parameter pairs with high, medium or low correlations from an inference result output by infer.

Parameter pairs with high correlation are those with an absolute correlation higher than high_cor_threshold. Parameter pairs with a low correlation are those with an absolute correlation lower than medium_cor_threshold. The remaining parameter pairs have a medium correlation.

Example

Assume fpm is the output of fit. The following is an example of how to use correlation_diagnostic:

fpmi = infer(fpm)
correlation_diagnostic(fpmi)

function covariate_select(
    model::PumasModel,
    population::Population,
    param::NamedTuple,
    approx::LikelihoodApproximation;
    control_param::Union{Tuple,Nothing},
    criterion = aic,
    method::CovariateSelection.T = CovariateSelection.Forward,
    fit_options::NamedTuple = (;),
    verbose::Bool = false,
) -> CovariateSelectionResults

Use a covariate selection procedure to find the best fit of model to population with starting values param and estimation method/likelihood approximation approx according to criterion by restricting parameters to zero for subsets of control_param traversed according to method.

The model should be written such that it includes all relevant covariates, possibly included with multiple different functional forms, such that restricting on a parameter listed in control_param to zero eliminates the effect of the covariate (see example below).

Currently, three methods are supported: CovariateSelection.Forward (the default), CovariateSelection.Backward, and CovariateSelection.All. The method CovariateSelection.Forward starts from the models where all parameters in control_param are restricted to zero. In the next step, a model is fit where all parameters in control_param but one are restricted to zero. If none of these models attain a smaller value of criterion then the model that restricts all parameters to zero is chosen. Otherwise, the model with the lowest criterion is chosen and the procedure continues by fitting all models with two estimated parameters where one of the estimated parameters must be the one corresponding to the best fit in the previous step. New steps are taken as long as the criterion is reduced. If the CovariateSelection.Backward method is chosen then the procedure starts with all parameters in control_param unrestricted and then one parameter is restricted at a time until there is no further improvement. The CovariateSelection.All method fits all possible 2^n models from the power set of control_param and selects the one with the best criterion value. This is the most thorough approach but scales exponentially with the number of control parameters.

Any callable object that returns a real value when a fitted model is passed can be used as criterion. It is always assumed that a smaller value is better. The default value is the Akaike Information Criterion aic but another commonly used criterion is the Bayesian Information Criterion bic which puts a higher penalty on the number of parameters compared to aic.

To avoid that logging messages from the calls to fit fill up the output window or default log, all logging messages produced by fit are by default suppressed. This might hide important issues with the model fitting so for it is possible to enable complete logging by passing verbose = true.

Example

model = @model begin
  @param begin
    β0 ∈ RealDomain()
    β1 ∈ RealDomain()
    β2 ∈ RealDomain()
    β3 ∈ RealDomain()
    β4 ∈ RealDomain()
  end
  @covariates x1 x2 x3 x4
  @pre μ = β0 + β1 * x1 + β2 * x2 + β3 * x3 + β4 * x4
  @derived begin
    y ~ Normal.(μ, 0.1)
  end
end

selection = covariate_select(
  model,
  population,
  init_params(model),
  NaivePooled();
  control_param = (:β1, :β2, :β3, :β4),
  method = CovariateSelection.Forward,
  criterion = bic,
  fit_options = (; constantcoef = (:β2,)),
)

crossvalidate(b::BayesMCMCResults, cv_method::ExactCrossvalidation; ensemblealg = EnsembleThreads())

Runs cross-validation by re-running the MCMC sampling while leaving out some data points each time. The first argument b is the output of the MCMC fit. The choice of how the data is split in each run is specified in the cv_method argument which is an instance of ExactCrossvalidation. Please refer to the documentation of ExactCrossvalidation by typing ?ExactCrossvalidation for details on the parameters of the method.

The cross-validation re-fits are parallelized across the splits of the data according to ensemblealg.

crossvalidate(b::BayesMCMCResults, cv_method::PSISCrossvalidation)

Performs cross-validation using the Pareto smoothed importance sampling (PSIS) approach by re-weighing the samples while leaving out some data points each time. The first argument b is the output of the MCMC fit. The choice of how the data is split in each run is specified in the cv_method argument which is an instance of PSISCrossvalidation.

discard(b::BayesMCMCResults; burnin = 0, ratio = 1.0)

Removes burnin samples from the beginning of each chain and then keeps only ratio * 100% of the remaining samples in each chain (i.e. thinning).

dosecontrol(fpm::AbstractFittedPumasModel)::Vector{ConstantInterpolationStructArray}

Return the dosecontrol parameters from fpm. The dosecontrol parameters are the variables bioav, lags, duration, rate if they're defined in the @dosecontrol block in the @model macro. The function returns a vector of covariate interpolation objects. Each of them can be evaluated at an arbitrary time point, but the interpolant is piecewise constant between the event times. Note, that in the model solution, subject simulation, and parameter fitting the dose control parameters can be time-varying, but they only affect the dosing according to their value at the beginning of each dose event. This means that the dose control parameters are only relevevant at these time points. Each of the covariate interpolation objects can be converted to a DataFrame by calling the DataFrame constructor on each element, or a DataFrame with all subjects can be obtained by calling the DataFrame constructor on the vector output by this function. This DataFrame will only contain the times that correspond to the start of each dose event for the reason discussed above.

The computation is parallelized across subjects.

dosecontrol(
    model::PumasModel,
    population::Population,
    param::NamedTuple,
    vrandeffs::AbstractVector{<:NamedTuple};
    kwargs...,
)

Return the dosecontrol parameters of the model (i.e., the variables bioav, lags, duration and rate if defined in the @dosecontrol block in the @model macro) for each subject in the population, evaluated with population parameters param and random effects vrandeffs (one named tuple of random effects for each subject).

The function returns a vector of covariate interpolation objects. Each of them can be evaluated at any time point, but the interpolant is piecewise constant between the event times. Each of the covariate interpolation objects can be converted to a DataFrame by calling the DataFrame constructor on each element, or a DataFrame with all subjects can be obtained by calling the DataFrame constructor on the vector output by this function.

Keyword Arguments

• ensemblealg::SciMLBase.EnsembleAlgorithm (default: EnsembleThreads()): mode of parallelization across subjects. Possible options are: EnsembleSerial(), EnsembleThreads() and EnsembleDistributed().

dosecontrol(
    model::PumasModel,
    subject::Subject,
    param::NamedTuple,
    randeffs::NamedTuple,
)

Return the dosecontrol parameters of the model (i.e., the variables bioav, lags, duration and rate if defined in the @dosecontrol block in the @model macro) of the subject, evaluated with population parameters param and random effects randeffs.

If the subject has observations, the individual coefficients are evaluated with population parameters param and random effects set to the empirical Bayes estimates (EBEs) given param, estimated with likelihood approximation method approx. If the subject does not have any observations, the individual coefficients are evaluated with population parameters param and random effects set to the mode of the domain-transformed prior given param.

The function returns a covariate interpolation object. It can be evaluated at any time point, but the interpolant is piecewise constant between the event times. Moreover, it can be converted to a DataFrame by calling the DataFrame constructor on the interpolation object output by this function.

dosenum(t, events) -> # of doses for each time

Creates an array that matches t in length which denotes the number of events that have occurred prior to the current point in time. If t is a scalar, outputs the number of events before that time point. Only dose events (evid 1 and 4) are counted.

dynamics_names(model)::Vector{String}

Return the names of the dynamic variables of model, i.e. the variables defined in the @dynamics block. For a model built from a ModelingToolkit system, these are the unknowns of the compiled system. The names match the dynamics columns of the DataFrames produced from simobs, predict, and inspect results.

elpd(r::CVResult)

Estimate of the expected log pointwise predictive density (ELPD) for the held-out data in the result r of crossvalidate.

The leave-one-out estimate of the ELPD is defined as

\[\widehat{\mathrm{elpd}}_{\mathrm{loo}} = \sum_{i=1}^{n} \log p(y_i \mid y_{-i}),\]

i.e. the sum of pointwise leave-one-out log predictive densities over the n held-out observations, following eq. 4 of ^{[Vehtari2017]}.

empirical_bayes(
  fpm::FittedPumasModel,
  data::Union{Subject, Population};
  [diffeq_options::NamedTuple],
)

Return empirical bayes estimates for a fitted model, fpm, when applied to data.

empirical_bayes(fpm::Pumas.FittedPumasModel)
empirical_bayes(fpm::Pumas.FittedPumasEMModel)
empirical_bayes(insp::Pumas.FittedPumasModelInspection)

Return sampled random effects or empirical bayes estimates from a fit or model inspection. If the model was estimated with the FO likelihood approximation methods the empirical bayes estimates will be obtained using the LaplaceI approximation. If either FOCE or LaplaceI was used the final empirical bayes estimates will be returned. If Pumas.SAEM was used to fit the empirical bayes estimates are obtained using the LaplaceI approximation.

empirical_bayes(
  model::PumasModel,
  data::Union{AbstractSubject, Population},
  param::NamedTuple,
  approx::LikelihoodApproximation;
  diffeq_options::NamedTuple = NamedTuple(),
)

Return sampled random effects or empirical bayes estimates. The empirical bayes estimates will be obtained using the LaplaceI approximation regardless of the value of the approx argument except for NaivePooled where an empty NamedTuple is returned.

eventnum(t, events) -> # of doses for each time

Creates an array that matches t in length which denotes the number of events that have occurred prior to the current point in time. If t is a scalar, outputs the number of events before that time point. The following events are counted: 1 (dose), 3 (reset), 4 (reset and dose).

expectation(
    g,
    ::MonteCarloExpectation,
    model::PumasModel,
    subject::Subject,
    dist::MvNormal;
    imaxiters::Int = 10000,
    diffeq_options::NamedTuple = (;),
    rng::Random.AbstractRNG = Random.default_rng(),
)

Compute the expected functional g of the predictions for the given subject with respect to the multivariate normal distribution dist of population parameters (truncated to the parameters' domains) and the subject-specific random effects of the Pumas model using Monte Carlo integration.

Typically, distribution dist is inferred, e.g., dist = MvNormal(pmi) where pmi is the output of infer. Population parameters are sampled using rejection sampling to ensure the parameter values are in the parameters' domains.

Keyword Arguments

• imaxiters: The number of Monte Carlo samples to use.
• diffeq_options: The options of the differential equation solver used.
• rng: The random number generator used in the sampling, which defaults to default random number generator in Julia.

expectation(
    g,
    ::MonteCarloExpectation,
    model::PumasModel,
    subject::Subject,
    param_dists::NamedTuple;
    imaxiters::Int = 10000,
    diffeq_options::NamedTuple = (;),
    rng::Random.AbstractRNG = Random.default_rng(),
)

Compute the expected functional g of the predictions for the given subject with respect to the distributions param_dists of population parameters and the subject-specific random effects of the Pumas model using Monte Carlo integration.

Keyword Arguments

• imaxiters: The number of Monte Carlo samples to use.
• diffeq_options: The options of the differential equation solver used.
• rng: The random number generator used in the sampling, which defaults to default random number generator in Julia.

findinfluential(m::AbstractPumasModel,
                data::AbstractPopulation,
                param::NamedTuple,
                approx::Union{MAP, LikelihoodApproximation};
                diffeq_options::NamedTuple = NamedTuple())
findinfluential(fpm::AbstractFittedPumasModel,
                diffeq_options::NamedTuple = NamedTuple())

Return a vector of the negative log-likelihood values.

It can be used before fitting a model in order to find the subjects that contribute the most to the log-likelihood value, or to identify subjects for which the evaluation of the log-likelihood fails. This can happen due to incorrect data, causing NaNs to show up.

It can also be used on a fitted model to identify the subjects with maximum contribution to the log-likelihood.

If used with MAP, the prior contribution will be ignored.

A DataFrame method is defined on the resulting object for easy viewing and plotting of the results.

fi = findinfluential(model, pop, parameters, FOCE())
fi_df = DataFrame(fi)

If the above results contain NaNs, the id column can be useful in finding problematic data points. To visualize the negative log-likelihood values across subjects, some options are available. A bar plot can be produced as follows:

using AlgebraOfGraphics
draw(
  data(fi_df) * mapping(:nll) * visual(BarPlot);
  axis = (;
    xlabel = "id",
    ylabel = "negative log-likelihood",
    xticks = 1:length(pop),
    xtickformat = x -> fi_df.id,
  ),
)

A histogram can be used to see the distribution:

using AlgebraOfGraphics
draw(
  data(fi_df) *
  mapping(:nll) *
  histogram() *
  visual(BarPlot);
  axis = (; xlabel = "negative loglikelihood"),

)

Alternatively, a density plot can be created:

using AlgebraOfGraphics
draw(
  data(fi_df) *
  mapping(:nll) *
  density() *
  visual(BarPlot);
  axis = (; xlabel = "negative loglikelihood"),
)

icoef(fpm::AbstractFittedPumasModel; obstimes)::Vector{ConstantInterpolationStructArray}

Return the individual coefficients from fpm. The individual coefficients are the variables defined in the @pre block in the @model macro. The function returns a vector of covariate interpolation objects defined on obstimes. If no obstimes are given, the time points used to define the covariates will be used. Each of them can be evaluated at any time point, but the interpolant is piecewise constant between the obstimes. Each of the covariate interpolation objects can be converted to a DataFrame by calling the DataFrame constructor on each element, or a DataFrame with all subjects can be obtained by calling the DataFrame constructor on the vector output by this function.

icoef(
    model::PumasModel,
    population::Population,
    param::NamedTuple,
    vrandeffs::AbstractVector{<:NamedTuple};
    kwargs...,
)

Return a vector of individual coefficients (i.e., the variables defined in the @pre block in the @model macro) for each subject in the population, evaluated with population parameters param and random effects vrandeffs (one named tuple of random effects for each subject).

The function returns a vector of covariate interpolation objects. Each of them can be evaluated at any time point, but the interpolant is piecewise constant between the interpolation time points. Each of the covariate interpolation objects can be converted to a DataFrame by calling the DataFrame constructor on each element, or a DataFrame with all subjects can be obtained by calling the DataFrame constructor on the vector output by this function.

Keyword Arguments

• obstimes::Union{AbstractVector{<:Real},Nothing} (default: nothing): Time points at which the individual coefficients are evaluated. If obstimes is nothing, for each subject the time points used to define the covariates are used.
• ensemblealg::SciMLBase.EnsembleAlgorithm (default: EnsembleThreads()): Mode of parallelization across subjects. Possible options are: EnsembleSerial(), EnsembleThreads() and EnsembleDistributed().

icoef(
    model::PumasModel,
    subject::Subject,
    param::NamedTuple,
    randeffs::NamedTuple;
    kwargs...,
)

Return the individual coefficients (i.e., the variables defined in the @pre block in the @model macro) of the subject, evaluated with population parameters param and random effects randeffs.

The function returns a covariate interpolation object. It can be evaluated at any time point, but the interpolant is piecewise constant between the interpolation time points. Moreover, it can be converted to a DataFrame by calling the DataFrame constructor on the interpolation object output by this function.

Keyword Arguments

• obstimes::Union{Vector{<:Real},Nothing} (default: nothing): Time points at which the individual coefficients are evaluated. If obstimes is nothing, for each subject the time points used to define the covariates are used.

infer(fpm, bts::Pumas.Bootstrap; level=0.95)

Perform bootstrapping by resampling the Subjects from the Population stored in fpm. The confidence intervals are calculated as the (1-level)/2 and (1+level)/2 quantiles of the estimated parameters. The number of samples used in the bootstrapping is bts.samples. bts.ensemblealg specifies the ensemblealg used here. If ensemblealg is EnsembleSerial(), a single thread will be used. If it is EnsembleThreads() (the default value), multiple threads will be used. If it is set to EnsembleDistributed() each fit will be sent to a worker in the worker pool. See the documentation of Bootstrap for more details on constructing an instance of Bootstrap.

infer(fpm, sir::SIR; proposal_covariance, level = 0.95, ensemblealg = EnsembleThreads()) -> FittedPumasModelInference

Perform sampling importance re-sampling (SIR) for the Subjects from the Population stored in the fitted Pumas model fpm, with the covariance matrix proposal_covariance of the initial proposal distribution. See SIR for the settings of the SIR procedure.

Rows and columns of proposal_covariance must be labelled with the names of the estimated parameters of fpm such that the covariance matrix is unambiguous, regardless of how the parameters are vectorized internally. Since proposal_covariance is supplied explicitly, the SIR procedure can be performed even if the computation of the variance-covariance matrix of fpm fails or yields a matrix that is not positive definite.

The confidence intervals are calculated as the (1-level)/2 and (1+level)/2 quantiles of the estimated parameters.

Keyword Arguments

• proposal_covariance: the covariance matrix of the initial proposal distribution, in the space of the parameter estimates. The initial proposal distribution is used in the space of that matrix, i.e., neither the proposal_scale nor the proposal_shape of sir applies to it (see SIR). It must be a DimensionalData.DimArray that is symmetric and positive definite over the free parameters, and its rows and columns must be labelled with the names of the estimated parameters. The names of the estimated parameters are the parameter entries of the rows of coeftable(fpm) for which constant is false: Constant parameters are not estimated and must be omitted.
• level::Real (default: 0.95): the confidence level of the confidence intervals.
• ensemblealg (default: EnsembleThreads()): the mode of parallelization. Possible options are EnsembleThreads() to use multi-threading and EnsembleSerial() to use a single thread.

Examples

Assume fpm is the output of fit and Σ is a covariance matrix of the estimated parameters.

using DimensionalData

sir = SIR(; samples = 1000, resamples = 200)

# Names of the estimated parameters, e.g., `["tvcl", "tvvc", "Ω₁,₁", "Ω₂,₁", "Ω₂,₂", "σ"]`
ct = coeftable(fpm)
labels = ct.parameter[.!ct.constant]

infer(fpm, sir; proposal_covariance = DimArray(Σ, (X(labels), Y(labels))))

See also SIR, infer, coeftable, vcov.

infer(fpm; kwargs...)

Estimate the variance-covariance matrix and return a struct used for inference based on the fitted model fpm.

Keyword Arguments

• level::Real (default: 0.95): the confidence level of the confidence intervals. The confidence intervals are calculated as the (1-level)/2 and (1+level)/2 quantiles of the estimated parameters.
• sandwich_estimator::Bool (default: true): whether to use the sandwich estimator of the variance-covariance matrix.
• rethrow_error::Bool (default: false): whether an error is thrown if the estimation of the variance-covariance matrix fails, or the exception is stored without throwing it.
• ensemblealg (default: ensemblealg that was used during model fitting): the mode of parallelization. Possible options are: EnsembleSerial(), EnsembleThreads(), EnsembleDistributed() and EnsembleSplitThreads().

See also vcov.

infer(fpm, alg::MarginalMCMC; kwargs...)

Infer the uncertainty of the parameters of the fitted model fpm by sampling from the marginal likelihood with the Markov Chain Monte Carlo (MCMC) algorithm alg.

Keyword Arguments

• level::Real (default: 0.95): the confidence level of the confidence intervals. The confidence intervals are calculated as the (1-level)/2 and (1+level)/2 quantiles of the sampled parameters.

infer(pmi::FittedPumasModelInference, sir::SIR; level = pmi.level, ensemblealg = EnsembleThreads()) -> FittedPumasModelInference

Perform sampling importance re-sampling (SIR) for the Subjects from the Population stored in the fitted Pumas model of the inference results pmi. See SIR for the settings of the SIR procedure.

The covariance matrix of the initial proposal distribution of the SIR procedure is derived from pmi: It is the variance-covariance matrix if pmi was computed with infer(fpm; ...), and the covariance matrix of the parameter samples if pmi was computed with Bootstrap, SIR, or MarginalMCMC. Hence the SIR procedure can be performed even if the computation of the robust variance-covariance matrix of the fitted Pumas model fails or yields a matrix that is not positive definite. In the latter case the initial proposal distribution is fitted to the parameter samples in the space that the proposal_scale and proposal_shape of sir select, whereas a variance-covariance matrix is used in its own space (see SIR).

The confidence intervals are calculated as the (1-level)/2 and (1+level)/2 quantiles of the estimated parameters.

Keyword Arguments

• level::Real (default: pmi.level): the confidence level of the confidence intervals.
• ensemblealg (default: EnsembleThreads()): the mode of parallelization. Possible options are EnsembleThreads() to use multi-threading and EnsembleSerial() to use a single thread. It is only used for the SIR procedure, hence the variance-covariance matrix of pmi can be computed with an ensemble algorithm that the SIR procedure does not support, e.g., infer(infer(fpm; ensemblealg = EnsembleDistributed()), sir).

Examples

Assume fpm is the output of fit.

sir = SIR(; samples = 1000, resamples = 200)

# Initial proposal based on the robust variance-covariance matrix of the fitted Pumas model
infer(infer(fpm), sir)

# Initial proposal based on the negative inverse Hessian
infer(infer(fpm; sandwich_estimator = false), sir)

# Initial proposal based on the covariance matrix of bootstrap estimates
infer(infer(fpm, Bootstrap(; samples = 200)), sir)

# Initial proposal based on the covariance matrix of marginal MCMC samples
infer(infer(fpm, MarginalMCMC(; approx = FOCE())), sir)

# Initial proposal based on the covariance matrix of the samples of a previous SIR procedure
infer(infer(fpm, sir), sir)

See also SIR, Bootstrap, MarginalMCMC, vcov.

init_params(model::PumasModel)

Create a parameter set populated with the initial parameters of PumasModel.

inspect(fpm::AbstractFittedPumasModel; wres_approx::LikelihoodApproximation, nsim::Int, rng)

Output a summary of the model predictions, residuals, Empirical Bayes estimates, and NPDEs (when requested).

Called on a fit output and allows the keyword argument wres_approx for approximation method to be used in residual calculation. The default value is the approximation method used for the marginal likelihood calculation in the fit that produced fpm, unless the estimation method does not natively perform likelihood approximation, e.g. expectation maximization algorithms. In this case, wres_approx must be explicitly specified if weighted residuals are to be calculated for models with Gaussian error models. The keyword nsim controls the number of times each subject is simulated for npde computation. A FittedPumasModelInspection object with pred, wres, ebes, ewres, and npdes is output.

The per-subject computations (predictions, weighted residuals, dose control parameters, NPDEs/EWRES, and individual coefficients) are parallelized across subjects.

lrtest(fpm_0, fpm_A[, approx])::LikelihoodRatioTest

Compute the likelihood ratio test statistic of the null hypothesis defined by the fitted Pumas Model fpm_0 against the alternative hypothesis defined by the fitted Pumas model fpm_A, using the likelihood approximation method approx.

The pvalue function can be used for extracting the p-value based on the asymptotic Χ²(k) distribution of the test statistic.

lrtest(loglikelihood_0::Real, loglikelihood_A::Real, dof_0::Integer, dof_A::Integer)::LikelihoodRatioTest

Compute the likelihood ratio test using number.

mean_of_LKJ(d::LKJCholesky)

Return the mean of the LKJ distribution with the same parameters as d.

postprocess(obs::SimulatedObservations)
postprocess(func::Function, obs::SimulatedObservations)
postprocess(obs::SimulatedPopulation; stat = identity)
postprocess(func::Function, obs::SimulatedPopulation; stat = identity)

Post-processes the output of simobs. The output of a simobs call stores the simulated observations but also all the intermediate values computed such as the parameters used, individual coefficients, dose control parameters, covariates, differential equation solution, etc. There are a number of post-processing operations you can do on the simulation output to compute various queries and summary statistics based on the simulations.

The postprocess function is a powerful tool that allows you to make various queries using the simulation results. There are multiple ways to use the postprocess function. The first way to use the postprocess function is to extract all the information stored in the simulation result in the form of a vector of named tuples. Each named tuple has all the intermediate values evaluated when simulating 1 run. Let sims be the output of any simobs operation.

generated = postprocess(sims)

generated is a vector of simulation runs. Hence, generated[i] is the result from the ith simulation run. Time-dependent variables are stored as a vector with 1 element for each time point where there is an observation. Alternatively, if obstimes was set instead when calling simobs, the time-dependent variables will be evaluated at the time points in obstimes instead.

The second way to use postprocess is by passing in a post-processing function. The post-processing function can be used to:

• Transform the simulated quantities, or
• Compare the simulated quantities to the observations.

We use the do syntax here which is short for passing in a function as the first argument to postprocess. The do syntax to pass in a post-processing function is:

postprocess(sims) do gen, obs
  ...
end

where gen is the named tuple of all generated quantities from 1 simulation run and obs is the named tuple of observations. Besides the model-generated quantities, gen also exposes the observation times as gen.t and the subject's dosage regimen as gen.events. For instance to query the ratio of simulated observations dv that are higher than the observed quantity dv at the respective time points, you can use:

postprocess(sims) do gen, obs
  sum(gen.dv .> obs.dv) / length(gen.dv)
end

gen.dv is the simulated vector of dv whose length is the same as the number of observation time points. obs.dv is the observation vector dv. gen.dv .> obs.dv returns a vector of true/false, with one element for each time point. The sum of this vector gives the number of time points where the simulation was higher than the observation. Dividing by the number of time points gives the ratio. When using postprocess in this way, the output is always a vector of the query results. In the query function body, you can choose to use only gen or only obs but the header must always have both gen and obs:

postprocess(sims) do gen, obs
  ...
end

The third way to use the postprocess function is to compute summary statistics of the simulated quantities or of functions of the simulated quantities. Summary statistics can be computed by passing a statistic function as the stat keyword argument. For example in order to estimate the probability that a simulated value is higher than an observation, you can use:

postprocess(sims, stat = mean) do gen, obs
  gen.dv .> obs.dv
end

This function will do 2 things:

• Concatenate the query results (e.g. gen.dv .> obs.dv) from all the simulation runs into a single vector.
• Compute the mean value of the combined vector.

In order to get the probability that the simulated quantity is higher than the observation for each time point, you can call the mean function externally as such:

generated = postprocess(sims) do gen, obs
  gen.dv .> obs.dv
end
mean(generated)

This returns a vector of probabilities of the same length as the number of time points without concatenating all the queries together.

To compute a summary statistic of all the generated quantity, you can also call:

postprocess(sims, stat = mean)

without specifying a post-processing function.

Besides mean, you can also use any of the following summary statistic functions in the same way:

• std for element-wise standard deviation
• var for element-wise variance
• cor for correlation between multiple quantities
• cov for covariance between multiple quantities

These functions can be passed in as the stat keyword argument to postprocess, e.g.,

generated = postprocess(sims, stat = std) do gen, obs
  ...
end

or

generated = postprocess(sims, stat = var) do gen, obs
  ...
end

The cor and cov statistics are unique in that they require a post-processing function that outputs a vector. For example, to estimate the correlation between the CL and Vc parameters in a 1-compartment model, you can use any of the following:

postprocess(s, stat = cor) do gen, obs
  [only(unique(gen.CL)), only(unique(gen.Vc))]
end

cov can be used instead of cor to compute the covariance matrix. The output of this operation is either a correlation or a covariance matrix.

probstable(fpm::FittedPumasModel)

Return a DataFrame with outcome probabilities of all discrete dependent variables.

read_pumas(df::AbstractDataFrame; kwargs...)

Import NMTRAN-formatted data from a data frame. Data stored in a file has to be read into a data frame first, e.g., with CSV.read.

Keyword Arguments

• observations::Vector{Symbol} (default: [:dv]): a vector of column names of dependent variables.

• covariates::Vector{Symbol} (default: Symbol[]): a vector of column names of covariates.

• id::Symbol (default: :id): the name of the column with the IDs of the individuals. Each individual should have a unique integer or string.

• time::Symbol (default: :time): the name of the column with the time corresponding to the row. Times have to be non-decreasing per ID. Several dose records may share a time. They take effect one after another in the order in which they appear in the data set, and each one acts on the compartment amounts the previous one left, so the row order decides the outcome.

• evid::Union{Symbol,Nothing} (default: nothing): the name of the column with event IDs, or nothing. Possible event IDs are:

  • 0: observation
  • 1: dose event
  • 2: other type event
  • 3: reset event (the amounts in each compartment are reset to zero and the on/off status of each compartment is reset to its initial status)
  • 4: reset and dose event

  The event ID defaults to 0 if the dose amount is 0 or missing, and 1 otherwise.

• amt::Symbol (default: :amt): the name of the column of dose amounts. If the event ID is specified and non-zero, the dose amount should be non-zero. The default dose amount is 0.

• addl::Symbol (default: :addl): the name of the column that indicates the number of repeated dose events. The number of additional doses defaults to 0.

• ii (default: :ii): the name of the column of inter-dose intervals. When the number of additional doses is specified and non-zero, this is the time to the next dose. For steady-state events with multiple infusions or bolus doses, this is the time between implied doses. The default inter-dose interval is 0. It is required to be non-zero for steady-state events with multiple infusions or bolus doses, and it is required to be zero for steady-state events with constant infusion.

• cmt::Symbol (default: :cmt): the name of the column with the compartment to be dosed. Compartments can be specified by integers, strings and Symbols. The default compartment is 1.

• rate::Symbol (default: :rate): the name of the column with the rate of administration. A rate=-2 allows the rate to be determined by Dose Control Parameters (DCP). Defaults to 0. Possible values are:

  • 0: instantaneous bolus dose
  • >0: infusion dose administered at a constant rate for a duration equal to amt/rate
  • -2: infusion rate or duration specified by the dose control parameters (see @dosecontrol)

• ss::Symbol (default: :ss): the name of the column that indicates whether a dose is a steady-state dose. A steady-state dose is defined as the result of having applied the same dose from the infinite past. Possible values of the steady-state indicator are:

  • 0: dose is not a steady state dose.
  • 1: dose is a steady state dose, and the compartment amounts are to be reset to the steady-state amounts resulting from the given dose. Compartment amounts resulting from prior dose event records are "zeroed out", and infusions in progress or pending additional doses are cancelled.
  • 2: dose is a steady state dose and the compartment amounts are to be set to the sum of the steady-state amounts resulting from the given dose plus the amounts at the event time were the steady-state dose not given.

  The default value is 0.

• route: deprecated and ignored. Route is no longer part of the dose representation; setting it warns. To keep route information, pass it as a covariate.

• mdv::Union{Symbol,Nothing} (default: nothing): the name of the column that indicates if observations are missing, or nothing.

• event_data::Bool (default: true): toggles assertions applicable to event data. More specifically, checks if the following columns are present in the DataFrame, either as the default values or as user-defined values: id, time, and amt. If no evid column is present, then a warning will be thrown, and evid is set to 1 when amt values are >0 or not missing, or evid is set to 0 when amt values are missing and observations are not missing. Otherwise, read_pumas will throw an error. If event_data = false, it is only checked that id and time columns are present, but the time does not have to be monotonically increasing.

• covariates_direction::Symbol (default: :left): the direction of covariate interpolation. Either :left (Last Observation Carried Forward, LOCF) (default), or :right (Next Observation Carried Backward, NOCB). Notice, that for models with occasion variables it is important to use :left for the correct behavior of the interpolation.

• check::Bool (default: true): toggles NMTRAN compliance check of the input data for event_data == false the only check is that there is an id column. More specifically, checks if the following columns are present in the DataFrame, either as the default values or as user-defined values: :id, :time, :amt, :cmt, :evid, :addl, :ii, and :ss. Additional checks are:

  • all variables in observations are numeric, i.e. Integer or AbstractFloat.
  • :amt column is numeric, i.e. Integer or AbstractFloat.
  • :cmt column is either a positive Integer, an AbstractString, or a Symbol.
  • :amt values must be missing or 0 when evid = 0; or >=0 otherwise.
  • all variables in observation are missing when evid = 1.
  • :ii column must be present if :ss is specified or present.
  • :ii values must be missing or 0 when evid = 0.
  • :ii column must be >0 if :addl values are >0, and vice-versa.
  • :addl column, if present, must be >=0 when evid = 1.
  • :evid values must be !=0 when :amt values are >0, or :addl and :ii values are >0.

• adjust_evid34::Bool (default: true): toggles adjustment of time vector for reset events (evid = 3 and evid = 4). If true (the default) then the time of the previous event is added to the time on record to ensure that the time vector is monotonically increasing.

If a column does not exist, its values are imputed to be the default values.

sample_params(
    mcmc::BayesMCMCResults;
    samples::Int,
    rng::AbstractRNG = default_rng(),
    replace::Bool = true,
)

Draw samples sets of population-level parameters from the mcmc results.

By default (replace = true), sampling is performed with replacement. If replace is set to false, parameters are sampled without replacement.

The random number generator used for sampling can be set with the rng keyword argument.

sample_params(pmi::FittedPumasModelInference; rng = Random.default_rng(), resample::Bool = false, samples::Int)

Create a parameter set populated with values sampled from the parameter uncertainty.

The mandatory samples keyword controls how many samples are generated. The rng keyword controls the pseudo-random number generator used to generate the samples. When resample is false, vcov(pmi) is calculated and used to sample parameter values from a multivariate normal distribution with mean coef(pmi) and covariance matrix vcov(pmi). When resample is true, random values will be drawn from the sampled coefficients of the infer procedure. This is only possible when the inference procedure was performed with the Bootstrap or SimulatedInference algorithms. If the inference was calculated as the asymptotic covariance matrix, the case resample = true will throw an error.

sample_params(model::PumasModel; rng = Pumas.default_rng())

Create a parameter set populated with values for the parameters of a PumasModel. Random values will be drawn from the prior distributions for parameters that have a prior. Otherwise, a deterministicitic value will be assigned to the parameter similar to init_params.

sample_randeffs(
  [rng::AbstractRNG=Random.default_rng(),]
  model::Union{PumasModel,PumasEMModel},
  [subject_or_population::Union{AbstractSubject,AbstractPopulation},]
  param::NamedTuple,
)

Generate a random set of random effects for model model using parameters param.

Optionally, a random number generator object rng can be passed as the first argument. The optional argument subject_or_population takes a subject or population for which the random effects are generated. If a population of subjects is provided, a vector of random effects is returned with one element for each subject within the population.

sample_randeffs(
    [rng::AbstractRNG = default_rng(),]
    mcmc::BayesMCMCResults,
    samples::Int;
    subject::Union{Int,Nothing} = nothing,
    replace::Bool = true,
)

Draw samples sets of subject-level parameters from the mcmc results.

By default (replace = true), sampling is performed with replacement. If replace is set to false, parameters are sampled without replacement.

By default (subject = nothing), subject-level parameters are returned for all subjects that were used to obtain the mcmc results, concatenated in a single vector. If subject is set to an integer, only the subject-level parameters of the subjectth subject are returned.

The random number generator used for sampling can be set with the rng argument.

samples_table(infer_result) -> DataFrame

Return a DataFrame of parameter samples from a non-asymptotic inference result returned by infer, with one row per sample and one column per parameter. Supported for Bootstrap, SIR, and MarginalMCMC results.

For Bootstrap results, an additional loglikelihood column reports the log-likelihood of the bootstrapped parameter estimates.

For SIR results, the returned data frame contains the SIR proposal samples across all iterations. Columns iteration and resampled specify the iteration number of each proposal sample and whether the proposal sample was resampled. The additional loglikelihood column reports the log-likelihood of the samples.

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

Simulate random observations from model for population with parameters param at obstimes (by default, use the times of the existing observations for each subject). If no randeffs is provided, then random ones are generated according to the distribution in the model.

Arguments

• model may either be a PumasModel or a PumasEMModel.

• population may either be a Population of Subjects or a single Subject.

• param should be either a single parameter set, in the form of a NamedTuple, or a vector of such parameter sets. If a vector then each of the parameter sets in that vector will be applied in turn. Example: (; tvCL=1., tvKa=2., tvV=1.)

• randeffs is an optional argument that, if used, should specify the random effects for each subject in population. Such random effects are specified by NamedTuples for PumasModels (e.g. (; tvCL=1., tvKa=2.)) and by Tuples for PumasEMModels (e.g. (1., 2.)). If population is a single Subject (without being enclosed in a vector) then a single such random effect specifier should be passed. If, however, population is a Population of multiple Subjects then randeffs should be a vector containing one such specifier for each Subject. The functions center_randeffs, zero_randeffs, and sample_randeffs are sometimes convenient for generating randeffs:

randeffs = zero_randeffs(model, population, param)
solve(model, population, param, randeffs)

If no randeffs is provided, then random ones are generated according to the distribution in the model.

• obstimes is a keyword argument where you can pass a vector of times at which to simulate observations. The default, nothing, ensures the use of the existing observation times for each Subject.

• ensemblealg is a keyword argument that allows you to choose between different modes of parallelization. Supported options are:

  • EnsembleSerial(): sequential evaluation.
  • EnsembleThreads(): multithreaded evaluation.
  • EnsembleDistributed(): sequential evaluation within multiple worker processes.
  • EnsembleSplitThreads(): multithreaded evaluation within multiple worker processes.

• diffeq_options is a keyword argument that takes a NamedTuple of options to pass on to the differential equation solver. The absolute tolerance, relative tolerance, and the maximum number of iterations of the numerical solvers for computing steady states can be specified with the ss_abstol, ss_reltol, and ss_maxiters options. By default, ss_abstol and ss_reltol are set to the absolute and relative tolerances of the differential equation solver, and ss_maxiters is set to 1000.

• rng is a keyword argument where you can specify which random number generator to use.

• simulate_error is a keyword argument that allows you to choose whether a sample (true) or the mean (false) of the predictive distribution's error model will be returned for each subject.

simobs(pmi::FittedPumasModelInference, population::Population = pmi.fpm.data, randeffs::Union{Nothing, AbstractVector{<:NamedTuple}} = nothing; samples::Int, rng = default_rng(), kwargs...,)

Simulates observations from the fitted model using a truncated multivariate normal distribution for the parameter values when a covariance estimate is available and using non-parametric sampling for simulated inference and bootstrapping. When a covariance estimate is available, the optimal parameter values are used as the mean and the variance-covariance estimate is used as the covariance matrix. Rejection sampling is used to avoid parameter values outside the parameter domains. Each sample uses a different parameter value. samples is the number of samples to sample. population is the population of subjects which defaults to the population associated the fitted model. randeffs can be set to a vector of named tuples, one for each sample. If randeffs is not specified (the default behaviour), it will be sampled from its distribution.

simobs(fpm::FittedPumasModel, [population::Population,] vcov::AbstractMatrix, randeffs::Union{Nothing, AbstractVector{<:NamedTuple}} = nothing; samples::Int, rng = default_rng(), kwargs...)

Simulates observations from the fitted model using a truncated multi-variate normal distribution for the parameter values. The optimal parameter values are used for the mean and the user supplied variance-covariance (vcov) is used as the covariance matrix. Rejection sampling is used to avoid parameter values outside the parameter domains. Each sample uses a different parameter value. samples is the number of samples to sample. population is the population of subjects which defaults to the population associated the fitted model so it's optional to pass. randeffs can be set to a vector of named tuples, one for each sample. If randeffs is not specified (the default behaviour), it will be sampled from its distribution.

The rows and columns of vcov are the estimated, i.e. non-constant, parameter entries in the order of the rows of coeftable(fpm), which is the layout that vcov returns.

simobs(model::PumasModel, population::Population; samples::Int = 10, simulate_error = true, rng::AbstractRNG = default_rng())

Simulates model predictions for population using the prior distributions for the population and subject-specific parameters. samples is the number of simulated predictions returned per subject. If simulate_error is false, the mean of the predictive distribution's error model will be returned, otherwise a sample from the predictive distribution's error model will be returned for each subject. rng is the pseudo-random number generator used in the sampling.

simobs(model::PumasModel, subject::Subject; samples::Int = 10, simulate_error = true, rng::AbstractRNG = default_rng())

Simulates model predictions for subject using the prior distributions for the population and subject-specific parameters. samples is the number of simulated predictions returned. If simulate_error is false, the mean of the predictive distribution's error model will be returned, otherwise a sample from the predictive distribution's error model will be returned for each subject. rng is the pseudo-random number generator used in the sampling.

simobstte(
  model::PumasModel,
  subject::Union{Subject,Population},
  param::NamedTuple,
  randeffs::Union{Vector{<:NamedTuple}, NamedTuple, Nothing}=nothing;
  minT::Real=0.0,
  maxT::Real,
  nT::Integer=10,
  repeated::Bool=false,
  rng::AbstractRNG = default_rng())

Simulate observations from a time-to-event model, i.e. one with a TimeToEvent dependent variable and return either a new Subject or Population with random event time stamps in the time vector. The observations vector contains the censoring information (see TimeToEvent for more details on the data coding).

The simobstte function first computes nT values of the survival probability from t=minT to maxT and then interpolates with a cubic spline to get a smooth survival function. The user must supply maxT. Given a survival function, it is possible to simulate from the distribution by using inverse cdf sampling. Instead of sampling a uniform variate to be compared with the survival function, we sample an exponential variate and compare with cumulative hazard. The two approaches are equivalent. The Roots package is then used for computing the root. If repeated=true then new event times are drawn until maxT has been reached.

successful_fits(vfpm::Vector{<:AbstractFittedPumasModel)
successful_fits(fpm::FittedPumasModelInference{<:Any, <:Bootstraps})
successful_fits(fpm::Bootstraps)

Returns a vector of all successful fits stored in a vector of fits or the results of Bootstrap inference. Can be compared to the input number of samples or length of vector to calculate the fraction of successful fits.

tad(t, events) -> time after most recent dose

Converts absolute time t (scalar or array) to relative time after the most recent dose. If t is earlier than all dose times, then the (negative) difference between t and the first dose time is returned instead. If there are no dose events, missing is returned.

truncated_latent(
    distribution::UnivariateDistribution,
    lower::Union{Real,Nothing},
    upper::Union{Real,Nothing},
)
truncated_latent(
    distribution::UnivariateDistribution;
    lower::Union{Real,Nothing} = nothing,
    upper::Union{Real,Nothing} = nothing,
)

Construct a truncated distribution based on distribution and the truncation limits lower and upper. When fitting a model with a truncated_latent dependent variable, the log probability of the truncated distribution is used in the likelihood calculation. However, when calculating the mean or other summary statistics, the underlying untruncated distribution is used.

Examples

When truncating a distribution it is necessary to provide lower and upper bounds. Keyword arguments can be used to specify the bounds. The default is to use nothing for both lower and upper bounds to indicate no truncation at either end. The keyword for can be used whenever the bound is fixed for all subjects and observations. If the cp is the modelled concentration and σ is the standard deviation of the error, then the following code will create a truncated normal distribution with a lower bound of 1e-3 similar to how the functions would be used in a @derived block.

julia> cp = [10.0, 4.0, 1.0]
3-element Vector{Float64}:
 10.0
  4.0
  1.0

julia> σ = 0.35
0.35

julia> @. truncated_latent(Normal(cp, σ); lower = 1e-3);

This can be used to set a lower truncation point, and upper truncation point (with the upper keyword) or both. However, if the truncation point depends on a covariate we cannot broadcast (@.) if the keyword is a vector. Then we can use the first version mentioned at the top where limits are given af positional arguments.

julia> cp = [10.0, 4.0, 1.0]
3-element Vector{Float64}:
 10.0
  4.0
  1.0

julia> lower = [1e-3, 1e-3, 1e-3]
3-element Vector{Float64}:
 0.001
 0.001
 0.001

julia> σ = 0.35
0.35

julia> @. truncated_latent(Normal(cp, σ), lower, nothing);

Similarly, the second positional argument can be nothing and the third argument can be a vector or all three arguments can be vectors if both lower and upper limits are relevant.

vpc(simulations::AbstractVector{<:Pumas.SimulatedPopulation}; kwargs...)

Compute the quantiles for VPC for a Pumas model with simulated confidence intervals around the empirical quantiles based on simulated populations simulations, typically obtained with simobs.

The function supports the same keyword arguments as vpc(::Pumas.AbstractFittedPumasModel; kwargs...).

vpc(simulations::Pumas.SimulatedPopulation; kwargs...)

Compute the quantiles for VPC for a Pumas model with simulated confidence intervals around the empirical quantiles based on a simulated population simulations, typically obtained with simobs.

The function supports the same keyword arguments as vpc(::Pumas.AbstractFittedPumasModel; kwargs...).

vpc(
  fpm::AbstractFittedPumasModel;
  samples::Integer = 499
  qreg_method = IP(),
  observations::Vector{Symbol} = [keys(fpm.data[1].observations)[1]],
  stratify_by::Vector{Symbol} = Symbol[],
  quantiles::NTuple{3,Float64}=(0.1, 0.5, 0.9),
  level::Real=0.95,
  ensemblealg=EnsembleThreads(),
  bandwidth=nothing,
  maxnumstrats=[6 for i in 1:length(stratify_by)],
  covariates::Vector{Symbol} = [:time],
  smooth::Bool = true,
  rng::AbstractRNG=default_rng(),
  nodes::Dict = Dict{Tuple{},Vector{Float64}}(),
  nnodes::Integer = 0,
  prediction_correction = false,
)

Compute the quantiles for VPC for a FittedPumasModel or FittedPumasEMModel with simulated confidence intervals around the empirical quantiles based on samples simulated populations.

The following keyword arguments are supported:

• samples: The number of simulated populations to generate, defaults to 499
• quantiles::NTuple{3,Float64}: A three-tuple of the quantiles for which the quantiles will be computed. The default is (0.1, 0.5, 0.9) which computes the 10th, 50th and 90th percentile.
• level::Real: Probability level to use for the simulated confidence intervals. The default is 0.95.
• observations::Vector{Symbol}: The name of the dependent variable to use for the VPCs. The default is the first dependent variable in the Population.
• stratify_by: The covariates to be used for stratification. Takes an Vector of the Symbols of the stratification covariates.
• ensemblealg: This is the parallelization choice.
• bandwidth: The kernel bandwidth in the quantile regression. If you are seeing NaNs or an error, increasing the bandwidth should help in most cases. With higher values of the bandwidth you will get more smoothened plots of the quantiles so it's a good idea to check with your data to pick the right bandwidth. For a ContinuousVPC, the default value is 2.0. For a DiscreteVPC, the default value is nothing to indicate automatic bandwidth tuning. For a CensoredVPC, the bandwidth should be given as a tuple where the first value applies to the ContinuousVPC and the second to the DiscreteVPC for the censoring variable. For example, (4.0, nothing) applies a bandwidth of 4.0 to the ContinuousVPC and automatically tunes the DiscreteVPC. Automatic tuning is only applicable to DiscreteVPCs at this moment.
• maxnumstrats: The maximum number of strata for each covariate. Takes an Vector with the number of strata for the corresponding covariate, passed in stratify_by. It defaults to 6 for each of the covariates.
• covariates: The independent variable for VPC, defaults to [:time].
• smooth: In case of discrete VPC used to determine whether to interpolate the dependent variable if independent variable is continuous, defaults to true.
• rng: A random number generator, uses the default_rng from Random as default.
• nodes: A Dict specifying the local regression nodes. If stratify_by argument has been specified, the dict should have tuples of values as keys corresponding to each stratum. Otherwise, the key should be an empty tuple. The values of the Dict should be the nodes of the local regressions. If the nodes argument is not passed, at most nnodes quantiles is used from each stratum.
• nnodes: An integer specifying the number of points to use for the local regressions for each stratum. The default value is 11 but the actual number of nodes will not exceed the number of unique covariate values in the stratum.
• qreg_method: Defaults to IP(). For most users the method used in quantile regression is not going to be of concern, but if you see large run times switching qreg_method to IP(true) should help in improving the performance with a tradeoff in the accuracy of the fitting.
• prediction_correction: Default to false. Set to true to enable prediction correction that multiplies all observations with the ratio between the median population prediction (from a quantile regression on the binning covariate) and each individual's population prediction.

wresiduals(fpm::AbstractFittedPumasModel, approx::LikelihoodApproximation; nsim=nothing)

Calculate the individual and population weighted residual.

Takes a fit result and an approximation method for the marginal likelihood calculation which defaults to the method used in the fit as positional arguments.

Using the FO approximation method corresponds to the Weighted Residuals (WRES) and FOCE corresponds to Conditional Weighted Residuals (CWRES).

If nsim is specified only the simulation based Expected Individual Weighted Residuals (EIWRES) are included in the output as individual residuals and population residuals are not computed.

The output is a SubjectResidual object that stores the population (wres) and individual (iwres) residuals along with the subject and approximation method (approx).

The computation is parallelized across subjects.

zero_randeffs(
  model::Union{PumasModel,PumasEMModel},
  [subject_or_population::Union{AbstractSubject,AbstractPopulation},]
  param::NamedTuple,
)

Return a set of zero random effects for model model using parameters param.

The optional argument subject_or_population takes a subject or population for which the random effects are returned. If a population of subjects is provided, a vector of random effects is returned with one element for each subject within the population.

ηshrinkage(
    fpm::AbstractFittedPumasModel;
    stratify_by::Union{Nothing,Vector{Symbol}} = nothing,
)

Calculate the η-shrinkage values of the fitted Pumas model fpm, and return them as a named tuple with shrinkage values for each random effect.

A vector of names of constant covariates can be passed to stratify_by to compute η-shrinkage values separately for each stratum of the population defined by the values of the specified covariates. By default (stratify_by = nothing), no stratification is performed.

Details

The shrinkage s of a normally distributed random effect η ~ Normal(μ, ω) is computed as the relative difference of the standard deviation of the empirical Bayes estimates (EBEs) ηᴱᴮᴱ of η and the standard deviation ω of the distribution of η:

s = 1 - std(ηᴱᴮᴱ) / ω

For multivariate normally distributed random effects η ~ MvNormal(μ, Ω), shrinkage is computed element-wise:

s = 1 .- std(ηᴱᴮᴱ) ./ sqrt.(diag(Ω))

ϵshrinkage(fpm::Union{FittedPumasModel,FittedPumasEMModel})

Calculate the ϵ-shrinkage.

The function returns a named tuple with shrinkage values for each derived variable with Gaussian error model. For derived variables with non-Gaussian error models, no shrinkage value is computed.

aic(fpm::AbstractFittedPumasModel)

Calculate the Akaike information criterion (AIC) of the fitted Pumas model fpm.

bic(fpm::AbstractFittedPumasModel)

Calculate the Bayesian information criterion (BIC) of the fitted Pumas model fpm.

coeftable(pmi::FittedPumasModelInference; level = pmi.level) -> DataFrame

Construct a DataFrame of parameter names, descriptions, estimates, standard error, and confidence interval with confidence level level. The default value of level is the one that was passed in infer. The description of each parameter comes from the model's @metadata block, and falls back to the parameter name when the model gives none.

coeftable(fpm::FittedPumasModel) -> DataFrame

Construct a DataFrame of parameter names, descriptions, and estimates from fpm. The description of each parameter comes from the model's @metadata block, and falls back to the parameter name when the model gives none.

coeftable(vfpm::Vector{<:FittedPumasModel}) -> DataFrame

Construct a DataFrame of parameter names, descriptions, and estimates and their standard deviation from vector of fitted single-subject models vfpm. The description of each parameter comes from the model's @metadata block, and falls back to the parameter name when the model gives none.

fit(
    fpm;
    reset_optim::Bool = false,
    verbose::Bool = true,
    checkidentification::Bool = true,
    ignore_numerical_error::Bool = true,
    [diffeq_options::NamedTuple,]
    [ensemblealg::SciMLBase.EnsembleAlgorithm,]
)

Start a new fit where an old one left off.

Both the fitted parameters and the fitted random effects from the fitted Pumas model fpm are used as starting values for the new fit.

The state of the optimization algorithm can be reset by setting reset_optim to true.

The other supported keyword arguments allow further customization of the fit in the same way as fit for non-fitted Pumas model. By default, the fit uses the differential equation solver settings diffeq_options and the parallelization ensemblealg of the optimization that were used for fitting fpm.

fit(model::PumasModel, data::Population, param::NamedTuple, alg::BayesMCMC; kwargs...)

Samples from the posterior distribution of the model's population and subject-specific parameters given the observations in data using the Hamiltonian Monte Carlo (HMC) based No-U-Turn Sampler (NUTS) algorithm. The MCMC sampler is initialized starting from param for the population parameters. The subject-specific parameters are also appropriately initialized. The Bayesian inference algorithm's options are passed in the alg struct. Run ?BayesMCMC for more details on the options that can be set.

Keyword arguments

• constantcoef (default: ()): A Tuple of Symbols of the parameters to keep constant during sampling, e.g. constantcoef = (:σ,). Their values are read from param. This can be used to sample from conditional posteriors with some parameters fixed to specific values.
• diffeq_options (default: (;)): A NamedTuple of options of the differential equation solver, e.g. (reltol = 1e-8, abstol = 1e-8).
• ensemblealg (default: EnsembleThreads()): The parallelization strategy. Can be set to EnsembleSerial() for serial sampling, EnsembleThreads() for multi-threaded sampling, EnsembleDistributed() for multi-processing sampling, or EnsembleSplitThreads() for multi-processing over chains and multi-threading over subjects.
• parallel_chains (default: true if enough threads/workers exist): If set to true, the chains will be sampled in parallel using either multi-threading or multi-processing depending on the value of ensemblealg.
• parallel_subjects (default: true if enough threads/workers exist): If set to true, the log probability computation will be parallelized over the subjects using either multi-threading or multi-processing depending on the value of ensemblealg.
• rng (default: Random.default_rng()): The random number generator used for sampling.
• verbose (default: true in interactive sessions, false otherwise): Controls printing of HMC sampling diagnostics (acceptance rate, tree depth, step size, divergent-transition flag, log posterior density) to the terminal.
• ignore_numerical_error (default: true): If set to false, numerical errors during sampling are raised. This is useful when debugging model issues.

fit(
    model::PumasModel,
    population::Population,
    param::NamedTuple,
    approx::Union{LikelihoodApproximation, MAP};
    init_randeffs::Union{Nothing,NamedTuple,Vector{<:NamedTuple}} = nothing,
    optim_alg = Optim.BFGS(linesearch = LineSearches.HagerZhang(; rho = 2.0), initial_stepnorm = 1.0),
    optim_options::NamedTuple = (;),
    optimize_fn = nothing,
    constantcoef::Tuple = (),
    ensemblealg::SciMLBase.EnsembleAlgorithm = EnsembleThreads(),
    diffeq_options::NamedTuple = (;),
    checkidentification::Bool = true,
    verbose::Bool = true,
    ignore_numerical_error::Bool = true,
  )

Fit the Pumas model model to the dataset population with starting values param using the estimation method approx. Currently supported values for the approx argument are FO, FOCE, LaplaceI, NaivePooled, and BayesMCMC. See the online documentation for more details about the different methods.

To control the optimization procedure for finding population parameters (fixed effects), use the optim_alg and optim_options keyword arguments. In previous versions of Pumas the argument optimize_fn was used, but is now discouraged and will be removed in a later version of Pumas. These options control the optimization all approx methods except BayesMCMC. The default optimization function uses the quasi-Newton routine BFGS method from the Optim package. It can be changed by setting the optim_alg to an algorithm implemented in Optim.jl as long as it does not use second order derivatives. Optimization specific options can be passed in using the optim_options keyword and has to be a NamedTuple with names and values that match the Optim.Options type. For example, the optimization trace can be disabled and the algorithm can be changed to L-BFGS by passing optim_alg=Optim.LBFGS(), optim_options = (; show_trace=false) to fit. See Optim for more details.

It is possible to fix one or more parameters of the fit by passing a Tuple of Symbols as the constantcoef argument with elements corresponding to the names of the fixed parameters, e.g. constantcoef=(:σ,).

When models include an @random block and fitting with NaivePooled is requested, it is required that the user sets all variability parameters to zero with constantcoef such that these can be ignored in the optimization, e.g. constantcoef = (:Ω,) while overwriting the corresponding values in param with (; init_params..., Ω = zeros(3, 3)).

Starting values of the random effects can be provided with init_randeffs. If init_randeffs is set to nothing (the default), then the random effects are initialized with center_randeffs(model, population, param) as central points of the subjects' random effect distributions for the initial parameter values param.

If init_randeffs is a NamedTuple, then these values are used as initial random effects for all subjects. Alternatively, init_randeffs can be set to a vector of NamedTuples with starting values for each subject.

Parallelization of the optimization is supported for most estimation methods via the ensemble interface of DifferentialEquations.jl. Currently, the supported options are:

• EnsembleThreads() (default): fit by using multiple threads.
• EnsembleSerial(): fit by using a single thread.
• EnsembleDistributed(): fit by using a single thread in multiple worker processes.
• EnsembleSplitThreads(): fit by using multiple threads in multiple worker processes.

The fit function will check if any gradients and throw an exception if any of the elements are exactly zero unless checkidentification is set to false.

Further keyword arguments can be passed via the diffeq_options argument. This allows for passing arguments to the differential equations solver such as alg, abstol, and reltol. The default values for these are AutoVern7(Rodas5P(autodiff=true)), 1e-12, and 1e-8 respectively. See the DifferentialEquations.jl documentation for more details.

The keyword verbose controls if info statements about initial evaluations of the loglikelihood function and gradient should be printed or not. Defaults to true.

Since the numerical optimization routine can sometimes visit extreme regions of the parameter space during it's exploration we automatically handle the situation where the input parameters result in errors when solving the dynamical system. This allows the algorithm to recover and continue in many situations that would have otherwise stalled early. Sometimes, it is useful to turn this error handling off when debugging a model fit, and this can be done by setting ignore_numerical_error = false.

fit(
  model::PumasEMModel,
  population::Population,
  param::NamedTuple,
  approx::Union{SAEM,LaplaceI};
  diffeq_options::NamedTuple = NamedTuple(),
  ensemblealg = EnsembleThreads(),
  rng = default_rng()
)

Fit the PumasEMModel to the dataset population using the initial values specified in param. The supported methods for the approx argument are currently SAEM and LaplaceI. See the online documentation for more details about these two methods.

When fitting with SAEM the variance terms of the random effects (Ω) and the dispersion parameters for the error model (σ) are initialized to the identity matrix or 1 as appropriate. They may also be specified. With SAEM, it is recommended to choose an init larger than the true values to facilitate exploration of the parameter space and avoid getting trapped in local optima early. Currently LaplaceI fits with a PumasEMModel require Ω to be diagonal.

Options for ensemblealg are EnsembleSerial() and EnsembleThreads() (the default). The fit will be parallel if ensemblealg == EnsembleThreads().

Log-Likelihood

Since SAEM is an expectation-maximization algorithm, the (log-)likelihood is not calculated during fitting. To approximate the log-likelihood of a fitted model, use likelihood approximation methods such as FOCE or LaplaceI:

loglikelihood(fpm, FOCE())
loglikelihood(fpm, LaplaceI())

loglikelihood(fpm::AbstractFittedPumasModel)

Compute the log-likelihood of a fitted Pumas model.

loglikelihood(
    model::AbstractPumasModel,
    population::AbstractPopulation,
    param::NamedTuple,
    approx::Union{LikelihoodApproximation, MAP};
    diffeq_options::NamedTuple = (;),
    ensemblealg::SciMLBase.EnsembleAlgorithm = EnsembleThreads(),
)

Compute the log-likelihood of a population for a Pumas model with parameters param, using the estimation method approx.

Keyword Arguments

• diffeq_options: Options that are passed to the differential equation solver.

• ensemblealg: Mode of parallelization across subjects. Supported options are:

  • EnsembleSerial(): sequential evaluation.
  • EnsembleThreads(): multithreaded evaluation.
  • EnsembleDistributed(): sequential evaluation within multiple worker processes.
  • EnsembleSplitThreads(): multithreaded evaluation within multiple worker processes.

loglikelihood(
    model::AbstractPumasModel,
    subject::AbstractSubject,
    param::NamedTuple,
    approx::Union{LikelihoodApproximation, MAP};
    diffeq_options::NamedTuple = (;),
)

Compute the log-likelihood of a subject for a Pumas model with parameters param, using the estimation method approx.

Keyword Arguments

• diffeq_options: Options that are passed to the differential equation solver.

loglikelihood(b::BayesMCMCResults; split_method, split_by)

Computes a 3D array of the in-sample loglikelihoods of all the data points/subsets using the output of the MCMC fit function. The output is a 3D array of size (ndata, nsamples, nchains) where ndata is the number of data points, nsamples is the number of MCMC samples per chain and nchains is the number of chains.

The way by which the data is split into points/subsets is determined using the split_method and the split_by keyword arguments. The split_method argument can be of any of the following types:

• LeaveK for a leave-K-out splitting
• KFold for a K-fold splitting
• LeaveFutureK for leaving K future points at a time while split_by can be of any of the following types:
• BySubject for leaving out subjects
• ByObservation for leaving out individual observations per subject

Please refer to the documentation of each of the above types for the parameters of each split method, e.g. by typing ?LeaveK or ?BySubject.

Examples:

Assume there are 5 subjects and 10 observations per subject. The following are some of the combinations in which the above inputs can be used:

• Leave-1-observation-out splitting, splitting 1 observation for all the subjects at a time. allsubjects = true means that the same observation index is removed for all the subjects, e.g. the 10th observation for all the subjects is used as a single subset, then the 9th observation is used for all the subjects is used as another subset, etc.

      split_method = LeaveK(K = 1)
      split_by = ByObservation(allsubjects = true)
      loglikelihood(b; split_method = split_method, split_by = split_by)

• Leave-2-future-observations-out splitting, splitting 2 future observations per subject such that no less than 4 past observations are remaining. allsubjects = false means that only the data of one subject at a time will be split. So observations 5 and 6 of subject 1 are the first data subset. The second subset is observations 7 and 8 of subject 1. The third subset is observations 9 and 10 of subject 1, etc. Then in the fourth to sixth runs, subject 2's data gets split. In the seventh to ninth runs, subject 3's data gets split, etc.

      split_method = LeaveFutureK(K = 2, minimum = 4)
      split_by = ByObservation(allsubjects = false)
      loglikelihood(b; split_method = split_method, split_by = split_by)

• Leave-1-subject-out splitting, marginalizing the subject-specific parameters using FOCE when computing the in-sample likelihoods. Using this method, each subject is a data point. The in-sample likelihoods computed are the marginal likelihoods of the subjects approximated using FOCE. FO and LaplaceI can also be used instead.

      split_method = LeaveK(K = 1)
      split_by = BySubject(marginal = FOCE())
      loglikelihood(b; split_method = split_method, split_by = split_by)

• Leave-1-future-subject-out splitting, where the in-sample likelihoods are conditional likelihoods computed using 0s for the subject-specific parameters instead of marginalizing. In this method, each subject after the third one is considered a data point, i.e. the 4th and 5th subjects are the only 2 points for which the in-sample likelihood is computed.

      split_method = LeaveFutureK(K = 1, minimum = 3)
      split_by = BySubject(marginal = nothing)
      loglikelihood(b; split_method = split_method, split_by = split_by)

loglikelihood(r::BayesianCVResult)

Computes the out-of-sample loglikelihoods from the result r of crossvalidate. The output is vector of matrices, e.g. loglik where loglik[i] is the loglikelihoods matrix of data point/set i which was not used in training. loglik[i][j, k] is the loglikelihood using the parameter values from the jth sample of the kth chain.

loglikelihood(fitted_model, approx::LikelihoodApproximation)

Compute the log-likelihood of a fitted Pumas model using the likelihood approximation method approx.

predict(
  model::AbstractPumasModel,
  population::Population,
  param::NamedTuple,
  [vrandeffs,];
  [obstimes::AbstractVector,
  diffeq_options::NamedTuple]
)::Union{SubjectPrediction,Vector{SubjectPrediction}}

Compute population and individual predictions for either the single subject or vector of subjects (Population) population based on model and the population parameters param. Optionally, random effects can be passed as a Vector of NamedTuples where the ith element is used for the individual predictions of the ith subject.

If omitted, the empirical bayes estimates given the input param will be calculated if the input subjects have observations. If any subject does not have any observations, the mode of the domain-transformed prior will be used.

If the optional obstimes argument is passed then the time points in obstimes and the time points of the observations of the subject are used for the predictions. Otherwise, only the time points of the observations of the subject are used for the predictions.

The function allows for extra keyword arguments to be passed on to the differential equations solver through the diffeq_options keyword. See the online documentation for more details.

predict(
  model::AbstractPumasModel,
  subject::Subject,
  param::NamedTuple,
  [randeffs,];
  [obstimes::AbstractVector,
  diffeq_options::NamedTuple]
)::Union{SubjectPrediction,Vector{SubjectPrediction}}

Compute population and individual predictions for the single subject based on model and the population parameters param. A NamedTuple of random effects, randeffs, can be omitted or provided by the user. If they are omitted, they will be estimated from the data in the subject.

If the optional obstimes argument is passed then the time points in obstimes and the time points of the observations of the subject are used for the predictions. Otherwise, only the time points of the observations of the subject are used for the predictions.

The function allows for extra keyword arguments to be passed on to the differential equations solver through the diffeq_options keyword. See the online documentation for more details.

predict(
  fpm::AbstractFittedPumasModel,
  [population::Union{Subject,Population}];
  [obstimes::AbstractVector]
)::Union{SubjectPrediction,Vector{SubjectPrediction}}

Compute population and individual predictions for the fitted model fpm. By default, the predictions are computed for the estimation data but the predictions can also be computed for user supplied data by passing either a single subject or a vector of subjects (Population) as the population argument.

If the optional obstimes argument is passed then the time points in obstimes and the time points of the observations of the subject are used for the predictions. Otherwise, only the time points of the observations of the subject are used for the predictions.

Any optional keyword arguments used when fitting fpm are reused when computing the predictions.

pvalue(t::LikelihoodRatioTest)::Real

Compute the p-value of the likelihood ratio test t based on the asymptotic Χ²(k) distribution of the test statistic.

stderror(f::AbstractFittedPumasModel) -> NamedTuple

Compute the standard errors of the population parameters and return the result as a NamedTuple matching the NamedTuple of population parameters.

vcov(fpm; kwargs...)

Estimate the variance-covariance matrix for the parameters of the fitted model fpm.

Keyword Arguments

• sandwich_estimator::Bool (default: true): whether to use the sandwich estimator of the variance-covariance matrix.
• rethrow_error::Bool (default: false): whether an error is thrown if the estimation of the variance-covariance matrix fails, or the exception is returned without throwing it.
• ensemblealg (default: ensemblealg that was used during model fitting): the mode of parallelization. Possible options are: EnsembleSerial(), EnsembleThreads(), EnsembleDistributed() and EnsembleSplitThreads().

The matrix is singular if the entries of a parameter are linearly dependent: those of a SimplexDomain parameter sum to one and those of a ZeroSumDomain parameter sum to zero, so the row and column of the last entry are minus the sum of those of the others. The number of free parameters is dof, not the size of the matrix.

See also infer.

sample(
    [rng::AbstractRNG = default_rng(),]
    mcmc::BayesMCMCResults,
    samples::Int;
    subject::Union{Int,Nothing} = nothing,
    replace::Bool = true,
)

Draw samples sets of population-level and subject-level parameters from the mcmc results.

Samples are returned as a tuple of a vector of population-level parameters and a vector of subject-level parameters.

By default (replace = true), sampling is performed with replacement. If replace is set to false, parameters are sampled without replacement.

By default (subject = nothing), subject-level parameters are returned for all subjects that were used to obtain the mcmc results. If subject is set to an integer, only the subject-level parameters of the subjectth subject are returned.

The random number generator used for sampling can be set with the rng argument.

@covariance

Allows specifying a block-diagonal covariance structure for the Ω parameters in an @emmodel. The default is fully diagonal.

For example, to specify that the Ω is block diagonal with 1x1, 2x2, and 2x2 blocks:

@emmodel begin
  @random begin
    Ka ~ 1 | LogNormal
    CL ~ 1 | LogNormal
    Vc ~ 1 | LogNormal
    Q ~ 1 | LogNormal
    Vp ~ 1 | LogNormal
  end
  @covariance 1 2 2
end

This would mean that both CL and Vc may be correlated with each other, and Q and Vp may also be, but all other pairwise correlations are zero. @covariance 5 would result in a dense matrix.

Currently, estimation with LaplaceI is only supported for models with diagonal covariance matrices.

@covariates

Defines the model's covariates, e.g. wt. Must be used in an @model block. For example:

@model begin
  @covariates wt
end

@delay(
    distribution::Distribution,
    cmt,
) -> Real

Compute the distributed delay contribution for distribution for all dose amounts in administered to cmt up until the current time point. This macro should be used in the @pre block to compute an absorption rate explicitly added to the ODE in the @dynamics block.

Examples

Weibull delay

The model below implements the Weibull absorption model of Hartmann, D., Gysel, D., Dubach, U. C., & Forgo, I. (1990)

mdl_weibull = @model begin
  @param begin
    θCb ∈ RealDomain(lower = 0.0)
    θVc ∈ RealDomain(lower = 0.0)
    θVm ∈ RealDomain(lower = 10.0)
    θKm ∈ RealDomain(lower = 0.0)

    θA ∈ RealDomain(lower = 0.0)
    θB ∈ RealDomain(lower = 0.0)

    θf ∈ RealDomain(lower = 0.0)

    σ ∈ RealDomain(lower = 0.0)
  end

  @pre begin
    Baseline = θCb * θVc
    Vc = θVc
    VVm = Vc * θVm
    VKm = Vc * θKm

    λ = θA^(-1 / θB)
    B = θB

    f = θf

    dose_density = @delay(Weibull(B, λ), Central)
  end

  @dosecontrol begin
    bioav = (; Central = 0.0)
  end

  @init begin
    Central = Baseline
  end

  @dynamics begin
    Central' =
      VVm * (Baseline / (VKm + Baseline) - Central / (VKm + Central)) +
      f * Vc * dose_density
  end

  @derived begin
    μ := @. Central ./ Vc
    y ~ @. Gamma(1 / σ^2, μ * σ^2)
  end
end

Gamma delay

The model below uses the delay function to implement a gamma distributed delay demonstrated in Savic, Jonker, Kerbusch, and Karlsson (2007).

mdl_gamma_delay = @model begin
  @param begin
    θn ∈ RealDomain(lower = 0.0)
    θmat ∈ RealDomain(lower = 0.0)
    θCL ∈ RealDomain(lower = 0.0)
    θVc ∈ RealDomain(lower = 0.0)

    σ ∈ RealDomain(lower = 0.0)
  end

  @pre begin
    kₑ = θCL / θVc
    Vc = θVc
    dose_density = @delay(Gamma(θn, θmat / θn), Central)
  end

  @dosecontrol begin
    bioav = (; Central = 0.0)
  end

  @dynamics begin
    Central' = dose_density - kₑ * Central
  end

  @derived begin
    μ := Central ./ Vc
    y ~ @. Gamma(1 / σ^2, μ * σ^2)
  end
end

@derived

Defines the error model of the dependent variables. Must be used in an @model block. For example:

@model begin
  @derived begin
    dv ~ @. Normal(1000 * Central / Vc, 1000 * Central / Vc * σ)
  end
end

Variables defined with ~ must define an error model for each observation time point. Variables defined with = are not used for model fitting but only returned when simulating. If they are vectors, they must provide exactly one value for each observation time point.

@dosecontrol

Define dose control parameters as a function of fixed and random effects. Options include bioav, duration, lags, and rate. Must be used in an @model block. For example:

@model begin
  @dosecontrol begin
    bioav = (Depot1 = max(0, θ[5]), Depot2 = clamp(1 - θ[5], 0.0, 1.0), Central = 1)
    lags = (Depot1 = 0, Depot2 = max(0.1, θ[6]), Central = 0)
  end
end

@dynamics

Defines the dynamic system. Must be used in an @model block.

Supports three forms:

1. A system of differential equations:

@model begin
  @dynamics begin
    Depot1' = -ka * Depot1
    Depot2' = ka * Depot1 - ka * Depot2
    Central' = ka * Depot2 - CL / Vc * Central
  end
end

1. A pre-defined differential equation model:

@model begin
  @dynamics Depots1Central1
end

1. Chemical reactions using the Catalyst.jl DSL:

@model begin
  @dynamics begin
    ka, Depot1 --> Depot2
    ka, Depot2 --> Central
    CL / Vc, Central --> 0
  end
end

This is equivalent to the first example above.

@emmodel

Define a PumasEMModel. It may have the following blocks:

• @param and @random: Defines fixed and random effects, e.g.

@random begin
  CL ~ 1 + wt | LogNormal
  θbioav ~ 1 | LogitNormal
end

Distributions specify the parameter's support. Only Normal (-∞,∞), LogNormal (0,∞), and LogitNormal (0,1) are currently supported.

These define 1 unconstrained parameter per term in the formula, for example CL = exp(μ + wt * β + η) where μ, β ∈ (-∞,∞) and wt is a covariate constant with respect to time. η is a random effect. The @param block is equivalent, with the exception that there is no random effect η. These variables are accessible in @pre, @dosecontrol, @dynamics, @post, and @error blocks.

• @covariates: Covariates available in the @pre, @dosecontrol, @dynamics, and @post blocks.
• @pre: Block evaluated before the @dynamics.
• @dosecontrol: Block specifying dose control parameters. Function of @param and @random variables only.
• @init: Block specifying initial values of dynamical variables.
• @dynamics: Dynamics, equivalent to the functionality in the @model macro.
• @post: Block evaluated after the @dynamics and before the @error.
• @error: Block describing the error model. Dispersion parameters are implicit.

@error

Defines the error model of the dependent variables in a PumasEMModel. Distributions are parameterized by the mean only, the dispersion parameter(s), if any, are left implicit. See the documentation for details on the available error models. The mean parameter(s) must be (a) symbol(s). Use an earlier block, e.g. @post to define any necessary transformations.

Must be used in an @emmodel block. For example:

@emmodel begin
  @random begin
    CL ~ 1 | LogNormal
    Vc ~ 1 | LogNormal
  end
  @dynamics Central1
  @post begin
    cp = Central / Vc # cp is calculated here
  end
  @error begin
    dv ~ CombinedNormal(cp) # cp is the mean of the error distribution
  end
end

@init

Defines initial conditions of the dynamic system. Must be used in an @model block. For example:

@model begin
  @init begin
    Depot1 = 0.0
    Depot2 = 0.0
    Central = 0.0
  end
end

@metadata

Allows including additional details of a model as metadata. These take the form of key/value pairs defined with = as follows, where the values can take on <:Any value:

using Unitful

@model begin
  @metadata begin
    desc = "A short description of the model."
    timeu = u"hr"
  end
end

In the above example, a textual description of the model is provided by the desc key. It must be written in plaintext and should preferably not take up more than a single line. The timeu key describes the time unit used in the model. The example uses a Unitful unit, which requires loading Unitful explicitly, but a String such as "hr" works as well. At present desc and timeu are the only metadata keys that are used.

@model

Defines a Pumas model, with the following possible blocks:

• @param: defines the model's fixed effects and corresponding domains, e.g. tvcl ∈ RealDomain(lower = 0).
• @random: defines the model's random effects and corresponding distribution, e.g. η ~ MvNormal(Ω).
• @covariates: Names of covariates used in the model.
• @pre: pre-processes variables for the dynamic system and statistical specification.
• @dosecontrol: defines dose control parameters as a function of fixed and random effects.
• @vars: define variables usable in other blocks. It is recommended to define them in @pre instead if not a function of dynamic variables.
• @init: defines initial conditions of the dynamic system.
• @dynamics: defines the dynamic system.
• @derived: defines the statistical model of the dependent variables.
• @observed: lists model information to be stored in the solution.
• @options specifies model options, including checklinear, cse, inplace, and subject_t0.

@observed

Lists model information to be stored in the solution. Must be used in an @model block. For example:

@model begin
  @observed begin
    conc = @. 1000 * Central / Vc
  end
end

Vector-valued variables must provide exactly one value for each observation time point. Other variables, e.g. scalars, are repeated for every time point in the tabular output.

@options

Specifies model options, including the following:

• checklinear::Bool: determines whether the solver should check if the system defined in the @dynamics block is linear. If the system is linear, setting this option to true (default) enables calculating the solution through matrix exponentials. If it is set to false or time (t) appears in the @pre block, this optimization is disabled. This option can be useful when the matrix exponential solver is not superior to general numerical integrators or for debugging purposes.

• cse::Bool: determines whether common subexpressions of the system defined in the @dynamics block are eliminated when the right-hand side of the system is generated. By default, this option is set to true and every distinct subexpression is evaluated only once per evaluation of the right-hand side, which can improve performance for systems that use the same expensive subexpression in multiple equations. Setting it to false disables the transformation, e.g. to check how it affects the results of a model or the time it takes to construct it. Since the option changes the order in which floating point operations are performed, solutions of the system can differ within the tolerances of the solver. It only applies to systems that are solved with general numerical integrators; subexpressions shared by the entries of the coefficient matrix of a linear system are always recomputed.

• inplace::Bool: controls whether the solver uses mutable operations on Arrays (true), which can improve performance for large systems. By default, it is set to true for systems containing more than 5 dynamic variables, and set to false otherwise.

• subject_t0::Bool: determines the starting time point from which the solver integrates the system. By default, all systems are solved from time 0 to the last time point in the subject, making subject_t0 false. If subject_t0 is set to true, the solver will begin integration from the first time point found within the subject.

• symbolic_jacobian::Bool: determines whether the Jacobian of the system defined in the @dynamics block should be computed symbolically. By default, this option is set to false and the Jacobian is computed by automatic differentiation. This option only applies to systems that are solved with general numerical integrators.

• ss_states: selects the states for which a steady state is sought under steady-state dosing (ss = 1 and ss = 2). By default a steady state is sought for all states of the system, which is not possible if a state does not admit one, such as a compartment accumulating the AUC. The option takes either the states to include, as in ss_states = [Depot, Central], or the states to exclude, as in ss_states = Not([AUC]), using the state names of the @dynamics block. The remaining states are held at the values they had before the steady-state dose. Naming a state that does not exist, naming a state twice, and excluding every state are errors. The option requires a dynamical system and is not supported for the pre-defined analytical models, which have a closed-form steady-state solution, nor for mixed analytical/numerical models.

Must be used in an @model block. For example (defaults are shown below):

@model begin
  @options begin
    checklinear = true
    cse = true
    inplace = length(odevars) < 5
    subject_t0 = false
    symbolic_jacobian = false
  end
end

ss_states has no entry above since its default, a steady state for all states, is what applies when the option is not set.

@param

Defines the model's fixed effects and corresponding domains, e.g. tvcl ∈ RealDomain(lower = 0). Must be used in an @model block. For example:

@model begin
  @param begin
    tvcl ∈ RealDomain(lower = 0)
    tvv ∈ RealDomain(lower = 0)
    Ω ∈ PDiagDomain(2)
    σ_prop ∈ RealDomain(lower = 0, init = 0.04)
  end
end

For PumasEMModels, the @param block uses a formula syntax for covariate relationships. A distribution from the normal family is used to indicate domain (currently only Normal, LogNormal, and LogitNormal are supported):

@emmodel begin
  @param begin
    CL ~ 1 + wt | LogNormal
    Vc ~ 1 | Normal
  end
end

this gives CL = exp(log(θ.CL[1]) + θ.CL[2] * wt) and VC = θ.VC, where θ is the parameter named tuple. If we only have 1 in the formula, then it is a scalar. If we have covariates, we have a tuple of parameters for variable.

@post

Allows for post-processing dynamical variables, to define any variables needed in the @error block.

Must be used in an @emmodel block. For example:

@emmodel begin
  @random begin
    CL ~ 1 | LogNormal
    Vc ~ 1 | LogNormal
  end
  @dynamics Central1
  @post begin
    cp = Central / Vc
  end
  @error begin
    dv ~ CombinedNormal(cp)
  end
end

@pre

Pre-processes variables for the dynamic system and statistical specification. Must be used in an @model block. For example:

@model begin
  @params begin
    tvcl ∈ RealDomain(lower = 0)
    tvv ∈ RealDomain(lower = 0)
    ωCL ∈ RealDomain(lower = 0)
    ωV ∈ RealDomain(lower = 0)
    σ_prop ∈ RealDomain(lower = 0, init = 0.04)
  end
  @random begin
    ηCL ~ Normal(ωCL)
    ηV ~ Normal(ωV)
  end
  @covariates wt
  @pre begin
    CL = tvcl * (wt / 70)^0.75 * exp(ηCL)
    Vc = tvv * (wt / 70) * exp(ηV)
    ka = CL / Vc
    Q = Vc
  end
end

@random

defines the model's random effects and corresponding distribution, e.g. η ~ MvNormal(Ω). Must be used in an @model block. For example:

@model begin
  @random begin
    η ~ MvNormal(Ω)
  end
end

For PumasEMModels, the @random block is equivalent to the @params block, except that each variable defined here has an additional random effect. For example:

@emmodel begin
  @random begin
    CL ~ 1 + wt | LogNormal
    Vc ~ 1 | Normal
  end
end

this gives CL = exp(log(θ.CL[1]) + θ.CL[2] * wt + η_CL) and VC = θ.VC + η_VC, where θ is the parameter named tuple. If we only have 1 in the formula, then it is a scalar. If we have covariates, we have a tuple of parameters for variable.

@vars

Define variables usable in other blocks. It is recommended to define them in @pre instead if not a function of dynamic variables. Must be used in an @model block. For example:

@model begin
  @vars begin
    conc = Central / Vc
  end
end

RandeffsInitialization

Strategy for initializing random effects at each iteration when fitting a model with FOCE or LaplaceI.

Available options are:

• RandeffsInitialization.Center
• RandeffsInitialization.PreviousOrth
• RandeffsInitialization.CenterOrPreviousOrth

Initialize random effects at central values (MCETA = 0 in NONMEM).

Initialize random effects at central values or orthogonalized values of the previous iteration, whatever yields the better conditional log-likelihood, and at central values in the initial iteration.

Initialize orthogonalized random effects at values of the previous iteration, and at central values in the initial iteration.

Experimental

DataFrame(
  fpmi::FittedPumasMixtureModelInspection;
  include_covariates = true,
  include_observations = true,
  include_events = true,
  include_icoef = true,
  include_dosecontrol = true,
  include_ebes = true,
  include_dynamics = true,
)

Create a long-format DataFrame from an inspect result on a fitted mixture model, with one row per (subject, time), following Monolix's conditional / marginal split for mixtures:

• individual quantities (*_ipred/*_iwres, random effects, individual coefficients, dynamics) are taken from each subject's most-probable (MAP) component, i.e. its subpopulation;
• population quantities (*_pred/*_wres) are the mixture marginal, prior-weighted Σₖ πₖ·qₖ over the components;
• *_npde/*_ewres are marginal predictive diagnostics (per subject) and kept as-is.

Both the per-subject responsibilities wᵢₖ (responsibilities, as probability₁ … probabilityₖ columns) and the MAP assignment (subpopulation) are included.

Binless(; bandwidth, qreg_method=IP(), default_bandwidth=nothing)

Quantile-regression smoothing with no bins.

bandwidth is required: a kernel width that suits one dataset usually suits no other, and a silent default does more harm than good. Pass a number, or a Dict keyed by stratum NamedTuple for per-stratum control with default_bandwidth covering the strata the Dict omits. For a discrete VPC, bandwidth = 0.0 asks for the cross-validated automatic choice.

Binned(; strategy = EqualWidth(), nbins)
Binned(strategy::ManualBreaks)

Binned VPC: summarize the observations per x bin (per-bin quantiles for a continuous VPC, per-bin proportions for a discrete one), plotted at the bin midpoints. strategy picks the bin boundaries (EqualWidth(), EqualCount(), or ManualBreaks); nbins sets the bin count for the computed strategies and must be omitted for ManualBreaks, whose breaks fix it.

CensoredVPCSpec(; y, x, pred=nothing, quantiles=(0.1, 0.5, 0.9))

Specify a censored VPC. y and x are required. Produces a continuous VPC panel (quantile regression with censored values imputed) and separate %BLQ and %ALQ panels (proportion of observations below/above the quantification limit over time with simulation CI). Censoring bounds are extracted from the model's Distributions.Censored distribution in @derived.

ContinuousVPCSpec(; y, x, pred=nothing, quantiles=(0.1, 0.5, 0.9))

Specify a continuous VPC. y and x are required.

DirectBins()

Binned VPC with one bin per unique x value: no smoothing and no boundary choice, each observed x value is summarized on its own and plotted at that value. Suited to designs with a small set of nominal times.

DiscreteVPCSpec(; y, x, count_vals=nothing, cumulative=:none)

Specify a discrete VPC (e.g. Bernoulli, Categorical, OrderedLogit outcomes). y and x are required. count_vals is nothing for Bernoulli and a vector of category values for Categorical/OrderedLogit. cumulative controls panel aggregation: :none (default) shows the proportion in each category, :leq the proportion at or below each category, :geq the proportion at or above it. A discrete VPC compares observed proportions against simulated ones, so the panels are proportions rather than probabilities.

Discrete VPCs do not support prediction correction, so there is no pred field (passing pred is an error).

EqualCount()

Bin-boundary strategy for Binned: nbins bins holding approximately equal numbers of observations. The edges are non-interpolating quantiles, so they are actual observed x values.

EqualWidth()

Bin-boundary strategy for Binned: nbins bins of equal x width spanning the observed range.

FiniteMarkov(transitions::Vector{Matrix{<:Real}}, initial::Vector{<:Real}}

Distribution-like structure for discrete-time finite Markov chains with time varying transition matrices transitions and initial value distribution initial. The elements are stored as logits, i.e. the transition probabilities are defined

\[P(X_t=j|X_{t-1}=i) = \begin{cases} \frac{\exp{\tau_{tij}}}{1 + \sum_{l=1}^{k-1} \exp{\tau_{til}}} & \text{if } j \in \{1, \ldots, k - 1\},\\ \frac{1}{1 + \sum_{l=1}^{k-1} \exp{\tau_{til}}} & \text{if } j = k, \end{cases}\]

where $T_t = (\tau_t)_{ij}$ is an $k \times k - 1$ matrix and the tth element of the transition vector. The vector initial contains the logits of the probabilities for the initial value with the last logit implicitly being zero similarly to the logits for the transition probabilities. The vector transitions should have the same length as the number of observations. As a consequence, the first transition matrix is not referenced.

The struct supports evaluation of the log-likelihood as well as simulation.

Examples

@vars begin
  conc = Central / V

  # Drug effects on transition probabilities (logit scale)
  drugep₁₂ = slope₁₂ * conc
  drugep₂₃ = slope₂₃ * conc

  # Transition probabilities with and without drug effect
  θde₁₂ = θ₁₂ + drugep₁₂
  θde₂₃ = θ₂₃ + drugep₂₃
end

@derived begin
  Pθ := [
    [
      θ₁₁ _θ₁₂  θ₁₃
      θ₂₁  θ₂₂ _θ₂₃
      θ₃₁  θ₃₂  θ₃₃
      -Inf -Inf -Inf
    ] for (_θ₁₂, _θ₂₃) in zip(θde₁₂, θde₂₃)
  ]

  DV ~ Pumas.Experimental.FiniteMarkov(Pθ, [θ₁₁, θ₁₂, θ₁₃])
end

ImportanceSampling(; proposal::Union{Prior,LaplacePosterior})

Monte Carlo sampler that draws random effects from proposal and reweights them to the conditional distribution of the random effects given the data. proposal can be Prior or LaplacePosterior.

Used as the sampler of MCEM.

LaplacePosterior(alg::Union{FOCE,LaplaceI}; dof::Real = 4)

Proposal that samples the random effects from a Laplace approximation to their posterior (conditional) distribution, based on the likelihood approximation method alg, which can be either FOCE or LaplaceI.

The proposal is a multivariate Student-t with dof degrees of freedom (dof = Inf gives a Gaussian proposal); heavier tails (lower dof) keep the variance of the importance-sampling weights finite (Geweke, 1989). By comparison, NONMEM's IMP method uses a Gaussian proposal by default (dof = Inf), whereas saemix and Monolix default to Student-t proposals with dof = 4 and dof = 5, respectively.

Used as the proposal of ImportanceSampling.

MCEM(; iterations::Int, samples::Int[, optim_alg, optim_options::NamedTuple, sampler::ImportanceSampling])

Monte Carlo expectation-maximization (MCEM) algorithm for maximum-likelihood estimation of population parameters in Pumas models.

Keyword Arguments

• iterations::Int: number of MCEM iterations.
• samples::Int: number of Monte Carlo samples of random effects for each subject in each iteration.
• optim_alg (default: Optim.NelderMead()): algorithm implemented in Optim.jl for optimizing the population parameters in each iteration ("M step").
• optim_options::NamedTuple (default: (;)): optimization-specific options with names and values that match the Optim.Options type of the optimization algorithm optim_alg. By default, options show_trace = false (optimization trace is not printed) and allow_f_increases = false (optimization is aborted if the objective function value increases) are used, unless specified otherwise in optim_options.
• sampler::ImportanceSampling (default: ImportanceSampling(; proposal = Prior())): Monte Carlo sampler for the random effects in the E-step. See ImportanceSampling and its proposals Prior and LaplacePosterior.

Log-Likelihood

Since MCEM is an expectation-maximization algorithm, the (log-)likelihood is not calculated during fitting. To approximate the log-likelihood of a Pumas model fpm_mcem fitted with MCEM, you can use likelihood approximation methods such as FOCE or LaplaceI:

loglikelihood(fpm_mcem, FO())
loglikelihood(fpm_mcem, FOCE())
loglikelihood(fpm_mcem, LaplaceI())

Post-Processing and Diagnostics

Pumas models fitted with MCEM do not yet support:

1. asymptotic inference and inference methods based on sampling importance re-sampling (SIR), bootstrapping and marginal MCMC using infer,
2. calculations of empirical Bayes estimates using empirical_bayes,
3. predictions without explicitly provided values of random effects using predict,
4. diagnostics and calculations of weighted residuals using inspect.

To work around these limitations, you can construct a fitted Pumas model with the parameter estimates obtained with MCEM that supports these features by performing a follow-up "estimation" with e.g. FOCE or LaplaceI for zero iterations:

fpm = fit(model, data, coef(fpm_mcem), FOCE(); optim_options = (; iterations = 0))
infer(fpm)
empirical_bayes(fpm)
predict(fpm)
inspect(fpm)

ManualBreaks(breaks)

Bin-boundary strategy for Binned: user-supplied bin edges. breaks must be strictly increasing, contain at least two points, and span the data (an x value outside the outermost breaks is an error). The break points determine the bin count (one fewer than the number of breaks), so construct the binning as Binned(ManualBreaks(breaks)) without nbins.

Prior()

Proposal that samples the random effects from their prior distribution.

Used as the proposal of ImportanceSampling.

PumasMixtureModel(models::Tuple{Vararg{PumasModel}}, mixprobs::Pair{Symbol, <:Union{SimplexDomain, Dirichlet}})

A mixture of PumasModels. The marginal likelihood for a subject is

$L(y|\theta) = \sum_k \pi_k L_k(y|\theta_k)$

where $\pi_k$ are the mixture probabilities (estimated as part of the parameter vector) and $\theta_k$ is the subset of parameters needed by submodel $k$.

Parameters are shared by name: if two submodels both reference tvcl, the same value is used for both. The combined parameter set is the union of all submodel parameters plus a user-specified mixture probability parameter constrained to the unit simplex.

Constructor

PumasMixtureModel(models, :name => SimplexDomain(N))
PumasMixtureModel(models, :name => Dirichlet(N, 1.0))

Constructs a mixture model from a tuple of PumasModels and a mixture probability specification. The specification is a Pair mapping a parameter name to either a SimplexDomain (flat prior) or a Dirichlet distribution (informative prior). The number of probabilities must match the number of models. All submodels that define a parameter with the same name must declare it identically, whether as a domain or as a prior distribution.

Random effects and diagnostics such as predict, empirical_bayes, wresiduals, icoef and dosecontrol are per-component quantities: they are given and returned as a Tuple with one entry per submodel, in the order the submodels were passed to the constructor. For a population the result is a Vector of such tuples, one per subject.

VEM

VEM stands for Variational Expectation Maximization. It uses gradient-based optimization to maximize a lower bound on the log marginal likelihood with an optional prior on the population parameters. VEM finds:

1. Point estimates of the population parameters, and
2. An approximation of the individual posterior distribution over the random effects of each subject, conditional on the population parameters.

There are 2 variants of VEM:

1. Presampling VEM (default): the base random variable of the variational family is pre-sampled at the start of each phase, and reused within the iterations of that phase. This variant uses a deterministic optimizer such as L-BFGS. Note that the seed still affects the samples used, so different seeds will lead to different results.
2. Resampling VEM: the base random variable of the variational family is re-sampled at each iteration. This variant uses a stochastic optimizer such as Adam.

The variant of VEM used can be controlled through the presample keyword argument.

Example fit

To fit a model using VEM, you can use:

fpm = fit(model, data, param, VEM())

Keyword arguments of the VEM algorithm struct

All of the keyword arguments of the VEM algorithm and their default values are described below.

• presample::Bool = true: whether to pre-sample the base random variable of the variational family before optimization starts; pre-sampled for each phase but reused within iterations of each phase. If false, they will be re-sampled at each iteration.
• nphases::Int = default_vem_nphases(presample): number of phases. Each one resets the optimizer and has its own learning rate and number of samples.
• nsamples::Int = 0: if not zero, it sets the number of samples to always be this value, but only through the default values of init_nsamples and nsamples_increment.
• init_nsamples::Int = default_vem_init_nsamples(presample, nsamples): number of samples of the base distribution of the variational family per Subject and iteration, during the first phase. Shared between Subjects (if share_samples is true) and/or iterations (if presample is true).
• nsamples_increment::Int = default_vem_nsamples_increment(presample, nsamples, init_nsamples): by how much to increment the number of samples from one phase to the next.
• share_samples::Bool = false: whether to share samples between all Subjects, as opposed to each having their own samples.
• prior::Bool = true: whether to consider the prior on the population parameters, if any, as defined in the @param section of a PumasModel.
• diagonal::Bool = true: whether to approximate the posterior on the random effects as having a diagonal covariance matrix (i.e. as all elements being independent).
• base_dist_std::Float64 = 0.1: standard deviation of the base distribution of the variational family.
• max_nsamples::Int = default_vem_max_nsamples(presample): maximum number of samples from the base distribution of the variational family allowed per phase.
• use_ebes::Bool = false: whether to use the empirical Bayes estimates of the random effects as the posterior approximation's mean after every phase.

Keyword arguments of fit

The optimizer, parallelism, and differential equation simulator options are passed as keyword arguments to fit:

• optim_alg = presample ? default_vem_optim() : default_vem_stoch_optim(): optimizer to use. Must be of type Optim.FirstOrderOptimizer if presample is true, and of type Optimisers.AbstractRule otherwise.
• optim_options::NamedTuple = (;): options of the optimizer.
• rng::AbstractRNG = Random.default_rng(): the pseudorandom number generator. By default, a TaskLocalRNG.
• diffeq_options::NamedTuple = (;): options of the differential equation simulator.
• ensemblealg = EnsembleThreads(): EnsembleSerial() for computing gradients using a single thread, or EnsembleThreads() for multithreading.
• constantcoef::Tuple = (): a Tuple of Symbols of the parameters to keep constant during the fit.

Objective function value

When calling fit, the negative of the lower bound is minimized, which is equivalent to maximizing the lower bound on the log marginal likelihood. This has the effect of indirectly maximizing the marginal likelihood. The objective function value displayed during optimization is the negative of the lower bound up to an omitted constant. The log-likelihood is not calculated during fitting.

Log-likelihood

To approximate the log-likelihood after a VEM fit, you can use one of the following:

loglikelihood(fpm, FO())
loglikelihood(fpm, FOCE())
loglikelihood(fpm, LaplaceI())

Confidence intervals and standard errors

Models fitted with VEM do not support inference methods yet, such as the asymptotic estimator, SIR and bootstrap. To get confidence intervals and standard errors after fitting a model with VEM, you can use the VEM fit as a starting point in a second fit using an algorithm that supports inference, e.g. FOCE() or LaplaceI(). Example:

fpm_vem = fit(model, data, param, VEM())
fpm_foce = fit(model, data, coef(fpm_vem), FOCE(), init_randeffs = empirical_bayes(fpm_vem))
infer(fpm_foce)

Diagnostics and weighted residuals

When calling inspect on a VEM-fitted Pumas model (fpm) with a:

1. Gaussian,
2. Censored latent Gaussian, or
3. Truncated latent Gaussian error model, you need to specify a likelihood approximation method to calculate the weighted residuals using the wres_approx keyword argument, e.g.

inspect(fpm; wres_approx = FOCE()) # uses FOCE()
inspect(fpm; wres_approx = LaplaceI()) # uses LaplaceI()

If the error model is not one of the above, the weighted residuals are not calculated and the wres_approx argument is ignored.

VIS(nsamples::Int)

A likelihood approximation method that uses importance sampling with a number of samples nsamples for each subject, using the variational approximation of the individual posterior distribution as a proposal distribution.

This can be used to evaluate the log-likelihood of a model fitted with VEM using importance sampling given the final approximations of the individual posteriors.

icl(fpm::FittedPumasMixtureModel) -> Float64

Integrated Completed Likelihood for choosing the number of mixture components. Uses the completed-data form ICL = -2·ℓc + dof·log(nobs), where the completed log-likelihood ℓc = Σᵢ log(π_ẑᵢ · Lᵢ,ẑᵢ) evaluates each subject under its most-probable component ẑᵢ = argmaxₖ wᵢₖ (subpopulation). Equivalently ICL = BIC - 2·Σᵢ log(maxₖ wᵢₖ), i.e. BIC plus a penalty for classification uncertainty. Lower is better; unlike BIC it favors mixtures whose components are well separated.

Introduced by Biernacki, Celeux & Govaert (2000); see McLachlan & Peel (2000), Finite Mixture Models, for the finite-mixture/classification framing.

impute_censored(obs_df, sim_dfs; id=:id, y, lloq=nothing, uloq=nothing, rng=Random.default_rng())

Replace censored observations in obs_df with draws from the simulated values in sim_dfs that fall within the censored region.

Censoring is inferred from y and the limits, and the bounds are inclusive: an y is below the limit of quantification (BLQ) when it is at or below the lower limit (y <= lloq) and above (ALQ) when it is at or above the upper limit (y >= uloq). A value exactly on a bound is therefore treated as censored. No explicit censoring flag is required. lloq/uloq can be a Real scalar or a Symbol column name (per-row thresholds); at least one must be provided.

impute_missings(obs_df, sim_dfs; id=:id, y, rng=Random.default_rng())

Replace missing y values in obs_df with draws from the simulated values in sim_dfs at the same subject and y position. Rows are matched by their order within each id group, so obs_df and every frame in sim_dfs must share the same per-subject row layout. Returns a copy of obs_df with the missing y values filled in.

npde(obs_df, sim_dfs; id=:id, y, x, lloq=nothing, uloq=nothing, rng)

Compute NPDE and EWRES from pre-simulated DataFrames, decorrelating against the empirical covariance of the replicates.

That covariance is a sum of nsim outer products of deviations from their own mean, so its rank is at most nsim - 1, and a subject with at least as many observations as replicates is rejected rather than silently regularized. For the model-based covariance, which stays well conditioned at modest replicate counts, pass the fit as the first argument: npde(fpm, obs_df, sim_dfs; ...).

If lloq/uloq are provided, censored observations (inferred from y and the limits) are imputed first. Returns an NPDEResult.

npde_tests(result::NPDEResult)

Run the standard NPDE test battery (matching the npde R package):

• Wilcoxon signed-rank: H₀ = NPDE symmetric around 0
• Variance test: H₀ = Var(NPDE) = 1 (chi-squared test)
• Shapiro-Wilk: H₀ = NPDE ~ Normal
• Global: Bonferroni-adjusted minimum p-value

Returns a DataFrame.

observed_vpc(vpctype, binning, df; stratify_by=Symbol[], transform=identity, ...)
observed_vpc(binning, fpm; y, x, stratify_by=Symbol[], ...)

Summarize the observed data, and resolve every choice the simulated pass will reuse.

This is where a VPC is configured. The binning strategy and its resolved bin boundaries or quantile-regression nodes, the per-stratum bandwidth, the strata, the y transform and the prediction-correction factors are all settled here and stored on the returned ObservedVPC. simulated_vpc reads them back rather than resolving anything of its own, so the band and the curve are always computed on the same grid: a difference between them can then only come from the model, which is what a VPC exists to show.

To re-bin or re-smooth, build a new ObservedVPC and summarize the same replicates against it. There is deliberately no way to pass a different binning to simulated_vpc.

The first form takes a prepared DataFrame and needs no model. The second builds that frame from a fit, and is the one that needs y and x, since a spec has to be resolved first.

Returns an ObservedVPC, or an ObservedCensoredVPC for a CensoredVPCSpec, whose continuous, censored_lower and censored_upper fields hold the quantile, %BLQ and %ALQ panels.

responsibilities(fpm::FittedPumasMixtureModel) -> Vector{Vector{Float64}}

Per-subject posterior component probabilities (responsibilities, a.k.a. IPmix): for subject i and component k,

wᵢₖ = πₖ exp(-nllᵢₖ) / Σⱼ πⱼ exp(-nllᵢⱼ),

where πₖ are the (population) mixture probabilities and nllᵢₖ is subject i's marginal negative log-likelihood under component k. Returns one probability vector (summing to 1) per subject. See also subpopulation.

simulated_vpc(obsvpc, sim_dfs; level=0.95)

Summarize each simulated replicate the same way the observed data was summarized, and reduce the replicates to a confidence band.

Every summarization setting comes from obsvpc: the binning, the resolved bin boundaries or nodes, the bandwidth, the strata, the transform and the prediction-correction factors. Nothing is resolved here, and there is deliberately no keyword to override any of it, because a band computed on a different grid than the observed curve would differ for reasons that have nothing to do with the model.

sim_dfs must reproduce the observed design: each replicate is matched to the observed frame on (id, time), and a replicate with rows the observed data does not have, or missing rows it does, is an error rather than a silently thinner band.

level sets the width of the interval around each simulated quantile. It is the only thing decided at this step, and vpc_summarize recomputes it from the stored per-replicate summaries without re-simulating.

Returns a SimulatedVPC, or a SimulatedCensoredVPC for a censored observed VPC.

subpopulation(fpm::FittedPumasMixtureModel) -> Vector{Int}

Most-probable component (MAP / NONMEM MIXEST) for each subject, i.e. argmaxₖ wᵢₖ of the responsibilities.

vpc(vpcspec::VPCSpec, obs_df, sim_dfs; binning, level=0.95, ...)

Run the VPC pipeline from pre-simulated DataFrames. The VPC type must be provided explicitly (e.g. ContinuousVPCSpec(; y=:conc, x=:time)).

level is a confidence level, as elsewhere in the package: the interval around each simulated quantile, which sets the shaded band width.

Returns a SimulatedVPC or a SimulatedCensoredVPC.

vpc_setup(derived_distribution; y, x, pred=nothing, quantiles=nothing)

Return the VPC specification for y against x, taking the outcome family from one distribution of the error model.

y and x are required. pred turns on prediction correction and applies to continuous and censored outcomes; quantiles chooses which quantile curves to estimate and does not apply to discrete outcomes. Passing either where it does not apply is an error, as is a bare censored(...), since a censored VPC needs the latent scale.

With no distribution at hand, construct the spec directly: ContinuousVPCSpec, DiscreteVPCSpec, CensoredVPCSpec.

vpc_summarize(sv::SimulatedVPC; level=0.95)

Recompute the simulated band from sv at a different confidence level, without simulating again.

Returns a DataFrame with :lower, :middle and :upper, keyed by the independent variable and the stratification columns, plus :τ for a continuous VPC or :category for a categorical one — the same layout as sv.summarized_vpc.

Only level can be changed this way. The binning, bandwidth, nodes, stratification and transform come from the ObservedVPC; to change one of those, build a new ObservedVPC and call simulated_vpc with the same replicates.

wider = vpc_summarize(sim; level = 0.99)

function GlobalSensitivity.gsa(
  func,
  model::PumasModel,
  population::Union{Population,Subject},
  params::NamedTuple,
  method::GlobalSensitivity.GSAMethod,
  p_range_low::NamedTuple = NamedTuple(),
  p_range_high::NamedTuple = NamedTuple();
  constantcoef::Tuple = (),
  rng::AbstractRNG = Random.default_rng(),
  diffeq_options::NamedTuple = NamedTuple(),
  ensemblealg::SciMLBase.EnsembleAlgorithm = EnsembleThreads(),
  samples = 1000,
  kwargs...,
)

Perform a global sensitivity analysis (GSA) of the population mean of quantities computed by the post-processing function func on the parameters of the Pumas model.

func is a post-processing function with the same calling convention as postprocess: it takes (gen, obs), where gen is a NamedTuple of all the generated quantities of one simulated subject and obs is the NamedTuple of that subject's observations, and it must return a NamedTuple of the quantities to run GSA on. GSA is run on every entry of the returned NamedTuple, named by its keys. The function is typically passed using do syntax. For a Population, the per-subject results are aggregated into a single value with the population mean.

Since gen holds everything postprocess provides, a quantity of interest can be computed without adding it to the model. An entry of the returned NamedTuple may be a vector, in which case GSA is run on each element separately. The quantities should be deterministic: variance from sampling an error model cannot be attributed to any parameter.

Arguments

• func: the post-processing function, see above.
• model: a PumasModel.
• population: a Population of Subjects or a single Subject.
• params: a named tuple of parameters.
• method: one of the GSA methods from GlobalSensitivity.jl such as Sobol(), Morris(), eFAST(), or RegressionGSA().
• p_range_low & p_range_high: the lower and upper bounds for the parameters you want to run the GSA on. For parameters without lower or upper bounds in p_range_low or p_range_high, respectively, bound(s) in the @param block are used. If bounds for a parameter are neither specified as arguments nor in the @param block, an error is thrown.

Keyword arguments

• constantcoef: A tuple of the names of the parameters that are not to be varied in the GSA.
• rng: The random number generator to use.
• samples: The number of samples to use for the GSA.
• diffeq_options: The options of the differential equation solver, passed on to simobs when computing the derived variables of interest.
• ensemblealg: The mode of parallelization of the simobs calls when computing the derived variables of interest. Possible options are: EnsembleSerial(), EnsembleThreads(), EnsembleDistributed() and EnsembleSplitThreads().
• kwargs...: The remaining keyword arguments are passed on to GlobalSensitivity.gsa, as documented in GlobalSensitivity.jl.

Examples

using Pumas, GlobalSensitivity

model = @model begin
  @param begin
    tvka ∈ RealDomain(lower = 0.1, upper = 3.0)
    tvcl ∈ RealDomain(lower = 0.5, upper = 10.0)
    tvvc ∈ RealDomain(lower = 10.0, upper = 40.0)
    σ ∈ RealDomain(lower = 0.0)
  end
  @covariates wt
  @pre begin
    Ka = tvka
    CL = tvcl * (wt / 70)^0.75
    Vc = tvvc * (wt / 70)
  end
  @dynamics begin
    Depot' = -Ka * Depot
    Central' = Ka * Depot - (CL / Vc) * Central
    AUC' = 1000 * Central / Vc
  end
  @derived begin
    dv ~ @. ProportionalNormal(1000 * Central / Vc, σ)
  end
end

# A single 100 mg dose, sampled over 24 hours, across a range of body weights
population = map(enumerate([55.0, 70.0, 85.0, 100.0])) do (i, wt)
  Subject(;
    id = i,
    time = [0, 0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 16, 24],
    covariates = (; wt),
    events = DosageRegimen(100),
  )
end

# σ = 0 so that the simulated `dv` is the deterministic model prediction
params = (; tvka = 1.5, tvcl = 1.1, tvvc = 20.0, σ = 0.0)

# How sensitive are cumulative exposure and the concentration 24 hours after the dose to
# absorption and disposition? `tvvc` is not given a range, so its `@param` bounds are used.
gsa(
  model,
  population,
  params,
  eFAST(),
  (; tvka = 0.1, tvcl = 0.5),
  (; tvka = 2.5, tvcl = 3.0);
  constantcoef = (:σ,),
  samples = 500,
) do gen, _
  return (; AUC = last(gen.AUC), C24 = last(gen.dv))
end

# output

eFAST Sensitivity Analysis

First Order Indices
2×4 DataFrame
 Row │ dv_name  tvka        tvcl      tvvc
     │ String   Float64     Float64   Float64
─────┼──────────────────────────────────────────
   1 │ AUC      0.00974594  0.630539  0.229083
   2 │ C24      0.0169295   0.88679   0.0308932

Total Order Indices
2×4 DataFrame
 Row │ dv_name  tvka       tvcl      tvvc
     │ String   Float64    Float64   Float64
─────┼─────────────────────────────────────────
   1 │ AUC      0.0163925  0.757388  0.351633
   2 │ C24      0.0243564  0.94935   0.0954693

References

For more details regarding the GSA methods, their options, and keyword arguments please consult the GlobalSensitivity.jl documentation.