Global Sensitivity Analysis

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.