API reference

NLMEEstimation.NLMEEstimation Module

julia

NLMEEstimation

Model-agnostic estimators for nonlinear mixed-effects (NLME) models.

A model is any Julia value for which the protocol functions are implemented (see nsubjects, nlme_layout, and per evaluation_style either conditional_loglikelihood or conditional_loglikelihoods). The package owns the random-effect prior η_i ~ N(0, Ω) and its parameterization (see NLMELayout); the model owns everything else.

Derivatives come from automatic differentiation: a model either implements the derivative methods itself or declares an AD backend through ad_backend. The package promises and tests reverse-mode AD with Enzyme. Second order is optional and only ever per-subject blocks (provides_eta_hessian).

Estimators:

empirical_bayes (alias posthoc): per-subject posterior modes (EBEs), by Newton on Hessian blocks or by gradient-only BFGS, at fixed population parameters.
focei: population estimation using first-order conditional curvature and a derivative-free outer optimizer.
saem: population estimation by MCMC stochastic approximation.

covariance_step gives standard errors of a fit's population coordinates, and marginal_loglikelihood its likelihood, OFV, AIC and BIC (for SAEM, by importance sampling or adaptive Gauss–Hermite quadrature). residual_diagnostics gives PRED, IPRED, IWRES, CWRES and NPDE.

Each of these functions takes progress = nothing: a Treebars node under which it builds its own labelled progress node, or a backend symbol such as :term for a new root. With nothing no node is created. Treebars.request_interrupt! on that node stops saem and focei early with their estimate so far; the other functions raise Treebars.ProgressInterrupted.

NLMEEstimation.AdaptiveGaussHermite Type

julia

AdaptiveGaussHermite(; nodes = 5, max_points = 100_000)

Adaptive Gauss–Hermite quadrature of each subject's marginal log-likelihood: nodes points per random-effect dimension, as a tensor product, centred at the subject's conditional mode and scaled by its inverse curvature (computed as for ImportanceSampling). nodes = 1 is the Laplace approximation. A subject with d random effects costs nodes^d batched evaluation points; a subject needing more than max_points is refused.

NLMEEstimation.BatchedEvaluation Type

julia

BatchedEvaluation()

The model evaluates a batch of subjects per call: it implements conditional_loglikelihoods. Per-subject entry points evaluate a batch of one; derivatives of a batch come from AD of the batch's summed log-density, which yields every subject's η-gradient at once because subjects are conditionally independent.

NLMEEstimation.CovarianceResult Type

julia

CovarianceResult

Result of covariance_step. point is the full concatenated population vector [θ; σ; ω] at which the step was taken, in the model's unconstrained coordinates and the package's log-Cholesky ω. coordinates is the same point in the coordinates the estimators fit: [θ; σ; free Ω coordinates]. It equals point unless the model declares fixed_omega; then a fully fixed family has no coordinates, and a partially fixed family contributes its estimated log standard deviations (named log sd(…)) and, when its correlation is estimated, its correlation coordinates m_jl = L_jl / L_jj of the Cholesky factor L (named m(…)). Every other vector and matrix field is indexed by the estimated coordinates free (positions in coordinates), named in names:

covariance, standard_errors: the estimated sampling covariance of those coordinates and the square roots of its diagonal, for method :rs (R⁻¹ S R⁻¹), :r (R⁻¹) or :s (S⁻¹).
R: minus the Hessian of the fitted objective, the summed marginal log-likelihood plus any population_logprior (objective is :ML or :MAP); nothing for method = :s, which does not compute it.
S: Σᵢ gᵢ gᵢᵀ over the subjects' marginal log-likelihood gradients gᵢ, the columns of subject_gradients. A population prior is not a subject term and enters R only.
gradient: the gradient of the fitted objective, close to zero at a converged estimate.
step: the central-difference step per coordinate; evaluations counts the marginal-approximation evaluations.

Use delta_method for natural-unit quantities.

NLMEEstimation.CovarianceStepError Type

julia

CovarianceStepError

Thrown by covariance_step when the requested covariance does not exist: R is not positive definite (methods :rs and :r) or S is singular (method :s, for example with fewer subjects than estimated coordinates). The fields R, S, gradient and names hold what was computed, for diagnosis.

NLMEEstimation.DualAveragingMALA Type

julia

DualAveragingMALA(; metric=:prior, initial_stepsize=0.1,
                    target_acceptance_rate=0.574,
                    gamma=0.05, t0=10.0,
                    min_stepsize=1e-5, max_stepsize=10.0,
                    curvature_floor=1e-4)

MALA with one independent dual-averaging step size per subject. metric=:prior uses the current Ω as proposal covariance. metric=:hessian obtains a fixed per-subject covariance from the initial regularized local posterior Hessian, falling back to Ω if the model has no Hessian. Gradients are evaluated in subject batches. The source gradient is reused after each accept or reject until a population M-step changes the target. Adaptation continues at every transition, including during population updates. This changing-target, continuously adapted mode is experimental; its transition is Metropolis corrected at each fixed population state, but a full SAEM convergence guarantee is not implied.

NLMEEstimation.EBEConvergenceError Type

julia

EBEConvergenceError(subjects, result)

Thrown by empirical_bayes when no mode was found for some subjects; result holds the partial EBEResult.

NLMEEstimation.EBEMethod Type

julia

EBEMethod

Per-subject mode search of empirical_bayes:

NewtonEBE: Newton on the η Hessian blocks, needs provides_eta_hessian.
QuasiNewtonEBE: BFGS on η gradients only.

default_ebe_method picks Newton when Hessian blocks are available.

NLMEEstimation.EBEOptions Type

julia

EBEOptions(; maxiter = 200, decrement_tol = 1e-14, gtol = 1e-8, max_backtracks = 30)

Settings of empirical_bayes.

maxiter: maximum number of iterations (accepted steps) per subject.
decrement_tol: NewtonEBE convergence, half the Newton decrement in log-density units.
gtol: QuasiNewtonEBE convergence, the standardized gradient ‖Lᵀg‖∞.
max_backtracks: maximum trust-region reductions (NewtonEBE) or step halvings (QuasiNewtonEBE) per iteration before a subject is reported as not converged.

NLMEEstimation.EBEResult Type

julia

EBEResult

Result of empirical_bayes. Column/entry k of every field belongs to subject subjects[k].

modes: neta × length(subjects) matrix of empirical Bayes estimates η̂_i.
hessians: the neta × neta Hessian in η of the complete-data log-density log p(y_i, η_i | θ, σ, Ω) at η̂_i per subject, or nothing when the model provides no second order (see ebe_covariances).
conditional_loglikelihood, complete_data_loglikelihood: the two log-densities at η̂_i.
converged, iterations: per-subject solver diagnostics.
gradient_norm: the standardized complete-data gradient ‖Lᵀg‖∞ at η̂_i.
method: the EBEMethod used.

NLMEEstimation.EvaluationStyle Type

julia

EvaluationStyle

How a model evaluates its conditional log-density, declared by evaluation_style:

PerSubjectEvaluation (default): the model implements conditional_loglikelihood for one subject at a time.
BatchedEvaluation: the model implements conditional_loglikelihoods for a batch of subjects in one pass, as joint-program backends (one log-density over all subjects) do.

NLMEEstimation.FOCEApproximation Type

julia

FOCEApproximation(; interaction = true, ebe_method = nothing,
                    ebe_options = EBEOptions(gtol = 1e-10))

The FOCE(I) approximation of each subject's marginal log-likelihood, the per-subject terms of focei_evaluate, as used by covariance_step. interaction = false gives FOCE without interaction. ebe_method = nothing selects default_ebe_method: Newton when the model provides η Hessian blocks, otherwise BFGS. The default mode tolerance is tighter than a fit's because finite differences of the objective amplify the error of the inner mode search. Requires GaussianObservations and an NLMELayout.

NLMEEstimation.FOCEEvaluation Type

julia

FOCEEvaluation

FOCE(I) log-likelihood at fixed population parameters. ebes.modes contains the subject EBEs, and curvatures[k] is the positive η curvature for subject subjects[k]. subject_loglikelihoods[k] is that subject's term of loglikelihood. logposterior adds the optional population log prior exactly once; ofv is minus twice logposterior (equal to minus twice the likelihood when objective == :ML).

NLMEEstimation.FOCEOptions Type

julia

FOCEOptions(; maxeval = 500, rhobeg = 0.2, rhoend = 1e-5,
             interaction = true, ebe_options = EBEOptions())

Settings for the first-order conditional estimation objective and its derivative-free population search. interaction = false omits the variance Jacobian term (FOCE without interaction). The radii are in the unconstrained (θ, σ, ω) coordinates.

NLMEEstimation.FOCEResult Type

julia

FOCEResult

Population estimate from focei. theta, sigma, and omega use the model's declared coordinates; omega_matrix is the estimated Ω on covariance scale. loglikelihood remains the approximate marginal likelihood; logposterior includes a model-provided population prior once, and ofv is minus twice that fitted objective. objective distinguishes :ML from :MAP. free lists the fitted positions in the search vector [θ; σ; free Ω coordinates], which is the concatenated (θ, σ, ω) unless the model declares fixed_omega (see CovarianceResult for those coordinates), and options the FOCEOptions of the fit; covariance_step reuses both.

NLMEEstimation.FixedOmega Type

julia

FixedOmega(Ω::AbstractMatrix)
FixedOmega(; sd = nothing, correlation = nothing)

A known value for all or part of one random-effect covariance family, declared through fixed_omega. FixedOmega(Ω) holds the whole family at the covariance matrix Ω, like nlmixr2's fix() on an ω block or NONMEM's OMEGA … FIX.

The keyword form fixes parts of a family in standard-deviation and correlation terms, Ω = D R D with D = Diagonal(sd):

sd is nothing (every margin estimated) or a vector with one entry per margin: a positive standard deviation for a fixed margin, and nothing or missing for an estimated one;
correlation is nothing (estimated) or the known correlation matrix R, symmetric positive definite with unit diagonal.

The estimators fit the remaining parts with the fixed parts held at their values. A family of size one has only its standard deviation.

NLMEEstimation.GaussianObservations Type

julia

GaussianObservations()

Observation-model declaration (see observation_model): the observations of subject i are independent Gaussians. The model then implements observations and observation_moments.

Contract: conditional_loglikelihood(model, i, θ, σ, η) equals Σ_j log N(y_ij | μ_ij, v_ij) with y_i = observations(model, i) and (μ_i, v_i) = observation_moments(model, i, θ, σ, η).

NLMEEstimation.ImportanceSampling Type

julia

ImportanceSampling(; draws = 2000, df = 4, scale = 1.0)

Importance-sampling estimate of each subject's marginal log-likelihood log ∫ p(y_i | θ, σ, η) N(η; 0, Ω) dη. The proposal for subject i is a multivariate Student t with df degrees of freedom (a positive integer, or Inf for a Gaussian), centred at the subject's conditional mode, with scale matrix scale² times the inverse curvature of the complete-data log-density there. The curvature comes from the η Hessian blocks when the model provides them, and otherwise from central differences of batched η gradients. Each draw costs one batched conditional evaluation over all subjects. The estimate reports its Monte Carlo standard error and per-subject effective sample sizes.

NLMEEstimation.MarginalLikelihoodResult Type

julia

MarginalLikelihoodResult

Result of marginal_loglikelihood. point is [θ; σ; ω]. loglikelihood is the marginal log-likelihood with every constant kept, and subject_loglikelihoods[k] is subject subjects[k]'s term. standard_error and subject_standard_errors are Monte Carlo standard errors. They are zero for the deterministic approximations, which report no approximation error. effective_sample_size holds the per-subject importance sampling effective sample sizes, or nothing.

logposterior adds a model population_logprior once (objective :MAP). ofv is -2 logposterior, as in FOCEResult, and equals -2 loglikelihood under :ML. nlmixr2's OBJF omits the Gaussian constant: it is ofv - nobs * log(2π). aic = -2 loglikelihood + 2 nparameters and bic = -2 loglikelihood + nparameters * log(nobs), where nobs counts the observations of the listed subjects. nobs, and so bic, is missing for a model that declares no observation_model unless it is passed explicitly. evaluations counts batched conditional evaluations (missing for FOCEApproximation).

NLMEEstimation.MuReferencing Type

julia

MuReferencing(theta_indices, design)

Mu-referencing declaration (see mu_referencing).

theta_indices: the θ coordinates that are mu-referenced.
design: one neta × length(theta_indices) matrix X_i per subject.

Contract: for every subject i, conditional_loglikelihood depends on θ[theta_indices] and η only through the individual linear predictor ϕ_i = X_i * θ[theta_indices] + η. Equivalently, it is unchanged under θ[theta_indices] .+= δ, η .-= X_i * δ. The model applies its own inverse link to ϕ_i. A zero row of X_i marks an η component without a mu-referenced typical value.

NLMEEstimation.NLMELayout Type

julia

NLMELayout(; ntheta, nsigma, eta_blocks)

Parameter layout of an NLME model, returned by nlme_layout.

ntheta: length of θ, the population (fixed-effect) coordinates.
nsigma: length of σ, the other scalar coordinates the conditional density depends on (residual-error parameters, for example).
eta_blocks: sizes of the mutually independent blocks of the per-subject random-effect vector η_i, in η order. sum(eta_blocks)is the η dimension (neta). A diagonal Ω isones(Int, neta); one dense Ω is[neta].

θ and σ are in the model's own unconstrained coordinates: every real vector of the right length is a valid input to conditional_loglikelihood.

Ω is owned by this package. Its unconstrained parameter vector ω is the concatenation, block by block, of the log-Cholesky coordinates of each block Ω_b = L_b L_bᵀ: the entries of the lower triangle of L_b in row-major order (L[1,1], L[2,1], L[2,2], L[3,1], …), with each diagonal entry stored as its logarithm. For a block of size one, the coordinate is log(sd). nomega gives the length of ω.

NLMEEstimation.NewtonEBE Type

julia

NewtonEBE()

Trust-region Newton on the complete-data Hessian blocks. Each step maximizes the local quadratic model of the complete-data log-density over a ball in the standardized coordinates L⁻¹η, with Ω = LLᵀ. Inside the ball that is the full Newton step when -H is positive definite; where -H is indefinite (a saddle or a non-concave region) the step follows the directions in which the density curves upward, out to the radius. The radius shrinks after a rejected step or a poor prediction and grows after a good step that reached it. A subject converges once gᵀ(-H)⁻¹g / 2 ≤ decrement_tol; the Newton step from there is taken as a final polishing step, so the returned mode is stationary to roundoff.

NLMEEstimation.NormalREPrior Type

julia

NormalREPrior(layout::NLMELayout, ω)

The random-effect prior η_i ~ N(0, Ω) with Ω = omega_matrix(layout, ω), with its Cholesky factor, precision and log-determinant computed once for repeated per-subject evaluation. This package owns the η prior; models never implement it.

NLMEEstimation.PerSubjectEvaluation Type

julia

PerSubjectEvaluation()

The model evaluates one subject per call: it implements conditional_loglikelihood. Batched entry points loop over subjects.

NLMEEstimation.ProtocolViolation Type

julia

ProtocolViolation(model_type, message)

Thrown by check_protocol when a model breaks the protocol contract.

NLMEEstimation.QuasiNewtonEBE Type

julia

QuasiNewtonEBE()

BFGS per subject on η gradients, starting from the prior covariance Ω as the inverse-Hessian approximation (rescaled to the observed curvature at the first update), with backtracking (halving) to the Armijo condition. A subject converges once its standardized gradient ‖Lᵀg‖∞ ≤ gtol, with Ω = LLᵀ (the gradient with respect to the standardized random effects L⁻¹η), and the complete-data density is concave there.

BFGS keeps a positive-definite curvature model, so on its own it cannot follow a direction in which the density curves upward: a saddle point passes the gradient test, and next to one the steps along that direction stay tiny. The solver therefore checks the complete-data Hessian block where BFGS meets upward curvature (a curvature pair with sᵀy ≤ 0, which BFGS skips) and, once no subject is still iterating, at every subject that passed the gradient test. The block is the model's own when it provides_eta_hessian (the one the result reports, so the check adds no model calls), and otherwise forward differences of batched gradients (neta batched gradient calls per check, not reported in the result). Where -H is not positive definite, the next step is NewtonEBE's trust-region step with a radius of one prior standard deviation, halved until the density increases, and BFGS restarts from Ω.

NLMEEstimation.RaggedEBEResult Type

julia

RaggedEBEResult

Empirical Bayes result for a RaggedNLMELayout. modes[k] and hessians[k] have the actual dimension of subjects[k], without padding.

NLMEEstimation.RaggedNLMELayout Type

julia

RaggedNLMELayout(; ntheta, nsigma, family_sizes, subject_families)

A nested random-effect layout with one shared covariance per family. family_sizes[f] is the dimension of family f; subject_families[i] lists the family of each successive block in subject i's η vector. A family may occur more than once in a subject (for example, once per lesion). No absent block is padded into η. The model's subject order must agree with subject_families.

NLMEEstimation.RaggedNormalREPrior Type

julia

RaggedNormalREPrior(layout, ω)

Shared Gaussian covariance and precision for each random-effect family of a RaggedNLMELayout. Repeated blocks use the same family parameters.

NLMEEstimation.RaggedSAEMResult Type

julia

RaggedSAEMResult

SAEM estimate for a RaggedNLMELayout. eta[k] is subject subjects[k]'s compact latent vector. omega_matrix[f] is the shared covariance of random-effect family f. objective is :MAP when the model supplies a population prior and :ML otherwise. iterations and interrupted are as in SAEMResult.

NLMEEstimation.ResidualDiagnostics Type

julia

ResidualDiagnostics

Per-observation residual diagnostics from residual_diagnostics. Every field below except subjects, modes and replicates holds one vector per subject in subjects, in observation order:

observations: y.
pred: population predictions, the observation means at η = 0.
ipred: individual predictions, the means at the subject's EBE η̂ (modes).
iwres: (y - ipred) / sqrt(v(η̂)), with v the observation variance.
cpred: the population mean of the first-order linearization around η̂, μ(η̂) - J η̂ with J = ∂μ/∂η at η̂.
cwres: conditional weighted residuals: U (y - cpred). Here U is the upper Cholesky factor of the inverse of the linearized covariance C = J Ω Jᵀ + diag(v(η̂)), so UᵀU = C⁻¹ and cwres is decorrelated within the subject.
cwres_per_observation: (y - cpred) ./ sqrt.(diag(C)), each observation standardized without decorrelation; this is what nlmixr2 reports as CWRES.
epred, npd, npde: simulation-based diagnostics over replicates datasets simulated at the population parameters. epred is the simulated mean, and npd the normal quantile of the prediction discrepancy, the fraction of simulated values below y. npde is the same after decorrelation by the Cholesky factor of the inverse simulated covariance, as for cwres. pde is the decorrelated discrepancy before the quantile. A discrepancy of 0 or 1 is replaced by 1/(2K) or 1 - 1/(2K) for K replicates. These fields are nothing when replicates = 0.

Under a correct model, cwres and npde are approximately independent standard normal values.

NLMEEstimation.SAEMOptions Type

julia

SAEMOptions(; iterations=300, burnin=150, mcmc_burnin=20,
             mh_rounds=1, rw_scale=0.5, adapt_iterations=150,
             step_exponent=0.8, numeric_maxiter=12, mstep_interval=1,
             numeric_mstep=:auto,
             omega_floor=1e-8)

Settings for saem. The first burnin stochastic-approximation iterations use step size one; subsequent sizes are (k-burnin)^(-step_exponent). mcmc_burnin runs kernel sweeps at the starting population parameters without updating them. With StandardSAEMMH, each round makes one prior-independence and one joint random-walk proposal for every subject, then one random-walk proposal per η coordinate, each evaluated in one batch. The joint random-walk scale adapts toward 0.3 acceptance and each coordinate's scale toward 0.4 for the first adapt_iterations sweeps; then they are frozen. With WarmupHMCNUTS, each round makes one fixed-kernel NUTS transition per subject. omega_floor is a covariance eigenvalue floor in the M-step. With few subjects for the mu-referenced coefficients, the covariance that maximizes an M-step's Gaussian term can be singular; its closed-form update moves towards the floor by at most 100 coordinate updates per M-step, so Ω may approach it over several iterations. numeric_mstep=:auto uses an online update for remaining σ coordinates when all θ coordinates have a closed-form MuReferencing update; otherwise it uses the historical running-mixture objective. Either numerical update also moves the mu-referenced θ, holding each sampled η fixed. :online explicitly applies the online update to every remaining θ and σ coordinate: it optimizes the current complete-data conditional likelihood, then moves toward that optimizer by the SA step size, retaining one sampled latent-effect matrix. Its use for a structural θ without a declared sufficient statistic is an approximation. :mixture optimizes the exact running conditional-density mixture and retains every sampled matrix. The default mstep_interval=1 updates population parameters every SA iteration; a larger interval is an approximate, lower-frequency variant in :online mode. With DualAveragingMALA, each round makes one batched MALA proposal and updates a separate subject step size throughout the run. Set mcmc_burnin=0 to omit a separate initial MCMC phase; burnin still controls the SA step sizes.

NLMEEstimation.SAEMResult Type

julia

SAEMResult

Population estimate from saem. eta holds the final subject-chain states, one column per entry of subjects. acceptance gives prior and joint random-walk and per-coordinate random-walk acceptance fractions (prior, random_walk, coordinates); rw_scale is the final joint proposal scale. These MH-only fields are missing for other kernels. kernel identifies the method and nuts_divergences counts its fixed-kernel transitions that diverged. For MALA, mala_stepsizes and mala_acceptance hold per-subject final values; mala_adaptation_steps counts updates per subject. numeric_mstep records the resolved numerical M-step mode (:online or :mixture). objective is :ML without a model population prior and :MAP otherwise. The result is a population point estimate, not a draw from a population posterior. iterations counts the completed SA iterations; interrupted is true when a Treebars interrupt request (see progress in saem) stopped the run before SAEMOptions.iterations, so the estimate is that of the last completed iteration. The acceptance and MALA fields are NaN when no kernel sweep ran.

NLMEEstimation.StandardSAEMMH Type

julia

StandardSAEMMH()

The standard SAEM per-subject kernel: per subject and MH round, one prior-independence proposal, one covariance-scaled random-walk Metropolis proposal, and one random-walk proposal per η coordinate, each coordinate with its own adapted scale (saemix's componentwise kernel). The coordinate sweep keeps a subject's chain moving when its posterior widths differ widely across coordinates, which a single joint scale cannot follow.

NLMEEstimation.WarmupHMCNUTS Type

julia

WarmupHMCNUTS(; warmup_draws=20, warmup_evaluations=100,
               stepsize_adaptation_limit=20, max_tree_depth=8,
               target_acceptance_rate=0.8,
               curvature_floor=1e-4, metric=:diagonal)

Warm up a separate NUTS chain for each subject with WarmupHMC. Its initial covariance comes from the inverse local negative log-posterior Hessian when available, and otherwise from Ω. Subsequent SAEM rounds use fixed-kernel in-memory NUTS transitions with the current population parameters. Warmup checkpoints are temporary and removed when saem returns. The default :diagonal metric uses Hessian diagonal curvatures; metric=:dense requests the inverse regularized Hessian as a full initial covariance. target_acceptance_rate controls WarmupHMC's step-size adaptation during the one-time subject warmup; it must be strictly between zero and one.

NLMEEstimation._observation_moments Method

julia

observation_moments(model, subjects, θ, σ, H)

One (mean, variance) pair per subject in subjects order. Column k of H is the random effect for subjects[k] (entry k, a compact vector, under a RaggedNLMELayout). A batched model implements this method for a joint-program pass.

NLMEEstimation._observation_moments_and_jacobians Method

julia

observation_moments_and_jacobians(model, subjects, θ, σ, H)
    -> Vector{Tuple{μ_i, v_i, Jμ_i, Jv_i}}

Batched observation moments and first-order η Jacobians. Result entry k belongs to subjects[k], with H[:, k] as its random effect (H[k] under a RaggedNLMELayout). Each Jacobian has length(observations(model, subjects)[k]) rows and one column per η component of that subject. A joint-program model provides this directly in one batched computation; a per-subject model uses the per-subject derivative method for each entry.

NLMEEstimation._observations_batch Method

julia

observations(model, subjects::AbstractVector{<:Integer})

Observed vectors in subjects order. A batched model implements this method to obtain all observations in one joint-program query.

NLMEEstimation.ad_backend Method

julia

ad_backend(model) -> Union{Nothing, ADTypes.AbstractADType}

Optional. The automatic-differentiation backend this package uses, through DifferentiationInterface, for every derivative method the model does not implement itself. The default nothing declares no backend; calling a missing derivative method is then an error.

What this package promises and tests is reverse-mode AD with Enzyme, ADTypes.AutoEnzyme(); second-order defaults are forward-over-reverse Hessian-vector products. Other backends are passed through untested.

NLMEEstimation.apply_fixed_omega Method

julia

apply_fixed_omega(model, ω) -> Vector

The Ω parameters ω with each family declared by fixed_omega moved onto its fixed values. The estimated parts of a constrained family keep their values in ω, in standard-deviation and correlation terms: an estimated margin keeps its standard deviation, and an estimated correlation keeps the correlation of ω. Unconstrained families are unchanged.

NLMEEstimation.check_protocol Method

julia

check_protocol(model, θ, σ, ω; subjects = 1:nsubjects(model), η = nothing,
               reference_backend = nothing, reference_hessian_backend = reference_backend,
               population_gradient = false, rtol = 1e-8) -> true

Check that model implements the protocol consistently at the given parameters, and throw a ProtocolViolation naming the first violation otherwise.

Checked for the batch subjects at random effects η (zeros by default; a neta × length(subjects) matrix, one vector for all subjects, or one vector per subject):

the layout matches the lengths of θ, σ and ω, and evaluation_style is an EvaluationStyle;
the batched log-densities are finite reals and equal the per-subject ones;
the batched gradients equal the per-subject gradients, and the values returned alongside derivatives equal the plain values;
subjects are separable: changing one subject's η leaves every other subject's log-density unchanged;
when provides_eta_hessian: the Hessian blocks are square, symmetric and equal the per-subject Hessians;
a declared mu_referencing trait has one design per subject of the right size, and the density is invariant under the mu-referencing shift;
a declared observation_model reproduces the conditional log-density.
a declared Gaussian model's batched moments and first-order η Jacobians have the promised per-subject shape and agree with its plain moments;
a declared fixed_omega is valid for the layout, and ω holds its fixed values.

With reference_backend (an AD backend), the model's η gradients are also compared against that backend, and its Hessian blocks (when provided) against reference_hessian_backend (for Enzyme, a forward-over-reverse DifferentiationInterface.SecondOrder). With population_gradient = true, the optional θ/σ gradient is checked as well (against reference_backend when given).

NLMEEstimation.complete_data_loglikelihood Method

julia

complete_data_loglikelihood(model, i, θ, σ, ω_or_prior, η) -> Real

The complete-data log-density log p(y_i, η_i | θ, σ, Ω) = conditional_loglikelihood(model, i, θ, σ, η) + re_logprior(prior, η). The fifth argument is either the Ω parameter vector ω or a precomputed NormalREPrior.

NLMEEstimation.complete_data_loglikelihood_gradient_and_hessian Method

julia

complete_data_loglikelihood_gradient_and_hessian(model, i, θ, σ, ω_or_prior, η)
    -> (ℓ, f, ∇η f, ∇²η f)

The conditional log-density ℓ, the complete-data log-density f, and the gradient and Hessian of f in η_i. The conditional part comes from conditional_loglikelihood_gradient_and_hessian (AD), the prior part in closed form.

NLMEEstimation.conditional_loglikelihood Method

julia

conditional_loglikelihood(model, i, θ, σ, η) -> Real

The conditional log-density log p(y_i | θ, σ, η_i) of subject i's observations, without the random-effect prior (which this package owns, see re_logprior). θ and σ are in the model's unconstrained coordinates.

Required under PerSubjectEvaluation; derived from conditional_loglikelihoods under BatchedEvaluation.

Return -Inf where the density is zero (a failed solve, for example); estimators treat non-finite values as infeasible and report them, never as success.

NLMEEstimation.conditional_loglikelihood_and_gradient Method

julia

conditional_loglikelihood_and_gradient(model, i, θ, σ, η) -> (ℓ, ∇η ℓ)

The conditional log-density of subject i and its gradient in η_i. Implement it, or let it derive from ad_backend (per-subject models) or from the batched gradient (batched models).

NLMEEstimation.conditional_loglikelihood_and_population_gradient Method

julia

conditional_loglikelihood_and_population_gradient(model, i, θ, σ, η) -> (ℓ, ∇θ ℓ, ∇σ ℓ)

Optional, or provided by ad_backend. The conditional log-density of subject i with its gradients in θ and σ at fixed η_i. Used by population-level estimators, not by empirical_bayes.

NLMEEstimation.conditional_loglikelihood_gradient_and_hessian Method

julia

conditional_loglikelihood_gradient_and_hessian(model, i, θ, σ, η) -> (ℓ, ∇η ℓ, ∇²η ℓ)

The conditional log-density of subject i with its η gradient and neta × neta η Hessian. Available when provides_eta_hessian.

NLMEEstimation.conditional_loglikelihood_hessian_blocks Method

julia

conditional_loglikelihood_hessian_blocks(model, subjects, θ, σ, H) -> Vector{Matrix}

The per-subject η Hessian blocks of a batch: entry k is ∇²η ℓ_i for i = subjects[k] at H[:, k]. Available when provides_eta_hessian. Under BatchedEvaluation the default takes neta Hessian-vector products of the batch's summed log-density, one per η component across all subjects at once.

NLMEEstimation.conditional_loglikelihoods Method

julia

conditional_loglikelihoods(model, subjects, θ, σ, H) -> Vector

The conditional log-densities of a batch of subjects: entry k is log p(y_i | θ, σ, η_i) for i = subjects[k] and η_i = H[:, k], with H of size neta × length(subjects) for NLMELayout. For RaggedNLMELayout, H[k] has length neta(layout, subjects[k]).

Required under BatchedEvaluation: a joint-program model places each column of H into its subject's random-effect slots (the other subjects' slots are irrelevant, by conditional independence) and returns the per-subject terms. Derived from conditional_loglikelihood under PerSubjectEvaluation.

NLMEEstimation.conditional_loglikelihoods_and_gradients Method

julia

conditional_loglikelihoods_and_gradients(model, subjects, θ, σ, H) -> (ℓ::Vector, G::Matrix)

Batched values and η-gradients: G[:, k] is ∇η ℓ_i for i = subjects[k] at H[:, k]. Under BatchedEvaluation the default takes the AD gradient of sum(conditional_loglikelihoods(model, subjects, θ, σ, H)) with respect to H, which by conditional independence is every subject's gradient in one pass; implement it directly to fuse that pass with the values.

NLMEEstimation.covariance_step Method

julia

covariance_step(model, fit::FOCEResult; method = :rs, step = 1e-3,
                approximation = FOCEApproximation(interaction = fit.options.interaction))
covariance_step(model, fit::Union{SAEMResult, RaggedSAEMResult}; method = :rs,
                step = 1e-3, approximation = <see below>, rng = Random.default_rng())
covariance_step(model, θ, σ, ω; subjects = 1:nsubjects(model), free = nothing,
                η0 = nothing, method = :rs, step = 1e-3,
                approximation = FOCEApproximation(), rng = Random.default_rng(),
                progress = nothing)
    -> CovarianceResult

The covariance step: the sampling covariance and standard errors of the estimated population coordinates (θ, σ, ω), in the unconstrained coordinates in which they are fitted.

R is minus the Hessian of the fitted objective and S the sum over subjects of the outer products of their marginal log-likelihood gradients. method = :rs (the default) returns the sandwich R⁻¹ S R⁻¹, which remains consistent when the model's random-effect or residual distribution is misspecified; :r returns the inverse observed information R⁻¹, and :s the inverse cross-product S⁻¹. With a population_logprior, R includes the prior and R⁻¹ is the Laplace approximation of the population posterior covariance at the MAP point.

Both matrices are central finite differences of the per-subject marginal log-likelihoods of approximation, which needs objective values only: no population derivatives of the model. Each coordinate j moves by step * max(1, |x_j|). :s takes 2p + 1 evaluations for p estimated coordinates; :r and :rs take 2p² + 1. Every evaluation warm-starts the subject modes from the same solution at the estimate.

For a FOCE(I) fit, the defaults reuse its subjects, its fitted coordinates (fit.free) and its interaction setting, so the step differentiates the fitted objective. For an SAEM fit, all (θ, σ, ω) are estimated. A model with GaussianObservations and an NLMELayout defaults to the FOCE(I) linearization at the SAEM estimate. Any other model, including every RaggedNLMELayout, defaults to ImportanceSampling. Its draws, from rng, are held fixed across the difference points, so the estimated marginal is a smooth function of the coordinates. AdaptiveGaussHermite is the deterministic alternative for few random effects per subject. For those two approximations, perturbations of ω and of mu-referenced θ reuse the conditional values and cost no model evaluations. free selects positions in the estimated coordinates [θ; σ; free Ω coordinates], which are the concatenated (θ, σ, ω) unless the model declares fixed_omega (see CovarianceResult). Fixed covariance parts carry no uncertainty and are never differentiated; ω must hold their values. Throws CovarianceStepError when the requested matrix is not positive definite.

progress (a Treebars node or backend symbol) gets a Covariance step node: the subject modes at the point, then a Marginal likelihoods counter over the 2p² + 1 (or 2p + 1) evaluations. An interrupt raises Treebars.ProgressInterrupted.

NLMEEstimation.default_ebe_method Method

julia

default_ebe_method(model) -> EBEMethod

NewtonEBE if provides_eta_hessian, else QuasiNewtonEBE.

NLMEEstimation.delta_method Method

julia

delta_method(result::CovarianceResult, f; jacobian = nothing)
    -> (estimate, covariance, standard_errors)

First-order (delta-method) uncertainty of natural-unit quantities f(x), a real vector function of the full population vector x = [θ; σ; ω] (result.point). jacobian(x) may supply ∂f/∂x as a length(f(x)) × length(x) matrix; otherwise the Jacobian is taken by central differences. Coordinates outside result.free, and covariance parts held by fixed_omega, carry no uncertainty; the Jacobian is taken with respect to result.coordinates, mapped through the fixed parts. For example, with nθ = length(θ) and nσ = length(σ), the covariance of exp.(θ) and of the Ω entries is

julia

delta_method(result, x -> exp.(x[1:nθ]))
delta_method(result, x -> omega_elements(layout, x[nθ+nσ+1:end]))

NLMEEstimation.ebe_covariances Method

julia

ebe_covariances(result::EBEResult) -> Vector{Matrix}

The Laplace (Gaussian) approximation of each subject's η posterior covariance, (-H_i)⁻¹, from the Hessians at the modes. Throws if the result carries no Hessians or some -H_i is not positive definite.

NLMEEstimation.empirical_bayes Method

julia

empirical_bayes(model, θ, σ, ω; subjects = 1:nsubjects(model), η0 = nothing,
                method = default_ebe_method(model), options = EBEOptions(),
                require_convergence = true, progress = nothing) -> EBEResult

Empirical Bayes estimates (posthoc): for fixed population parameters θ, σ and Ω parameters ω, the per-subject mode

julia

η̂_i = argmax_η  log p(y_i | θ, σ, η) + log N(η | 0, Ω).

All subjects are solved in lockstep: each iteration makes one batched gradient call (conditional_loglikelihoods_and_gradients) for every subject still iterating, plus batched value calls for its trial steps, so a joint-program model pays O(N) per iteration rather than O(N) per subject. Models are never called concurrently.

η0 is the starting point: nothing (zeros), a neta × length(subjects) matrix, one vector per subject, or an earlier EBEResult covering these subjects (warm start). method is NewtonEBE or QuasiNewtonEBE.

Throws an EBEConvergenceError if some subject's mode is not found, unless require_convergence = false, in which case the result records the failure in converged. A non-finite complete-data log-density at the starting point is always an error. A trial point at which the model returns -Inf or throws a DomainError, OverflowError or PosDefException (a point outside its domain) is rejected, and a shorter step is tried.

Evaluating the complete-data log-density f can round off by many ulps of f when f sums terms much larger than itself, such as observations far from zero with a small residual SD. Near a mode, a full step whose predicted increase is below √eps · max(1, |f|) is therefore accepted unless f decreases by more than that. The convergence tests use gradients and Hessian blocks, which stay accurate there.

progress is a Treebars node (or a backend symbol such as :term) under which an Empirical Bayes (Newton) or (BFGS) node reports the solver's iteration and how many subjects have converged. Treebars.request_interrupt! on it raises Treebars.ProgressInterrupted at the next iteration.

NLMEEstimation.evaluation_style Method

julia

evaluation_style(model) -> EvaluationStyle

Optional. Defaults to PerSubjectEvaluation.

NLMEEstimation.fixed_omega Method

julia

fixed_omega(model) -> nothing, or family => FixedOmega pairs

Optional declaration that some random-effect covariance families of model are known, wholly or in part. Return nothing (the default) when every family is estimated, or a collection of family => FixedOmega(…) pairs (a vector of pairs, a single pair, or a dictionary). A family is an entry of the layout's eta_blocks (NLMELayout) or family_sizes (RaggedNLMELayout), by position.

focei and saem estimate only the free parts and keep the fixed parts at their values, so the fitted omega and omega_matrix report them. Their starting ω0 must already hold the fixed values; apply_fixed_omega builds such a vector. empirical_bayes evaluates at whatever ω it is given, and check_protocol verifies that its ω holds the declared values.

NLMEEstimation.focei Method

julia

focei(model, θ0, σ0, ω0; subjects = 1:nsubjects(model), free = nothing,
      options = FOCEOptions(), progress = nothing) -> FOCEResult

Estimate population parameters by minimizing the FOCE(I) OFV with NLopt's BOBYQA, a quadratic-interpolation trust-region algorithm. The model's θ and σ and the package's log-Cholesky ω are already unconstrained. Each objective call warm-starts subject EBEs from the preceding successful call. No population derivatives or η Hessians are taken. free may select one-based positions in the concatenated (θ, σ, ω) vector; other coordinates stay at their starting values. With a population prior, the outer objective is the MAP point's negative twice log posterior; the result still reports its likelihood separately. Covariance families declared by fixed_omega keep their fixed parts. The search covers the estimated standard deviations and correlation of a partially fixed family; a fully fixed family is not searched. free may not name a position inside a constrained family's ω segment. ω0 must hold the fixed values (apply_fixed_omega).

progress (a Treebars node or backend symbol) gets a FOCE(I) node counting objective calls against maxeval, with the best OFV so far as its message and each call's mode search as a transient child. Treebars.request_interrupt! on it stops the search at the current objective call and returns the best point so far, with converged = false and optimizer_info.status == :FORCED_STOP. An interrupt during the evaluation at θ0, σ0, ω0, before any search, raises Treebars.ProgressInterrupted.

NLMEEstimation.focei_evaluate Method

julia

focei_evaluate(model, θ, σ, ω; subjects = 1:nsubjects(model), η0 = nothing,
               interaction = true, ebe_method = QuasiNewtonEBE(),
               ebe_options = EBEOptions(), progress = nothing)

Evaluate the first-order conditional (with interaction by default) marginal log-likelihood. The EBEs maximize the exact complete-data log-density; their Laplace curvature uses only the first η derivatives of the Gaussian observation means and variances. Supply η0 as for empirical_bayes to warm-start the mode search; ebe_method and ebe_options configure it. A model's population_logprior, if present, is added once after the subject likelihood terms; it does not alter the subject EBEs. progress (a Treebars node or backend symbol) gets a FOCE(I) objective node with the mode search's node under it; an interrupt raises Treebars.ProgressInterrupted.

NLMEEstimation.hessian_ad_backend Method

julia

hessian_ad_backend(model) -> Union{Nothing, ADTypes.AbstractADType}

Optional. The backend used for default second-order derivatives (η Hessian blocks). Defaults to ad_backend. Return nothing to declare that the model provides no second order, so estimators use their gradient-only paths.

NLMEEstimation.marginal_loglikelihood Method

julia

marginal_loglikelihood(model, fit; approximation = ImportanceSampling(),
                       rng = Random.default_rng(), nobs = nothing, progress = nothing)
marginal_loglikelihood(model, θ, σ, ω; subjects = 1:nsubjects(model), η0 = nothing,
                       approximation = ImportanceSampling(), rng = Random.default_rng(),
                       nparameters = nothing, nobs = nothing, progress = nothing)
    -> MarginalLikelihoodResult

The marginal log-likelihood Σᵢ log ∫ p(y_i | θ, σ, η) N(η; 0, Ω) dη at a fit's estimate or at given population coordinates, with OFV, AIC and BIC. This gives SAEM fits a likelihood-based objective. approximation is ImportanceSampling, AdaptiveGaussHermite or FOCEApproximation. The first two need only conditional values and η gradients, so they apply to any likelihood and to ragged layouts. Subject modes warm-start from the fit's EBEs or final chain states (η0). nparameters, the AIC/BIC parameter count, defaults to the estimated coordinates: θ, σ and the Ω coordinates that fixed_omega leaves free, or fit.free for a FOCE(I) fit. nobs overrides the observation count used by BIC. progress (a Treebars node or backend symbol) gets a Marginal likelihood node with the subject modes and, while they are evaluated, a counter over the sampling or quadrature points; an interrupt raises Treebars.ProgressInterrupted.

NLMEEstimation.mu_referencing Method

julia

mu_referencing(model) -> Union{Nothing, MuReferencing}

Optional trait. Declares which θ coordinates enter the model as typical values of the random effects, ϕ_i = X_i θ_μ + η_i. SAEM's closed-form M-step needs it. Default: nothing (not declared).

NLMEEstimation.neta Method

julia

neta(layout::NLMELayout) -> Int

Dimension of the per-subject random-effect vector η_i.

NLMEEstimation.nlme_layout Function

julia

nlme_layout(model) -> NLMELayout or RaggedNLMELayout

Required. The parameter layout: lengths of θ and σ and the η block structure. See NLMELayout for the coordinate conventions.

NLMEEstimation.nomega Method

julia

nomega(layout::NLMELayout) -> Int

Length of the package-owned Ω parameter vector ω (see NLMELayout).

NLMEEstimation.nsubjects Function

julia

nsubjects(model) -> Int

Required. Number of subjects (independent units carrying their own η_i).

NLMEEstimation.observation_model Method

julia

observation_model(model) -> Union{Nothing, GaussianObservations}

Optional trait. Declares the observation-model structure FOCE(I) linearizes (first order: the Jacobians of the predicted means and variances in η). Default: nothing.

NLMEEstimation.observation_moments Function

julia

observation_moments(model, i, θ, σ, η) -> (mean::AbstractVector, variance::AbstractVector)

Predicted means and variances of subject i's observations, differentiable in η; required when observation_model is declared.

NLMEEstimation.observation_moments_and_jacobians Method

julia

observation_moments_and_jacobians(model, i, θ, σ, η) -> (μ, v, Jμ, Jv)

One subject's observation means, variances, and their Jacobians in η. Jacobian rows follow observation order. The default derives them with reverse-mode AD from observation_moments; a model may implement this directly.

NLMEEstimation.observations Function

julia

observations(model, i) -> AbstractVector

Subject i's observations; required when observation_model is declared.

NLMEEstimation.omega_cholesky Method

julia

omega_cholesky(layout, ω) -> LowerTriangular

Block-diagonal lower-triangular Cholesky factor L of Ω = L Lᵀ, built from the unconstrained Ω parameters ω (see NLMELayout). Generic in the element type of ω, so it can be differentiated.

NLMEEstimation.omega_elements Method

julia

omega_elements(layout, ω) -> Vector

The entries of Ω on the covariance scale, in one-to-one correspondence with the log-Cholesky coordinates ω: block by block (family by family for a RaggedNLMELayout), the lower triangle in row-major order. Diagonal positions are variances and the others covariances. Pair it with delta_method for standard errors of Ω.

NLMEEstimation.omega_matrix Method

julia

omega_matrix(layout, ω) -> Matrix

The random-effect covariance Ω = L Lᵀ for the unconstrained Ω parameters ω.

NLMEEstimation.omega_parameters Method

julia

omega_parameters(layout, Ω::AbstractMatrix) -> Vector{Float64}

Inverse of omega_matrix: the unconstrained Ω parameters ω of a covariance matrix Ω. Throws if Ω has nonzero entries outside the layout's η blocks or if a block is not positive definite.

NLMEEstimation.parameter_names Method

julia

parameter_names(model) -> (theta = Vector{String}, sigma = Vector{String}, eta = Vector{String})

Optional. Human-readable coordinate names for reporting. Defaults to θ1, θ2, …, σ1, …, η1, ….

NLMEEstimation.population_logprior Method

julia

population_logprior(model, θ, σ, ω) -> Union{Nothing, Real}

Optional log prior on the population coordinates, evaluated once per fit, not once per subject or retained SAEM draw. Return nothing (the default) for maximum likelihood. A model with a prior returns its log density with respect to the unconstrained (θ, σ, ω) coordinates. If the original prior is defined on transformed parameters, the model must include the transformation's log absolute Jacobian here. ω is the package's log-Cholesky coordinate vector. The random-effect prior η ~ N(0, Ω) is already owned by this package and must not be repeated here.

NLMEEstimation.population_logprior_and_gradient Method

julia

population_logprior_and_gradient(model, θ, σ, ω)
    -> (logprior, ∇θ logprior, ∇σ logprior, ∇ω logprior)

Optional gradient for a model's population_logprior. The default uses ad_backend(model) when a prior is present and returns zero gradients when there is none. Models may implement this method directly, including when their conditional-likelihood gradient uses a different backend.

NLMEEstimation.population_prior_provenance Method

julia

population_prior_provenance(model)

Optional, model-supplied description of the population prior. Returned unchanged in SAEM and FOCE results when a prior is present; defaults to nothing.

NLMEEstimation.posthoc Function

julia

posthoc(model, θ, σ, ω; kwargs...)

Alias of empirical_bayes.

NLMEEstimation.provides_eta_hessian Method

julia

provides_eta_hessian(model) -> Bool

Optional. Whether η Hessian blocks are available (conditional_loglikelihood_hessian_blocks). Defaults to whether a hessian_ad_backend is declared; a model implementing the Hessian methods itself returns true. Second order is only ever requested as per-subject neta × neta blocks, never as a joint matrix.

NLMEEstimation.re_logprior Method

julia

re_logprior(prior::NormalREPrior, η) -> Real
re_logprior(model, ω, η) -> Real

The random-effect prior log-density log N(η | 0, Ω).

NLMEEstimation.re_logprior Method

julia

re_logprior(prior::RaggedNormalREPrior, i, η)

Log density for subject i's compact η vector, summing one Normal density per listed family block.

NLMEEstimation.re_logprior_gradient_and_hessian Method

julia

re_logprior_gradient_and_hessian(prior::NormalREPrior, η) -> (lp, ∇η lp, ∇²η lp)

The prior log-density with its gradient -Ω⁻¹η and Hessian -Ω⁻¹ in η, in closed form (exact).

NLMEEstimation.re_prior Method

julia

re_prior(model, ω) -> NormalREPrior or RaggedNormalREPrior

The η prior of model at Ω parameters ω.

NLMEEstimation.residual_diagnostics Method

julia

residual_diagnostics(model, fit; replicates = 1000, rng = Random.default_rng(),
                     progress = nothing)
residual_diagnostics(model, θ, σ, ω; subjects = 1:nsubjects(model), η0 = nothing,
                     replicates = 1000, rng = Random.default_rng(),
                     ebe_options = EBEOptions(), progress = nothing) -> ResidualDiagnostics

PRED, IPRED, IWRES, CWRES and NPDE for a model with GaussianObservations, at a fit's estimate or at given population coordinates. Subject modes (EBEs) warm-start from the fit's EBEs or final chain states, or from η0. Both NLMELayout and RaggedNLMELayout are supported.

CWRES follows Hooker, Staatz and Karlsson (2007): the observations of subject i are linearized around its EBE, giving mean μ(η̂) - J η̂ and covariance J Ω Jᵀ + diag(v(η̂)), and the residual is decorrelated by the Cholesky factor of that covariance's inverse. This is the approximation the FOCE(I) objective of focei_evaluate uses, with the residual variance at the EBE.

NPDE follows Brendel, Comets, Laffont, Laveille and Mentré (2006) and Comets, Brendel and Mentré (2008). replicates datasets are simulated at the population parameters (η ~ N(0, Ω), then y ~ N(μ(η), v(η))), each costing one batched observation_moments call. The observed and simulated vectors are decorrelated with the Cholesky factor of the inverse of each subject's simulated covariance. replicates must exceed every subject's number of observations; replicates = 0 skips the simulation.

progress (a Treebars node or backend symbol) gets a Residual diagnostics node with the subject modes and a Simulated replicates counter; an interrupt raises Treebars.ProgressInterrupted.

NLMEEstimation.saem Method

julia

saem(model, θ0, σ0, ω0; rng=Random.default_rng(), subjects=nothing,
     η0=nothing, options=SAEMOptions(), kernel=StandardSAEMMH(),
     progress=nothing) -> SAEMResult

Estimate population parameters using MCMC stochastic approximation. Subject chains persist across iterations. The Gaussian random-effect prior supplies a closed-form covariance M-step; MuReferencing also supplies the generalized-least-squares fixed-effect M-step. Remaining θ and σ coordinates use a first-order AD numerical M-step, which holds each chain's η fixed and so also moves the mu-referenced θ: a σ coordinate that trades off against a random effect then moves together with it. Its quasi-Newton search starts from a difference estimate of the objective's curvature, carried from one M-step to the next. By default, σ coordinates are updated online only when all structural θ coordinates have a MuReferencing sufficient-statistic update. Otherwise the historical running conditional-density objective is used and retains one latent-effect matrix per SA iteration. Set numeric_mstep=:online to force an online update for undeclared structural θ as a model-specific approximation. SAEMResult.numeric_mstep records the selected mode. kernel selects the standard lockstep MH proposals, batched continuously adapted MALA proposals, or per-subject in-memory WarmupHMC NUTS transitions. When population_logprior is supplied, the M-step maximizes the same objectives plus the prior, added once. The mu-referenced θ and Ω move on from their closed-form update to the maximizer with the prior, a search that does not evaluate the model, and the numerical update of the remaining θ and σ includes the prior; so a prior adds no model evaluations to the step. The model supplies a population-prior gradient or an AD backend from which to derive it. Covariance families declared by fixed_omega keep their fixed parts: a fully fixed family is not updated, and a partially fixed family maximizes the expected Gaussian random-effect density over its estimated standard deviations and correlation. ω0 must hold the fixed values (apply_fixed_omega); omega_floor applies to estimated parts only. A model rejects a point outside its domain by returning a log density of -Inf or by throwing a DomainError, OverflowError or PosDefException. The MH and MALA kernels then reject that subject's proposal alone, and the numerical and MAP M-steps halve their line-search step. The point each M-step starts from must evaluate finitely.

progress is a Treebars node (or a backend symbol such as :term) under which a SAEM node shows the run's phases: NUTS warmup (one transient WarmupHMC tree per subject), MCMC burn-in and Stochastic approximation, each counting its iterations. The running message gives the stage (M-step while it runs), the SA step size γ, the chains' summed conditional log-density at the current population estimate (MH and NUTS) and the kernel's acceptance, divergences or mean step size. Treebars.request_interrupt! on that node stops the run at the next iteration boundary: saem returns the estimate of the last completed iteration with interrupted = true and iterations counting the completed SA iterations.

NLMEEstimation.saem_population_gradient_style Method

julia

saem_population_gradient_style(model[, numeric_mode]) -> EvaluationStyle

Optional. Selects how SAEM computes the conditional-likelihood value and population-parameter gradient for its numerical M-step. Defaults to evaluation_style. A model with BatchedEvaluation can return PerSubjectEvaluation here when it implements an efficient conditional_loglikelihood_and_population_gradient for one subject; SAEM sums those subject values and gradients. This choice affects only that M-step gradient. Likelihood, η-gradient, and η-Hessian evaluation continue to use evaluation_style(model).

The optional second argument is Val(:online) or Val(:mixture), the resolved SAEM numerical M-step mode. Override it when one gradient route is efficient for online updates but costly for the retained-draw mixture. By default both modes use the one-argument style, preserving existing model methods.

NLMEEstimation.ReferenceModels Module

julia

NLMEEstimation.ReferenceModels

Small NLME models implementing the protocol, for tests, examples and cross-checks of other implementations of the same model. Not part of the protocol.

NLMEEstimation.ReferenceModels.AdditiveError Type

julia

AdditiveError()

Additive residual error: σ = [log a], variance a².

NLMEEstimation.ReferenceModels.CombinedError Type

julia

CombinedError()

Combined additive and proportional residual error: σ = [log a, log b], variance a² + b² C² at predicted concentration C.

NLMEEstimation.ReferenceModels.OneCompartmentOral Type

julia

OneCompartmentOral(subjects; backend, error_model = CombinedError(),
                   eta_blocks = [1, 1, 1], reference_weight = 70.0)

One-compartment model with first-order absorption and elimination after a single oral dose, solved analytically. With θ = [log ka, log CL, log V, β_WT]:

julia

log ka_i = θ[1] + η_i[1]
log CL_i = θ[2] + β_WT · log(WT_i / reference_weight) + η_i[2]
log V_i  = θ[3] + η_i[3]
C_i(t)   = dose_i · ka_i / (V_i (ka_i − ke_i)) · (exp(−ke_i t) − exp(−ka_i t)),  ke_i = CL_i / V_i
y_ij     ~ N(C_i(t_ij), residual variance)

Implements the full protocol, including the MuReferencing and GaussianObservations traits. Derivatives come from backend (an ADTypes backend whose package is loaded, for example ADTypes.AutoEnzyme() with Enzyme). eta_blocks sets the Ω block structure ([3] for a dense Ω).

NLMEEstimation.ReferenceModels.PKSubject Type

julia

PKSubject(dose, weight, times, concentrations)

One subject of a single-oral-dose study: dose amount, body weight, sampling times and observed concentrations.

NLMEEstimation.ReferenceModels.ResidualError Type

Residual-error model of OneCompartmentOral.

NLMEEstimation.ReferenceModels.simulate_one_compartment_oral Method

julia

simulate_one_compartment_oral(rng, θ, σ, Ω; error_model = CombinedError(),
    nsubjects = 20, dose = 320.0, weight_range = (50.0, 90.0),
    times = [0.25, 0.5, 1, 2, 3.5, 5, 7, 9, 12, 24], reference_weight = 70.0)
    -> (subjects::Vector{PKSubject}, η::Vector{Vector{Float64}})

Simulate a study from OneCompartmentOral with η_i ~ N(0, Ω) and body weights uniform on weight_range. Returns the subjects and the true random effects.