Simulation and Estimation Diagnostics

Model diagnostics are used to check if the model fits the data well. Diagnostics for continuous data are different from for discrete data.

The inspect Function

The main function to call when looking at diagnostics from a fit is inspect with the signature:

Pumas.inspect — Function
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.

The function outputs a summary of the model predictions, residuals, Empirical Bayes estimates, and Normalised Prediction Distribution Errors (NPDEs) and Expected Weighted RESiduals (EWRES) (when requested). It is 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. The keyword nsim controls the number of times each subject is simulated for the npde and ewres computations. rng is a random number generator that can be used to seed any random elements such as those in npde and ewres computations.

inspect returns a FittedPumasModelInspection object with the

  • model predictions
  • residuals
  • empirical Bayes estimates
  • individual coefficients
  • dose control parameters
  • NPDEs and EWRES with nsim set

Example syntax

my_inspect = inspect(my_fit)
my_inspect_df = DataFrame(my_inspect)

The result of inspect(my_fit) is a FittedPumasModelInspection object, which can then be converted into tabular form using the DataFrame call. This is useful if one is interested to export the results into a spreadsheet, or use the resulting DataFrame for plotting.

Weighted Residuals

In Pumas, you can calculate individual and population weighted residuals using wresiduals as well as expected individual weighted residuals with the signature

Pumas.wresiduals — Function
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.

Takes a fit result, an approximation method for the marginal likelihood calculation which defaults to the method used in the fit and the number of simulations with the keyword argument nsim. If nsim is specified only the Expected Simulation based Individual Weighted Residuals (EIWRES) is included in the output as individual residual and population residual is not computed. Using the FO approximation method corresponds to the WRES and while FOCE corresponds to CWRES. The output is a SubjectResidual object that stores the population (wres) and individual (iwres) residuals along with the subject and approximation method (approx).

Both forms weight the residual by the dispersion of the error model: the individual weighted residual divides by the standard deviation of the error distribution of the observation, and the population weighted residual weights by a first-order approximation of the marginal covariance of the subject's observations. Both therefore require an error model that is described by a mean and a variance, and weighted residuals are computed only for dependent variables whose error model is Gaussian.

The error model of a dependent variable is the distribution it is written with in the @derived block, so whether it qualifies is decided per dependent variable. The following qualify:

The following do not:

  • distributions that are not Gaussian, even when they describe a transformed Gaussian variable: y ~ @. LogNormal(log(μ), σ) does not qualify, whereas modelling the log-transformed observations as logy ~ @. Normal(log(μ), σ) does;
  • truncated(Normal(μ, σ); lower, upper) and censored(Normal(μ, σ); lower, upper), which truncate or censor the observed variable itself rather than a latent variable;
  • the discrete error models, and TimeToEvent.

A model that mixes Gaussian and non-Gaussian error models therefore gets weighted residuals for the Gaussian dependent variables only, and if none of them qualifies the result is empty rather than an error.

Predictions

To calculate predictions from a fit use the predict function with the signature:

StatsAPI.predict — Function
predict(d::NCA.DoseLinearityPowerModel, xs::Vector; level = nothing, interval = :confidence)

Predict the response of a power model d for the given vector xs of doses. The level argument sets the confidence interval level for the prediction. If level is not provided, the level used to fit the model is used. The interval argument can be set to :confidence or :prediction. The :confidence interval is the default and provides a confidence interval for the mean response. The :prediction interval provides a prediction interval for a new response.

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

For dependent variables with a Categorical, OrderedLogit or Binomial error model, the DataFrame of a prediction contains one column of outcome probabilities per category. The columns are named after the dependent variable with the category as a subscript, so a painord variable with four categories gives painord₁_pred to painord₄_pred, and likewise for _ipred. The categories are the values the dependent variable can take, which means they start at 1 for Categorical and OrderedLogit and at 0 for Binomial. Bernoulli variables are not expanded, since their prediction already is the success probability.

The scalar painord_pred and painord_ipred columns are kept alongside them, but they hold the mean of the error distribution, which generally is not a meaningful prediction for categorical data. They are retained for backwards compatibility for Categorical and Binomial only; for OrderedLogit they contain missing. Use the per-category probabilities instead.

Model Metrics

Akaike information criterion:

StatsAPI.aic — Function
aic(fpm::AbstractFittedPumasModel)

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

Bayesian information criterion:

StatsAPI.bic — Function
bic(fpm::AbstractFittedPumasModel)

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

Shrinkage

Pumas.ηshrinkage — Function
η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(Ω))
Pumas.ϵshrinkage — Function
ϵ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.

If the observations for an individual is informative, η-shrinkage will be low. If, on the other hand, the observations are non-informative, the individual estimates (EBEs) will shrink towards the population mean. Hence, the variance of the EBEs will also shrink and η-shrinkage will be high.

η-shrinkage is reported for random effects that are Normal or MvNormal and that do not depend on covariates. Random effects that do not meet these conditions are left out of the result and a warning naming them is shown.

Passing stratify_by computes the shrinkage separately within each stratum, and returns one set of shrinkage values per stratum rather than a single one:

ηshrinkage(my_fit; stratify_by = [:sex])

The covariates used for stratification have to be constant within each subject.

When data gets sparse the IWRES distribution also shrink towards it’s mean (zero). This is often called overfitting or ϵ-shrinkage.