Inference

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, 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, 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(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(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.

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.

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(fit::DissolutionFit) -> DataFrame

One row per formulation: formulation, then one column per parameter of the fitted model in field order (finf, mdt, b, int, tlag for a Weibull). The parameters in fit.constant were held fixed.

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)

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.

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.

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.

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.

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.

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.