nlmixr2est 7.1.0

Breaking changes

• saemControl(lbfgsLmm=, lbfgsFactr=, lbfgsPgtol=, lbfgsMaxIter=) have been removed and now error as unused arguments. They were announced in 7.0.2 as controlling a bounded L-BFGS-B refinement of the fixed-effect-only (phi0) parameters of a general log-likelihood model, but no such refinement was ever implemented: the options were validated and stored and then read by nothing. That phi0 step is optimized by the bounded derivative-free routine (bobyqa, or stats::optimize for a single parameter), which honors the ini-block bounds and takes no L-BFGS-B settings. Passing any of the four never changed a fit, so removing them changes no result.

Changed defaults

• est="saem" now uses nu=c(4,4,4) when nu is not set and a residual error component depends on an eta (for example a <- add.sd * exp(eta.sd) with cp ~ add(a) or cp ~ add(a) + dnorm()). At the default c(2,2,2) those chains mixed too slowly and the fit stopped short of the optimum; an explicit saemControl(nu=) is used as given.

• est="saem" now refines a population theta that carries no random effect with newuoa over all such thetas at once, under a budget of 25 objective evaluations per iteration (saemControl(nonMuThetaOpt="newuoa", nonMuThetaMaxEval=25)), rather than with sweeps of golden-section coordinate descent. That refinement is where such a model spends most of its time – each of its evaluations re-solves every subject and chain – and the sweeps solved it far more precisely than a stochastic-approximation step that then moves a fraction of the way there can use. Fits of models that have a non-mu theta will report slightly different estimates. Measured by the FOCEi conditional objective at each run’s converged estimates: on the nimotuzumab target-mediated model 1.6x faster at an indistinguishable objective (143.53 vs 143.48), and on the mavoglurant PBPK model 1.9x faster at a clearly better one (1977.0 vs 2055.8), where the cheaper refinement escapes a poor additive-error basin the old default settles into. Pass saemControl(nonMuThetaOpt="optimize") for the previous behavior.

New features

• est="saem" now records its MCMC mixing diagnostics on the fit instead of discarding them. saem computed an acceptance rate, used it for nothing and threw it away, so a chain that had stopped moving looked exactly like one exploring properly. One row per iteration:

  • $mcmcAccept, pooled acceptance per kernel (prior, rw, coord);
  • $mcmcAcceptCol, kernel 3’s acceptance PER SAMPLED PARAMETER. Kernel 3 is Metropolis-within-Gibbs – it proposes one coordinate at a time – so its acceptance is already a per-column quantity, and pooling it hides one parameter behaving differently from the rest. A column near 1 while the others sit near their target is a coordinate whose proposals the likelihood is not rejecting, whose chain is then exploring the prior rather than the posterior, and whose omega therefore measures the wandering;
  • $mcmcStuck, the fraction of SUBJECTS that accepted nothing that iteration. A healthy-looking pooled 0.3 is equally consistent with everyone at 0.3 and with half the population never moving at all;
  • $mcmcPhiSd and $mcmcPhiAcf, the pooled SD and the lag-1 autocorrelation of each sampled parameter. An autocorrelation of 1.0 means the chain did not move.

  Diagnostics only: they accumulate counters and read the sampled parameters, and change no fitted value. Asserted against pinned estimates.

• setCov() is now an S3 generic dispatched on the covariance method, so other packages (for example SIR or bootstrap) can add a method with setCov.<method>(); setCovAllMethods() lists them. A method’s own options are changed with control=, holding only covariance options: rsControl() for "r,s"/"r"/"s", saControl() for "sa" and impCovControl() for "imp". Each covariance records the options it was computed with (fit$env$covOptions), and setCov() reuses a cached covariance only when the requested options match, recomputing otherwise.

• The cache key of a covariance method’s options is the new generic setCovOptions(), so a method whose covariance depends on more than its control (for example SIR, seeded from another covariance on the fit) can put that state in the key.

• setCov(fit) <- value installs an already-computed covariance: a matrix (as "user", or setCov(fit, "name") <- matrix), or any result with a setCovValue() method, which also records the options that produced it and can keep the result itself on the fit.

• The fit print separates the other calculated covariances with ;, since "r,s" contains a comma.

• est="imp", "impmap" and "qrpem" with nIter=0 evaluate the fit at the supplied parameters: one E-step and no M-step (like NONMEM EONLY=1), with the importance-sampling objective in $impObj. $runInfo notes the E-step-only run and where the starting etas came from (#1091, #1092).

• setOfv(fit, "imp") and setOfv(fit, "impmap") add an importance-sampling objective to any fit through that E-step-only run at the fit’s estimates.

• est="vae" groups near-interchangeable covariates into colinearity clusters, controlled by the new vaeControl(covSelectColinearCut=) (default 0.9). A cluster never restricts what may be selected. It does two things: the covariate M-step keeps the previous iteration’s choice unless a cluster mate beats it by a full covariate’s L0 cost, so the selection stops chattering between columns the design cannot tell apart; and the mates that came within that margin are reported in the fit’s $covNearTie. Clusters are a coarsening of the mutual-exclusion groups, so two shapes of one covariate never cluster together.

• vaeCovariates() reports the same clustering in a new cluster column and takes the threshold as colinearCut.

• est="vae" now refines covariate attribution across correlated latent dimensions. Each dim’s covariate search only sees the other dims through a frozen Gauss-Seidel offset, so it cannot notice that a covariate on one dim would be better explained on a correlated one. A new pass groups the dims by the empirical correlation of the posterior means, scores joint moves through an exact group-restricted GLS, and only writes back a move that strictly improves the group. It is gated on a correlated omega: with a diagonal omega the objective is separable and each per-dim search is already exact, so the pass reports the correlated dims in $runInfo (advising you to declare the omega block) instead of running. Controlled by vaeControl(covSelectPhiCor=, covSelectPhiJoin=, covSelectPhiLeave=, covSelectPhiMaxDim=); the counters and the sticky pair adjacency are reported in the fit’s $vae.

• est="saem" now fits residual error components that are modeled rather than estimated directly, such as a <- add.sd*exp(eta.sd); cp ~ add(a) or a <- add.sd + WT*cov.sd; cp ~ add(a). These endpoints are fit as the equivalent cp ~ add(a) + dnorm() log-likelihood, and $runInfo notes the promotion.

• est="saem" estimates every theta without an eta that informs a general likelihood (dnorm(), t(), cauchy(), the discrete and continuous densities, and ll()), including a boxCox()/yeoJohnson() lambda, through a temporary mu-referenced eta on the scale of its range: exp() for a positive parameter such as a standard deviation or degrees of freedom, expit() for a probability, additive when unbounded. The theta is reported as its back-transformed theta + mean(eta), the temporary eta is removed from the fit, and $runInfo lists the thetas that received one. These parameters were previously left near their initial values.

• Added a native analytical outer Hessian for fast Gaussian FOCE/FOCE+/FOCEI/AGQ fits, using the existing sensitivity pool. Fast nlminb fits used it automatically.

• Added optional full conditional inner curvature for fast Gaussian FOCEI via innerHessian="conditional", used by inner trust and n1qn1’s warm="calc" seed. The FOCEI marginal objective was unchanged.

• Added foceiControl(detHessian="conditional"): the objective’s Laplace log-determinant uses the full conditional Hessian (the observed information at the conditional mode) instead of FOCEI’s Gauss-Newton expected information. The analytic outer gradient carries the matching third-order term, probed the way the analytic outer Hessian already does (2 solves per eta); the analytic covariance and outer Hessian fall back to finite differences under it. On theo_sd this objective sits 0.4 above the exact (AGQ) marginal -2 log-likelihood where FOCEI sits 1.3 below, at about 1.75x the FOCEI fit time (#1068).

• Evaluated the conditional inner value, gradient and full curvature jointly in one pooled sensitivity solve, including M2/M3/M4 censoring.

• Added foceiControl(outerOpt="trust"), a trust-region Newton outer optimizer (RcppTrust) driven by the analytical outer Hessian. outerTrustHessian= selects the curvature – the analytical Hessian under fast=TRUE, a damped BFGS update, or a finite difference of the outer gradient – with outerTrustRinit/outerTrustRmax, outerTrustFterm/outerTrustMterm, outerTrustRelStep and outerTrustRestarts controlling the region, its tolerances and the step handed to the Hessian. Because the solver’s own convergence test is satisfied by a collapsing trust region, the reported point is checked with its Newton decrement and the region re-entered when it is not stationary. Measured on one model only (theo_sd, a fast FOCEi fit of the one-compartment ODE): 116.807191 against outerOpt="nlminb"’s 116.808709, at comparable cost once the model cache is warm.

• Added est="flaplace", "mflaplace", "iflaplace", "fagq", "mfagq" and "ifagq" – the Laplace and adaptive-quadrature methods (plus their mu-referenced "lin"/"irls" variants) run with the full conditional curvature (fast=TRUE, innerHessian="conditional", detHessian="conditional"): both the inner optimizer’s curvature and the objective’s Laplace log-determinant, which is what Full names. AGQ places its nodes with that same curvature, so the quadrature and the determinant cannot disagree. They report as Full Laplace/Full AGQ, and require Gaussian endpoints. The analytic outer gradient carries the determinant’s third-order term for Laplace; under AGQ it declines to the finite-difference gradient, which still needs the node spread’s derivative.

• Both shapes of a focei covariance are named and cached, so setCov() can swap between them. foceiControl(covFull=) decides whether fit$cov is the structural-theta block or the full theta + residual sigma + Omega matrix; the full shape is now reported as "r,s (full)", "r (full)", "s (full)" or "analytic (full)", and the theta-only shape keeps the unqualified name. A fit computes both, so the shape it does not install is kept in fit$covList and setCov() reinstalls it directly rather than recomputing it. On the finite-difference path the two are different estimators – "s" inverts the theta block of the cross-product while "s (full)" takes the theta block of the full inverse, which also carries the Omega estimation uncertainty – so their standard errors differ; on the analytic path the assembly is always full and the theta standard errors agree. A default focei fit now reports "r,s (full)" where it reported "r,s".

• impmapControl(proposal=) selects the importance-sampling proposal family for est="imp", "impmap" and "qrpem": "normal" and "t" as df already reached, plus "laplace" (a spherical multivariate Laplace, whose exponential tail dominates the joint target’s so the importance weights are bounded by construction) and "mixture" (a defensive scale mixture about the same mode, set by propMixScale/propMixWeight). The default "auto" resolves to the historical df behaviour, so existing fits are unchanged. The resolved family is reported in fit$env$impProposal and per subject in fit$env$impPropInd.

• impmapControl(qrScramble = "owen" | "lms") scrambles the QRPEM Sobol point set (qr = TRUE). The Cranley-Patterson shift randomizes the set but leaves the correlation structure between the sequence’s high-order dimensions intact, so it helps least on the many-random-effect models that need it most; a nested uniform (Owen) or linear matrix scramble permutes the digits and breaks it. Scrambling replaces the shift rather than composing with it, and its key is derived arithmetically from impSeed and the (iteration, subject, dimension) indices, so it draws nothing from the RNG and the fit stays reproducible and independent of the thread count. Defaults to "none".

• impmapControl(nBurn=, burnFreezeOmega=) add burn-in EM iterations to the importance-sampling family (est="imp", "impmap", "qrpem"). They run before the nIter budget rather than out of it, let the gamma and auto controllers settle, and can hold Omega at its starting value while the structural and residual-error thetas update. Convergence is not tested until the whole trailing nConvWindow lies past the burn-in, so a frozen Omega cannot be mistaken for a settled one. Both default off.

• saemControl(iovMethod = "twoLevel") estimates inter-occasion variability the way the rest of saem estimates a variance. The shared pre-processing rewrite that every estimation method uses carries the occasion magnitude as a population parameter multiplying per-occasion unit-variance etas, which makes it non-mu-referenced – so saem had to estimate a variance through its fixed-effect-only path, a stochastic sampled mean over draws whose pseudo-variance is deliberately annealed followed by a bounded direct optimization, while every other variance component gets a closed-form M-step. In practice the estimate collapsed toward zero. "twoLevel" writes the occasion term out as a second variance component instead, following Panhard and Samson (2009): one zero-mean eta per occasion level entering additively, with the per-occasion variances constrained equal so they estimate the single inter-occasion variance the model declares. A model the newer handling does not cover (more than one occasion variable, a correlated occasion term, an occasion parameter that is not mu-referenced) falls back to the shared rewrite and says so in the fit’s $runInfo. Either way the fit presents the same: $omega split into $id and $occ, the iov.x ~ v | occ row restored in $ui, and an $iov table of the per-occasion deviations. $cov carries one row for the occasion variance, om.iov.x, rather than one per occasion level under an internal name – the per-occasion columns estimate a single variance, so they are contracted by averaging (the delta method for Psi = mean(v_1, ..., v_K)).

  iovMethod = "collapsed" is a third, opt-in setting that additionally uses the paper’s own sampler: one parameter per occasion carrying mu + b_i + c_ik together under a compound-symmetric prior, rather than sampling b_i and c_ik separately. It targets the same estimates – the shared mean and the compound-symmetric block are both imposed exactly, not by projection – and differs only in how the chain mixes. Note the Gaussian-quadrature objective is unreliable for it when the inter-occasion variance is much smaller than the between-subject one, because the prior is then nearly degenerate along b and an axis-aligned grid covers that badly; compare fits on the estimates rather than on objf.

  On the paper’s own design (n = 24, 1000 replicates) "collapsed" recovers the inter-occasion variances less accurately than the default – mean absolute relative bias 7.5% against 4.2%, and -16.1% on one component where the default gets -6.8% – which is why it is opt-in. It does run somewhat faster.

  "twoLevel" is now the default, so saem IOV estimates change. They were biased toward zero, badly. On the paper’s own simulation design (n = 24, inter-occasion variances of 0.0025/0.01/0.01, 1000 replicates each) the shared rewrite recovers them with -94.7%/-89.0%/-43.6% relative bias – it collapses them toward zero, and the variance it does not capture instead inflates the residual error by 25%. The two-level handling gets -2.2%/-6.8%/-3.6%, less biased than the figures Panhard and Samson report for their own implementation (-8.7%/-10.9%/-5.4%), and leaves the residual error 2.9% high. iovMethod = "theta" restores the old behavior.

• Inter-occasion variability can now be correlated. iov.cl + iov.v ~ c(0.1, 0.03, 0.2) | occ estimates the occasion covariance instead of erroring with “correlated inter-occasion random effects are not supported”.

  This adds foceiControl(iovMethod=), one of "auto" (default), "theta" or "omega":

  • "theta" is the long-standing expansion – one magnitude theta per occasion parameter, with unit-variance per-occasion etas fixed to it. That shape provably cannot carry a correlation between two occasion parameters, so it still refuses one.
  • "omega" fixes the magnitude theta at one and estimates the per-occasion eta blocks instead: occasion one is the block, and each later occasion repeats it. That is NONMEM’s $OMEGA BLOCK(n) SAME, and the correlation lives in the estimated block.
  • "auto" picks "omega" when the occasion block has any off-diagonal element and "theta" otherwise, per LEVEL of variability – a correlation on occ does not change how an unrelated diagonal occ2 is expanded.
  • "auto" only reaches for "omega" on an estimation method that honours the repeated block – the FOCEi family. saem, the variational (vae, fbvi, emvi), nonparametric (npag, …) and importance-sampling (imp, impmap, qrpem) methods estimate omega elsewhere and still refuse a correlated occasion block, since they would otherwise estimate each occasion independently and report only the first.

  The two expansions are the same statistical model, verified exactly: with the occasion variance at 1 (where the two parameterizations coincide) they agree on the objective to machine precision, and evaluated at matched random effects they agree to ~2e-8 at any variance. They are not interchangeable in practice, though – "theta" presents unit-scale etas to FOCEi’s inner optimizer, which converges the inner problem better when the occasion variance is far from 1 (on theo_sd with a variance of 0.1, "theta" reaches an inner optimum 0.059 lower and "omega" stops short of it). That is why "auto" keeps "theta" for a diagonal block and reaches for "omega" only when a correlation makes it necessary.

  iovXform parameterizes the "theta" magnitude and is inert under "omega", where the magnitude is fixed at one; asking for a non-"sd" value under "omega" now says so once.

  Analytic covariance falls back to finite differences for a repeated (SAME) omega block: its IOV special case reads the magnitude theta as the occasion standard deviation, which is 1 in this mode, and would have overwritten the estimated per-occasion variances rather than merely being conservative.

• foceiControl(innerOpt=) adds a trust-region Newton inner (per-subject eta) optimizer backed by the RcppTrust package, "trust". Unlike n1qn1, which gets an approximate Hessian only once as a warm-start seed, the trust-region step is supplied a fresh exact Gauss-Newton+Omega^-1 Hessian every iteration. Each eta is scaled by sqrt(diag(Omega)), with the trust-region radius derived from the eta confidence region (foceiControl(trustConf=), default 0.975); trustRinit/trustRmax override the derived radius directly. est="vae" does not use foceiControl()’s inner loop and is unaffected.

• foceiControl(innerOpt=) gains "auto", and it is now the default for every FOCEi-family method (focei/foce/foi/fo, and impmap’s MAP inner problem). It picks "n1qn1" for a generalized log-likelihood endpoint (dnorm(), ll(), dpois(), …) and "trust" for everything else, which is where each is faster: such an endpoint has no Gauss-Newton inner Hessian, so "trust" must finite-difference one (2*neta inner solves) at every trial point where "n1qn1" builds it once as a warm-start seed. On a 1-compartment oral model (120 subjects) as a dnorm() endpoint "trust" took 26.5s against n1qn1’s 16.3s, while on the same model with a normal endpoint it took 1.1s against 4.6s. "trust" and "n1qn1" remain selectable explicitly.

  This changes the exact numeric result of every FOCEi-family fit that does not pin innerOpt= explicitly (typically by a few objf units at most; see inst/benchmarks/results/ for a broad benchmark against n1qn1), though usually faster and comparably accurate. Pin foceiControl(innerOpt="n1qn1") to keep exact bit-for-bit reproducibility with prior releases.

• New nlm-family estimation method, est="trust" (trustControl()), a trust-region Newton optimizer for the population theta vector backed by the RcppTrust package – unrelated to foceiControl(innerOpt="trust") above, which optimizes per-subject eta instead. Unlike every other nlm-family method (nlm/nlminb/bobyqa/newuoa/uobyqa/n1qn1/ lbfgsb3c/optim), whose optimization loop lives in R and calls back into C++ once per iteration, trust‘s entire loop runs inside a single C++ call – RcppTrust needs no R API, so there is no per-iteration R round-trip. It optimizes in the same scaled-parameter space every nlm-family method (bobyqa included) already uses, and supplies a full gradient and a full Hessian every iteration. By default (trustControl(hessianMethod="fd")) this Hessian is a fresh finite-difference-of-the-gradient every outer iteration (there is no analytic outer-theta Hessian in this package, so this costs roughly ntheta extra full population-gradient solves per outer iteration – the price of true Newton-trust behavior). hessianMethod can instead build the Hessian as a quasi-Newton update from consecutive outer iterations’ gradients (already computed regardless of hessianMethod, so these add no extra evaluations): "bfgs" (damped BFGS, always positive definite), "sr1" (Symmetric Rank-1, not forced positive definite), or "bofill" (Bofill’s SR1/Powell-Symmetric-Broyden blend, the standard Berny/transition-state-search Hessian update) – see ?trustControl for full references. trust is unbounded, like n1qn1/nlm. trustControl(fterm=, mterm=), the Newton loop’s own function-value and predicted-decrease convergence tolerances, default to 10^(-sigdig-2) – two orders tighter than every other nlm-family tolerance here uses, matching foceiControl(trustFterm=, trustMterm=) (the analogous tolerance for the OTHER RcppTrust-backed solve in this package, the per-subject eta problem inside FOCEi) rather than the plain 10^(-sigdig) bobyqaControl() and foceiControl()’s own epsilon use.

• foceiControl(hessianMethod=) extends the same idea to FOCEi’s INNER (per-subject eta) problem: for a non-normal-endpoint model (any distribution other than norm), the per-subject inner Hessian has no Gaussian Gauss-Newton shortcut and falls back to a finite difference of the gradient every innerOpt="trust" Newton step (calcEtaHessian()). "fd" (default) keeps this original finite difference; "bfgs"/"sr1"/ "bofill" build the Hessian instead as a quasi-Newton update from consecutive Newton steps’ already-computed gradients (no extra evaluations) – the same three update formulas trustControl() above uses. Since this loop runs per subject, per Newton step, per outer iteration, avoiding a fresh finite difference at every one compounds into a much larger speedup than the outer-theta case: bfgs/sr1/bofill ran roughly 2.7-19x faster than "fd" on this package’s own small benchmark (a Poisson and a general ll() model, inst/benchmarks/benchmark-focei-hessian-method.R). "fd" stays the default: unlike the outer-theta case, this inner Hessian’s log-determinant is added directly into the reported objective (LikInner2()), and on a real one-compartment PK model fit as a general ll() endpoint every quasi-Newton option converged to the same wrong parameter estimate with a worse reported objective than "fd"’s correct answer – see the Bug fixes section below. Has no effect on normal-endpoint models. Only meaningful with innerOpt="trust": asking for "bfgs"/"sr1"/"bofill" under any other inner optimizer is now an error rather than a silent no-op, so a request that could not be honored is not lost. Since innerOpt="auto" (the default) sends exactly the non-normal-endpoint models this applies to to "n1qn1", using it means pinning innerOpt="trust" as well.

• impmapControl()/impControl() gain combSens (default TRUE): when est="impmap"/"imp"/"qrpem" has non-mu (structural or residual-error) thetas to estimate, combSens=TRUE carries their sensitivity columns on the INNER model itself instead of a second, dedicated model, and the E-step’s own per-sample inner solve now supplies the M-step’s Newton step directly (no second solve) whenever sir=FALSE (the default) – roughly halving the ODE solving the M-step’s theta gradient costs. Pass combSens=FALSE for the previous two-model behavior.

• A pure-linear matExp() model now solves natively through rxode2’s matrix-exponential driver (rxControl(method="indLin")) under SAEM instead of being flattened to an equivalent d/dt() ODE first. SAEM has no analytic-sensitivity consumer of the state derivatives, so native solving is all that changes; focei/nlm/nls are unaffected, and a model with an indLin() forcing term (e.g. Michaelis-Menten) still flattens (issue #859).

• A linCmt() parameter driven by both an eta and a time-varying covariate (for example cl <- tcl*(wt/70)^0.75*exp(eta.cl) with wt changing over a subject’s records) now gets an exact FOCEi-family eta gradient. The analytic linCmt() sensitivity reconstructs each row’s carried state as if the parameter had been constant over the subject, so a covariate that changes it between rows silently conflated the interval sensitivities (objective-function and converged-eta differences against the equivalent ODE model). The generated inner model now carries the exact sensitivity across rows through rxode2’s linCmtB() carry sentinels for every eligible (parameter, eta) pair; models without such a pair generate identical code. foceiControl(linCmtSensCarry=) opts out ("none"). Data with steady-state (ss > 0) or evid = 2 records fall back to the previous gradient with a note in $runInfo, and "linear" covariate interpolation on such a covariate is an error (a linCmt() model evaluates each interval at its row-end covariate value, so only a piecewise-constant interpolation is representable). The carry also covers the two ways an eta reaches the state through an event: a modeled f() whose d(ln F)/d(eta) depends on a covariate and a modeled alag() on a time-varying kernel each get a per-row jump contribution (#920‘s row-local terms are exact only while the parameters are constant), and every linCmt() parameterization (trans) is handled, with the observation-scaling term taken from rxode2’s own micro-constant translation. Data whose doses enter another compartment than the modified one, or an infusion with an alag() / covariate f() channel, fall back like ss records do. A generalized ll() endpoint (or any prediction that wraps the linCmt() value in a larger expression) is carried too: the concentration is read back once as rx_lcConc_, the carry supplies its eta sensitivity and symengine the outer chain rule, including any eta dependence the likelihood has with the concentration held fixed (#1004). The same carry also serves the population methods’ theta gradients (#1003): a theta on a covariate-driven linCmt() parameter gets the carried score in the nlm family and nls (nlmControl(linCmtSensCarry=) / nlsControl(linCmtSensCarry=)), with the concentration factored out of the wrapped log-likelihood; SAEM’s linearized FIM needs no change (it perturbs phi and re-solves values, which is exact under a time-varying covariate). Requires an rxode2 with the carry sentinels (the event channels need its which1 = -8 pin); older versions keep the previous behavior. The candidate detection itself is memoized by the focei model digest and persisted as a sidecar in rxode2’s cache directory (rxCreateCache()), so repeated fits – and, with a persistent cache, fresh sessions – skip the symbolic pass. The focei model-cache bundle itself now stores each generated model as its rxode2::rxNorm() text instead of a serialized model object (about 1 kb instead of 1 Mb; rebuilding from the text hits rxode2’s compiled-model cache); bundles written by an earlier version still load.

• focei/foce/agq/laplace/nlm now compute a matExp() model’s eta/ theta sensitivities natively via rxode2::rxSensMatExp(), instead of flattening the model to an equivalent d/dt() ODE first and differentiating that. A pure-linear matExp() model always takes this path; a model with an indLin() forcing term (e.g. Michaelis-Menten) takes it under focei and focep, and falls back to the ODE flatten (unchanged prior behavior) under foce (non-interaction), the mu-referenced/IRLS family (mfocei/ ifocei/mfoce/ifoce), and nlm – those combinations’ gradient/ covariance machinery is not yet compatible with the native forcing sensitivities and is tracked separately (issue #860; follow-up work in #861/#862). foceiControl(fast=TRUE) is automatically downgraded to fast=FALSE for a matExp() model taking the native path, since the analytic outer-gradient/covariance model (foceiCovAnalytic.R) is still ODE-flattened.

• A modeled alag() or f() on a linCmt() compartment now gets an exact FOCEi/FOCE eta gradient instead of a silently incomplete one. The structural linCmt() Jacobian only covers p1/v1/ka/…; the moving-boundary (dose-time) contribution of a modeled alag() is now added via rxode2’s linCmtB(which1=-3) (nlmixr2/rxode2#1235), and the bioavailability contribution via the exact d(pred)/dF = pred/F identity (issue #920).

  Both corrections require every dose reaching the linear system to share the same alag()/f() (rxode2/rxode2#1237); a model declaring more than one is left as before. This cannot be checked for a regimen that doses an unlagged/unscaled compartment alongside the lagged/scaled one (a common design for estimating f() from paired IV+oral data) – that combination returns a biased, not obviously wrong, gradient. An infused dose into the lagged/scaled compartment also cannot be checked, and returns NA (rxode2/rxode2#1236). foceiControl(eventSens="fd") opts out of both.

• A prior distribution given in the ini({}) block is no longer silently ignored. nlmixr2Est() now refuses any prior the estimation method cannot use before dispatching, so a model carrying one fails with an explanation naming the parameter and the est= instead of being fit to something other than what it says.

  A method declares what it supports with an attribute on itself, so it can opt in as it gains support:

  attr(nlmixr2Est.myMethod, "nlmixr2Priors") <- "general"

  The levels are "none" (the default when the attribute is absent), "theta" (population parameters only), "general" (everything the shared kernel supports, including a prior on an omega element under any convention), "nwpri" (NONMEM’s own $PRIOR NWPRI omega convention) and "tnpri" (Monolix’s/NONMEM’s own-estimation joint-normal convention, including a normal prior directly on an omega element) – see ?nlmixr2Est for what each accepts – and "all". Because the check happens in the generic, methods registered by other packages – babelmixr2’s nonmem, monolix, saemix and the rest – are covered without any change of their own.

• est="focei" and every method in its family (foce, focep, fo, foi, the mu-referenced mfoce*/IRLS ifoce* variants, laplace, agq and their quadrature/*f fast-path siblings) now honours a prior on a population parameter AND on an omega element – under any of the kernel’s three conventions ("general", NONMEM’s "nwpri", Monolix’s/NONMEM’s-own-estimation "tnpri"), auto-detected from what the model’s own ini({}) actually wrote – added to the objective as -2*log p(theta, omega) (nlmixr2/rxode2#1270, issue #929, issue #931) – declared nlmixr2Priors = "general".

  The convention is auto-detected by default (foceiControl(priorMethod= "auto")), but can be forced with foceiControl(priorMethod="general"/ "nwpri"/"tnpri"). Forcing a convention the model’s priors are not representable under (e.g. priorMethod="tnpri" on an invWishart() prior) errors before any estimation starts, naming the parameter and which method it needs instead.

  foceiControl(fast=TRUE)’s analytic outer gradient has a real d/dtheta log p(theta) term (a straight fold into the same natural-scale accumulator the outer FD substitution already uses) and a real d/d(chol(Omega^-1)) log p(omega) term (chain-ruled through the SAME estimation-scale derivative data – d.omegaInv/tr.28 from the model’s rxSymInvCholEnv handle – FOCEi’s own, non-prior omega gradient already relies on), so a prior no longer downgrades it: this package uses symbolic/analytic derivatives throughout, not finite differences, wherever one is available. (A prior referencing a theta a mu-referenced family profiles out of the outer problem entirely is the one case that first term cannot attribute; that specific fit declines to finite differences instead of silently under-counting, detected once at setup, not per evaluation.) covMethod="analytic" is still downgraded to a finite-difference covariance, since the analytic Hessian has no prior term yet – a finite difference of the (now prior-inclusive) objective picks the prior term up automatically.

  Fixed along the way: .nlmixr2FitUpdateParams() rebuilt a fit’s omega rows of iniDf from the raw Omega matrix whenever an mixed-effects model estimate was pinned back onto the model (piping, .setOfvFo()’s post-fit re-entry for setOfv()/addCwres()), and that rebuild had no way to carry the prior column lotri itself does not know about – a prior on an omega element silently vanished the first time a prior-carrying fit was re-entered, which is every fit’s own finalize step. rxUiPriors(fit$ui) now still reports it afterward.

• est="imp", est="impmap" and est="qrpem" now honour a prior on a population parameter AND on an omega element, declared nlmixr2Priors = "general" individually on each of the three (issue #932). Their shared M-step is an importance-sampling EM, not FOCEi’s outer optimizer, so the objective already picking up the prior (via #931’s plumbing) was not enough on its own – the estimates it reported would otherwise still be the maximum-likelihood ones. Each M-step update now folds in the prior’s own score/curvature before taking its step:

  • the non-mu structural/residual-error Newton step, the mu-referenced covariate regression (updateMuGroups()), and the plain mu-intercept mean-shift each fold in an FD-Hessian one-step Newton correction – exact for a Gaussian prior (a quadratic log-density has no Taylor truncation error), a reasonable one-step approximation otherwise (Cauchy, multiNormal()).
  • the Omega EM moment-average update gets the EXACT joint posterior mode for a conjugate invWishart() term (NONMEM’s own "nwpri" convention or the textbook "general" one), and a one-step Fisher-scoring (One-Step-Late) correction, reusing the same Abar construction FOCEi’s own omega-prior gradient already computes, for a normal prior directly on an omega element ("tnpri") or a multiNormal() block mixing omega with theta.

• A prior may now be placed directly on a single omega COVARIANCE (off-diagonal) element – prior(eta.cl, eta.v) ~ dnorm(0, 0.1) on a model with a correlated BSV block – a marginal, independent prior on that one cell, distinct from a whole-block invWishart()/ multiNormal() prior. NONMEM has no direct mechanism for this specific marginal form; it fills the same ergonomic gap this package’s om.<eta> shorthand already fills for one individual variance. This is entirely a new capability in the shared upstream kernel (rxode2/lotri, see their own NEWS) – FOCEi’s foceiPriorOmegaGradAdd() and imp/impmap/qrpem’s impPriorOmegaCorrect() already operated on the full gradOmega/Omega matrices generically, so both pick this up with no nlmixr2est source changes at all, confirmed by new tests in test-focei-prior.R and test-imp-prior.R that pass unmodified against the upgraded dependencies.

• est="saem" gained controls for the cost of the nonMuTheta="regress" refinement, which estimates population thetas that carry no random effect: saemControl(nonMuThetaOpt=, nonMuThetaSweeps=, nonMuThetaMaxEval=, nonMuThetaTol=, nonMuThetaEvery=). That refinement runs every iteration of the second half of the fit, and when the non-mu thetas are structural (they drive the ODE) each of its objective evaluations is a full re-solve of every subject and chain, so it can cost more than the rest of the algorithm put together – on the nimotuzumab target-mediated example and the mavoglurant PBPK example (both with five non-mu thetas) it is around 60% of the run time. nonMuThetaOpt="newuoa" (the new default, see Changed defaults) and nonMuThetaOpt="nelderMead" run one clamped multivariate optimization over all free phi0 coordinates under a fixed evaluation budget (nonMuThetaMaxEval); both see the coupling between those coordinates, which the previous coordinate descent (nonMuThetaOpt="optimize", still available) cannot. nonMuThetaEvery additionally runs the refinement only every k-th iteration.

• An estimated transform-both-sides lambda (boxCox()/yeoJohnson()) now carries a real theta-sensitivity column instead of a silent zero (#949). The conditional depends on lambda through both sides of the residual h(y; lambda) - h(f; lambda): the prediction side now comes from the sensitivity model (the direct partial is taken for residual-error thetas too, not hard-coded to zero, so rx_pred_’s rxTBS() is differentiated), and the DV side from a new d(lambda)/d(theta) output multiplied by the analytic d(h(y; lambda))/d(lambda) where the DV transform is applied. The censored (M2/M3/M4) score picks up the matching DV and LIMIT partials. The column agrees with central differences to ~1e-9 relative on Box-Cox and Yeo-Johnson fixtures where it was previously identically zero – the failure mode that made an estimated lambda an imp/advi M-step no-op and gave gradient-based callers a wrong direction.

• Dose-handling (alag()/f()/dur()/rate()) theta sensitivities are no longer silently zero (#946). The theta-sensitivity model is now compiled with rxode2’s analytic event (“jump”) sensitivities (following the control’s eventSens, the same source the inner model uses), and its solves run under its own event shape (OdeSwapEsBatch(odeSlotThetaSens) brackets the batch; the shape is a process global whose installer calls into R, so it cannot be swapped inside the parallel region). The inner batch in the conditional-likelihood C API likewise installs the inner model’s shape, so an eta entering dose handling gets its jump too. An estimated alag theta’s derivative now agrees with central differences (~1e-5 relative, ODE-tolerance-limited) where it previously came back identically zero – the failure mode that made an alag-estimating imp/advi theta update a no-op and broke gradient-based samplers.

• New FOCEi conditional-likelihood C API (#937): a plain-C, non-throwing, gradient-returning entry-point table over the foceiLikLoad()-ed problem, exposed to downstream packages through _nlmixr2est_foceiPtrs() / .nlmixr2estFoceiPtrs() with the caller side in inst/include/nlmixr2estFoceiPtr.h – the same external-pointer idiom as the likelihood-contribution registry. Seven entries: an ABI version, dimensions plus capability/hazard flags (so focep/fo/finite-difference etas/mixtures are refused at load rather than sampled wrongly), a return-code setTheta, the batched per-subject conditional log p(y_i | eta_i) with its d/d(eta) (subject-parallel, deterministic: the sticky solve-tolerance relaxation is reset every call so the value is a pure function of (theta, eta)), a global Omega^-1 conditioning knob, and the theta-sensitivity index/score pair for d/d(theta) of the conditional at fixed eta. The eta gradient is assembled inside nlmixr2est as Omega^-1 eta - fInd->lp – the same identity the ADVI outer gradient uses – so no caller ever reconstructs the sign convention. Built for the Stan linkage in nlmixr2/nlmixr2stan.

• foceiLikLoad() gains three arguments for external callers (#939): scale="natural" pins the parameter scaling to the identity so foceiLikRun()’s theta is directly comparable with ui$iniDf$est (no more re-implementing nlmUnscalePar() downstream); thetaSens=TRUE builds and wires the theta-sensitivity model (the d(f)/d(theta) / d(V)/d(theta) forward sensitivities the imp/advi engines use), reported on the handle as thetaSens/thetaSensIdx; and est= names the estimation method whose capability attributes the pre-process hooks consult, turning “the bounds survive preprocessing” from an accident of focei’s attributes into a guarantee the caller can request.

• nlmUnscalePar() is now exported (#940). External engines that drive the nlm-family objective (e.g. babelmixr2’s FME-based methods) previously had to reach it with get("nlmUnscalePar", envir = asNamespace("nlmixr2est")).

• Requires rxode2 (>= 5.1.7). The compatibility layer that also let this package build and run against 5.1.5 has been removed, so the event-sensitivity shape swap and the CMT re-basing of the shared solve pool always go through rxode2’s C API instead of writing its structures by field.

• est="npag" / est="npb" now support a hand-written general likelihood (ll()) properly. A model whose ll() is written as the exact normal log-density now agrees with the equivalent add() model to the known 0.5*log(2*pi) per observation, at every grid size. Requires rxode2 5.1.7 for the safeLog=2 log-domain mode.

• foceiControl(fast = TRUE) now keeps the analytic outer gradient when a subject’s augmented sensitivity solve fails, instead of sending the whole gradient evaluation to finite differences. Such a subject is finite-differenced on its own and folded into the otherwise-analytic sum, covering the omega directions as well as theta and sigma. The per-subject machinery was present but could never apply: a failed solve left the subject with no observation count, and the assembly declined the whole evaluation before reaching the substitution. fit$env$nOuterFdInd counts the substitutions applied, so “one subject was finite-differenced and the rest stayed analytic” is now distinguishable from “the gradient declined” – the two previously looked identical.

• foceiControl(fast = TRUE) now uses the analytic outer gradient for general-likelihood models with more than one endpoint, which previously fell back to finite differences. It was gated off as unverifiable, but what did not verify was the objective below rather than the gradient; against central differences of the corrected objective it agrees to 8e-3 relative.

• saemControl(revisitUninformativeEtas=) (default FALSE) re-runs the uninformative-eta test at the end of burn-in and replaces the verdict reached at the initial estimates. The test asks whether perturbing an eta moves that subject’s prediction, and is otherwise only run once, before the fit – so the initial estimates decide, for the whole fit, which etas saem may sample. The second test reuses the fit’s own model evaluation, so it adds a few solves at one iteration and leaves the random number stream alone: where it changes no verdict the fit is identical. It is off by default because the two verdicts only disagree when theta moved a long way during burn-in, which usually means it has not settled, and the second verdict can freeze an eta for the rest of the fit.

• est="saem" now recognizes t()/cauchy() residual-error endpoints as general-likelihood models, the same way a literal ll() endpoint already was, instead of erroring ("t isn't supported yet" / "Distribution not supported"). rxode2’s own FOCEi line generator already reduces these to the same shape (rx_pred_ an explicit log-density, rx_r_ ~ 0), so saem now dispatches through the same path.

• est="saem" now recognizes a general-likelihood model with any number of endpoints, including a genuine mix of norm and general-likelihood endpoints in the same fit (e.g. one add() condition alongside a t(), cauchy(), or literal ll() condition), instead of silently scoring the whole fit as if every endpoint were normally distributed. A norm condition mixed with a general-likelihood one is transparently expressed as the equivalent llikNorm() general likelihood (the same normal log-density, sharing the same variance-formula machinery), mirroring how FOCEi already handles this mix – no distribution()-family limitation remains.

Bug fixes

Estimation

• The standalone analytic-covariance entry point no longer installs a covariance that is not positive definite. An outer optimizer that stops short of a local minimum leaves an observed information with a negative eigenvalue, which inverts to negative variances and NaN standard errors – and once installed as the fit’s $cov that is what every later getVarCov() returns. The fit’s existing covariance is now kept and a warning says why. The live covMethod="analytic" seam, setCov() and the saem installer already guarded this; the standalone entry was the one that did not.

• Corrected objective scaling of the fast outer gradient.

• Corrected FOCE curvature’s row stride when solves included bookkeeping rows.

• Aligned FOCE+ objectives and derivatives with the live-variance ETA score root.

• Included M2/M3/M4 censoring in the analytical AGQ gradient and the FOCE/FOCE+/FOCEI/AGQ outer Hessian with censOption="gauss".

• A focei-family fit now reports whether its inner solves actually converged. fit$env$nTrustInner breaks the innerOpt="trust" per-subject Newton solves down by outcome (calls, error, notConverged, solverFail, newtonGate, warmRetry, radiusRetry, nudge, failed) and fit$env$nInnerRerank gains noGood and dropped. Only the call count existed before, and only through an internal accessor, so a fit whose inner solves were all failing was indistinguishable from one where they all converged.

• The inner restart candidates are no longer chosen from without regard to whether the attempt that produced them succeeded. A failed attempt’s eta could win the marginal re-rank over a converged one – its Laplace log|H| term is measured at a point the inner objective never descended to – which mceta >= 1 makes more likely, since the extra starting points are what produce a mixed candidate set. A failed attempt is now used only when no succeeded one is available.

• mceta >= 1 no longer fails a subject whose eta=0 floor pass produced nothing. The “did this pass find anything” check is per-pass, but the trust arm and the n1qn1 restart cascade both treated it as “did this subject find anything” and returned a failed inner solve, discarding a converged candidate an earlier starting point had already produced.

• innerOpt="trust" re-solves in place before falling back to its eta-nudge restarts. When the inner Newton step still fits inside the current trust radius, the solve stopped on its own step-size criterion rather than on the model’s remaining decrease, and the cascade threw that point away to restart from a nudge fill. On a 300-subject fit swept away from the true parameters, subjects whose whole cascade is exhausted drop by up to half and the total number of inner solves drops with them; a fit where nothing trips the check is unchanged.

• A fit no longer reports the PREVIOUS fit’s censoring. $censInformation is built from a process-global flag recording which censoring methods (M2/M3/M4) a fit used, and that flag was cleared only after being read at the end of a fit – so a fit that set it without reaching there left it set, and the next fit in the session reported the leftover. An uncensored model fit right after a censored one said “M2, M3 and M4 censoring”. It is now also cleared when a fit starts.

• est="agq"/"laplace"/"focei" with the default innerOpt="trust" and a finite-difference outer gradient no longer stop early at a much worse objective. trustFterm/trustMterm, the inner per-subject Newton solve’s convergence tolerances, defaulted to the plain 10^(-sigdig) used by the other tolerances here, but the inner solve is the function the outer problem differentiates: its stopping tolerance is the objective’s noise floor, and a finite-difference outer gradient cannot resolve a step below it. On an nAGQ=2 theo_sd fit the outer search was misled to an objective of 134.46 against 118.52, taking 6.8s against 1.5s – a noisy gradient both lengthens the search and ends it worse. They now default to 10^(-sigdig-2), two orders tighter, the same relationship lbfgsFactr already uses. Every fit with innerOpt="trust" and no explicit trustFterm/trustMterm changes numerically as a result; generalized-likelihood models are unaffected, since innerOpt="auto" sends those to "n1qn1". Isolated on this package’s own 100+ model benchmark corpus (inst/benchmarks/results/trust-inner-benchmark.md), the tighter tolerance is a real median accuracy improvement (|objf diff| 0.0115 -> 0.004) at a small median speed cost (1.10x -> 1.07x faster than "n1qn1") – not a free win on every model.

• A generalized log-likelihood model (dnorm(), ll(), dpois(), …) fit with foceiControl(innerOpt="trust") no longer converges to badly biased population parameters on oral or two-compartment models. Such an endpoint has no Gauss-Newton inner Hessian, so the per-subject eta Hessian is a finite difference of the analytic eta gradient; the Shi (2021) step search is told that gradient’s noise floor is rxControl(atolSens=), but the gradient comes out of a sensitivity solve rtolSens governs as well. Understating the noise inflates the search’s ratio test, which then shrinks the step until the difference is taken inside the noise, with only the absolute shi21hMin (1e-4) to stop it. n1qn1 absorbed that – it uses this Hessian only as a warm-start seed and then corrects it by its own quasi-Newton updates as it iterates, so a noisy seed is transient – but trust re-derives it as its trust-region model Hessian at every trial point, with nothing to correct it, and adds its log-determinant to the reported objective, so the noise steered both the step and the number being minimized: on a 1-compartment oral model (120 subjects) fit as a dnorm() endpoint, Vc came out 90.6 against a plain focei 66.4, with eta variance shrinkage 56/52/34% against 8/10/13% and the omegas left at their starting values. The step is now floored at a fraction of each eta’s own SD (foceiControl(hessEtaStepMin=), default 0.05), recovering Vc 65.5 at unchanged solve tolerances; tightening atolSens or rtolSens instead also fixed it, at 2-4x the runtime. The floor applies to every inner optimizer – the step a finite difference needs is a property of the problem, not of who consumes the Hessian – though only trust was visibly broken without it. Normal endpoints, which use the Gauss-Newton inner Hessian, never reach this code.

• foceiControl(mceta = ) is no longer discarded for a model whose etas are all mu-referenced. Any non-default setting was reset to -2 on the grounds that “the initial etas are all exactly zero, so the search has nothing to explore” – true only of the first inner solve, since every later one starts from the previous iteration’s mode, and the mceta > 0 candidates are draws from omega, which are not zero. This is what made mceta = 10 return an objective bit-identical to mceta = -2 on such a model. On theo_sd’s one-compartment model (every eta mu-referenced), mceta = 5 now starts 420 of its 1008 inner solves from a draw and reaches a different objective than mceta = -2.

• foceiControl(mceta = n) (n > 0) now actually uses the extra starting etas, and at fixed parameters can no longer give a worse objective than mceta = 0. The candidate set included the carried “last eta” – the previous outer iteration’s converged conditional mode – whose inner objective is essentially always the lowest, so it won for every subject and mceta = n returned an objective bit-identical to the keep-last behavior of mceta = -1/-2. The candidates are now eta=0 plus the n-1 draws from omega.

  Because a candidate is ranked by the objective at its starting point, which does not order the points the inner optimization converges to, a subject that starts from a draw now also solves from eta=0 and keeps whichever converged lower – without that, a draw that merely looked better could converge worse (measured on a fixed-omega inverse-CDF model at fixed parameters: mceta = 2 gave -2251.0 against mceta = 0’s -2302.5, and now gives -2338.8). That floor solve wraps the whole inner-optimizer dispatch rather than living inside one arm of it, so it holds for whichever inner optimizer is configured.

  The comparison is made on the objective the fit REPORTS – the marginal one, which carries the Laplace log|H| term – and not on the inner joint density the optimizer minimizes. The two order candidates differently often enough (on the model above, 23 of the 44 subjects that had a choice) that ranking on the inner objective alone handed the fit the worse candidate. Restarts are ranked only when there is more than one candidate, so an inner solve that never restarts pays nothing for it.

  Finally, mceta = 1 means eta=0 rather than being a silent no-op, a non-finite eta=0 no longer pins the search, and the draws are made once per fit instead of at every objective evaluation, so the objective is the same function at every evaluation. The fit records which candidate each inner solve started from in $env$nMcetaStart, and how often the two orderings disagreed in $env$nInnerRerank.

• focei now checks rxode2’s per-subject event counts before it sizes the per-subject blocks it strides with them (gVid, ga/gc, gB, gcH*, llikObsFull). When rxode2 and nlmixr2est are built against different solve layouts those counts are read from the wrong bytes, and the setup sized megabytes of storage from garbage – reported as “dataset too large”, as an R_Calloc failure, or as a segfault, depending on what the mis-read bytes held. A negative count, or a dose or evid=2 count larger than the subject’s own record count, now stops the fit with that as the reason instead (#1039).

• Re-fitting a completed importance-sampling fit from the fit object – nlmixr2(fit, est = "imp" | "impmap" | "qrpem") – failed outright with unused argument: 'impMuThetaIdx', .... The control is re-validated by do.call(impmapControl, ctl), which forwards anything it does not recognise to foceiControl(), and the four per-model M-step index maps stamped on a fit’s runtime control are arguments of neither. They are now carried through the round-trip.

  With that fixed, est also now wins over any field another method’s est stamped on its control. est="imp" stamps mapIter = 0 (never re-center) and est="qrpem" stamps qr = TRUE, sir = TRUE, so a re-fit would otherwise have run a different algorithm than the one asked for – re-fitting an imp fit as "qrpem" drew plain Monte-Carlo samples and still reported QRPEM, and a qrpem fit re-fit as "imp" kept quasi-random sampling on. Values the user wrote themselves are untouched: the rule is keyed on the est field a completed fit carries, not on the value.

• impmapControl(mapIter=) was accepted and then ignored: the kernel re-centered the importance-sampling proposal at each subject’s MAP mode on every EM iteration regardless of the value. It now sets the MAP-assist period – 1 (the default, and the previous behavior) re-centers every iteration, k > 1 every kth, and 0 not at all after the startup MAP pass. Affects est="impmap" and est="qrpem"; est="imp" never re-centers and is unchanged.

• With two or more occasion parameters on one level, fit$iov$<level> had NA for every occasion (and the fit warned “NAs introduced by coercion”). The occasion number was parsed out of the eta names by stripping the last occasion parameter’s rx.<name>. prefix from a column built out of the first one’s names, so nothing matched. This affected the existing IOV expansion too, not just the new one.

• Two occasion parameters whose names share a prefix (iov.v and iov.v2) had their per-occasion columns mixed together in fit$iov$<level> and fit$shrink, again leaving every occasion NA. The columns are named rx.<param>.<occ> and were selected by substring, so iov.v also matched rx.iov.v2.1. Also present on the existing IOV expansion.

• saem’s Gaussian-quadrature objective no longer silently attempts a grid it cannot finish. The grid is nnodesGq^nphi1 whole-population solves and nphi1 grows by one per occasion level for every IOV parameter, so three occasion parameters over two occasions already asks for 3^9 = 19683 solves and four over three asks for 3^16, which is 43 million. The node count is now stepped down until the grid fits a budget (getOption("nlmixr2.saemGqMaxNodes", 50000)), and the fit says which count it used in $runInfo; one node is the Laplace approximation, which the progress message already names.

• saem now gives the same answer every time for a model whose residual error needs the internal optimizer – anything richer than a pure add() or pure prop() endpoint, so combined1()/combined2(), pow(), and the transform-both-sides variants. Repeating one fit in a session, same data and same seed, could return a different residual estimate each time (measured: add.sd 0.0859, 0.0877, then 0.1807 for three identical calls). The residual step caches the transformed predictions and observations it scores, and the cache was keyed on the ADDRESS of the buffers holding them. Those buffers are rebuilt every M-step, so the allocator routinely handed back the same address with different contents and the cache reported itself still valid – the optimizer then scored later iterations against the FIRST iteration’s predictions, which still carry all of the between-subject variability as error. The residual therefore stayed near its starting value and the variance components shrank to compensate. The cache is now invalidated whenever the data behind it is rewritten.

• saem with IOV and a general log-likelihood (ll()) endpoint no longer returns an astronomically large objective function (#1000). Past half the iterations, saem refines its non-mu-referenced (phi0) parameters with a bounded optimizer while holding the ODE states fixed, on the grounds that a general-likelihood phi0 parameter is a likelihood standard deviation the solve never sees. The IOV magnitude is a phi0 parameter that does drive the structural model, so with the states held fixed the objective was exactly constant in it and the optimizer ran to its upper bound: the IOV magnitude grew geometrically (past 1e17), and the reported objective followed. Whether the states can be held fixed is now measured rather than assumed, and a phi0 parameter that drives the solve is refined inside a local trust region. A failed solve reaching the Gaussian-quadrature objective is also scored as a bad solve rather than as an extremely good log-density.

• A second focei fit of a model whose dosing depends on an eta (f(depot) <- exp(eta.f), alag(), dur(), rate()) no longer silently returns the wrong answer. The compiled model bundle is cached as model TEXT and rehydrated with rxode2(), but the event-sensitivity mode was not stored with it – so the FIRST fit of a model in a session built its sensitivity models with eventSens = "jump" (rxode2’s analytic dosing-parameter sensitivities) while every LATER fit rebuilt them without it. Nothing errored: the dose-parameter sensitivity was simply zero, so those etas never left their initial values and the objective, the etas and the analytic gradient were all wrong (etaF pinned at 0; the analytic outer gradient exactly 0 in the f/alag directions). The mode is now recorded when the model is built and replayed when the bundle is rehydrated, and the cache key carries a format marker so entries written before it are not reused.

• Fitting the same data twice in one session gives the same answer again (#1020). On a machine with at least twice as many threads as the problem has subjects, the objective function and the estimates differed between two identical nlmixr2(..., est = "focei") calls, and the optimizer frequently stopped at (or beside) the initial estimates. The cause is in rxode2, where two bugs met on the path FOCEi takes: sortIds() ordered subjects by their measured solve time, so the order was a function of wall-clock timing, and several per-individual drivers read ind_solve()’s subject id as a position in the reordered rx->ordId – so once the order stopped being the identity the wrong individual was integrated. The fixes are in rxode2; this package gains the regression test.

• IOV (iov.x ~ v | OCC) no longer copies an unrelated parameter’s prior onto the parameters the expansion creates, and now carries the prior the user declared. .uiApplyIov() builds the IOV magnitude theta and the per-occasion etas by copying an existing iniDf row as a template, and did not clear the template’s prior. So the magnitude theta silently inherited the FIRST theta’s prior (an estimation method with prior support sampled it against a distribution belonging to another parameter), a prior(iov.x) written on the occasion eta was dropped with the row the rewrite deletes, a fix()ed IOV parameter was refused outright (“a prior given for fixed parameter(s)”), and so was any model whose first eta carried a prior – every per-occasion rx.<iov>.<occ> eta inherited it. The magnitude theta now carries prior(iov.x) (on the iovXform scale, "sd" by default); the copied rows carry no prior otherwise.

• Several IOV parameters on ONE occasion variable (iov.cl ~ 0.1 | occ; iov.v ~ 0.04 | occ) work. The occasion variable was visited once per parameter riding it, duplicating every magnitude theta, and fix was read from a vector over the whole occasion rather than from each parameter’s own row, so the model errored with “replacement has 2 rows, data has 1”.

• The IOV parameter restored onto a finished fit keeps its own prior. .uiFinalizeIov() rebuilds the user’s iov.x ~ v | occ row from a template copied from the first remaining eta and restored eight fields from the original but not prior, so fit$ui$iniDf reported the FIRST eta’s prior on every IOV parameter – the same template-copy mistake as above, on the way back out.

• An occasion parameter with two variance declarations (iov.cl ~ 0.1 | occ; iov.cl ~ 0.15 | occ) is named in the error. rxode2 does build that ui, so every per-parameter field the rewrite read was a vector and it died on “replacement has 2 rows, data has 1” without saying which parameter was at fault.

• Correlated inter-occasion random effects (iov.cl + iov.v ~ c(...) | occ) are refused with an explanatory error. The expansion gives each occasion parameter its own magnitude theta and unit-variance etas, which cannot represent a correlation between two of them; the off-diagonal row was treated as one more occasion parameter named (iov.cl,iov.v), and the model died in rxRename() with unexpected '='.

• Fixed a heap overflow in the FOCEi theta-reset path. The buffers it saves and copies back on a restart each had their length re-derived from a second copy of the allocation’s formula, and every copy had fallen behind the layout it described – most damagingly the eta block, which claimed nall^2 where sum(nobs_i^2) had been allocated, so a reset read past the end of the block and wrote the overshoot back. Depending on what followed it in the heap, a fit that reset its thetas could return truncated state, corrupt an unrelated allocation, or abort the R process outright (test-matexp.R did the last of these). Each length now comes from the allocation itself, and a restore whose saved length does not match the current one errors instead of copying.

• est="saem" no longer prints solve(): system is singular; attempting approx solution on every iteration of the second half of the fit when the model has more than one population theta without a random effect. The nonMuTheta="regress" refinement (the default) wrote its result back with one least-squares against the whole phi0 design, but that design is block structured – each coefficient belongs to exactly one phi0 theta – and with no phi0 covariate every one of its columns is the same intercept column, so the normal equations are exactly rank deficient. The back-solve is now done per phi0 theta against its own columns. The rank-deficient solution also filled the coefficient matrix off-structure, which is why a fix()ed non-mu theta of a general log-likelihood model did not reproduce its fixed value; that is fixed with it. Estimates are otherwise unchanged.

• Fixed the analytic-gradient methods (est="foceif" and the rest of the *f family, or foceiControl(fast=TRUE)) stopping short of the optimum, after roughly as many outer evaluations as the model has parameters regardless of the data. Their outer optimizer is lbfgsb3c, whose lbfgsFactr rule tests the objective reduction of a single step rather than stationarity, so at the previous default of 10^-sigdig the fit stopped as soon as one step was small. foceiControl(lbfgsFactr=) now defaults two orders tighter, 10^(-sigdig-2)/.Machine$double.eps. On a 2-compartment oral fit this moves the objective from 8.63 above the derivative-free reference to 0.10 below it, at 13 outer evaluations instead of 6 (the reference needs 96); a 1-compartment fit and a sparse pheno_sd fit improve similarly. Fits that were stopping early now run further and report a lower objective, which can change model ranking. The other sigdig-derived optimizer tolerances are unchanged – tightening nlminb’s rel.tol/x.tol measured no benefit.

  Note that this does not make the analytic-gradient methods interchangeable with the derivative-free default on every surface. On models where the two optimizer classes converge to different points – measured on a correlated-eta (block omega) model and a Michaelis-Menten model – the gap is unaffected by any tolerance, so compare candidate models within one estimation method rather than across focei and foceif.

• Fixed est="saem" scoring a general log-likelihood endpoint (ll()) as a Gaussian observation in both its objective function and its standard errors. Such an endpoint estimates no residual error, so the residual step never runs and the placeholder values it starts from survive: every log-density was scored as a normal mean with standard deviation 10 + |ll|. The reported objf/logLik/AIC/BIC and every standard error were meaningless, on the default path – covMethod="sa" cannot be computed for these models and already fell back to the linearized Fisher information, which is where the defect lives. An ll() row now contributes its own log-density to the objective, with the log(2*pi) normalizer applied only to normally-distributed rows, and contributes the observed information of that log-density to the covariance. On an exponential time-to-event model with a closed-form marginal likelihood the reported -2LL goes from 1124 to within 1e-3 of the exact 1392; against the Gaussian twin of a one-compartment model (an ll() written as the exact normal log-density versus the equivalent add() model) the standard error ratios go from 4.0-80.5 to 0.99-1.02, and the two objective function values now agree outright rather than up to a constant.

• est="saem"’s mu-referenced (population, eta-carrying) theta update for a general-likelihood (ll()) endpoint is no longer a plain stochastic-approximation recursion over the MCMC-sampled phi, which carried no curvature information and could converge to the wrong basin. A direct bobyqa optimization of the exact joint log-likelihood (mirroring the existing nonMuTheta="regress" mechanism) now drives these thetas, sharing FOCEi’s own inner (eta-sensitivity) model through the shared ODE solve pool. saemControl(phi1Hessian=) (default FALSE) optionally adds a Laplace log|H| correction on top; measured to not be what fixes convergence for a Gaussian or exponential-TTE ll() model, so it stays opt-in. A near-Gaussian t()/cauchy() endpoint can still diverge under this step (nlmixr2/nlmixr2est#999); tracked separately.

• Fixed est="saem" applying the transform-both-sides log-Jacobian with the wrong sign in its Gaussian-quadrature likelihood, so the reported objf/logLik/AIC/BIC for an lnorm(), boxCox(), yeoJohnson(), logitNorm() or probitNorm() endpoint was off by 4*sum(log|dt/dy|). The quadrature builds the likelihood of the transformed observations, so the Jacobian has to be added to reach the likelihood of the original data – the convention FOCEi already uses. On a one-compartment lnorm() fit whose exact marginal likelihood is a one-dimensional integral, the reported -2LL goes from -214.37 to 113.86 against an exact 113.86. add(), prop() and ll() endpoints have a zero Jacobian and are unaffected.

• The saem quadrature likelihood now accumulates in the log domain. A transformed endpoint with many observations per subject makes the per-subject log-density large enough that the old exp() accumulation overflowed, and the whole fit reported an infinite objective function.

• saemControl(nsdGq=) is now honored when the likelihood is calculated with the fit; it was read under a name the control never stored, so any value other than the default was silently ignored.

• saemControl(covMethod=) "sa" and "fim" now say plainly that they do not apply to a general log-likelihood endpoint and use the linearized Fisher information, instead of reporting that the covariance “could not be computed”. The stochastic-approximation covariance phase is also skipped for such a model rather than run and discarded.

• saem’s uninformative-eta detection (saemControl(handleUninformativeEtas=TRUE), the default) could freeze an eta that the data does inform. The test asks whether perturbing an eta moves the prediction, and it is run once, at the initial estimates; when those are poor enough that the prediction underflows at the observed times, nothing moves and the eta is frozen at its mu for the whole fit. The verdict is now only taken when the subject’s largest prediction is finite and itself above the tolerance. On a warfarin fit started from k=1/h (true value near 0.02/h) this froze the volume eta for 19 of 32 subjects, biasing the population estimates and shrinking that eta’s variance about fourfold.

• Fixed est="saem" erroring with missing value where TRUE/FALSE needed while finalizing a fit with an unfixed logitNorm(), probitNorm(), propF(), powF() or powF2() residual sd (#915). .getSaemTheta() copied the saem-estimated residual sd back into ui$iniDf$est for the add()/lnorm()/prop()/pow()/pow2()-family endpoints but not for these, so the value stayed NA and the FOCEi scaleC setup that reenters to build the fit table hit it during an if() test.

  A second, independent gap in the same error types was found alongside it: rxUiGet.saemAres/saemBres/saemCres also never matched logitNorm()/probitNorm()/propF()/powF()/powF2(), so the SAEM kernel silently started every such fit from the hardcoded fallback (10/1) instead of the residual sd given in ini().

• Fixed foceiControl(fast=TRUE) FOCE fits (est="foce", "mfoce", "ifoce") discarding a usable analytic outer gradient and paying for a full finite-difference gradient instead. FOCE freezes the residual variance, so its mode is not the inner problem’s and an inner Newton has to find it; that Newton demanded a score below foceiControl(foceEbeTol=) (1e-9) even though the score is computed from the ODE solve and cannot be driven below the solve’s own noise. On the reported model it reached |S| = 1.5e-9 and then threw the whole gradient away over a Newton decrement – the objective the point still had left to give – of 2e-15. A stalled subject is now accepted at its best iterate when that decrement is negligible, and the iterate itself is kept rather than wherever the exhausted line search stopped. mfoce reaches a pure analytic gradient on the reported model, and the 3-ETA theophylline fit’s finite-difference fallbacks drop from 16 to 1 while mfoce and ifoce – which solve the same problem – now agree with each other instead of landing 6.7 objective-function units apart. A genuinely unconverged mode still declines.

• Fixed foceiControl(gradTrim=) lower gradient clamp testing g < gradTrim instead of g < -gradTrim. Since the branch above it had already caught everything over +gradTrim, every remaining component was replaced by -gradTrim, so a small positive gradient could reach the optimizer as a large negative one. Only reachable with a finite gradTrim; the default Inf skips these branches.

• Fixed the outer finite-difference gradient corrupting any component whose forward difference falls below foceiControl(gradCalcCentralSmall=). The confirming central difference overwrote the objective at theta+delta with the forward gradient before using it, so it returned roughly -objective/(2*h) rather than a derivative, and that value was left unclamped by gradTrim. The confirmation now also keeps the gradient it started from when its own solve fails, rather than replacing it with a non-finite value that resets the fit – the same rescue the two gradTrim recomputations were missing.

• Fixed the outer finite-difference gradient returning a sign-reversed or stale derivative when a central-difference term came back non-finite. The one-sided rescue used an objective that is never filled in on the central path, and on the path switched to central by foceiControl(gradCalcCentralLarge=) it read the previous parameter’s perturbed objective.

• Fixed covMethod="analytic" ignoring its own solve tolerance whenever the shared ODE solve pool was available, solving at the fit’s much looser tolerance instead of foceiControl(covSolveTol=) (or, unset, a value tightened from sigdig). The augmented solves are differenced twice to recover a 3rd-order tensor, so their error is the standard errors’ error: they carried the fit’s instead, and because the pool needs fast = TRUE, foceiControl(fast=) changed the standard errors (1.7e-2 relative on a 5-ETA 2-compartment model at the default sigdig). Both routes now agree exactly. Set covSolveTol to trade accuracy back for the slightly larger covariance step.

• Fixed foceiControl(covSolveTol=) being dropped part-way through the covariance step. Once the analytic route had restored the fit’s ODE solve – which it does whether it succeeded or declined – the finite-difference covariance work after that point ran at the fit’s tolerance again, because rebuilding the solve resets the tolerances along with it.

• Fixed a subject whose pooled augmented solve failed being scored into covMethod="analytic" as zeros – no prediction and no sensitivity – instead of sending the covariance to its fallback. The zero fill was written for the R outer gradient, which replaced such a subject’s column by a finite difference; that gradient is gone, and zeros are finite, so nothing downstream noticed. Such a population now falls back to the unpooled solve.

• Fixed the pooled covMethod="analytic" solve running single-threaded. It coerced rxControl(cores = 0) – the default, meaning “use rxode2’s thread setting” – to a literal 1, so its loop over subjects never went parallel, while the rxSolve route it replaced passed the 0 through and did. A 5-ETA 2-compartment covariance goes from 1.04s to 0.71s.

• Fixed est="npag" / est="npb" reporting a log-likelihood above its analytic maximum for a model with a hand-written general likelihood, with the residual parameters driven out of domain – including to a negative standard deviation. On a 2-endpoint PK/PD fit the log-likelihood read +2364 to +2831 where the data bounds it at -155. Three causes: the residual step scored every row with a Gaussian extended-least-squares form even where rx_pred_ is a log-density (the same defect class as #838, in a function that fix did not touch); the moment warm start took a moment of a log-density; and rxode2’s safeLog turned log(negative) into a large finite value, so an invalid negative SD was rewarded by about +36 per observation rather than rejected.

• Fixed npag’s reported objective being inflated whenever a residual variance collapsed, general likelihood or not. likInner0 floored the variance r to 1 for the err^2/r term but took log() of the unfloored value, so the two terms disagreed – and the disagreement paid +18.02 per affected observation.

• The nonparametric engines no longer accept an evaluation the inner problem refused. npEvalCondLik discarded likInner0’s NA return, and because the per-observation likelihoods are initialized only once, a rejected evaluation summed a finite blend of two different parameter vectors.

• The per-subject finite differences are taken at the fit’s ODE tolerance. A subject only reaches this path after its tolerance was loosened to try to rescue the failed solve, so it would otherwise be differenced on a looser function than the analytic terms it is summed with – and a retry part-way through could leave the two legs of one central difference at different tolerances.

• Fixed the outer finite-difference step store being allocated per optimizer parameter but indexed by full-theta position, so a fit with any fixed parameter wrote past the end of it.

• Fixed the objective function for a model that has a general-likelihood endpoint (ll(), pois(), binom(), …) alongside any other endpoint. Each observation’s distribution was read one row before the model had been evaluated for that row, so a subject’s FIRST observation was scored as normal: its log-density was treated as a prediction of DV against a variance forced to 1. On a two-endpoint warfarin model the objective read 11,463,666 where the correct value is 53,697, and the conditional estimates were shifted with it. This affected such fits at any foceiControl(fast=) setting. Models with a single endpoint, and models whose endpoints are all Gaussian, are unchanged.

• Fixed the est="imp"/"impmap"/"qrpem" theta score for endpoints with different DV transforms, e.g. an lnorm() PK endpoint alongside an add() PD one. The M-step read each observation’s transform and distribution without evaluating the model for that row, and the theta-sensitivity model did not emit rx_yj_/rx_lambda_ at all, so every observation was scored with one arbitrary endpoint’s transform – on a 2-endpoint PK/PD fit that put tka at -47.6 and the residual sigma at 2.8e4 where FOCEI gives 0.53 and 0.11. The two now agree to 1e-3. Models with a single endpoint, or whose endpoints share a transform, are unchanged.

• foceiControl(fo=TRUE) now rejects a general-likelihood endpoint or a censored observation wherever it appears in a subject, not only on that subject’s last observation. Both guards tested the last row’s value, so a subject whose final observation was Gaussian and uncensored slipped past them and was fit with an objective FO does not support.

• optExpression=FALSE (and sumProd=TRUE) are now honored for the last model each fit builds: the EBE / Llik EBE models of the nlm-family estimators (bobyqa, newuoa, n1qn1, nlm, nlminb, optim, uobyqa, lbfgsb3c, nls) and nlme, plus the saem predOnly model. The control was removed from the model while finalizing the fit, before those models were built, so the build read the defaults instead (issue #864). This matters because optExpression=FALSE is the workaround for a delay-differential model whose past() duration is an expression.

• Every nlm-family method’s default scaleType="nlmixr2" scale constant no longer explodes to its scaleCmax ceiling for a parameter whose starting gradient happens to be genuinely near zero – for example a bounded/ upper_exp-transformed theta whose sensitivity is tiny at the model’s default initial estimate (issue #994). The derivative-based formula scaleC[i] = |gradTo/gradient_i| had no guard analogous to FOCEi’s own scaleCtheta/.guardScaleC() safeguard, so a near-zero denominator was clamped to a scale constant up to 100000x too large, permanently distorting every later scaled gradient/Hessian entry for that dimension and derailing est="trust"’s Newton-based step decisions in particular. Each element is now guarded with the same .guardScaleC() band FOCEi already uses, falling back to the transform-aware ui$scaleCtheta value when out of band.

• est="trust" on a model mixing linCmt() with an ODE (e.g. an effect compartment) no longer converges to a badly wrong objective (issue #996). trust was left out of .linCmtOdeEstFamily (R/preProcessLinCmtOde.R) when it was added, so the linCmt()-to-ODE translation this list exists to trigger – needed because the extra theta-sensitivity states this family adds push linCmt()’s compartments past the numbers the data was translated against, the same problem #286 fixed for FOCEi – never ran for it. 5 of 7 population thetas had an exactly-zero starting gradient as a result. "trust" is now included in that list.

• nlmControl(normType=)/trustControl(normType=)/etc.’s "mean", "std", and "len" normalizations (with the default scaleType="nlmixr2") now compute their mean/standard-deviation/length constants from every estimated parameter instead of silently dropping the last one (issue #995). scaleSetup()’s (src/scale.h) per-normType setup loop used for (unsigned int k = scale->npars-1; k--;), which tests the pre-decrement k for truthiness before the body runs, so the body itself never executed with k=npars-1 – the top-indexed parameter’s value was excluded from the running mean/variance/sum-of-squares, and its scale constant was never reset for recomputation. "rescale"/"rescale2" (the default) were unaffected (their min/max accumulator is separately seeded with the top parameter before the loop runs).

• foceiControl(hessianMethod=)’s default is "fd" again (reverting a mid-development flip to "sr1" made for consistency with trustControl()’s own outer-theta default). A quasi-Newton inner Hessian feeds its log-determinant directly into the reported per-subject Laplace objective (LikInner2()), unlike the outer-theta case where the Hessian only affects the step; on a real one-compartment IV bolus PK model fit as a general ll()/dnorm() endpoint, "bfgs"/"sr1"/"bofill" all converged to the same wrong Vc (about 90 against a simulated 70 and a plain focei fit’s ~67-68) with a worse reported objective than "fd"’s correct answer – the inaccurate Hessian misled the outer search into a worse point it reported as better.

  trustControl(hessianMethod=)‘s own outer-theta default was briefly reverted to "fd" too pending confirmation, since that option’s earlier benchmark (showing "sr1" about as accurate and faster) predates the fixes for issues #994 and #996 above, both of which distorted several of its benchmark models’ results identically regardless of hessianMethod (a wrong raw gradient upstream of Hessian construction, in both cases). Re-run after both fixes, "sr1"/"bofill" track "fd" closely (median |objective diff| vs bobyqa, across the same 23-model corpus: 1.53/1.55 vs "fd"’s 1.55) and "bfgs" if anything tracks it slightly better (0.43) – confirming the earlier small accuracy gap was at least partly noise from those bugs, not a genuine Hessian-construction difference for this outer problem, which does not have the inner problem’s specific failure mode (nlmTrustObjfun()’s reported value is the plain log-likelihood, set before the Hessian is touched, so a less-accurate hessianMethod here can only cost step quality, not bias the reported number). trustControl()’s default is therefore "sr1" again – faster, with no demonstrated accuracy cost for this problem.

• foceiControl() gains trustFterm/trustMterm: innerOpt="trust"’s own function-value and predicted-decrease convergence tolerances, independent of epsilon (which is also "n1qn1"’s unrelated “precision of estimate” tolerance – tightening epsilon to tune one inner optimizer used to tighten the other too, even on a fit that never uses it). Default 10^(-sigdig), the same formula every other tolerance in this control uses.

• est="focei" with the default innerOpt="trust" no longer silently returns a fit whose outer bobyqa search never actually explored. On a real 2-compartment IV infusion steady-state model, the default outer rhobeg=0.2 made bobyqa’s initial quadratic model collapse: bobyqa reported normal convergence, but 4 of 5 population parameters never moved from their starting values at all. rhobeg=0.25 or larger escaped it and matched a SAEM reference fit closely (root cause not understood – n1qn1, FOCEi’s other inner optimizer, was never stuck at the same rhobeg on the same model). .bobyqa() now detects this directly (the final point never left its own starting exploration radius) and retries with a wider rhobeg, capped at 0.3, with a $runInfo warning when it engages – a pragmatic safeguard against the specific symptom rather than a fix for whatever makes bobyqa’s model collapse at that one radius.

Mixture models

• A mixture proportion now binds to the mix() component it was written for, whatever order ini() declares the proportions in. mix(a, p1, b, p2, c) means p1 is component 1’s share, but the theta slots were collected with which(names(theta) %in% mixProbs), which returns ascending theta positions – ini() order. Every consumer reads that index as “component m’s slot” (op_focei.mixProb[m], .getMixFromLog(), the back-transform that writes the estimates back), so declaring the proportions in a different order than mix() uses them ran component 1 on p2’s value and reported each component’s proportion under the other’s name. Measured at zero iterations with ini({p2 <- 0.20; p1 <- 0.70}): $mixProbabilities came back 0.2 0.7 0.1 instead of 0.7 0.2 0.1. A model whose ini() order already matched its mix() order – which is the usual way to write one – was never affected.

• Mixture proportions now get standard errors, under every covariance method that supports them: covMethod="r", "s", "r,s" and "imp". They were forced out of the covariance entirely (skipCov), so p1 reported SE = NA no matter what was asked for. This is the NONMEM 7 Technical Guide’s own construction, eq. (7.51)-(7.54): the mixture parameters’ information is the outer product of the same per-subject scores g_ia that give the gradient, and nlmixr2est’s S matrix already IS that outer product, so the block drops in once the per-subject score exists. The per-subject mixture score is taken analytically rather than by finite difference, so it costs no extra solves.

  The reported covariance is rotated onto the probability scale with the FULL mexpit Jacobian J = diag(p) - p p' – so $cov, the SE, the %RSE and the CI for p1 all sit on the same scale as the estimate, and the proportions’ cross-covariances (with each other and with the structural thetas) are carried across rather than dropped. Same principle as covFull reporting Omega on the natural variance scale instead of chol(solve(omega)).

  The confidence interval is taken on the logit scale, expit(logit(p) +/- z*SE/(p(1-p))), from the SE actually reported. The generic symmetric est +/- z*SE interval walks out of (0, 1) for a proportion (a fit reported p1 = 0.648 (-0.045, 1.34)) and was built from the covariance before the rotation, so it did not agree with the SE beside it.

  That rotation carries a factor of p(1-p), so a proportion sitting near 0 or 1 gets an SE that shrinks toward zero. It is the right delta-method answer but it reads as certainty, when in fact the symmetric Wald interval has stopped being meaningful there, so the fit’s $runInfo now says when a proportion is at a boundary.

• The S matrix is no longer singular for mixture models, so covMethod="r,s" stops silently degrading to "r". foceiS() built each subject’s score by finite-differencing component 0’s likelihood instead of the marginal log(sum_m p_m L_im), so a parameter entering only another component got an exactly-zero score. Measured on a 3-component fit, two of six diagonal entries came out at 1.9e-18 and 6.6e-10, the S matrix was reported non-positive-definite, and the sandwich was dropped. Each perturbation now re-optimizes every component before combining, in the components-serial / subjects-parallel order the rest of the mixture code uses.

• covMethod="imp" builds its Monte-Carlo proposals for every mixture component rather than component 0 alone, so the importance-sampling objective it differences is the mixture marginal. Without this the proportions’ directions were exactly flat and their SEs came back as 0.

• covMethod="analytic" now declines a mixture model and falls back, instead of reporting a single-component observed information as if it were the mixture’s. The augmented sensitivity model differentiates one component’s conditional likelihood and has no mixture-proportion block at all.

• est="saem" mixture fits can report a mixture-proportion SE as well. saem leaves the proportions out of the parameter vector its kernel converges (they are updated by a separate EM step), so its linFim/fim/sa covariance has no mixture rows to extend; the (7.51) block is appended from the fit’s own posterior responsibilities. The cross terms (7.52)-(7.54) are NOT formed – they need per-subject scores for the other parameters on the same footing, which that covariance does not expose – so the block is uncorrelated with the structural parameters and its SEs are mildly optimistic.

  The block is only reported when the fit is actually at the mixture’s fixed point, judged by the score statistic s' I^-1 s (which is on a chi-square scale, so unlike a tolerance on mean(r) - p it does not loosen as the number of subjects grows). An information matrix describes the precision of a maximum-likelihood estimate, and away from that point it is a confident-looking number attached to something that is not one; instead the fit’s $runInfo says the SE was skipped. est="saem" currently lands far from it (nlmixr2/nlmixr2est#1058), so in practice this declines today and will start reporting once that is fixed.

• est="focei" estimates the mixture proportions under a gradient-based outerOpt. mixGrad() supplies an analytic value that short-circuits the finite difference in numericGrad(), and it chained the per-subject responsibility sum through the diagonal of the mexpit Jacobian (dmexpit() returns only the diagonal) while dropping the -2 of the objective scale. That flipped the sign at every number of components, so the line search rejected the first step and the proportions never moved off their initial values; from three components up the magnitude was wrong too. The full Jacobian collapses to -2 * sum_i (r_il - pi_l), which is what it now uses – checked against a central difference of the objective to eight significant digits at three components, and equal to a literal evaluation of the NONMEM 7 Technical Guide’s own equations (1.194) and (1.197) for the mixture-proportion gradient. Fits left on the default derivative-free outerOpt="bobyqa" never reached this code and are unchanged, as are the reported standard errors (a mixture proportion is skipCov).

• fit$etaMat no longer carries the mixnum column that $eta gains for a mixture fit. Every consumer that hands it back as foceiControl(etaMat=) compared neta + 1 columns against the model’s neta and stopped with “The etaMat must have the same number of ETAs (cols) as the model” – so $cov, addCwres(), the FO objective and re-fitting a fit were all failing for every mixture fit, in most cases inside a try() that swallowed it.

• A mixture fit’s $ui now carries the mixture probability on the probability scale rather than the mlogit scale it is estimated on. fullTheta is left on the estimation scale and only $theta was back-transformed, so $ui reported p1 = -0.847 where fixef() reported 0.3, and re-fitting the fit failed its own ini() validation (“the probabilities in a mixture must sum to a number between 0 and 1”).

• est="vae" can fit a mixture (mix()) model. It could not before: the two defects above stopped every such fit during assembly, so only the training loop had ever run with a mixture.

• est="vae" estimates the mixture proportions. They were read once from ini() and never updated, so the reported proportion was whatever the model started at. They are now estimated on the mlogit scale through their own analytic gradient – the same closed form focei’s corrected mixGrad() uses, each component’s summed responsibility against its expected count – consumed by the same Adam loop that trains the encoder. The gradient agrees with finite differences to 1e-4 at two components and at three. On simulated data with a 3:1 split started from 0.5, the fitted proportion comes back at 0.74.

• est="vae" no longer treats a mixture proportion as an ordinary structural theta. p1 appears inside mix(a, p1, b), so it passed the “does this theta appear in a model expression” filter and became a nonMuTheta regression parameter, moved by bobyqa against the (-Inf, Inf) bounds it carries in iniDf – writing values like p1 <- 10.6 back into ini(). est="npag" already excluded them for the same reason.

• est="vae"’s mixture objective is the marginal -2LL. It exponentiated -obj/2 and negated the result twice, so it marginalized a square root likelihood and carried a spurious -log of the winning component’s proportion. Both errors cancel exactly when there is one component, and again when the components are identical under a uniform proportion, so no existing test could see it. The prediction-model M-step scored candidate thetas the same wrong way, and the two must agree or the M-step optimizes a different function than the reported ELBO.

• est="vae"’s encoder now characterizes every subject-component pair, each at its own random effect, the way the other mixture methods do. It produced a single posterior per subject which was then reused for every component, so the components were compared at a random effect that had been fitted to none of them. The component enters the encoder at its head, alongside the covariates, so a mixture model’s encoder head is nMix inputs wider. Everything reported per subject – the random effect, the predictions, the residual variance – comes from the selected component.

• est="vae"’s encoder is trained on the gradient of the objective it optimizes. The objective was the marginal over components but the gradient handed to the encoder was the single best component’s; it is now the responsibility-weighted mean. Predictions and the residual variance stay on the best component, since the closed-form residual step is a moment estimator.

• A subject whose components all fail to solve no longer improves est="vae"’s objective. It contributed nothing at all; it is now charged the same bad-solve penalty the focei mixture likelihood charges.

• est="vae" reads the mixture proportion on the scale the inner problem reads it on. The prepared theta vector held the raw ini() probability while the inner problem passes that slot through mexpit(), so p1 = 0.3 was used as mexpit(0.3) = 0.574.

• A mixture fit’s table reported mixest, mixnum and the result of mix() itself as 0 for every row, and PRED/IPRED were computed from those zeros – silently wrong predictions rather than an error (#1041). The prediction model is built through symengine, which expands the mix() call away, and rxode2 then no longer read the model as a mixture at all, so the per-individual component never reached the solve. Fixed in rxode2 (nlmixr2/rxode2#1358); this release stops working around it.

• A theta that saem does not estimate as a parameter of its own, declared before a mu-referenced population parameter, shifted every eta after it onto the wrong parameter in the model saem solves for its table. The eta to theta map is an index into the SAEM estimation parameter vector, and it was used to subscript the model text built in ini() order; the two differ by exactly those thetas. This is not mixture-specific – a mu-referenced COVARIATE parameter is dropped from that vector as well, so an ordinary covariate model that declares tcl.wt before tv had the volume’s eta land on tcl.wt and the volume get none, and its table lost the volume’s between-subject variability entirely (v came back constant, and IPRED with it). Now paired by name.

• The mixest/mixnum iCov handed to the table step is rejected by rxode2 when its ID is a factor, which it always was: it is built as an integer and output creation re-levels every ID in the fit environment afterwards. The whole table step was then dropped and the fit came back without a table.

• A rejected iCov no longer takes the table step down with it – the retry without it now covers the rxode2 messages that can actually be raised, and the retry says so in the fit’s $runInfo rather than quietly handing back a table whose mixture columns are all 0 (which is what happens on an rxode2 without nlmixr2/rxode2#1358).

• The post-hoc correction of the mixest/mixnum/mixunif output columns is removed. It only ever fired for columns literally named me, mn and mu, so a model that named them anything else kept the zeros; and now that the solve is right it was corrupting correct values – it wrote the per-subject component into the column holding mixnum, which is the component count. A model that reads mixunif in an expanded prediction model gets the supplied component back rather than a fabricated 1/nMix; use mixest.

• est="saem" now reports a mixture proportion that agrees with the fit’s own posterior responsibilities (sum_i (r_i - p) == 0) and with what the data identifies. Three things were wrong: the shared parameter-table hook mlogit back-transformed saem’s proportions, which are already on the natural scale, so the reported value was expit(p) and could not equal the responsibilities it was averaged from; mixProbMethod="regress" classified each subject at that subject’s phiM draw, whose fixed-effect-only columns carry a search variance of 1 rather than a real BSV, so the draw swamped the between-component signal and misclassified a fifth of the subjects even with 8-fold separated components; and an eta shared by every mixture component was marked as owned by whichever component mentioned it last, sending shared-eta mixtures down the split-ETA code paths. Every mixProbMethod now reports the proportion at the score-zero point: one exact M-step at the final responsibilities, rather than a value still carrying the annealing or Dirichlet-style shrinkage that stabilizes the trajectory.

Crashes and stability

• An over-parameterized est="saem" fit no longer dies with “nearest PD calculation failed” after the last iteration (#923). A singular final Omega was already projected to the nearest positive-definite matrix before the residual/table step, but that projection itself errors on the fully degenerate cases – an all-zero, non-finite, or negative-definite Omega – which is exactly what an over-parameterized model produces. Those now fall back to a floored diagonal so the completed run is returned as a fit, with a $runInfo note saying the Omega was singular; the reported Omega is left as estimated. A collapsed or non-finite saem Omega is also reported in $runInfo on its own, and a failure while assembling the fit object retries once without the table step rather than throwing the finished run away. When nmNearPD() does error it now says which degenerate case it hit.
• The shared solve pool’s lhs-width probe could segfault instead of declining. It verifies a model by calling that model’s calc_lhs, and generated calc_lhs dereferences per-subject pointers that are bound by a solve, not by building the pool – so an inner problem that had been set up but had not solved yet crashed inside the check written to make a mismatched model fall back safely. The probe now binds the subject itself, and additionally verifies that the pool holds the model’s states and that its parameter layout matches the one the pool’s parameter vector was filled with (calc_lhs reads it by index, so a same-width model in a different order mis-reads). Fits are unchanged; .odeSwapInfo() reports the new npars/parLayoutOk columns and the probeIniN/probeDenyN counters. The lhs column map, which is installed separately from the model it describes, is now checked against that model’s width at the pooled entries as well.
• est="npag" / est="npb" now exclude an observation or a residual parameter whose endpoint cannot be determined from the residual moment warm start, instead of attributing it to the first endpoint. Both the observation CMT lookup and the residual parameter’s condition lookup resolved “no match” to endpoint 0, which is indistinguishable from the correct answer for a single-endpoint model; on a multi-endpoint model that warm-started (and, where there is a lone scale per endpoint, estimated) one endpoint’s residual SD from another endpoint’s residuals. Anything dropped this way is now reported in $runInfo rather than being silent.

Other

• An omega block declaring one of its covariances at exactly 0 (for example eta.ka + eta.cl + eta.v ~ c(0.1, 0.01, 0.1, 0, 0.01, 0.1)) no longer aborts the fit with theta has to have N elements. The block’s cholesky factor is dense, so that 0 cannot be held; it is now estimated from ~0, as a 0 element of a NONMEM $OMEGA BLOCK is, and $runInfo names the random effects involved. The same applies to a correlated pair that is not adjacent in eta order (eta.a with eta.c, eta.b between them), which used to be refused as well. The FOCEi family, est="vae" and the general-likelihood inner driver all took the same route (#1079, rxode2#1365). A covariance declared at 0 in a two-eta block is unchanged: it leaves the two etas uncorrelated, as it always has.

• foceiControl(warm="save") now restarts the n1qn1 inner problem from the curvature the subject’s previous inner solve left, as it was always meant to. It reconstructed that Hessian from a buffer it had just zeroed, so n1qn1 was handed an all-zero factorization and self-initialized on every inner solve – the option reused nothing since FOCEi was first imported (#1043). A single-eta model was additionally unseedable because the one-by-one case multiplied the factorization back out as a zero matrix. With mceta sampling the eta=0 floor pass now gets that same seed rather than self-initializing, so it stays the run mceta=0 would have made. The previous self-initialized behavior is available as the new foceiControl(warm="none"), and warm="save" reuse is reported in the fit’s $nWarmSave.

• nlmixr2() names a model the way rxode2() does, through rxode2::rxModelNameFromExpr(): a symbol keeps its name, a call becomes its text or the name its rxModelName() method gives, and an anonymous model function is unnamed (NULL). It used to name the model with as.character(substitute(object)), one element per part of the call, so $modelName was "function" for an anonymous model function and "readModelDb" for nlmixr2(readModelDb("PK_1cmt")). The internal .nlmixr2objectNameAssign() that carried that name is removed.

• saemControl(nu=) was silently ignored. saem rebuilds its control from a plain list, and that path copied nBurn, nEm and nmc out of mcmc but only validated nu, so every fit ran the default c(2, 2, 2) kernel iterations whatever was asked for. Measured, saemControl(nu = c(10, 11, 12)) reached the estimator as c(2, 2, 2) and fitted identically to the default.

• nIter=0 in impmapControl(), impControl() and qrpemControl() crashed R at fit time; it now runs an E-step-only evaluation (#1091).

• est="saem" and est="npb" seed their random draws sequentially, one seed per individual per step, computed from the iteration so any draw’s seed is known from its position. The phi1 and phi0 SAEM MCMC blocks used to draw from the same seed, as could neighboring mixture components. Seeded fits give different, still reproducible, results.

• predict(fit, newdata) solves with the fit’s own rxControl, the same options nlmixr2(fit, data, "predict") uses; options passed to predict() still override it. The two gave different predictions for an ODE model.

• addNpde() and vpcSim() give the same result on every call and at any thread count, even when an earlier rxSetSeed() left rxode2’s own seed sequence in force.

• The SAEM Gaussian-quadrature objective caps its grid at 25 nodes per dimension instead of crashing R, and saemControl(nnodesGq=) rejects a value above 25.

• est="saem" estimated a boxCox()/yeoJohnson() lambda without the transform’s log-Jacobian, both in the closed-form residual step and with a general likelihood (dnorm(), t(), cauchy(), also covering lnorm(), logitNorm() and probitNorm()), so lambda and the residual SD were fit to the wrong density; the reported objective of a general-likelihood fit also used the starting lambda.

• est="ifocei" and est="mfocei" fit models with inter-occasion variability. They were the only two FOCEi family methods that did not declare IOV support, so the occasion parameters were never expanded and the fit errored with default method not implemented for type 'list' (#1083).

• est="flaplace", "mflaplace", "iflaplace", "fagq", "mfagq" and "ifagq" fit a correlated inter-occasion block. These are thin delegates to the base Laplace/AGQ methods, but were missing from the list of methods that honour a repeated (SAME) omega block, so they refused a model their own base method fits (#1083).

• est="saem" fits with a dnorm(), t() or cauchy() endpoint reported the log-density instead of the prediction as PRED/IPRED (and the residuals derived from them) in the fit table (#1084).

• A focei inner ETA solve that has spent every etaNudge/etaNudge2 restart and still failed now falls back on draws from Omega (foceiControl(etaRestart=), 4 by default, 0 to disable). Every nudge sets every ETA to the same constant, which explores poorly once the inner problem has more than one basin; the draws are starting points from the distribution the ETAs come from. Measured on #1044’s model at a displaced parameter set, subjects that ended with every attempt spent fell from 45 to 25 of 300 and the objective from 687874 to 139103. The draws are taken once per fit from rxode2’s seeded engine, so the objective stays a function of theta alone, and they are read only after a solve has already failed – a fit whose inner solves converge is bit-identical with the fallback on and off. This applies to innerOpt="trust", which reports a convergence verdict per solve; n1qn1 reports none, so its own restart cascade is unchanged (#1044).

• saem reports a fix()ed eta variance as the value it was fixed at. The reported omega was snapshotted before the fixed values were restored, so it carried the M-step’s unconstrained estimate instead – fix(0.3) came back as 0.318 while the fit itself correctly sampled with 0.3 (#1073).

• covMethod="analytic" for FOCE and foce="foce+" no longer carries the inner solver’s residual score into the observed information. The FOCE kernel uses the general total-derivative form, whose last term is Phi_eta . eta_ab; it evaluated Phi_eta in full, but the FOCE inner problem zeroes S_FOCE = Omega^-1 eta + sum(q0 a) by construction, so only the interaction remainder Phi_f - q0 belongs there. The rest was the inner tolerance multiplied by a term that is not small. The FOCE assemblers also re-solved each subject’s EBE before building R, where the FOCEI assembler uses the fit’s own; they now agree. On an additive model, where the two methods must coincide, the foce+ and FOCEI observed informations agreed to 1.2e-2 at the same EBEs and now agree to 1.1e-13 (#1056).

• est="focep" (foce="foce+") converges at the default sigdig. FOCE+ polishes the inner optimizer’s eta onto the truncated-score root it defines its EBE by, and that polish stops once the score reaches the noise floor the solve tolerance buys – near 1e-3 at the default sigdig=3, well above the 1e-9 it asked for. A subject that stalled there reported its likelihood as NA, which the outer search read as a cliff: on theo_sd it stopped 4.8 objective units high (121.560 against 116.804) with the omegas barely off their starting values, and the observed information at that point had a negative eigenvalue so covMethod="analytic" refused to install it. The polish is now best-effort – the eta it was handed is the inner optimizer’s own answer and is always usable, so every exit path keeps the best point found and evaluates the likelihood there (#1069).

• saem refuses a model whose random effect has no population parameter of its own – added to none, or sharing one with another random effect – naming the random effects, instead of fitting the model without them and then failing with “subscript out of bounds” while assembling the reported omega at the end of the run (#1047).

• The analytic outer gradient takes Omega^-1 and its estimation-scale derivatives from the native Cholesky map the objective already maintains (verified against the R handle on first use, which stays the fallback), so a gradient evaluation no longer calls the rxSymInvChol closure.

• Swapping a peer model into the shared FOCEi solve pool (the augmented outer-gradient model for every analytic gradient, the outer Hessian’s probes, the AGQ node model) installs its event-sensitivity shape from the registry with a C call. Only the first swap of a model goes through R (rxEventSensLoadModel()); the shape it installed is read back and reused. Before, every swap re-derived it in R, which kept every batch boundary on the R side of the derivative passes. .odeSwapInfo() reports esInstallC/esInstallR.

• A model containing mtime() can be fit again, with every estimation method. etTrans() materializes the modeled times as EVID 10-99 records (TIME=0, AMT=NA) and $dataSav persisted them, so re-translating it for each estimation solve rejected them as doses with a missing amount. They are now dropped, and since rxode2::rxS() keeps only the mtime variable name, the mtime() declarations are re-emitted into every generated model so the modeled times survive estimation and the mtime variable stays defined. The extra records the solve regenerates are also kept out of the fit’s output table and out of the nlme objective (#919).

• An mtime() whose time depends on an estimated parameter (a boundary that moves with a theta or an eta) now contributes to the sensitivities. The declaration is loaded into symengine as an ordinary assignment, so the switch time is differentiated like the same branch written in place (ifelse(t < exp(tsw), ...)); it used to reach symengine as a free symbol, which made every derivative through the boundary zero – the EBEs for an eta the boundary depends on stayed pinned at their initial values, and foceiControl(fast=TRUE) reported a gradient of exactly 0 for the theta.

• A model that declares mtime(v) and also assigns v as an ordinary variable is refused instead of silently using the declared value. The declaration is re-emitted at the top of every generated model, so a later mtime() reading v got the declared value where rxode2 gives it the reassigned one.

• fit$cor returns NULL instead of erroring when the fit has no covariance (covMethod=""), matching fit$cov (#1038).

• Fit accessors, the control getters, setCov() and print() look items up only in the fit environment itself. A fit reloaded by nlmixr2save is parented on the global environment, so fit$cov resolved to stats::cov (which broke print()), fit$ranef to nlme::ranef, and $mixNum, $mixList or $parHist to a variable of that name in the user’s workspace (#1038).

• print() on a fit shows the fixed-parameter correlation line again. It was gated on exists("cor", fit$env), which is never true for a fit that has not been through a save/load round trip, so a strong theta correlation was never reported (#1038).

• A focei fit of a model whose dosing depends on an eta – f(), alag(), rate(), dur() – now says so when rxode2’s analytic event (“jump”) sensitivities cannot be installed, instead of silently returning a fit whose dosing etas never left their initial values while their omegas stayed finite. The event-sensitivity mode already rides with the cached model bundle; it now also survives a second deflate/inflate round trip, and the model bundle records which etas enter a dosing expression so the fit can tell a model that needs the jumps from one that does not (#1016).

• saemControl(covMethod="sa"|"fim") keeps its own Omega standard errors on a diagonal-Omega model. The variance parameters the analytic FIM cannot cover (a non-additive endpoint’s residual error) are still taken from the linearized FIM, but that splice replaced the WHOLE variance block, so a near-singular residual pair propagated into the Omega rows – on a two-endpoint add()+prop() / add() model om.eta.ka reported an SE of 560 against an estimate of 1.1, where the analytic FIM gives 0.51. A declared Omega block still takes the whole block from the linearized FIM, because the Louis score only ever sees the diagonal of Omega (#1022).

• saemControl(covMethod="sa"|"fim") no longer reports a standard error for a fix()ed additive residual error. The SAEM kernel fills an endpoint’s residual slot whether or not the value is estimated, so a fixed add.sd came back with a covariance row – printed as a back-transformed 95% interval on a value the fit never estimated – and the remaining parameters were given the marginal instead of the conditional information. Found while fixing #1022.

• foceiControl(fast=TRUE) no longer converges to a non-stationary point when an external likelihood contribution is registered (nlmixrRegisterLikContrib, e.g. from nlmixr2nn). The analytic outer gradient re-derives d(objective)/d(theta) from the model sensitivities alone, so a contributor that adds a theta-dependent log-likelihood term left the objective right and the gradient wrong. It now declines to the finite-difference gradient for a contributor that writes back llik or d(LL)/d(eta); a pure observer changes nothing and stays on the fast path (#1051).

• covMethod="analytic" judges positive definiteness on the full theta + sigma

  • Omega matrix even when only the structural-theta block is installed. A submatrix of the inverse of an indefinite observed information can look positive definite on its own, so a covFull = FALSE fit (or setCov(fit, "analytic")) could install a covariance derived from a point that is not a local minimum (#1055).

• A focei fit reports standard errors that match its own covariance again. .foceiInstallFdFullCov() replaces $cov with the full theta+omega matrix after the C++ step has already derived popDf$SE from the native theta-only covariance it discards, so parFixedDf$SE described a matrix the fit no longer held and a setCov() round trip silently changed the reported SEs. The parameter table is now refreshed from the covariance actually installed (nlmixr2extra#125).

• addCwres() works on a fit that already reports the focei (or foce) objective function and has no CWRES, instead of stopping with “objective function ‘FOCEi’ already present”. setOfv(fit, "focei") adds that objective function row without the residual columns, so setOfv(fit, "focei") followed by addCwres(fit) – and any estimation method that reports the focei objective function of its own, such as the nlmixr2bayes methods run with ofv="focei" – had no way to add CWRES at all. addCwres() now adds the residual columns and leaves the objective function row the fit already carries alone. The table step also calculates CWRES for any fit whose own objective function is already the focei one, since there is no way to add them afterwards.

• The analytic outer Hessian no longer drops to a finite-difference Hessian the first time a 3rd-order probe fails to solve. The probes run tighter than the fit (1e-12), which is where they fail, so the probe tolerance is now loosened a rung at a time (foceiControl(outerOdeRecalcFactor=), up to outerMaxOdeRecalc rungs, never looser than the fit’s own tolerance) and retried before falling back – the escalation the inner problem already takes on a bad solve (maxOdeRecalc). The rungs taken are reported as the fit’s $nHessTolRelax.

• foceiControl(fast=) no longer moves a maxOuterIterations = 0 fit’s ETAs. That fit evaluates the analytic outer gradient once so .foceiGradDirect() has something to report, and the evaluation ran an extra inner optimization pass per subject before the final one; under the default warm start (mceta < 0 keeps the last eta) the inner solve converges only to its own tolerance, so the extra pass shifted the reported ETAs and everything derived from them. On theo_sd at sigdig = 4 the covMethod = "analytic" standard errors differed by 1.8e-4 relative between fast = TRUE and fast = FALSE; they now agree to 2.7e-12.

• The FOCEi family (focei, foce, focep, laplace, agq, posthoc, and their i/m prefixed variants) now applies the M2/M3/M4 censoring correction to a t(), cauchy() or dnorm() endpoint instead of silently scoring a censored row with its ordinary, uncensored density (#992, completing #979 which covered the nlm family only). Such an endpoint is compiled to a scalar log-density, which hid the location, scale and degrees of freedom the correction needs; the inner model now carries them (together with d(f)/d(eta) and d(R)/d(eta), so the inner eta gradient is corrected as well, not just the objective). A censored t()/cauchy() fit under these methods changes its objective and its parameter estimates, and no longer emits the “censoring ignored” note in $runInfo. An ar() endpoint is excluded: its reported scale is the marginal, not the conditional, one.

• foceiControl(fast=TRUE) now downgrades to fast=FALSE for a CENSORED log-likelihood endpoint, as the documentation already said it did. Both the augmented outer-gradient model and the exact second-order inner Hessian differentiate the uncensored log-density, which the M2/M3/M4 correction replaces for a censored row; the outer gradient already refused such a fit at run time, but the inner Hessian did not.

• The persisted FOCEi model cache (rxUiGet.foceiModelCache()) now keys on the nlmixr2est version. That cache lives in rxode2’s user cache directory when rxode2::rxCreateCache() has been run, so it survives an upgrade: any release that changes the generated inner-model text (the censoring columns above, for one) was silently ignored for a model already cached there.

• est="impmap"’s inner Hessian (impGetHessian), which builds the importance-sampling proposal, could read a stale cached linCmtB() Jacobian on a linCmt() model with a non-mu structural theta (a theta with no random effect, e.g. ka fixed but V estimated on log scale). Several compiled peer models share one solve pool (odeSwap); linCmtB() caches its Jacobian in a field gated by rx->ndiff, a process-global that odeSwapSolveInd() never restored per peer, so a solve could read a Jacobian built for a DIFFERENT peer’s structural-parameter set. The resulting proposal was artificially wide, which masked a real heavy tail (Pareto k-hat) as healthy rather than repairing it. odeSwapSolveInd() now restores each peer’s own ndiff before every solve. That restore is itself a write to a field on the single shared solve struct, so impGetHessian’s parallel per-subject loop (a subject that falls back to the doFD/pred path can pick a different peer, and so a different ndiff, than a subject still on the plain inner path, concurrently) now serializes that write-then-solve window whenever the fit has a peer that needs it, rather than risk one subject’s solve reading another’s in-flight ndiff.

• The FOCEi family (focei, foce, foi, fo, posthoc, agq, laplace) no longer refuses an ordinary population dataset with “dataset too large for this mixture model configuration” (#1010). The per-subject residual variance block gVid was sized as nall * (nMix + 1) squared, but it holds one nobs_i x nobs_i matrix per subject, so it only needs sum(nobs_i^2). nall counts dose (and evid=2) records as well, and (sum x)^2 exceeds sum(x^2) by roughly the number of subjects, so the request was inflated by several orders of magnitude; the > 65535 guard added to turn the resulting 32-bit overflow into a clean error was therefore rejecting fits that need well under a megabyte. The block is now sized from the observation counts it is actually indexed by, and the guard is replaced by an overflow check on the allocation the setup really makes. The per-subject offset accumulators in that setup are now size_t as well, so a subject with more than 46,340 observations no longer wraps its own nobs_i^2 stride.

• FOCEi’s geta block (the per-subject current eta vector) was allocated neta doubles inside the same setup, but it is indexed once per subject – it needs (neta + 1) * nsub like the ten per-subject eta arrays that follow it. Every subject past the first therefore wrote its eta through storage belonging to gtryEta (the trial eta of the eta-reset/nudge path). Both arrays now get their own block.

• A FOCEi fit of a model with no random effects (neta == 0) sized its per-subject thetaGrad block with neta * nsub, which is zero on exactly that path, so the block had no storage and the per-record log-likelihood array started at the same address. The gradient writes ran past the allocation whenever the parameter count times the subject count exceeded the event-record count.

• A fit’s per-record llikObs came back reversed by subject. Its per-subject offsets were a running total taken in the FOCEi setup loop’s backwards order, while the array is handed to R as one contiguous block in record order, so the first subject was given the tail of the buffer. With unequal observation counts the subject boundaries did not line up either. The offsets now follow the subject order.

• est="saem" with saemControl(nMix > 1) (mixture SAEM) inverted the propT()/powT() (transformed-basis) vs. plain prop()/pow() (raw-basis) proportional/power error term in the two MSAEM-only E-step helpers, mixObsLoss() (softmax mixture-component responsibility weights) and mixNaiveClassify() (chain-init classification) (#982). Every other E-step call site picks the boxCox/yeoJohnson-transformed prediction when propT()/powT() is used and the raw prediction otherwise; these two passed the pair in the opposite order, so a propT()/powT() model was scored against the raw prediction and a plain prop()/pow() model against the transformed one, biasing mixture-component assignment for any nMix > 1 fit with a transformed residual model. nMix == 1 fits are unaffected.

• Every NLM-family method (nlm, bobyqa, newuoa, uobyqa, n1qn1, lbfgsb3c, optim, nlminb) silently mis-scored any M2/M3/M4-censored normal endpoint, whether or not ar() was present (#976). NLM forces every normal endpoint through a log-likelihood (dnorm) path so a single scalar objective can be emitted; that path always set rx_r_ ~ 0 (a full log-density has no separate variance to report), but the censoring correction reads rx_r_ as a real variance, so a hardcoded zero corrupted the correction for every censored observation. rx_rll_ (the standard deviation actually used to build the log-density) is emitted immediately beforehand in the same branch, so it is now squared back into rx_r_ instead of being discarded.

• Every NLM-family method also silently ignored M2/M3/M4 censoring on a t()/cauchy() endpoint entirely – a censored row was scored with its ordinary (uncensored) density (#979). Added a Student-t/Cauchy CDF-based correction (doCensT1(); cauchy is Student-t with nu=1, so one function covers both) alongside the existing normal one. FOCEi/FOCE/AGQ/Laplace and every other generalized-likelihood distribution (pois, binom, beta, and so on) still silently ignore censoring for now, but a fit now warns when that combination is used instead of staying silent.

• The rx_r_ fix above (.fixCensRNuLine(), R/focei.R) rebuilt rx_r_ by re-inlining rx_rll_’s defining expression, which for a transformed propT()/powT() error model contains the symbol rx_pred_ – a symbol already overwritten with the scalar log-likelihood by that point in the same branch, silently corrupting the variance for any censored propT()/powT() endpoint (found by an independent Antigravity review). Fixed by referencing the already-computed rx_rll_ variable instead of its expression.

• An ar() endpoint’s censoring correction now uses a self-consistent marginal (not the exact AR(1)-conditional) mean/variance for M2/M3/M4 scoring – a real improvement over the previous corrupted state, but still an approximation for ar() specifically; tracked as #1001.

• est="saem" with a pow() residual error model (rmPow/rmAddPow/ rmPowLam/rmAddPowLam) never applied the estimated power exponent in the E-step’s MCMC-acceptance likelihood – it always scored proposals as if the exponent were 1, no matter what the M-step estimated (#972). The M-step’s own objective did use the exponent, so the two steps disagreed about what model they were fitting: the power estimate collapsed toward 0 while the proportional-SD estimate inflated to compensate. The per-observation combined-error-SD builder now applies the current power exponent (refreshed from the M-step every iteration, not just read once at setup), matching the M-step’s combined1/combined2 formulas.

• est="saem" scored a censored (M2/M3/M4) observation with the wrong sign, the wrong scale, and the untransformed DV, so a censored row’s contribution to the chain’s acceptance could drive the fit away from, rather than toward, the true parameters (#876). The residual-error M-step also still counted a censored row’s recorded LOQ/limit as if it had been measured, biasing the residual SD low relative to focei on the same data (#916). Both are fixed: the E-step now scores a censored row the way focei’s inner likelihood does, and the M-step now simulates each censored row’s value from the truncated normal implied by the current fit (data augmentation, Samson/Lavielle/Mentre 2006) before building the residual sum of squares. The truncated-normal draw itself uses the same Botev

  (2015) minimax-tilting algorithm CWRES’s censored-observation simulation already relies on (censResid.h’s truncnorm(), via rxode2’s rxRmvn), ported to the seeded per-thread engine this file already uses everywhere else, rather than a plain inverse-CDF draw – which loses precision once the truncation bounds are a few SDs from the mean, the regime a BQL row’s bound often sits in.

• est="npag"/est="npb"’s residual-error moment (npResidMoments()) counted a censored (M3/M4) observation’s recorded LOQ/limit as if it had been measured, the same bug shape as #916 but in the nonparametric methods (#978). For the common configuration – one add()/prop() scale per endpoint, no regressor theta – that biased moment is installed as the final residual-error estimate with no further optimizer correction, so a censored row could badly distort it (an extreme recorded LOQ was measured to inflate est="npag"’s proportional SD from 0.03 to 14.7, and est="npb"’s from 0.03 to 8.1). A censored row’s DV is now excluded from the moment entirely, matching how the function already excludes an endpoint with no defined prediction to take a moment of. Excluding the row still leaves the moment a biased-low estimate of the true residual variance whenever the data really is censored (the variance of a truncated normal is always less than the untruncated one), so an endpoint that had any row excluded this way no longer takes the direct-install fast path either – it is refined against the already censoring-aware ELS objective instead, the same way a regressor theta already was.

• foceiControl(fast=TRUE) could converge to a different fit than fast=FALSE when a transform-both-sides (lnorm/boxCox) endpoint’s untransformed prediction was non-positive at some observation – for example a depot model’s TIME==0 row, where the central compartment is exactly zero (#867). rxode2 floors such a row’s transformed prediction at a constant, so the objective is locally flat there, but the analytic (“fast”) outer gradient kept differentiating the unclamped expression – not the gradient of the objective actually being minimized. Such a subject’s analytic gradient contribution is now detected and replaced with the same per-subject finite difference already used for a failed augmented solve, which differences the objective rxode2 actually evaluates.

• est="saem" scored a censored (M2/M3/M4) row on an ar() endpoint against the marginal normal distribution instead of the AR(1) conditional one that its uncensored neighbors already used (#918). arDYFhyp whitened a discarded copy of the prediction/SD to build the uncensored loss, then handed the censored-loss calculation the original marginal values. The whitened prediction/SD are now kept and passed through, so a BQL row after an AR-active observation is scored consistently with the rest of its chain.

• Post-estimation machinery no longer refuses a model whose ini({}) declares prior distributions (#938). Two parts:

  • The "output" and "posthoc" pseudo-methods now declare nlmixr2Priors = "all". Neither estimates anything – they evaluate an already-specified model and build its tables – so there is no prior they could silently ignore.

  • The internal zero-iteration est="focei" re-entries behind setOfv(), addCwres() and the impmap objective recompute now run with the prior gate bypassed (scoped, restored on exit). By the time they run, the priors were already accepted or refused by the estimation method that produced the fit; refusing again only broke post-processing. A user-initiated est="focei" on a prior-carrying model is still refused.

  This was latent while no estimation method declared prior support; it would have broken assembling a finished fit for the first method that does.

• The Hessian-based covariance for every NLM-family method (nlm, bobyqa, newuoa, uobyqa, n1qn1, lbfgsb3c, optim, nlminb, nls) was 4x too large (SEs 2x too large), for both built-in error models and custom ll() likelihoods. .nlmFinalizeList() converted the Hessian to a covariance using the R = 0.5*Hessian, cov = 2*solve(R) convention that is only correct when the Hessian is of a -2*log-likelihood objective (as FOCEI/FOCE/SAEM/ Laplace/AGQ use). NLM-family methods instead optimize a plain -log-likelihood, so their Hessian already is the observed information; halving then doubling compounded into a 4x inflated covariance. Point estimates, the objective value, and the log-likelihood were unaffected.

• est="saem" with covMethod="sa"/"fim" reported a nonsense (~1e-6) SE for a theta with no random effect (a “phi0” parameter, e.g. a covariate-free structural parameter), and on a model mixing mu-referenced and non-mu-referenced thetas could attribute the Fisher information rows to the wrong parameter names entirely (#906). The stochastic-approximation kernel orders its Fisher information [mu-referenced thetas][non-mu-referenced thetas], not iniDf/model order, and a non-mu-referenced theta’s mu information is a pseudo-variance the algorithm decays toward 0 (it is a fixed effect carried as a degenerate random effect), which blows up to a near-zero SE when inverted in place. The Fisher information block is now read out in its actual order, and non-mu-referenced theta rows are dropped before inverting and their SE spliced in from the linearized FIM instead. The kernel keeps a Fisher information row for a fix()ed theta too, so that drop is computed against the FIM’s raw row order rather than a fixed-filtered one (a fixed theta ahead of the dropped row previously shifted every later position and dropped the wrong one); a model whose Fisher information order cannot be verified (a mu-referencing covariate, or an old cached fit) now refuses "sa"/"fim" and falls back to the linearized FIM instead of reporting from an unverified order.

• est="saem" on a boxCox()/yeoJohnson() model (estimated or fixed()) fit at the identity transform instead of the declared one (#914). Two compounding defects: the kernel’s working lambda was never refreshed from the M-step’s lres, so every _powerD() call used the transform’s initial value (1) all fit long, and transMat reported that same wrong value instead of the fitted/fixed lambda; and the pure additive-plus-lambda and proportional/power-plus-lambda residual models never zeroed the error component they do not use (bres/ares), leaving it at its nonzero default and corrupting the residual SD g = ares + bres*|f|. Together these biased theta, BSV, and the residual SD (the BSV would collapse and the residual SD would inflate to absorb the untransformed data), matching neither the declared model nor FOCEi’s fit of the same data. The unzeroed component defect also affected the plain (non-boxCox/yeoJohnson) pure power error model pow() with no add()/prop() term, which shares the same fix.

• SAEM’s analytic Fisher information (covMethod="fim"/"sa") had exactly one residual slot no matter how many endpoints a model declared, so a multi-endpoint fit’s residual score/Hessian always came from whichever endpoint the internal loop happened to process last, divided by the first endpoint’s residual variance (#893). Because that slot is coupled to the fixed-effect/BSV block through the full Fisher information matrix, this corrupted the reported theta and Omega standard errors for any multi-endpoint fim/sa fit, not just the (previously unreported) residual SE. The analytic FIM now carries one residual slot per endpoint; a pure additive endpoint gets a real entry, and any other endpoint’s slot is held at exactly zero and dropped before the matrix is inverted, falling back to the linearized FIM for that endpoint’s residual SE as before.

• calc.COV()‘s covFull residual-variance block (saemix func_FIM.R blocB) is now masked to the endpoint each residual parameter belongs to (#904). For a multi-endpoint SAEM model with separate residuals per endpoint, every residual parameter’s dVi/d(param) previously spanned all endpoints’ observation rows instead of only its own, so the reported residual standard errors were wrong. Single-endpoint models, including combined add()+prop(), were unaffected.

• est="saem"’s E-step (the simulated chain and mixture responsibilities) now honors saemControl(addProp=) instead of always forming the combined-error SD as a + b*|f| (combined1) (#912). The M-step objective already branched on addProp, so a combined2 endpoint (a + b*f combined as sqrt(a^2+b^2*f^2), the default) was simulated under the wrong SD: the chain targeted a different posterior than the one being estimated. Only addProp="combined2" (or model-declared combined2()) fits with both an additive and a proportional/power term move; combined1 fits are bit-for-bit unchanged since that branch’s formula did not change.

• est="saem"’s M-step objective for a plain add()+pow() endpoint (no boxCox()/yeoJohnson()) formed the combined2 residual SD as a^2+b^2*f^(2*pw) and used it directly in place of the SD, missing the sqrt() every sibling combined objective (add()+prop(), and add()+pow()+boxCox()/yeoJohnson()) applies. Found while auditing the addProp branches for the E-step fix above; only a plain add()+pow() endpoint under the default addProp="combined2" was affected.

• "indLin" is no longer excluded from the ODE-method fallback candidates a post-fit table/residual solve tries when the fit’s own ODE method is neither "dop853", "liblsoda", nor "lsoda" (#858). rxode2/#1183-#1185 restored indLin()/matrix-exponential correctness, which was the reason for the exclusion.

Internal

• The SAEM predOnly model (used for residuals, tables and the covariance step) no longer emits a THETA/ETA alias assignment that exactly duplicates one the mu-reference replacement block already emitted. These were trailing dead stores that rxode2 repeated in both dydt and calc_lhs; the emitted model, its solve column layout, and every value it produces are unchanged.

• getBaseSimModelFit() for the focei family (focei, foce, focep, fo, foi, posthoc) no longer does three times the work for the same answer. The method built a predOnly-based simulation model expression and then discarded it, and called getBaseSimModelFit.default() twice – once with the result thrown away – so lowering a focei fit to a simulation model lowered it three times, one of those through a rxNorm() of the focei predOnly model. These methods are now aliases of the default, which is what they already amounted to.

• rxode2::rxSolve() on a fit no longer re-derives the model on every call (nlmixr2/rxode2#1289). Each call used to lower the fit to an rxode2 simulation model and re-run the pre-process hooks to build $simInfo; for an ODE model that was most of the ~0.1 s per call, and it grew process memory by a couple of MB per call that neither gc() nor rxode2::rxUnloadAll() gave back, so simulating from a fit in a loop eventually exhausted memory. The lowered simulation model is now cached (keyed on the fitted model itself, so a piped or refit model gets its own; set options(nlmixr2.simModelCache = FALSE) to disable), and $simInfo is only derived when the simulation actually uses the model’s uncertainty – which a plain rxSolve(fit, events) does not. On the issue’s reprex (one-compartment ODE fit of theo_sd) repeated rxSolve(fit, ev) went from 0.106 s and +2.0 MB per call to 0.008 s and no measurable growth; the solved results are unchanged, seed for seed.

• A covariate whose value is carried on the model in rxode2::rxForcedPars() is no longer required to be a column of the data. Such a covariate is supplied by the model itself, so demanding it from the data rejected a well-specified fit and forced the caller to add a placeholder column. This is the same rule rxode2 applies when resolving solve parameters, and it lets a model own parameters the user never sees – for example neural-network weights.

• ini() on a fit now calls rxode2::.iniHandleLine() rather than the rxode2::.iniHandleFixOrUnfix() alias for it. They are the same function; this was the last caller of the old name anywhere in the ecosystem, so rxode2 can now drop it (nlmixr2/rxode2#1250).

• est="npb"’s Gibbs sampler (support-point/stick-breaking/mixture-proportion draws) and the shared npbSampleMixProbs() mixture Dirichlet step now draw from rxode2’s per-thread threefry engine instead of R’s own RNG (R::rnorm/R::unif_rand/R::rbeta/R::rgamma, seeded via an R-level set.seed() call). A distribution the engine does not cover directly (Beta, Gamma) is drawn by inverse-CDF from a threefry uniform, the same technique already used for est="impmap"’s chi-square proposal scale. This is the convention every other estimation method already follows, and it means npbControl(seed=) reproducibility no longer depends on R’s ambient RNG state; a fit’s exact draws (and so its reported values, given the same seed) change as a result.

nlmixr2est 7.0.2

Breaking changes

• est="vae": naming a covariate in vaeControl(shapes=) now also limits the search to it. The list form gained a fixCov element defaulting to TRUE, so shapes = list(WT = "power") searches WT and nothing else, where previously it searched every covariate with WT restricted to "power". Add fixCov = FALSE to restore the old meaning. Excluded covariates are listed in $runInfo. A character vector (shapes = c("power", "lin")) names no covariate and is unaffected.

• Dropped the qs2 dependency (and with it stringfish, which no longer loads against RcppParallel >= 6.0.0): the focei model disk cache now uses RDS files and compressed fit components use base R serialization (rxode2::rxGetDefaultSerialize(), “bzip2” by default). Old fits holding qs2-serialized components can still be read when the qs2 package is installed; otherwise accessing them warns and returns NULL. Requires rxode2 (>= 5.1.5) for rxDeserialize().

New features

• The mu-referenced FOCEi family is experimental. est = "mfocei", "ifocei", "mfoce", "ifoce", "mfocep", "ifocep", "magq", "iagq", "mlaplace", "ilaplace" and their fast=TRUE siblings ("mfoceif" and relatives) are research methods. They are not validated to the standard of the established estimation methods, their results should not be relied on without independent checking, and their interface and defaults may change or be withdrawn in a future release without a deprecation cycle.

  • The same machinery is reachable from the ordinary methods with foceiControl(muModel=) ("lin" or "irls", default "none"), which is where it will continue to live.

• The default sigdig is now 3 (was 4) for every estimation method except est="nls". sigdig drives the ODE solver tolerances as rtol = 10^-sigdig and atol = 10^(-sigdig-3), so the default solve is now rtol = 1e-3, atol = 1e-6 – what most open-source ODE solvers default to, and still tighter than the precision the optimizer targets. Fits are faster. Pass sigdig = 4 to any control function to restore the previous tolerances.

  • est="nls" keeps sigdig = 4: its Levenberg-Marquardt step is sensitive to solver noise, and it already requests a solve three orders tighter than the optimizer target.

  • The optimizer tolerances that are tuned values rather than the 10^-sigdig formula (est="nlm", est="nlme") stay anchored at sigdig = 4, so at the new default they also sit one order looser.

  • Printed parameter tables now show 3 significant digits rather than 4. sigdigTable follows sigdig when it is not set explicitly, and that coupling is deliberate: a fit converged to about 3 digits should not report 4. Set sigdigTable = 4 to keep the previous output.

• Importance-sampling EM (est="imp" / "impmap" / "qrpem"): the proposal density is adapted per subject rather than by one global setting, and a diagnostic is reported that can tell when it matters.

  • The proposal scale is adapted from the second iteration on. The first iteration normalizes its weights against the starting mode, which is not yet a meaningful reference, so its coverage statistic reads far worse than the truth and would otherwise inflate the proposal for the whole fit.

  • fit$env$impPsisK gives a Pareto k-hat per subject – the tail index of that subject’s importance weights. k > 0.7 means those weights have infinite variance and that subject’s contribution is untrustworthy. This is worth checking because the two statistics already reported cannot detect the problem: xi (NONMEM’s IACCEPT quantity) and the Kish effective sample size are both means over samples drawn from the proposal, so neither sees a tail the proposal rarely visits. On plain theophylline, two of twelve subjects have k-hat of 2.60 and 1.28 while xi reads ~0.97 and the effective-sample fraction ~0.99 for those same subjects.

  • impmapControl(df=) switches the proposal from a multivariate normal to a multivariate t (NONMEM DF). This is the remedy for a bad k-hat, because it changes the proposal’s tails rather than its width, and tail weight is what decides whether the weights are well behaved. More samples does not help – boosting a failing subject tenfold moved its k-hat from 0.76 to 3.28 – whereas df = 20 cleared every failing subject for 0.25% of the effective sample size.

  • impmapControl(isample=) additionally accepts one count per subject.

  • impmapControl(gammaMethod=) selects how the proposal scale is adapted: one shared value, or per subject two-sided toward iaccept on that subject’s own xi (NONMEM’s rule). "auto", the default, uses the per-subject law only for models that are not transformably normal, since gamma = 1 is already efficient when the individual posterior is close to Gaussian.

  • impmapControl(auto=) is NONMEM’s AUTO=1: choose df, isample and iaccept per subject. It defaults to TRUE. It escalates df only for subjects whose k-hat says they need it, leaving the rest on the cheaper Gaussian, and shifts sample budget from data-rich subjects to difficult ones. Measured on theophylline against a high-accuracy reference, it takes the worst k-hat from 2.44 to 0.49 and improves Omega accuracy about 20%, for about 19% more Monte-Carlo noise on the objective; infinite-variance weights are a correctness problem whose error is unbounded in the worst case, while the added noise is bounded and measurable. Set auto = FALSE for the un-adapted behaviour, which is the better choice when fit$env$impPsisK is already comfortably below 0.7 everywhere and the tightest possible objective is wanted.

  Note NONMEM does not publish the values its AUTO=1 uses; only the nobs < neta trigger and IACCEPT ~ 0.2 are documented. The concrete numbers here (df = 30, the k-hat thresholds, the sample-budget rule) are nlmixr2’s own, tuned on the measurements above.

• Importance-sampling EM: the covMethod="imp" covariance is now evaluated at the proposal the fit actually converged on, rather than at the control’s initial gamma with a Gaussian proposal.

• Importance-sampling EM: $runInfo now names which sampling-efficiency statistic a fit is reporting, and states that xi and the Kish effective-sample fraction are not comparable with each other.

• foceiControl() gains outerMaxOdeRecalc, outerOdeRecalcFactor and outerStickyRecalcN, which loosen ODE tolerances and retry the analytic outer (augmented sensitivity) solve for a single subject that fails at the requested tolerance. Previously one subject’s failed augmented solve dropped the whole gradient to finite differences; now that subject can still contribute an analytic gradient, which is generally more accurate than the FD approximation. The loosening is per subject, so it is safe under the parallel outer solve, and it is tracked separately from the inner problem’s maxOdeRecalc/odeRecalcFactor/stickyRecalcN – a fit may loosen one and not the other, and the warning names whichever applied.

• est="vae": vaeControl(shapes=) list elements are now dispatched individually, so the covariate-named and list(var=, covar=, shapes=) forms can be mixed in one list. A named element is exact shorthand for the covariate-wide rule, and a shape value of TRUE means “eligible, default shapes” – which is how a categorical covariate is named, since it takes no parameterization:

  vaeControl(shapes = list(list(var = "cl", covar = "WT", shapes = "power"),
                           SEX = TRUE))

  This is also how a covariate is restricted to particular parameters without writing the effect into the model: a var+covar rule makes only that pair eligible.

• The est="vae" automatic covariate search gained a "hockey" shape, a two-armed piecewise-linear relationship knotted at the covariate’s centering value and written as

  ka <- exp(tka + beta.tka.WT.hockey.low * (WT <  70.5) * (WT - 70.5)
                + beta.tka.WT.hockey.hi  * (WT >= 70.5) * (WT - 70.5)
                + eta.ka)

  It is continuous at the knot, so the structural theta keeps its meaning as the parameter value there. Both arms enter or neither does, and hockey competes with the covariate’s other shapes for the same slot, so a parameter never carries two parameterizations of one covariate. It costs two coefficients against a linear shape’s one, so BICc only takes it when the bend earns its keep. "hockey" is part of the default shapes=; name shapes= without it to opt out. A covariate with fewer than catCutoff of the subjects on one side of the knot is skipped, with a note in $runInfo – reachable only with a covCenter= override, since the median splits the subjects in half.

  A hockey stick you write yourself already worked and is unchanged: each arm is independently a mu2 reference, so vaeControl(pinCovariates=TRUE) (the default) keeps your model text and coefficient names exactly as written.

• L0Learn moved from Suggests to Imports. The covariate search already errored rather than fall back when it needed L0Learn and the package was absent, so it was effectively required; making that explicit removes the failure mode.

• The covariate coefficients est="vae" injects after covariate selection are now named with . separators instead of _: beta.tka.WT.lin rather than beta_tka_WT_lin. This matches the separator the rest of nlmixr2 uses for generated and conventional parameter names (eta.cl, add.sd, prop.sd). A categorical coefficient is built from the covariate and its level directly (beta.tka.SEX.M), so the separator is consistent there too. Coefficients you write yourself are untouched – with vaeControl(pinCovariates=TRUE) (the default) the model keeps your names exactly as written.

• The est="vae" automatic covariate search now explores several parameterizations (“shapes”) of each covariate rather than the single hard-coded log(cov/mean) form. vaeControl(shapes=) takes the same vocabulary as nlmixr2scm::runSCM() – "power" (beta*log(COV/ctr)), "lin" (beta*(COV - ctr)), "log" (beta*log(COV)), "identity" (beta*COV) – plus a new "center" (beta*(COV/ctr)). At most one shape of a covariate may enter a given parameter, as in a stepwise covariate search. shapes= also accepts a list named by covariate, or a list of list(var=, covar=, shapes=) items restricting a single parameter/covariate pair; which covariates are searched is still governed by pinCovariates. Because the selection objective is a least-squares fit with a free intercept, "power" and "log" span the same model, as do "lin", "identity" and "center"; the search chooses between the two families and shapes= chooses how the winner is written back, with the coefficient and the structural parameter adjusted together so the prediction is unchanged.

• est="vae" gains vaeControl(covCenterType=) ("median", the new default, or "mean"), vaeControl(covCenter=) for per-covariate centering values such as c(WT = 70), and vaeControl(catCutoff=).

• The est="vae" covariate search now considers factor and character data columns, which were previously dropped without comment. Each becomes a set of 0/1 indicators against the most frequent level per subject, with levels held by fewer than catCutoff (default 5%) of subjects lumped into that reference. Several levels of one factor may enter a parameter together; only alternate shapes of one covariate are mutually exclusive.

• vaeCovariates() now returns one row per candidate search column, adding raw, shape, level and group columns, and takes the same shapes, covCenterType, covCenter and catCutoff arguments as the fit.

  Together these change the default est="vae" covariate search: more candidate forms are considered and centering moves from the mean to the median, so selected covariates and estimates may differ from 7.0.1. Setting vaeControl(shapes="power", covCenterType="mean") reproduces the previous search.

• vaeControl(covSelectMaxExact=) is now measured in bits of feasible-support space (sum over covariates of log2(1 + shapes tried)) rather than a plain candidate count, so the exact branch-and-bound keeps the same worst-case node budget whether a covariate carries one shape or several. With a single shape per covariate the setting means exactly what it did before. The default stays 17: re-measuring with tools/benchVaeCovSelect.R puts the exact-vs-L0Learn crossover at roughly 16 bits in BOTH regimes (one shape per covariate and two), which is what makes a single threshold in these units meaningful.

• The variational inference method previously called est="advi" is now two methods, est="emvi" (variational EM) and est="fbvi" (full Bayes), sharing a shared control – emviControl() with fbviControl() as its thin wrapper, the way impmapControl()/impControl() already work (was adviControl()). The old name was wrong on both halves: there is no automatic differentiation in the implementation (the gradients come from the FOCEi forward sensitivities), and the default mode was never the published algorithm but a variational-EM hybrid. The two modes were previously selected by pointEstimate=, which is kept but now defaults to whichever the chosen est implies; est wins over a contradicting value and says so. covMethod="advi" is likewise now covMethod="vi". est="advi" never appeared in a released version, so no deprecation shim is provided.

• est="vae" and est="emvi" now estimate the omega off-diagonals of a correlated random-effect block (eta.cl + eta.v ~ c(0.1, 0.01, 0.1)), like saem and the focei family. Both previously kept only the variances and reported the ini correlation unchanged. The estimated block appears in fit$omega and in the updated model’s ini(). Only the declared off-diagonals are estimated – a diagonal model is unchanged, and est="fbvi" (full Bayes) errors on a correlated block rather than silently dropping it.

• est="vae" gains vaeControl(covSelectMethod=) and vaeControl(covSelectMaxExact=), which make covariate selection practical on large candidate sets. The exact branch-and-bound blows up past a few dozen covariates (a single 30-covariate latent dimension takes ~43s, and the M-step runs one per dimension per iteration). With the suggested L0Learn package installed, covSelectMethod="auto" (the default) has L0Learn propose candidate supports for any latent dimension holding at least covSelectMaxExact (default 17, the measured wall-clock crossover) candidates, counted after pinCovariates trimming. Those are candidates only: each is scored with the same exact RSS/omega + penalty*|S| objective, the same OLS and the same tie-break the branch-and-bound uses, then improved by an add/drop/swap local search – so L0Learn’s own objective and scaling cannot shift a selection. Below the threshold the search stays exact and unchanged. When the exact search would be impractical but L0Learn is not installed, the fit errors rather than run it silently; covSelectMaxExact = Inf forces the exact branch-and-bound everywhere. A fit that used the approximate search says so in $runInfo and records it in fit$vae$covSelectMethodUsed.

• est="vae" gains vaeControl(nonMuTheta="grad"), which estimates a structural population theta with no random effect using the exact analytic outer gradient (the machinery behind foceiControl(fast=TRUE)) rather than the bounded bobyqa regression nonMuTheta="regress" uses: one augmented sensitivity solve per M-step replaces the derivative-free sweep. Both modes target the same (full outer) objective, so this is an optimizer change: on theo_sd with a non-mu-referenced tv it reaches a slightly better objective than "regress" and lands within 0.0005 of the FOCEi maximum-likelihood value. It is chosen for that accuracy, not for speed – it runs slower than "regress" (1.47x with one non-mu theta, 1.13x with three). It covers a conditionally Gaussian model and a single non-Gaussian (ll()/generalized) endpoint, which differentiates the log-density directly. A model outside analytic scope (linCmt(), IOV, fo, a multi-endpoint or censored ll() model) reverts to "regress" with a note in $runInfo.

• est="vae" residOptimize="twoStage" now applies to a log-likelihood (ll()) or generalized endpoint. Stage two eligibility was “the parameter has a slot in the error-parameter vector”, and such a model has none, so stage two never ran and "twoStage" silently behaved like the experimental joint "optimize" solve. Eligibility is now decided per parameter – an error parameter (as before), OR a parameter no d/dt() right-hand side, initial condition or dosing modifier can reach – so a theta read only by the log-density is optimized in its own frozen-ODE block as intended. A multi-endpoint model with one Gaussian and one ll() endpoint gets both its error parameter and its log-density-only theta into stage two.

• The est="vae" ELBO now includes the transform-both-sides Jacobian, so a model with lnorm()/boxCox()/yeoJohnson() reports its objective on the DV scale – matching what est="focei" already does – instead of the transformed scale. No effect on a model without a both-sides transform.

• est="vae"’s non-mu theta M-step (both nonMuTheta="regress" and "grad") now optimizes the FULL outer objective – the Laplace determinant, 0.5*log|Omega^-1| and the transform-both-sides Jacobian – rather than the joint likelihood at frozen encoder etas. Every mu-referenced theta is held at its current M-step value, so the two modes now differ only in optimizer (exact analytic gradient vs derivative-free bobyqa) and are directly comparable. On theo_sd with a non-mu tv this moves "regress" from 3.4175 to 3.4324 against a FOCEi maximum-likelihood value of 3.4299.

• Fixed est="vae" diverging when a structural theta with no random effect had no ini() bounds. With infinite bounds nothing constrained the non-mu theta M-step, and a parameter whose likelihood is flat in one direction ran away (an unbounded tv on theo_sd reached ~1e68). An unbounded such theta now falls back to a generous finite window around its ini() estimate, chosen wide enough not to bind at a sane optimum; a user ini() bound still wins. The unbounded model now converges to the same value as the bounded one (tv 3.4324 for nonMuTheta="regress", 3.4294 for "grad", against a FOCEi maximum-likelihood value of 3.4293).

• est="vae" gains vaeControl(residRhoend=), the convergence tolerance of the bounded optimizer that estimates the residual parameters (defaults to rhoend). Worth setting separately because that step runs with the ODE frozen, so tightening it is far cheaper than tightening rhoend, which also tightens the structural regression.

• est="vae" gains an experimental vaeControl(residOptimize="twoStage"), which estimates the residual-error parameters by block coordinate descent: the non-mu-referenced structural thetas first (driven by dv - f), then the residual parameters alone against the extended least-squares objective over the cached (y, f) pairs, needing no ODE re-solve. It is the only path that can estimate an error model with no closed form, and it beats the moment estimator on both additive (131.79 vs 131.81) and combined (121.03 vs 122.47) theo_sd fits. It is now the DEFAULT, so an est="vae" fit with a residual-error parameter changes; residOptimize="moment" restores the previous estimator. It is also the only path that estimates an error model with no closed form. pow() and lnorm() residuals were previously classified “other” and left SILENTLY at their ini() values – on theo_sd, pow(prop.err, pw) returned 0.300/0.800 unchanged (objective 154.4 against 134.8 estimated) and lnorm(add.err) returned 0.500 unchanged (objective 26163 against 849). A transform-both-sides boxCox()/yeoJohnson() lambda was frozen the same way and is now estimated too, bounded to (-2, 2) (boxCox 181.6 -> 43.1, yeoJohnson 131.8 -> 108.4, boxCox 181.6 -> -29.2 on theo_sd). Residual scale parameters are also floored strictly above zero, since the likelihood’s zero-variance floor (r == 0 -> r = 1) would otherwise make a collapsed residual look attractive to the optimizer.

• est="vae" gains vaeControl(sigma0Interp=) for how sigma0 becomes the encoder’s initial posterior spread. "sd" (default) makes the initial posterior SD sigma0, as documented; "reference" makes it sigma0 squared, reproducing the reference implementation (which documents sigma0 as a standard deviation, so its squaring appears unintended).

• est="vae"’s encoder is now conditioned on the covariates, as in Rohleff et al. (2025), which concatenates them to the LSTM’s final hidden state before the head that emits the posterior (torch.cat((hidden[-1], covariates), dim=1)). The covariates were previously not passed to the encoder at all, so the approximate posterior could not express a covariate relationship and the covariate M-step had a weaker signal to read off the posterior means. Fixing it moves the neonatal case study’s covariate estimates close to the reference’s (kin ~ GA 3.51 against its 3.45, previously 2.45) and removes a spurious effect. This changes the results of any est="vae" fit on a model with covariates.

• est="vae" gains vaeControl(gammaSeries=), selecting the decaying step-size series used in the smoothing phase: "reference" (default) 1/(iter - gammaIter), the textbook Kuhn-Lavielle series the reference uses, or "saem" 1/(1 + iter - gammaIter), the continuation form saemControl() uses (its decay starts at 1/2 rather than repeating a gain of 1).

• est="vae" aligns three more details with Rohleff et al. (2025): the smoothing gain is now 1/(iter - gammaIter) (it was 1/(1 + iter - gammaIter), smoothing a step harder than the reference throughout the tail); new vaeControl(omegaUpdate="suffStat") (default) forms the population variances from the EMA sufficient statistics and assigns them instead of blending them a second time at the M-step gain (omega only – the residual error is still smoothed on the SD scale, a documented remaining difference); and new vaeControl(inputScale="reference") (default) computes the encoder-input centering/scaling across the whole padded observation matrix as the reference does, rather than over the observed values only – on a ragged dataset the two differ materially (neonatal SD 1582 vs 506). omegaUpdate="blend" and inputScale="observed" restore the previous behavior.

• est="vae" covariate selection now regresses the SAEM sufficient statistic (an exponential moving average of the posterior means) rather than the current posterior means, matching Rohleff et al. (2025); vaeControl(covSelectSmooth=) restores the previous behavior. The effect is small in practice, since the M-step gain is 1 until gammaIter.

• est="vae" gains vaeControl(mStepObjective=), selecting the objective the M-step for a structural theta with no random effect is optimized against: "outer" (default) uses the full FOCEi outer objective (the frozen-eta joint likelihood plus the Laplace determinant, 0.5*log|Omega^-1| and the transform Jacobian), while "elbo" reproduces the plain variational bound of Rohleff et al. (2025). The default is a deliberate deviation from the reference: the Laplace term is what makes an analytic gradient available for those parameters (the gradient differentiates the marginal likelihood), so under "elbo" nonMuTheta="grad" is downgraded to "regress" with a note in $runInfo. The deviation is confined to that M-step – it does not touch the encoder, the ELBO training step or the covariate-selection criterion – so a model whose structural parameters are all mu-referenced fits identically under either setting.

• est="vae" gains vaeControl(pinCovariates=) (default TRUE) to respect the covariates already written in the model. When the model declares covariate effects, the automatic BICc covariate search is restricted to those covariate/parameter pairs – it may still drop a declared covariate, but never adds one on a parameter the model did not specify – and the original model is updated with the estimates, writing a dropped covariate’s coefficient as 0. A declared covariate that cannot be searched (time-varying, or a raw-linear form on a continuous covariate) is estimated in place by the regress M-step. With pinCovariates=FALSE a model’s declared covariates are estimated in place and the search is turned off; with no declared covariates the full search runs. Each case is noted in $runInfo. (Time-varying covariates are still reported as excluded from the search regardless of the setting.)

• est="vae" now honors mu2/mu3 (algebraic/centered) covariate references, like saem and the mu-focei family, via vaeControl(muRefCovAlg=) (default TRUE). A centered covariate such as wt.cl*(WT/70) or wt.cl*log(WT/70) is evaluated into an internal linear nlmixrMuDerCov# column – the centering is carried by the mu2/mu3 data rather than re-applied by the VAE covariate search – so it can be pinned and selected like any other covariate; the original expression is restored in the reported model.

• The est="vae" covariate search no longer adds its own centering on top of the model’s. A pinned covariate is searched at its MODEL value (the centering the model specifies – typically already applied by mu2/mu3 referencing – is retained), so the structural theta is the model’s intercept directly. A 0/1 indicator covariate (e.g. SEXF) is never centered, since it is already in its natural parameterization; other categorical covariates remain mean-centered and continuous ones remain log(cov/mean).

• The optimization sigdig now sets both the ODE solver tolerances and every estimation method’s optimizer convergence tolerance with one consistent formula, so the optimizer converges to exactly the precision the solve supports. The ODE rtol exponent IS sigdig and atol sits three orders below – rtol = 10^-sigdig, atol = 10^(-sigdig-3) – the same for every solver (stiff, non-stiff, auto-switching); the sensitivity (atolSens/rtolSens) solves match the main solve (the outer gradient and covariance are built from them, so a looser sensitivity tolerance would degrade analytic gradient/covariance accuracy), while steady-state (ssAtol/ssRtol) solves run one order looser. Every optimizer’s convergence tolerance is 10^-sigdig to match (n1qn1 epsilon; bobyqa/newuoa/uobyqa rhoend; nlminb rel.tol/x.tol; lbfgsb3c/optim factr as 10^-sigdig/eps; the FOCEi outer optimizer; saem’s inner residual tol; the standalone nlm and optim). At the default sigdig = 4 this is ODE atol = 1e-7, rtol = 1e-4 and optimizer tolerance 1e-4 (previously a symmetric ODE 5e-7 with optimizer 1e-5). sigdig is routed through all of focei/foce/fo/laplace, saem, emvi/fbvi, vae, nlme, nls, and the nlm family. est="nls" keeps a tighter ODE (three orders below the shared target) because its Levenberg-Marquardt step is sensitive to solver noise. An explicit atol/rtol passed through rxControl still overrides the sigdig-derived value.

• The default sigdig is now 4 for every estimation method. The FOCE family (foce/fo/foi/focep), agq/laplace, impmap, posthoc, and the mu-referenced / IRLS variants previously defaulted to sigdig = 3; with sigdig now driving the ODE tolerances, that inconsistency solved those methods a decimal looser than focei. A single default keeps every method at rtol = 1e-4.

• Added sugar aliases for the optim() methods so est = "neldermead", "bfgs", "cg", "lbfgsb", "sann" and "brent" stand in for est = "optim" with optimControl(method = ...). Any other optimControl() options still apply; the alias only sets the method (and its bounded/unbounded handling, so "brent"/"lbfgsb" honor bounds).

• The inner bounded-bobyqa optimizer that fits the residual-error thetas in est="npag", est="npb" and the est="vae" regress M-step now takes a configurable rhoend (final trust-region radius) via npagControl(rhoend=), npbControl(rhoend=) and vaeControl(rhoend=), threaded to the C++ engine. It defaults to 1e-4, matching the optimizer convergence tolerance 10^(-sigdig) at the default sigdig=4; vaeControl derives it from sigdig when set (npag/npb have no sigdig, so they use the fixed default). (est="saem" already routes its inner tolerance through saemControl(tol=).)

• FOCEi guards each theta’s scaling constant per transform, keeping the derivative-based scaleC where it is well-behaved and falling back only in that transform’s singular / out-of-range region. Each parameter keeps 1/|init| (linear/additive), 1 (log-normal), or its transform-specific formula while the value stays inside a band tailored to that transform (the linear band is foceiControl(scaleCband=), default c(0.1, 10)). Outside the band it falls back to the parameter’s native magnitude |init| (NONMEM7 Appendix K, eq 15.2); for a bounded transform (logit/expit/probit/probitInv), if |init| is also out of range it uses the geometric middle of the band. This fixes the singular cases that froze or destabilized the fit – 1/|init| blowing up for a small covariate initial estimate (and the issue-641 large-additive case, whose special handling this subsumes), log() at init 1, logit at the interval midpoint, factorial/gamma at a digamma zero – while leaving the well-scaled common case, and its results, unchanged.

• The bounded-transform (logit/expit/probit/probitInv) scaleC band is now built from each parameter’s OWN low and high bound instead of a fixed cutoff. The derivative-based scaleC factors as N * M, where N is a per-parameter scale using the distance to each bound ((x-low)(hi-x)/(hi-low) for logit/probit, E/(hi-low) for expit/probitInv) and M is a bounds-invariant factor that carries the singularity. Guarding scaleC to N * [lo, hi] applies the same dimensionless band at every bound, so logit(x, 0, 1) and logit(x, 1, 100) are guarded identically at equal fractional position. Previously a wide interval (e.g. logit(x, 1, 100)) had its healthy large scaleC clipped by the fixed c(1e-4, 10) band and slammed to the midpoint; (0, 1) results are unchanged.

• Fixed FOCEi scaleC for a gamma()-transformed population parameter: rxode2 reports it as curEval="lgammafn", which the scaling setup did not recognize, so it silently received the linear 1/|init| default instead of its 1/digamma scaling.

• The FOCEi family nudges a structural population parameter (theta) initialized at exactly 0 off zero before estimation, controlled by foceiControl(zeroTheta=) (default 0.001), since a zero initial estimate has no native scale to scale by. +zeroTheta is used when within the parameter’s bounds, otherwise -zeroTheta; if neither is within the bounds it errors. Fixed parameters (including those fixed at 0) are left untouched. Residual error parameters are also left untouched: they carry their own scaleC, so an error sd set to exactly 0 still disables that component and a combined error model reduces to the smaller model as before.

• foceiControl() gains shi21hMax and shi21hMin (defaults 2.0 and 1e-4), the upper and lower bounds on the adaptive shi21 finite-difference step used for FOCEi gradients (both the inner eta and, when shi21maxOuter != 0, the outer theta/covariate finite differences). A larger upper bound lets the gradient of a flat, small-magnitude parameter clear the ODE-solver noise floor. The NLM family keeps its own fixed bounds.

• The imp / impmap / qrpem importance-sampling family is faster: the theta-score M-step, the Monte-Carlo covariance (covMethod="imp", the default) and the per-subject proposal build in the E-step are now parallelized over subjects, using the cores set in the control’s rxControl (defaulting to rxode2::getRxThreads()), joining the already-threaded E-step weight loop. All are bit-identical to the single-threaded run at any thread count. This also fixes a latent bug in the theta-sensitivity M-step where d(V)/d(theta) was read from an under-sized per-thread lhs buffer, so a residual-error variance that depends on a structural parameter now contributes the correct M-step gradient.

• est="vae" now runs multi-threaded. The per-subject encoder forward pass and the exact branch-and-bound covariate M-step (previously serial, dominating the EM and covariate-selection phases) are parallelized over the cores set in vaeControl(rxControl=rxode2::rxControl(cores=)) (defaulting to rxode2::getRxThreads()), joining the already-threaded decoder solve. The encoder forward pass and the covariate branch-and-bound are bit-identical to the single-threaded run. The encoder backward (gradient) pass is also parallelized by default (vaeControl(parEncoderBackward=TRUE)); its cross-subject sum cannot be reduced in parallel bit-identically, so it is deterministic for a fixed cores but differs slightly from the serial path. A note is added to the fit’s $runInfo when it is active. For bit-identical, fully reproducible results set options(nlmixr2.identical=TRUE) (flips the default to serial) or vaeControl(parEncoderBackward=FALSE).

• The SAEM Louis stochastic-approximation FIM (covMethod="sa") and the importance-sampling Monte-Carlo observed information (covMethod="imp") are no longer tied to est="saem"/est="imp". They can now be requested as the covMethod of any mixed-effects estimation method (computed post-fit at the converged estimates) and switched onto any completed fit with setCov(fit, "sa") / setCov(fit, "imp"). (The population-only NLM family has no random effects, so sa/imp do not apply there.)

• Several estimation families changed their default covMethod now that any covariance can be applied to any mixed-effects method:

  • the FOCEI family (focei/foce/laplace/agq) now defaults to the "r,s" sandwich (was "analytic");
  • est="vae" now defaults to "r,s" (was "analytic");
  • est="nlme" now keeps nlme’s own covariance ("nlme") by default;
  • the nonparametric family (npag/npb) now defaults to the importance-sampling covariance ("imp"). est="saem" ("sa"), the importance-sampling family ("imp"), the NLM family ("r"/optimizer Hessian), est="emvi"/est="fbvi" ("vi") and fo/foi (no covariance) keep their previous defaults.

• vaeControl(bnbStrategy=) selects the frontier discipline for the exact branch-and-bound covariate selection in est="vae": "lifo" (default, the existing last-in-first-out depth-first search), "fifo" (first-in-first-out) or "lc" (least cost / best-first). The solver is exact, so the selected covariates are identical for every strategy; only the search order differs.

• est="vae" can now estimate structural population parameters that have no random effect (are not mu-referenced). Previously such a theta was frozen at its ini() value because the VAE only estimates parameters in the latent space. vaeControl(nonMuTheta=) selects the treatment: "regress" (default, matching saemControl(nonMuTheta=)) injects no eta and estimates each such theta directly by a bounded bobyqa regression against the FOCEi inner likelihood every M-step (bounds from the ini() lower/upper, blended with the M-step gain), recovering a no-random-effect population parameter without adding a spurious random effect. The eta-injection alternatives estimate it as theta + mean(eta) (the temporary eta is dropped from the output model): "eta" estimates the injected omega and the typical value; "fix" holds both the injected omega AND the typical-value theta fixed at their ini() values (nothing about the parameter is estimated, so it does not appear in the iteration table); "none" keeps the old freeze behavior. A $runInfo note lists which parameters were converted.

• The analytic observed-information covariance is now the preferred covMethod across the mixed-model estimation methods, falling back to each method’s previous default when a model is out of analytic scope:

  • est="saem" keeps the stochastic-approximation FIM ("sa") as the default covMethod, now followed by "analytic" and "linFim". covMethod="analytic" computes the FOCEI analytic covariance at the converged SAEM estimates and falls back to the linearized FIM ("linFim") with a message when out of scope or not positive definite; the "linFim" covariance stays selectable via setCov(fit, "linFim").
  • est="nlme" gains a covMethod argument (c("analytic", "r,s", "r", "s", "nlme", ""), default "analytic") that recomputes the covariance at the converged nlme estimates; "nlme" keeps nlme’s own standard errors (also available via setCov(fit, "nlme")).
  • est="npag"/"npb" (and their m/i variants), which previously reported no covariance, now compute one post-fit at the converged estimates (default "analytic" with the finite-difference fallback chain).
  • est="imp"/"impmap"/"qrpem" gain a covMethod argument (c("imp", "analytic", "r,s", "r", "s", ""), default "imp"). "imp" is the Monte-Carlo importance-sampling covariance that the old impCov=TRUE selected (the impCov argument is removed); the other tokens compute the post-fit FOCEI covariance.
  • est="emvi"/est="fbvi" keep their variational covariance ("vi") as the default but now honors an explicit covMethod (e.g. "analytic") without overwriting it with the variational covariance.
  • setCov()/getVarCov() accept covMethod="analytic" post-fit.

• A general FOCE-family per-subject log-likelihood can now be built from an rxode2 UI model and used outside of a fit, for MCMC/SAMBA-style algorithms (issue #414). foceiLikLoad() compiles the inner model and sets up the problem (including the data) in memory, foceiLikRun() evaluates the individual log-likelihood at a supplied population parameter vector and eta matrix – in parallel per subject – and foceiLikUnload() frees it. The likelihood type may be "focei" (with interaction), "focep" (FOCE+) or "foce" (NONMEM-style), and foceiLikRun(type=) selects the individual joint density log p(y_i, eta_i) ("joint", the default) or the conditional data log-likelihood log p(y_i | eta_i) alone ("cond"). Only one likelihood system may be loaded at a time; loading a second errors until the first is unloaded.

• fit$etaCI returns per-subject confidence intervals for each individual’s eta, complementing the existing fit$etaSE and fit$etaRSE. The intervals are the empirical-Bayes estimate plus/minus a normal quantile times the eta standard error, using the fit’s ci level (default 0.95). Like etaSE, it requires CWRES in the fit (add with addCwres() for non-focei methods).

• est="agq" now supports the analytic outer gradient (agqControl(fast=TRUE)), which was previously available only to the FOCEi family. The AGQ objective is the FOCEi objective with one term swapped – l(etahat) becomes log(sum_k a_k) over the quadrature nodes, while the log det, Omega and tbs terms are unchanged – so its gradient reuses the same sensitivity solve and adds the node terms plus tr(Ht^-1 dHt/dp) for the node placement. As with FOCEi this replaces the finite-difference outer gradient, so it is exact rather than a difference approximation and costs one augmented solve instead of one extra solve per parameter. The quadrature nodes solve a cheaper 1st-order model than the eta-hat point needs (they never read the 2nd-order block), which is where most of the node cost goes once the grid grows. Requires interaction=TRUE; a fit that cannot use it falls back to finite differences rather than failing.

• FOCEi fast=TRUE (and the *f wrappers) now handle general log-likelihood (ll()) and generalized (Poisson, binomial, …) endpoints analytically, where they previously fell back to finite differences. For such an endpoint the per-observation prediction is the log-density, so the inner Hessian is the exact H = Omega^-1 - sum d2(logLik)/deta2 assembled from a second-order sensitivity model at the empirical Bayes estimate. Both the objective’s log|H| and the Almquist outer gradient use it; the gradient’s parameter derivative of H comes from a batched central finite difference of the analytic second-order sensitivities (no third-order tensor). This is markedly faster than the finite-difference outer gradient for models with many subjects. Endpoints outside the analytic gradient’s scope (multiple endpoints, censored observations, linCmt(), IOV, nAGQ>1, or a bounded parameter transform) fall back to the finite-difference gradient, and a model whose second-order expansion is unsupported keeps the finite-difference inner Hessian – all transparently, rather than failing.

• covType="analytic" now covers est="agq" as well (it previously declined for nAGQ > 1 and fell back to the finite-difference covariance). The AGQ observed information is the FOCEi one with the same single term swapped, so the log det half is reused unchanged and only the data half becomes an expectation over the quadrature nodes plus a covariance between their score contributions. At nAGQ=1 it reduces to the FOCEi observed information exactly, and the FOCEi and Laplace results are unchanged. Validated against a finite-difference oracle (tight ODE tolerance, Richardson extrapolation): the AGQ standard errors agree to that oracle’s own noise floor. As with the gradient, a model outside its scope – a general or multi-endpoint residual variance, censoring, IOV, a finite agqLow/agqHi clamp, cholSECov=TRUE, or interaction=FALSE – reports why and keeps the finite-difference covariance.

• Requesting an unsupported est= method (e.g. a typo) now prints the available estimation methods grouped by category (Linearized, Integral approximation, Stochastic EM, Nonparametric, Machine learning, Optimizer (NLM family)) with a short description of each, instead of a single flat list. Calling nlmixr2() with no arguments prints the same grouped list (and invisibly returns it). The new nlmixr2AllEstType() returns the same information as a data frame, and each built-in method carries type and description attributes (e.g. attr(nlmixr2Est.focei, "type")) that third-party methods can set to join the list.

• est="npag"/est="npb" now PIN the current ODE solve during the residual-error (err) parameter optimization instead of re-integrating. Those parameters do not change the prediction f, so each subject’s states are cached at its posterior etas and the ODE is frozen (op_focei.freezeOde) while only r is recomputed – for a mixture the frozen recompute reuses each component’s cached states rather than re-solving them. A structural regressor (which does move f, including an estimated per-component clearance) still re-solves. Results are unchanged.

• SAEM mixture models now fix per-subject membership by default (saemControl(mixProbMethod="regress"), the new default): each subject is hard-classified to its best component once, held fixed, and the soft-EM responsibility step is skipped (reusing the existing responsibility-weighted machinery via a 0/1 mixWeights). This avoids the soft-EM collapse (a component running away to a degenerate value) and is lower-bias on both well-separated and overlapping component evaluations; on heavily overlapping components it can be higher-variance, so the previous soft-EM behavior remains available with mixProbMethod="regularized". Because membership is fixed, the S-step solves each subject once under its own component (a per-subject mixest regressor) instead of running one MCMC chain per component – roughly an nMix-fold reduction in ODE solves per iteration. Split-ETA mixtures (a separate eta per component, which start symmetric and must differentiate during the fit) automatically fall back to soft-EM (regularized).

• SAEM warm-starts its residual-error parameters from the observed per-endpoint moments at the initial predictions (additive SD from sqrt(mean(err^2)), proportional SD from sqrt(mean((err/f)^2))), the same moment estimate est="npag"/est="npb" use – saemControl(residWarmStart=TRUE), the default. Because SAEM forms this at the unconverged population prediction, the proportional moment excludes near-zero predictions (where between-subject variability dominates) and the warm-started value is clamped to a sane multiple of the ini value. Set residWarmStart=FALSE to start from the ini residual values. For mixture models the warm-start is disabled (the poor population initial fit would inflate the residual and stop the components from separating).

• The proportional residual moment used to warm-start est="npag"/est="npb" (and now SAEM) guards against a near-zero prediction: the ratio is err / (abs(f) <= 1e-6 ? 1 : f), so an f at (or near) zero no longer blows up the proportional moment.

• SAEM now estimates population theta parameters that have no associated random effect (the SAEM phi0 fixed effects) by a bounded direct optimization of the observation likelihood each iteration – saemControl(nonMuTheta="regress"), now the DEFAULT – keeping them as plain directly-estimated regressors instead of stochastic phi0 draws with a shrinking variance. The optimization uses robust coordinate descent within a local trust region, honoring each theta’s ini-block bounds, and holds phi0 fixed once the optimizer owns it. On a simulated model with three no-eta thetas (ka, V, a Hill power) this recovered them far more accurately than the old handling (e.g. the absorption theta RMSE dropped ~16x), at some extra runtime (the objective re-solves the ODE). The previous behavior is available with saemControl(nonMuTheta="eta"). For mixture models this falls back to the stochastic phi0 block (the direct optimizer cannot partition a per-component structural theta by subject membership).

• est="npag"/est="npb" now ESTIMATE a mixture (mix()) model’s component structural parameters (e.g. a per-subpopulation clearance) instead of holding them at their initial values. The residual/regressor step optimizes them against the exact mixture negative log-likelihood -sum_i log(sum_m a_m exp(cll_m)) (NONMEM7 eq 1.182), marginalizing over the components with the current proportions a_m (which the proportion update step moves); each per-component conditional log-likelihood carries the -0.5*log(r) penalty, so the additive residual does not collapse. Verified: a two-subpopulation clearance model recovers both component clearances and the mixing proportion, with a non-zero additive SD.

• The per-endpoint residual moment warm start now attributes each observation to its endpoint via a new rxode2 accessor (getIndCmt, reading the CMT time-varying covariate), so a multi-endpoint model warm-starts each endpoint’s residual from its own moment. Requires the matching rxode2 (function-pointer table index 82).

• est="npag"/est="npb" now estimate the residual-error parameters with EXTENDED LEAST SQUARES at the individual predictions instead of the marginal likelihood. The marginal likelihood over a flexible nonparametric support rewards a vanishing residual (each support point can then fit its subjects arbitrarily well), so the residual could drift toward zero. The residual step now minimizes the exact conditional normal negative log-likelihood sum_obs(0.5*(f-dv)^2/r + 0.5*log(r) + 0.5*log(2*pi)) at the posterior-mean etas (equivalently extended least squares – same minimizer) – the 0.5*log(r) term penalizes r -> 0, giving the saem/focei residual (e.g. theophylline add.sd ~ 0.73, prop.sd ~ 0.15) rather than a collapsed one. Each variance-scale parameter is warm-started (and, for a single scale per endpoint, set) from the saem-style per-endpoint moment: an additive SD from sqrt(mean(err^2)), a proportional SD from sqrt(mean((err/f)^2)), both on the transform-both-sides scale (so lognormal / box-cox are handled on the transformed residual). A non-mu structural “regressor” is optimized in the same step, with the posterior-mean etas re-derived per candidate so the eta grid cannot stale-absorb the structural shift (this identifies it, e.g. recovering theophylline’s clearance from a deliberately-wrong start). After the residual + regressor thetas converge, a final adaptive-grid pass re-optimizes the support with those thetas held CONSTANT, so the support is the nonparametric MLE of the mixing distribution for the fitted residual and the D(F) global-optimality certificate is restored (~0). npb runs the same residual/regressor step inside its sampler. (A mix() model’s structural component parameters are held at their initial values – the ELS step is not mixture-aware; the components are handled by the mixture marginalization and proportion update.)

• est="npag" now picks the initial grid size automatically from the model’s dimensionality when npagControl(points=) is not supplied: max(2028, 512 * n_eta) (2028 is the Pmetrics NPAG default, which covers a low-dimensional model but grows sparse and can collapse in high dimensions). Theophylline (3 etas) resolves to 2028 (matching Pmetrics); warfarin (8 etas) to 4096. Supply points to override.

• est="npag" is more robust on high-dimensional models (many etas), validated by a golden comparison against Pmetrics NPAG on the Warfarin PK/PD model (transit absorption + Emax turnover, 8 parameters): the per-cycle Psi build is now per-row log-sum-exp normalized on the non-gamma path too, so a hard subject’s conditional density cannot underflow a whole row to zero (which aborted condensation); the Burke interior-point solve ridges the Newton matrix and retries instead of aborting when it is ill-conditioned; and npagControl() exposes gridWidth and gridBounds ("auto"/"ini"/"both") so a bounded, high-dimensional grid can be focused on the plausible region (an unbounded box collapses the support). These are numerically transparent for well-conditioned fits (the normalization restores the exact objective; Burke weights are scale-invariant).

• est="npag"/est="npb" now estimate non-mu structural fixed-effect parameters (a theta with no eta, e.g. ke <- exp(tke), which npag’s grid otherwise does not cover – it covers only mu-referenced and residual/likelihood parameters). By default they are optimized as “regressors” in the residual step: the bounded bobyqa moves them alongside the residual parameters, re-solving the ODE per candidate (they feed the states, so the ODE freeze is turned off for that step). This identifies them sharply – e.g. recovering theophylline’s clearance from a deliberately-wrong start, and a bimodal mixture proportion (p1 = 0.70) that the grid alternative recovered only weakly. The opt-in npagControl(muExpand=TRUE) instead uses the saem-style mu-expansion: inject a pseudo-eta (ke <- exp(tke + eta.tke)), grid-estimate, and recover it as a fixed effect at finalization (support-mean folded into the theta, injected random effect collapsed; the injected eta carries a FIXED omega, excluded from the free omega objective like IOV, so it also works in mixture models). residOptimize="none" holds the structural regressors together with the residual parameters. (A non-mu-referenced ETA – an eta with no paired theta – needs neither: the npag box already covers every eta, so it is a grid dimension estimated as a pure random effect.)

• est="npag" now supports generalized (non-normal) / user-ll() likelihoods. The nonparametric objective sums the inner per-observation llikObs, which for a non-normal endpoint is exactly the user’s log-likelihood, so the objective is already correct; the residual/likelihood parameters (e.g. a Student-t’s degrees of freedom, iniDf$err non-NA) are estimated with the same frozen-ODE bounded step as the residual parameters. Freezing the ODE during that step is valid only when every optimized parameter feeds the post-solve f/r alone (err-tagged) – if a non-err parameter ever enters the optimized set the step re-solves instead. gamma is forced off (a non-normal endpoint has r == 1). A non-mu-referenced structural fixed-effect parameter cannot be placed on the grid and is held at its initial value, reported in the fit’s $runInfo. est="npb" handles non-normal endpoints too (the Gibbs sweep sums the same llikObs).

• est="npb" now runs the residual/regressor optimization (previously it held the residual-error and non-mu structural “regressor” thetas at their initial values and only sampled the mixing distribution). With the sampled mixing distribution held fixed, the same bounded bobyqa step npag uses fits the residual thetas (add/prop/ lnorm/lambda/ar) and any structural regressor – recovering, e.g., theophylline’s clearance from a deliberately-wrong start. npbControl(residOptimize=) selects it: "alternate" (default) re-fits during burn-in and then holds the thetas fixed for the sampling phase (so every collected draw shares the converged residual scale), "final" fits once at the converged draw, "none" holds them at their initial values. Unlike npag, npb does not optimize the assay-error multiplier (gamma) – the residual thetas are fit directly.

• est="npag" and est="npb" now support mixture (sub-population) mix() models. Each subject is split into per-component pseudo-subjects and the conditional likelihood is marginalized over the components using the mixture proportions (p(y_i | phi) = sum_m mixProb_m * p(y_i | phi, component m)). npag updates an estimated proportion each cycle by an EM step (support points and weights held fixed); npb samples the proportions inside the blocked Gibbs sweep – each subject draws a component from its posterior responsibility and the proportions are drawn from Dirichlet(1 + component counts), with the posterior-mean proportions reported in $env$npbMixProb. A fix()ed proportion is held at its ini value in both.

• est="npb" now supports multiple independent chains (npbControl(nchains=)): the stick-breaking Gibbs sampler runs once per chain (seed offset per chain), the posterior-mean draws are pooled, and a Gelman-Rubin R-hat per eta is reported in $env$npbRhat (~1 at convergence; > ~1.1 flags non-convergence).

• est="npb" is faster: the two per-sweep loops that re-solve the ODE serially (the support-location Metropolis-Hastings step, and the mixture-proportion Gibbs step for mix() models) now solve their per-subject conditional likelihoods in parallel over subjects, matching the already-parallel Psi build. The proposal and accept/reject draws stay serial in their original order, so a fixed-seed fit is bit-for-bit identical regardless of thread count.

• est="npag" is faster: it no longer does a redundant full conditional-density build at the first cycle (the degeneracy check now reads the working build’s per-subject maxima), and the one-time D(F) global-optimality scan is smaller by default and configurable via npagControl(dfScan=) (-1 auto, 0 to skip the certificate, or an explicit scan size). Neither change affects the fitted support, Omega, thetas, or objective.

• npagControl(cores=) and npbControl(cores=) set the number of threads used for the parallel per-subject conditional-likelihood solves. The default (NULL) uses the current rxode2 thread count (rxode2::getRxThreads()); an integer sets the thread count for the fit and restores it afterwards.

• est="saem" now fits general log-likelihood (ll() ~ expr) models the saemix way (the model returns the per-observation loglik; the standard MCMC kernels use -ll as the observation loss). The solve event data keeps DV when the model references it (previously dropped, so the likelihood solve errored “parameter(s) required for solving: DV”); the fixed-effect-only (phi0) parameters are optimized with the bounded bobyqa honoring the ini-block bounds (so a likelihood SD stays non-negative). Normal-endpoint saem is unchanged.

• Nonparametric engines (cont.): est="npag" optimizes the residual parameters with the bounded minqa::bobyqa, honoring the ini-block lower/upper bounds of each residual parameter (e.g. an additive SD stays >= 0, an AR correlation in (-1,1)). An unbounded optimizer could wander into an invalid region, so newuoa / nelder-mead are no longer used for the residual step (the residType control is removed).

• SAEM general log-likelihood: the fixed-effect-only (phi0) refinement step (saemix “ind.fix10”, distribution=general) is now optimized with the same derivative-free optimizers as the residual step (nelder-mead / newuoa, selected by type) instead of L-BFGS-B – the model emits no analytic d(ll)/d(phi0), so the previous finite-difference-gradient L-BFGS was pure overhead. phi0 does not enter the ODE, so the states are solved once and held fixed while phi0 is optimized (ODE-freeze), each evaluation recomputing only the log-likelihood. The SAEM-side L-BFGS plumbing (phi0 gradient, trampolines, lbfgs* config) is removed; FOCEI’s outerOpt="lbfgsb" is unaffected.

• Nonparametric engines (cont.): the npag residual-parameter optimization now freezes the ODE states – the inner likelihood solves each (support point, subject) once and re-evaluates only the output f/r for each candidate residual theta, skipping the (costly) re-integration. Results are identical to the full re-solve; on a combined-error theo fit it is ~35% faster, and much more for models with expensive ODEs. Exposed as a general freezeOde option on the inner likelihood (off by default, so all other engines are bit-identical).

• Nonparametric engines (cont.): est="npag" now estimates the residual-error parameters generally. A single variance-scale parameter (pure additive or proportional) is handled by the fast gamma up/down search folded into that theta; anything else – combined additive+proportional (the add/prop ratio a single gamma cannot recover), multiple endpoints (each add.sd/prop.sd), and transform (boxCox/yeoJohnson lambda) or autocorrelation (ar) parameters – is optimized against the nonparametric objective with the support points and weights held fixed, using the same optimizers as SAEM (residType: "newuoa" default, or "nelder-mead"), with gamma as a warm start. The residOptimize control selects "alternate" (default, every cycle), "final" (once at the converged support), or "none" (hold at ini). On a simulated two-endpoint model npag recovers add.sd1=0.20 and add.sd2=1.47 (truth 0.20 / 1.50), matching FOCEI, where a single global gamma had forced them equal; on simulated AR(1) data (true ar1.cor=0.6) it recovers ~0.54 from a 0 start where gradient FOCEI stalls at the ar1.cor=0 saddle. The reported residual reflects the estimate. Note: because the support distribution is flexible it can absorb additive residual scatter, so the additive term of a combined error model may be smaller than a parametric fit (documented in ?npagControl).

• Nonparametric engines (cont.): the npag/npb conditional likelihood now folds in the transform-both-sides (dTBS) per-observation Jacobian, so lnorm, boxCox, and yeoJohnson residual models are handled correctly and lambda-type transform parameters are estimable. Proportional and combined additive + proportional error are supported, and the global-optimality certificate D(F) is now evaluated at the fitted gamma (so it reaches ~0 for proportional/combined models with gamma optimization on). A model whose transform link sees a non-positive prediction (e.g. lnorm at an observation where the prediction is

  0) now raises a clear error instead of an Armadillo empty-matrix crash.

• Nonparametric engines (cont.): the npag/npb engines now support fixed parameters. Fixed population thetas (including fixed residual parameters such as add.sd <- fix(0.7)) are held at their ini value. Fixed-Omega etas – for example a fixed inter-occasion variance iov.ka ~ fix(0.05) | occ – remain support-point dimensions but keep their variance held at the fixed value instead of being estimated, so IOV models fit.

• Nonparametric engines (cont.): est="npag" now reports the global-optimality certificate D(F) ($env$npagDF; ~0 certifies the nonparametric maximum likelihood), records a per-cycle parameter-history trace through the shared scale.h printer ($parHistData), and installs the reported Omega masked by the model’s sparsity so correlated-eta models keep their off-diagonal terms. AR(1) and other transform-both-sides / structured residual models are supported (the residual enters as f + sqrt(r)*eps, so any structure carried in r flows through the conditional likelihood).

• Validation: a bimodal-recovery test confirms est="npag" recovers a two-subpopulation (fast/slow absorption) parameter distribution – both modes carry substantial weight and the recovered cluster means land near the simulated truth – the defining nonparametric capability a single-mode parametric random-effect model cannot reproduce.

• est="npb" (nonparametric Bayes) is now a usable engine: a truncated stick-breaking Dirichlet-process mixture sampled by a blocked Metropolis-within-Gibbs sampler (cluster assignments, stick weights, MH support locations). It reuses the same conditional-likelihood primitive as npag and returns a nlmixr2FitData with the posterior mixing distribution ($env$npbSupport/npbWeights), per-subject posterior-mean etas, and posterior draws of the population mean (npbMeanDraws) for Bayesian credible intervals. npbControl() exposes points (truncation K), alpha, burnin, nsamp, propSd, and seed. (Gelman-Rubin multi-chain convergence is a follow-up.)

• est="npag" is now a usable engine: it returns a standard nlmixr2FitData object with the nonparametric population summary (mean + variance mapped to the reported theta/Omega), per-subject posterior-mean etas, and the discrete support-point distribution attached to the fit ($env$npagSupport, npagWeights, npagPosteriorEta, npagGamma, npagNspp). npagControl() exposes points, cycles, and gammaOptimize. (The reported Omega uses the support-point variances; correlated-Omega models and the global-optimality certificate are follow-ups.)

• Nonparametric engines (cont.): added the residual-error magnitude (gamma) optimization inside the NPAG cycle (per-cycle up/down search). Gamma scales the residual variance inside the FOCEi inner likelihood, so censoring (BLQ/ALQ via the M3 censored likelihood – the normal tail probability below/above the limit) and transform-both-sides are handled correctly at the scaled error. The objective uses a log-sum-exp row normalization for numerical stability. Generalized (non-normal) likelihoods are not supported and are rejected with an error. Note: the npag/npb objective is the nonparametric marginal log-likelihood and is NOT comparable to NONMEM/FOCEI -2LL.

• Nonparametric engines (cont.): assembled the NPAG adaptive-grid cycle (Yamada Alg 1) – Sobol grid, Psi, Burke IPM, weight/QR condensation, adaptive-grid expansion (npExpandGrid), and the eps/F convergence controller. Runs end-to-end on Theophylline (exposed as npagCycle_ ahead of the full fit-object wiring).

• Nonparametric engines (cont.): added the Sobol initial grid (npSobolGrid), weight-threshold and QR rank-revealing condensation (npCondenseWeights / npCondenseQR), and the eta-space support-point box (.npEtaBox, control-selectable via gridBounds/gridWidth).

• Nonparametric engines (cont.): added the conditional-likelihood primitive (npEvalCondLik) and the parallel Psi-matrix builder (npBuildPsi), reusing the FOCEi inner solve so residual-error models, transform-both-sides and censoring carry over unchanged.

• Scaffolding for two native nonparametric estimation engines, est="npag" (nonparametric adaptive grid) and est="npb" (nonparametric Bayes), plus their mu-referenced sugar variants mnpag/inpag and mnpb/inpb (OLS and IRLS covariate M-step). Both reuse the FOCEI inner likelihood machinery; the estimation loop runs in C++. The algorithm itself is added in subsequent releases (the drivers currently report that estimation is not yet implemented).

• Fix the covariance matrix ($cov) of a bounded-parameter fit run with an unbounded method (e.g. saem): the internal rxBoundedTr.<name> name leaked into $cov and the back-transform Jacobian was not applied to it, so the reported standard errors were on the internal (transformed) scale. $cov (and the stashed full theta+Omega covariance) are now renamed to the original parameter names and Jacobian-corrected; Omega and residual terms are untransformed so they pass through unchanged.

• The nlm parameter-history machinery can now be driven by an external optimizer. nlmerSolveGrad() gains a record argument that logs the evaluation’s population parameter estimate (the per-subject mean of the phi columns) into the resident scale, and nlmGetParHist() is now exported so an externally-optimized engine (e.g. babelmixr2’s nlmer, driven by lme4::nlmer) can recover the accumulated parameter history before .nlmFreeEnv(). A new optional showOfv field in the nlm solve control hides the objective column for these engines (they record parameters only).

New estimation methods

• est = "emvi" and est = "fbvi" (emviControl() / fbviControl()): variational inference in the style of Kucukelbir et al. (2017), mean-field or block full-rank family. emvi is variational EM (variational posterior over the etas, population parameters point-estimated by an M-step); fbvi adds the population vector to the variational posterior under flat priors. Neither is the published ADVI algorithm and neither is named for it: the gradient comes from the FOCEi forward sensitivities rather than automatic differentiation, and even fbvi carries omega as per-eta log-variances rather than freely. The whole optimization runs in one C++ call, reproducibly and independent of the thread count.

• est = "impmap" and est = "imp" (impmapControl() / impControl()): importance-sampling EM in the style of NONMEM METHOD=IMP, with the E-step proposal at each subject’s MAP mode (impmap) or running conditional mean (imp). Supports mu-referenced, mixture, bounded and fix()ed models; the reported objective is a FOCEi evaluation at the EM estimate. Quasi-random (Sobol) importance sampling (qr=, Leary & Dunlavey 2012) and SIR M-step acceleration (sir=) are available and stay thread-count independent.

• est = "qrpem" (qrpemControl()): sugar for the impmap EM with qr = TRUE and sir = TRUE.

• Mu-referenced FOCEI family: mfocei/ifocei, mfoce/ifoce, mfocep/ifocep, magq/iagq, mlaplace/ilaplace (with matching *Control() functions). Mu-referenced population and covariate-coefficient thetas are profiled out of the outer optimizer by an in-C++ OLS (m*) or IRLS (i*) regression; bounded mu parameters are regression-updated with a clamped step. New foceiControl() options muModel, muRefCovAlg, muModelTol, muModelMaxCycles, muModelClampRetries.

• focep/mfocep/ifocep: the foce/mfoce/ifoce methods with foce = "foce+" forced.

• *f convenience methods (focef, foceif, focepf, agqf and the mu/irls variants mfocef/mfoceif/mfocepf/magqf and ifocef/ifoceif/ifocepf/iagqf): the base method with foceiControl(fast = TRUE) as the default.

FOCEI / FOCE

• foceiControl(fast = TRUE): analytic FOCEI/FOCE outer gradient from Almquist

  (2015) sensitivity equations, solved for all subjects in one threaded rxode2 solve; out-of-scope models fall back to finite differences. Covers censored M2/M3/M4, an estimated boxCox/yeoJohnson lambda, matExp()/indLin(), foce+, modeled dosing (f()/lag()/rate()/dur()), and mu-referenced covariate reuse. Under fast the outer optimizer defaults to lbfgsb3c and mceta defaults to -2 (Eq-48 warm-start of the next inner problem).

• covMethod = "analytic" (folding in the old covType): exact analytic observed-information covariance for FOCEI/FOCE matching NONMEM $COV MATRIX=R, covering additive/proportional/combined error, censored M2/M3/M4 (censOption = "gauss"), estimated lambda, foce+, matExp()/indLin(), and mu-referenced/covariate parameters; out-of-scope fits fall back to the finite-difference sandwich.

• covFull = TRUE (now the default) reports the full theta + residual + Omega covariance for both the analytic and finite-difference methods, with Omega rows named by the random effect (om.eta.cl / cov.eta.cl.eta.v). covMethod = "r,s" is a true sandwich solve(Rfull) %*% Sfull %*% solve(Rfull), "s" is solve(Sfull), "r" is solve(Rfull); covFull = FALSE keeps the theta-only shape.

• foceiControl(foce = c("nonmem", "foce+")): "nonmem" (default) freezes the FOCE residual variance at the eta = 0 prediction to match NONMEM; "foce+" keeps the live conditional variance.

• foceiControl(censOption = c("gauss", "laplace")): censored (M2/M3/M4/BLQ) inner-Hessian treatment; "gauss" (default) matches common tools, "laplace" uses the exact censored second derivative.

• foceiControl(warm = c("calc", "save")): "calc" (default) warm-starts each n1qn1 inner problem from the eta Hessian recalculated at the current theta.

• Residual (error-model) parameters are now included in the focei-family covariance (only fixed, IOV and mixture-probability thetas skip).

• Mixture (mix()) support for focei/foce/fo/foi.

SAEM

• saemControl(covFull = TRUE) (default): full theta + residual + Omega covariance from the linearized FIM. New covMethod = "sa" (stochastic-approximation Fisher information, Kuhn & Lavielle 2005). parHistData records off-diagonal Omega block covariances.

• saem fits general log-likelihood endpoints (ll(name) ~ <expr>, e.g. time-to-event); fixed-effect-only parameters are refined by a bounded derivative-free optimization honoring the ini-block bounds. (This entry originally said the refinement used bounded L-BFGS-B and that saemControl() gained lbfgsLmm/lbfgsFactr/lbfgsPgtol/lbfgsMaxIter; that was never true, and those options are removed in 7.0.3.)

matExp() / indLin()

• Matrix-exponential / inductive-linearization models estimate with the focei, nlm and SAEM families, matching the equivalent ODE model; compartments are ordered source-first from the k_<from>_<to> graph so default dosing is placed correctly.

Output and utilities

• The nonparametric eta-space outputs now carry the eta names: for est = "npag" the support-point matrix (fit$env$npagSupport) and posterior eta matrix (fit$env$npagPosteriorEta) get eta column names; for est = "npb" the same two matrices plus the posterior mean draws (fit$env$npbMeanDraws) get eta column names, and the per-eta R-hat vector (fit$env$npbRhat) gets eta row names.

• est = "npb" now prints its per-sweep iteration history through the shared iteration printer (like every other method) and stores it on the fit as parHistData; the sampler’s results are unchanged (bit-identical).

• The importance-sampling (covMethod = "imp") covariance step now shows a progress bar over its finite-difference evaluations, like the focei covariance step (shown when iteration printing is on).

• New vaeCovariates() returns the covariates est = "vae" would explore.

• New formatMinWidth() for shorter $parFixed display; $parFixed is rebuilt with data.frame operations (#346, #516).

• All estimators share one iteration printer (iterPrintControl(), src/scale.h) with a common row layout; analytic gradients are tracked as their own parHist type and the fit header reports the gradient and mu-model used, e.g. (outer: lbfgsb3c; grad: analytic; mu: irls).

• est = "vae" training runs entirely in C++ (vaeTrainCpp_) and reparameterizes the inner problem in place, so fits are substantially faster.

Changed defaults

• outerOpt = "nlminb" for the finite-difference methods ("lbfgsb3c" when fast = TRUE), sigdig = 4, mceta = -2, censOption = "gauss", foce = "nonmem", covMethod = "analytic".

Bug fixes

Estimation

• covMethod="analytic" now works for models with an estimated boxCox() or yeoJohnson() lambda, which previously always fell back to the finite-difference covariance. The augmented model emits a residual-variance sensitivity for every sigma parameter including lambda, while the shared gradient/covariance model drops only the non-lambda sigma directions; the extra column widened the per-subject sensitivities past the covariance buffers and the assembly errored. Only the dropped directions are restored now.

• The analytic covariance says why it declined. Errors raised while it is assembled were caught and reported as the generic “not available for this model”, which is indistinguishable from a genuine out-of-scope model; they are now reported as analytic err<n>: <message> in $runInfo, where <n> identifies the entry point. A dozen internal bail-outs that returned silently now name their reason too.

• The FOCEi-family objective function is now reproducible, and no longer depends on how the ETAs were reached. The inner problem uses finite-difference steps (etahf/etahr for the ETA gradient, etahh for the FD Hessian) that are searched once per subject and then reused, so whichever call came first fixed them – during optimization that is the warm-start Hessian (foceiControl(warm="calc")) or an early inner iterate, at an ETA that is not the one being reported. All three are now re-searched at the reported ETAs before the final objective is computed. Two consequences:

  • Repeating a fit now gives the same objective function value, and the same value regardless of the number of threads. It previously varied between runs of the same model on the same data, and differed between a threaded and a single-threaded run.

  • Objective function values change, most visibly for models with a non-normal endpoint (ll(), dnorm(), t(), cauchy(), count or ordinal), which difference the whole inner Hessian. A fit evaluated at supplied ETAs (etaMat=, maxInnerIterations=0) and the same fit optimized to those ETAs now agree exactly, where before they could differ by more than 100 objective units on an 8-ETA model.

• The mu-referenced methods (est="mfocei", "ifocei", "mfoce", "ifoce", "mfocep", "ifocep", "magq", "iagq", "mlaplace", "ilaplace") no longer discard a control belonging to another method in the FOCEi family. Each *Control() replaces its class rather than appending, so a foceControl(), focepControl(), agqControl() or laplaceControl() was treated as invalid and silently replaced with defaults – sigdig, covMethod, fast, the tolerances and the iteration caps were all dropped, reported only as a note in the fit output. Such a control is now converted and the settings are kept.

  • The conversion keeps the METHOD’s identity. A setting is carried over only when it differs from the defaults of the control it came from, so a foceControl() cannot quietly run est="mfocei" as FOCE, nor est="magq" as FOCE in place of the quadrature – while a deliberate agqControl(nAGQ=5) or foceiControl(interaction=FALSE) is still honored.

• The ETA-drift theta reset (foceiControl(resetThetaP=), resetThetaFinalP=) now defaults to OFF. It re-centered a mu-referenced theta by the mean ETA and restarted the fit, but when the ETAs cannot re-center – every omega fixed, or a model whose misfit the ETAs must absorb – the shift did not stick and the reset repeated until the restart cap errored the fit out (“Maximum number of theta resets (10) exceeded”). Where it did converge it reached a worse optimum than leaving it off. Set resetThetaP= to restore the old behavior.

• Fixed a theta-reset restart reporting the PREVIOUS attempt’s objective function. The restart reuses the fit environment, and the objective was only computed when the environment did not already carry one, so a restarted fit could report an objective (and the OBJF/AIC/BIC/log-likelihood derived from it) belonging to the aborted attempt rather than to its own parameters.

• The nlm family (est="nlm", "nlminb", …), est="nls" and the importance-sampling EM sensitivity model now honor the covariate interpolation declared in the model (nocb(), linear(), midpoint()). Their gradient and prediction models were generated without those lines, so they always used the default locf() interpolation.

• Fitting many models in one R session uses far less memory. Each compiled model retained a source reference back to the session it was built in, and compiled models are kept for the life of the session, so the retained state grew with every model fitted. A compiled model now retains well under a megabyte instead of tens of megabytes.

• foceiControl(fast=TRUE) now solves its augmented outer-gradient model in the shared FOCEi solve pool (single-endpoint models), sized for the augmented model and with that model’s event (“jump”) sensitivities installed for the batch. This makes the analytic gradient exact for modeled dosing (f()/ lag()), which previously crashed or fell back to finite differences on that path; multiple-endpoint models keep the previous rxSolve route.

• est="vae" with nonMuTheta="grad" solved its augmented outer-gradient model through rxode2::rxSolve on every M-step iteration instead of the shared FOCEi solve pool. The pooled and fallback routes are numerically equivalent, so this cost time rather than accuracy.

• The analytic outer gradient could silently degrade to finite differences. vaeOuterSolve_() returned R_NilValue from a List-returning function, which builds an empty list rather than NULL, so every refusal and every failed augmented solve looked to the caller like a successful solve that returned nothing. Affects est="vae" with nonMuTheta="grad" and any caller sharing that path.

• foceiControl(fast=TRUE) now computes the analytic outer gradient entirely in C++ for est="foce"/"focep", est="agq" and general-likelihood (ll()) endpoints, as est="focei" already did. Those three shapes previously returned to R on every gradient evaluation to rebuild the fit’s etas, omega and setup as R objects; besides the cost, that let R run between the augmented solve and the assembly, where it could disturb the shared solve pool.

• foceiControl(fast=TRUE) fell back to finite differences for every model with no d/dt() – a purely algebraic ll()/generalized endpoint such as a Poisson or logistic regression. Such a model has no ODE state sensitivities and needs none (its prediction derivatives are plain symbolic ones), but the augmented sensitivity model refused to build on the empty expansion, and the pooled solve additionally required a non-zero ODE state count. Both are fixed, so these models now get the analytic gradient; measured against central differences of the objective, agreement is within 6e-7 relative.

• foceiControl(fast=TRUE) no longer returns to R for the outer gradient at all. The R implementation it used to fall through to has been removed: it was a second copy of the same mathematics that had to be kept in step by hand, and reaching it rebuilt the fit’s etas, omega and setup as R objects on every gradient evaluation. A model the analytic gradient cannot handle now goes straight to finite differences, as before, just without the intervening attempt. est="vae" with nonMuTheta="grad" evaluates the same C++ gradient.

• The est="nlm" family (nlm, nlminb, bobyqa, nls and relatives) solved its prediction model without compacting the shared solve pool to that model’s own state count. The pool is sized for the larger sensitivity model, so the predictions were read back at the wrong stride whenever the two differ. No current result changes – for the models covered by the tests the two size the same, so no compaction was needed – but the mismatch is removed rather than left latent.

• foceiControl(fast=TRUE) now uses the analytic outer gradient for multiple-endpoint models, which previously took the slower finite-difference route. Enabling this needed a fix: rxode2 normalizes CMT inside each compiled model by subtracting that model’s own sensitivity-compartment count, which is right for a standalone solve but means peers of different sensitivity depth cannot share one translated event table. Pooled, the inner model resolved every observation to no endpoint at all, so its prediction, residual variance and eta sensitivities evaluated to zero – the conditional estimates collapsed toward zero and DV was silently log-transformed. The shared solve pool now re-bases the CMT covariate for whichever model is reading. Single-endpoint models were never affected.

  • General-likelihood models (ll(), and named distributions such as pois() / binom()) with more than one endpoint likewise use the finite-difference gradient, with a message saying so. Single-endpoint models of that kind are unaffected and use the analytic gradient (nlmixr2/nlmixr2est#838).

• The FOCE EBE Newton convergence tolerance is no longer derived from sigdig; it is fixed at 1e-9, the value it shipped with, and foceiControl(foceEbeTol=) overrides it. Deriving it made the analytic FOCE gradient available or not depending on the requested digits.

• FOCEi: the inner eta-reset / eta-nudge machinery could make the objective function depend on the optimizer’s history rather than on theta alone, so the same theta could return values hundreds of objective-function units apart. With a derivative-free outer optimizer (the default bobyqa) this corrupts the interpolation model and the fit stalls, oscillates, and can exit “normally” at a point worse than one it already visited. Fixed by:

  • Making the inner restart cascade monotone: each nudge restart is now a candidate and the best eta found is the one kept. Previously every n1qn1 restart overwrote the previous result, so the last restart won even when it was worse.
  • Repairing the if (!tryAgain) re-check guards in that cascade, which were unreachable (always evaluated inside if (tryAgain)). Once the first nudge fired, every remaining restart ran unconditionally and the eta was then zeroed regardless of the result.
  • Making the standardized-eta reset per component. A single eta in its tail previously zeroed the subject’s entire eta vector, discarding every converged EBE that subject had.
  • Fixing eta1SD, which was computed as 1/sqrt(etaS) where etaS is Welford’s sum of squared deviations rather than the variance. It is now divided by n - 1, and a zero/non-finite variance disables that criterion for the component instead of producing Inf (which made it fire for every nonzero eta).

• The per-subject “did this ODE solve fail” check now scans only the part of the solve buffer that the subject’s solve actually wrote. When a method sizes the shared solve buffer for a larger model and runs the inner solves compacted against it (est="impmap", "imp", "qrpem", "advi", "emvi", "fbvi", est="vae" with nonMuTheta="grad", and foceiControl(fast=TRUE) with a general ll() endpoint), the check read past that point into slots holding stale values left by an earlier, wider solve of the same reused buffer. A stale NaN/Inf there was reported as a failed solve that had not happened, needlessly loosening ODE tolerances and, once the retry budget was spent, latching the loosened tolerance for the rest of the fit. Objective values for those methods may change slightly as a result.

• Fixed est="vae" freezing a declared covariate effect when the covariate reaches its coefficient’s model line only through an intermediate variable (e.g. wt70 <- WT/70; ka <- exp(lka + beta*log(wt70) + eta.ka)). The coefficient was mis-classified as a plain non-mu-referenced structural theta: frozen at its initial value under nonMuTheta="none" and, under nonMuTheta="eta"/"fix", an eta was injected into the mu-referenced expression, erroring the fit (“2+ single population parameters in a single mu-referenced expression”). Covariate-coefficient detection now reads rxode2’s own mu2RefCovariateReplaceDataFrame (the same table .uiModifyForCovs folds into an nlmixrMuDerCov# column), which already recognizes the coefficient through the intermediate, so the declared effect is estimated in every nonMuTheta mode (issue #801).

• Fixed est="vae" with vaeControl(nonMuTheta="grad") silently discarding every update to a residual-error parameter. An error parameter’s live value is the internal a vector, and the theta slot is rebuilt from it on each evaluation, so the gradient M-step’s theta-only write was overwritten before it was read (the "regress" path already wrote both). The residual was left near its starting value – on theo_sd, add.sd converged to 1.70 against 0.80 for "regress", with an objective ~86 units worse – while the structural theta still looked correct. The gradient step now writes the error parameter back to a, and "grad" reaches a slightly better objective than "regress".

• est="vae" with vaeControl(nonMuTheta="grad") now warm-starts a residual parameter from the closed-form moment estimate on its first gradient step, as the "regress" path already did. While the regress optimizer owns the error parameters the closed-form M-step leaves them alone, so a residual held its ini() value for the whole KL warmup and the gradient steps had to reach the optimum from there – a residual started far from it never arrived, and the result got worse the longer klWarmup was (on theo_sd starting add.sd at 3.0: 1.99 at klWarmup=50 and 2.50 at 150, against 0.80 for "regress").

• est="nlme" now honors sigdig for the ODE solver tolerances. A reversed condition made nlmeControl() fall back to atol=rtol=1e-4 whenever sigdig was set (i.e. always, since it defaults to 4) and only pass sigdig through when it was NULL; the tolerances are now derived from sigdig like every other method.

• Fixed the FOCEi scaleC band guard corrupting est="vae" covariate selection. The guard only rescues a genuinely-computed derivative-based scaling constant (> 0) now; an uninitialized scaleC of exactly 0 is left for the usual min/max clamp instead of being overwritten with |init|. The overwrite had broken VAE covariate discovery on theophylline (no covariates selected, betas collapsed to 0).

• est="vae" with covariateSelection=FALSE now estimates the covariate coefficients written into the model – both linear (beta*WT) and transformed (beta*log(WT/70)) effects – rather than holding them at their ini() value. They are fit in place by the regress M-step regardless of nonMuTheta (previously fixed under nonMuTheta="none" and errored under "fix"/"eta"); a coefficient set with ini(... ~ fix()) still stays fixed.

• est="impmap" now estimates the non-mu structural and residual-error thetas of a general (custom ll()) likelihood model. For such an endpoint rx_pred_ is the log-likelihood itself and rx_r_ is 0, so the Gauss-Newton M-step skipped every observation (V<=0) and left those thetas frozen at their initial values; the M-step now uses the analytic d(ll)/d(theta) directly (empirical-Fisher information), so a raw ll() fit recovers the same parameters as the equivalent add() model.

• est="npag"/est="npb" no longer error with unused argument: 'dfScan' when the post-fit importance-sampling covariance is recomputed (the dfScan field leaked into the down-converted foceiControl).

• est="npag"/est="npb" with a transform-both-sides (lnorm/log/box-Cox) endpoint whose model prediction is non-positive at some observation (e.g. a pre-dose observation where the structural prediction is 0) now records a note in the fit’s $runInfo instead of silently fitting the rxode2-floored value with no indication.

• est="vae" with nonMuTheta="regress" now shows the regressed non-mu-referenced thetas in the iteration table and parameter history. The M-step bobyqa regression already estimated them, but they were omitted from the printed parameter walk (only the latent-space thetas, omega, and residual error were shown), so their progress was invisible; they are now appended to each row with the correct back-transform.

• est="vae" covariate selection no longer silently selects nothing at 32 candidate covariates. The best-subset step enumerated all 2^nCov subsets, which is undefined behavior at nCov = 32 (1u << 32 wraps to 1, so only the empty model was ever tried) and intractable well before that. It now uses an exact branch-and-bound over the same L0/BIC objective, returning the identical optimum while scaling to a few dozen covariates. The selection penalty now also follows the reference implementation’s warmup ramp, tunable via vaeControl(covSelectAlpha=) (default 2, ramped to 1 over klWarmup iterations); ramp iterations are labeled CovSel ramp in the iteration table.

• est="vae" no longer errors with replacement has 0 rows on data that has no AMT column (dose-free datasets such as the neonate weight data); such rows are now treated as observations (EVID = 0).

• est="saem" no longer dies with argument is of length zero when building the SAEM model list. Some rxode2 versions omit the ar column from a model’s predDf, and the SAEM autocorrelation helpers indexed that column directly; they now fall back to the iniDf (err == "ar") representation when the column is absent.

• A mu-referenced or method-variant FOCEi fit (ifocei, mfocei, foce, focep, agq, laplace, and the *f fast variants such as ifoceif) that needed to restart – for example after a zero/bad-gradient theta reset – died with focei$control must be a focei control object. These controls are all built by foceiControl() and then reclassed to their own class, so they do not carry "foceiControl" in their class vector, and the restart-path environment check rejected them even though the fit had been set up from a valid control. The check now recognises the whole FOCEi control family.

• Models that combine linCmt() with ODEs (for example a solved PK driving an effect-compartment ODE) now estimate correctly with the FOCEi and nlm families; the linear compartments are solved as ODEs for those methods. Previously the sensitivity compartments those methods add (one per eta for FOCEi, one per theta for nlm) shifted depot/central past the compartment numbers the data was translated against, so the dose silently landed in a sensitivity compartment, every prediction came back 0 and the objective function was meaningless. Since the model is then no longer mixing a solved system with ODEs, these fits now warn (recorded in fit$runInfo) that the analytic linCmt() could not be used. est="saem" was never affected, keeps the analytic linCmt() and does not warn, as do linCmt() models with no other ODE (#286).

• est="saem" no longer estimates a fix()ed theta that has no eta attached to it; such a parameter now stays at its initial estimate, as it already did for the FOCEi family. The direct phi0 optimization (nonMuTheta="regress", and general-likelihood models) takes over phi0 partway through the fit and skips the update that restores fixed values, so a fixed non-mu-referenced theta drifted off its initial estimate. Estimates of non-fixed parameters are unchanged.

• foceiControl(freezeResidGrad=TRUE) (the default) no longer makes a fit die with “maximum number of theta resets (10) exceeded”. The base solve that caches the states/EBEs for the frozen gradient ran without the gradient flag set, so an ETA-drift theta reset raised inside a gradient restarted the whole fit – on every gradient, until the reset limit tripped (#641).

• A model that combines an inter-occasion variability (IOV) term with a zero inter-individual variability eta on another parameter (for example eta.ka ~ 0 alongside iov.cl ~ 0.1 | occ) no longer fails with “initial ‘omega’ matrix inverse is non-positive definite”. With IOV present the omega is a per-condition list, so the zero-eta detector could not read the eta names and left the zero eta in the matrix, making it singular; the zero eta is now detected and removed as usual. Restoring the original model after such a fit also no longer errors for est="saem" (including table=list(cwres=TRUE)), where the IOV eta is re-expressed as per-occasion id-level etas (#627).

• est="saem" no longer collapses subjects that combine two dosing episodes with overlapping clock times separated by an evid=4 reset – for example a crossover where an IV arm and a depot (f(depot)) arm share the same times. SAEM solves each subject in the ODE solver’s internal time-sorted order, which relocated the reset ahead of the first episode’s observations and merged the two episodes into one trajectory; SAEM then reported a nearly constant PRED and a grossly inflated residual (focei/posthoc already handled this correctly). The reset episodes are now offset internally so the solve times increase within a subject, matching rxSolve()/focei; predictions are unchanged because only time-since-reset matters (#455).

• The est="fo"/est="foi" linearization pass returned an intermediate fit object with an empty control, so .updateParFixed() silently fell back to default table settings (ci/sigdigTable) instead of the fit’s control (#517). The FO/FOI fit now carries its control, and an intermediate fit without a method-specific nmObjGetControl surfaces its stored control rather than returning NULL.

• est="nlme" now accepts the common print control alias, so nlmixr2(..., "nlme", list(print=0)) no longer errors with unused argument: 'print'. nlme prints through its own verbose option, so print maps to it (print=0 runs quietly, any positive value is verbose); an explicit verbose is still honored when print is not supplied.

• FOCEi/FOCE models with a trigonometric term whose argument is a compound expression divided by something (for example a sinusoidal enterohepatic-cycle release sin(2 * 3.14 * (time - mtime1) / period)) no longer fail to build with “too few arguments to function ‘sin’”. The fix is in rxode2’s rxFromSE() (which was dropping the whole argument, emitting sin()); a regression test is added here (nlmixr2/nlmixr2est#513).

• FOCEi now estimates a population parameter that is initialized at exactly 0 (e.g. a covariate effect or an additive term) instead of leaving it frozen at its starting value. The default scaling constant is 1/|initPar|, which is Inf when initPar is 0; it clamped to scaleCmax and made the parameter effectively unoptimizable. getScaleC() now falls back to unit scaling when the initial estimate is 0.

• A single-subject / fixed-effect (“N of 1”) model – one whose only random effects are fixed to zero, which are dropped before estimation – now gives an actionable error when a method that requires random effects (fo, foi, saem, nlme) is used, pointing to methods that can fit it (focei, foce, or a population method such as nlminb, bobyqa or nls). The error also keeps the user’s original model name instead of reporting the internal .mod (issue #493).

• A focei model whose predictions do not depend on any random effect (for example y ~ dpois(rate) where rate is a fixed population parameter rather than a model-predicted value) no longer reports the generic “Aborted calculation” message. The underlying cause is raised directly with guidance on linking each endpoint’s distribution parameter to an eta-varying model quantity (#515).

• est="saem"’s “mis-match in nbr endpoints in model & in data” error is now actionable: it reports the number of endpoints in the model versus the data, lists the observation compartments found in the data, and points the user to check that the CMT/DVID values match the number of model endpoints (error terms). This is the common case of a dataset with extra DVID levels that the model has no matching endpoint for (issue #579).

• est="emvi"/est="fbvi" now reject a mixture (mix()) model up front with a clear message (rxode2::assertRxUiNoMix) instead of running a wrong fit that ignored the mixture structure and then failed late in the output tables with a cryptic “the probabilities in a mixture must sum to a number between 0 and 1, they sum to: 0”.

Estimation and convergence

• A FOCEI fit that hits a theta reset and then restarts no longer aborts with Assertion on 'fitEnv$etaObj$ID' failed: Must be of type 'integer', not 'factor'. The restart re-validated the previous attempt’s etaObf, whose ID column is a factor of the original subject IDs; it is now coerced back to an integer so a genuinely non-converging fit reports its real reason instead of this spurious assertion (#470).

• Fixed the est = "agq" quadrature node scaling. The adaptive Gauss-Hermite nodes were placed without the change-of-variable factor, so increasing nAGQ did not converge to the marginal likelihood – it converged to a wrong value (still better than Laplace, so the objective looked reasonable). The nodes are Gauss-Hermite for the e^{-x^2} kernel while the integral has an e^{-z'z/2} kernel, so they belong at sqrt(2) * chol(Ht)^-1 * x with an exp(x'x) untilt. With the fix the objective converges to the exact marginal likelihood as nAGQ grows. Every nAGQ > 1 objective value (and any standard errors derived from it) changes; focei/foce/fo/laplace are unaffected.

• The analytic covariance (covType = "analytic") now falls back to finite differences under cholSECov = TRUE: the covariance step re-factors the eta Hessian with the generalized Cholesky, which for a non-positive-definite Ht differs from the chol() the analytic observed information assumes.

• Fixed the fast = TRUE analytic gradient for models whose residual variance depends on the prediction (prop, add+prop, combined1, pow, add+pow): a determinant chain-rule aliasing injected a spurious term.

• Fixed the fast = TRUE analytic gradient/covariance for a random effect shared across parameters, enabled sensitivity reuse for a covariate on an eta-less parameter, and fixed the gradient never being used live (it read finalize-only state and silently fell back to finite differences).

• Fixed the FOCE (interaction = FALSE) objective and empirical-Bayes estimates: the residual variance is now supplied at the eta = 0 prediction, so ODE and linCmt() FOCE agree and match the NONMEM reference.

• Bounded the Shi (2021) finite-difference step so a curvature-free search can no longer corrupt the shared solver state.

• Fixed muModel = "lin"/"irls" erroring with two or more covariate expressions (#711) and the user-fixed covariate-coefficient regression bias.

• Fixed impmapControl(impSeed = ) being ignored.

• FOCEI now updates additive mu-referenced population parameters with large-magnitude initial estimates (#641).

• FOCEI theta resets now keep every reset population parameter inside its bounds instead of restarting the optimization out of range, and stop with an informative error when a parameter’s bounds are infeasible (#454).

Covariance and standard errors

• setCov(fit, "analytic") no longer silently installs (and mislabels) the "r,s" finite-difference covariance when the analytic covariance cannot be computed for the model; the fit’s covariance is left unchanged instead.

• fit$etaSE columns are now labeled se(<eta>) (matching fit$etaRSE’s rse(<eta>)%); the label was previously applied to a matrix’s names() (a no-op) so the columns came back as bare eta names.

• covMethod = "r"/"s"/"r,s" standard errors were inflated by a constant factor (sqrt(2) for "r", 2 for "s") from using 2*R^-1/4*S^-1; they now match NONMEM $COV (#666).

• A bounded-parameter fit under an unbounded method (e.g. saem) leaked the internal rxBoundedTr.<name> into $cov without the back-transform Jacobian; $cov is now renamed to the original parameters and Jacobian-corrected.

• The analytic FOCE/foce+ covariance no longer falls out of bounds (from dropped eta = 0 solve slots) to the finite-difference Hessian; the general (f,R) covariance reports covMethod = "analytic" (was "r"), and foceiCovAnalytic()/getVarCov() reproduce it instead of falling back.

• Fixed a segfault in the analytic covariance for out-of-scope models (the augmented build freed the fit’s solve before the finite-difference fallback ran), and the sign of the M2 upper-tail term in the censored inner gradient.

• The mu-referenced/irls FOCEI-family fits (mfocei/ifocei/…) now report Condition#(Cov)/Condition#(Cor) in $objDf; the post-fit covariance install skipped them because the fit tables were rendered before the full-model covariance was recomputed.

• Converting a fit to a different covariance (setCov(), getVarCov()) now refreshes Condition#(Cov)/Condition#(Cor) and the eigen diagnostics from the newly installed covariance instead of leaving the previous method’s values in place.

• SAEM covMethod = "fim" adds the mu-block Hessian (was indefinite / NaN SEs), and "fim"/"sa" report off-diagonal Omega and combined residual SEs. Fixed covMethod = "linFim" and the SAEM covariance erroring for a single population/covariate parameter, and cov2cor for a one-nonzero-diagonal Omega.

Crashes and stability

• Fixed a Windows heap-corruption segfault at more than one core (rxode2 saw every worker as thread 0); the inner loops now pass the real thread id.

• Fixed a segfault in est = "vae" (thread count capped at the solve’s core count) and in nlmSetup on the first estimator call of a session.

• Fixed FOCEi aborting with Cube::slice(): index out of bounds when mceta >= 1 and maxInnerIterations == 0, and a heap-buffer overflow / wrong back-transform in SAEM Box-Cox residual models.

• A non-positive-definite Omega is projected to the nearest PD matrix (SAEM mid-run, with a fit$runInfo warning; and the sym-inv-chol setup for a degenerate fit) so residual/table diagnostics still run; NPDE with a degenerate simulated covariance sets the subject’s NPDE to NA instead of aborting.

• Fixed a segfault when a dataset has no observed subject at all (every subject is a placeholder with no EVID==0 row, as in an aggregate-data output eval such as babelmixr2/admixr2). The no-observation-subject drop now keeps the rows when there is no observed subject to fall back to, and foceiSetup_ no longer reads an empty id vector out of bounds. .nlmSetupEnv() also now supplies a default iterPrintControl when an external caller omits it, instead of erroring with Index out of bounds: [index='iterPrintControl'].

Output, tables, and printing

• vpcSimExpand() no longer merges the entire observed dataset into the simulation when a requested extra column is missing: a dropped filter result meant an unknown column (e.g. a misspelled stratify in vpcPlot()) spliced every observed column into the simulation, and valid columns dragged the rest of the observed data along with them (colliding with the simulation’s own, e.g. time.x/time.y). Only the requested columns are merged now, and a column found in neither the simulation nor the data warns and is ignored (#830).

• For models without etas, the BSV(SD) and Shrink(SD)% columns are no longer added to $parFixed and $parFixedDf; they were always blank for these models (#355).

• Model-defined variables (e.g. ka, cl, v, tad, dosenum, and any user-added line such as WT.OUT <- WT) are now included in the output table whether or not cwres is requested. Previously tableControl(cwres=FALSE) dropped these columns while cwres=TRUE (the default) kept them, so the same model produced different output columns depending on the residual request (#497).

• A zero-fixed eta (e.g. bsva ~ 0) is again restored into the fitted model’s ini()/model() blocks when the estimation makes a nested nlmixr2() call (e.g. adding the focei objective or CWRES), so fit |> ini(bsva ~ 0.1) works; the nested call used to wipe the restore info held in a global (#741).

• augPred() now works on a focei fit whose model has a zero-fixed eta that appears in the prediction (e.g. eta.v ~ 0 used in both the ODE and the residual), instead of erroring with parameter(s) are required for solving: eta.v; the simulation model drops the zero eta consistently with saem (#514).

• laplace/agq family fits label their $objDf row Laplace/AGQ<n> (matching $ofvType) instead of FOCEi; previously the default interaction=TRUE made the interaction label win over the quadrature one. The quadrature objective stays the active one after CWRES; setOfv(fit, "focei") (and addCwres()) now evaluate the true focei objective on a quadrature fit instead of re-labeling its quadrature value.

• Restored the Function Val. objective column and the $parFixed shrinkage coloring; periodic headers now repeat only the column labels.

• $parFixed honors a user sigdig/ci for fits with literally-fixed parameters.

• Literally-fixed population parameters now report their back-transformed value (exp/expit/probitInv) in the Back-transformed column instead of the raw log/logit-scale estimate.

• augPred() now keeps the fit’s original subject ids: the returned id factor carries the actual (character/factor) ids from the fit instead of the internal integer re-numbering (#450).

• vpcSim(fit, pred=TRUE) (and hence VPC plots with a pred line) now works for models with IOV. With IOV the fit’s omega is a list of matrices (id plus one per occasion level), which the pred path treated as a single matrix and errored with invalid 'times' argument; the population prediction now zeros every random effect across all omega levels (#629).

• fit$time again attributes model build/compile to setup/configure (and the nlm family times setup/optimize) instead of other.

• Aggregated ODE-solve warnings report the real subject id; parHistData shows mixture-probability parameters on the natural scale and fit$mixList returns all components; iteration printing labels the estimation phase (Burn in/ KL anneal/EM/Smooth for vae, SA/EM for saem).

• fast = TRUE with a linCmt() model downgrades to fast = FALSE with a message instead of silently falling back per gradient call.

• est = "vae" with automatic covariate selection now reports the selected covariate coefficients (beta_<par>_<cov>) in $parFixed/$parFixedDf instead of dropping them when a population parameter is fixed, and the covariate-bearing mu-parameters back-transform (exp) instead of printing on the raw log scale.

• est = "vae" no longer errors with cannot find parameter 'NA' when a structural (mu-referenced) parameter is fixed with fix(); its random effect is kept (variance estimated) with the fixed value carried in the model.

Data handling

• SAEM no longer errors with No data with ID for a dose-only subject; observation-less subjects are dropped before estimation and re-inserted into the output with a population PRED and NA individual columns, like FOCEi (#687).

• FOCEi no longer errors with 'names' attribute [n] must be the same length as the vector [m] when a subject’s records are all removed during data translation (e.g. every TIME is NA). Such a subject vanishes from the processed data entirely rather than losing only its observations, so it is now detected and dropped from the subject index alongside observation-less subjects (#606).

• Fixed nlmControl() listing eventSens/sensMethod twice. The “initial ETAs were nudged” warning fires only when a nudge actually happened, and a non-default mceta on a fully mu-referenced model falls back to the default with a warning. saemControl(covMethod = "") (skip covariance) no longer errors.

Other

• nlmixr2fix() now actually repairs serialized fit components: it previously tested the component name (not the object) for rawness, so the repair loop never ran, and a successful qs2 read was discarded.

• Fixed $parFixed reporting an uninitialized-memory denormal (e.g. 9.4e-323) as a residual-error parameter’s SE/%RSE for SAEM fits (#816). The finalization filled theta SEs positionally from a covariance that does not span the residual thetas, reading past the end of its diagonal; the SE fill now maps by the covariance dimnames. Post-fit covariance installs also refresh the displayed $parFixed (previously only $parFixedDf was updated), so the residual SE, %RSE, and confidence interval now carry sqrt(diag(fit$cov)); a theta with no covariance row gets a blank SE instead of garbage.

• A non-default confidence level (e.g. saemControl(ci=0.8)) is now honored when a covariance install refreshes $parFixed. The refresh read ci from the model rather than the fit’s control, so it fell back to 0.95: the column was labeled Back-transformed(95%CI) over an 80% interval, and any interval it recomputed used the wrong level.

Internal

• Removed an unreachable duplicate missingTable default assignment in nlmixr2Est0() (issue #385); the earlier default already fixes the value, so the second block could never run. No change to fit results.

• Removed the last bare Rf_error call from the C++ sources (issue #632): the Rcpp::compileAttributes() output now emits the parenthesized (Rf_error) form, and the internal rxError macro was switched to (Rf_error) as well, so the package no longer trips Rcpp’s upcoming Rf_error deprecation warning (RcppCore/Rcpp#1247). The C .Call entry-point validators keep their justified Rf_errorcall uses.

• Consolidated data preparation and the nlm-family control/fit functions, and the analytic-covariance augmented model now uses rxode2’s chunked rxOptExpr(); no change to fit results. The test suite runs a single testthat worker on CI/CRAN and parallel elsewhere, with within-solve threads capped to 2 on CRAN.

nlmixr2est 6.1.0

• Added focei, foce, foi, fo mixture support in nlmixr2est

• Fix focei mixture models with llik residual distributions erroring when a model had exactly one mixture probability parameter

• Fix fit$mixList returning only the first mixture component

• parHistData Back-Transformed rows now show mixture probability parameters on the natural probability scale (0, 1) instead of the raw mlogit estimation scale.

• Fix issue 641: FOCEI now updates additive mu-referenced population parameters whose initial estimates are large in magnitude. Previously a missing branch in .foceiOptEnvSetupScaleC() let scaleC fall through to the C++ default of 1/|init|, which mapped unit steps in scaled space to negligible steps in unscaled space and effectively pinned such parameters at their initial value (e.g. tvemax <- -40 with no transform).

• When model estimation fails, all errors raised during the run are now collected and reported together, instead of only the last error. This is supported by a new collectErr argument to the internal .collectWarn() helper, which captures errors alongside warnings and returns them in the error element of its result list. As a result, errors hidden by on.exit({rxode2::rxProgressAbort()}) handlers (such as the “Aborted calculation” message reported in issue 607) no longer mask the underlying cause; both the inner stop message and any follow-up error from on.exit are now reported to the user. parameters on the natural probability scale instead of the raw mlogit scale. parameters on the natural probability scale

• Hardened mixture-model (mix()) estimation: clearer errors for est="nlme" and invalid initial probabilities, warnings for underflowing/collapsing mixture probabilities, and a fix for the SAEM omega-diagonal floor being raised outside mixture fits

• Fix segfault in nlmSetup on the first estimator call of a fresh R session for pooled estimators

• Guard against null pointer arithmetic in inner.cpp

• Use OpenMP threading for S matrix calculation

• Use OpenMP threading while calculating NPDEs

nlmixr2est 6.0.1

• Fix LTO violation as requested by CRAN by adding -DARMA_DONT_USE_OPENMP to PKG_CXXFLAGS in src/Makevars.in

• Require rxode2 5.1.2 which has the fixed M1-san issues observed here.

nlmixr2est 6.0.0

• focei, foce, fo, laplace, and agq have all been successfully made thread safe and parallelized (for a single CPU). The default tolerance relaxation for difficult to solve ODEs has been changed to per individual instead of for the entire population (which is a breaking change, so major release). This should allow more precision for a majority of the subjects in the optimization process.

• Add predict(fit, level="ipred"), predict(fit, level="individual") or predict(fit, level=1) to predict individual fits (with possibly a new dataset).

• Change test files to .rds files

• Drop magrittr %>% in favor of |>.

• Breaking change: Minimum R version increased from 4.0 to 4.1.0. This change is required to support the native pipe operator |>. Users on R < 4.1.0 will need to upgrade R to install this version of nlmixr2est.

• Bug fixes for deparsing nlmixr2 control objects

• nlm and related pooled methods now run in parallel (based on ID)

• Tests are optimized to reduce redundant fits and run in parallel.

• nlm (and related pooled optimizers: bobyqa, newuoa, uobyqa, n1qn1, lbfgsb3c, optim, nlminb) now support the same censoring behavior (M2/M3/M4) as FOCEI and SAEM. The $censInformation field is populated for these fits in the same way as FOCEI/SAEM.

• agqControl() and laplaceControl() now have rxUiDeparse() methods so they can be saved better in packages like nlmixr2save and shinyMixR.

• Added new outerOpt; methods to focei and related methods (agq, laplace, foce, fo, foi): “uobyqa” and “newuoa”.

• saem and other methods now respect bounds by default by internally adding the appropriate transform and then applying the back-transformation just before returning.

  For parameters that are mu-referenced, this breaks mu-referencing. When it breaks mu-referencing there is a warning issued. The best practice is still to have unbounded parameters with mu-referencing.

  If you want to ignore this behavior you may use control=list(boundedTransform=FALSE) or for saem control=saemControl(boundedTransform=FALSE)

• The mu referencing covariate procedure was made less fragile to support mu referencing in conjunction with iov and bounded parameter transformations.

• Add some bench-marking capabilities and small speed fixes for focei/saem

nlmixr2est 5.0.0

• Remove qs and change to qs2. This breaks backward compatibility.

• Default to non-compressed nlmixr2 objects

nlmixr2est 4.1.1

• Request nlmixr2est’s pre-processing hooks for augPred(), vpcSim() and $simInfo, which fixes augPred in cases where etas=0 are used in nlmixr2 (#587)

• Fix scale.h so that scaleType="none" does not also require scaleTo=0

• Request Armadillo 15 with the special flag in the new RcppArmadillo

• Fix focei without etas (and without log-likelihood normal) to run ELS (See #590).

• Change the IOV implementation (#596):

  • Now shows estimates as CV% or sd without shrinkage calculation.
  • Allow different forms of iov estimation, controlled by iovXform.
  • Retains the iov parameter(s) in the output data.frame.
  • With iov, the $omega shows a list of variability by the conditioning variable(s).
  • fit$iov will show the IOV deviations by the conditioning variables(s) with the exception of id
  • IOV models can be used in other estimation methods and inherits the ETA values.

• Added $etaMat method for nlmixr2 fits to give the value that needs to be passed between each estimation method (related to iov #596)

nlmixr2est 4.1.0

• Updated inferring the estimation method from the control object. Requires the control object to have a class of length one and match the estimation method. For example foceiControl() would assume that the estimation method is related to focei.

• Changed Rstudio completion to not evaluate (in case it gets turned on for data.frames) (See #568)

• Turned on data completion for items like $fitMergeInner

• Breaking change: Changed the estimation method posthoc to add tables and calculate the covariance by default. It is now a method with it’s own control, posthocControl(). As previously the default is not to include the interaction term (but you can turn it on with posthocControl(interaction=TRUE)).

• Added foceControl(), foControl() and foiControl() for the foce, fo and foi methods, respectively. They try to convert the related control structures to the correct control structure for the estimation method.

• Added iov support for focei, foce, and saem (#614)

• Added new estimation method agq which uses adaptive Gauss-Hermite Quadrature to fit a nonlinear-mixed effect model. In this method, you can choose the number of quadrature points to estimate the likelihood, with higher numbers giving more accurate likelihoods. The AGQ implementation in nlmixr2est allows you to specify the number of quadrature points via the agqControl() function, and supports both single and multiple subject models. This method is particularly useful for models where accurate likelihood estimation is critical.

• Also added a laplace method which is the same as agq with 1 node (and is numerically the same as focei, foce or log-likelihood focei/laplace, etc), but uses the agq routine.

• Fixed saem mu-reference display by not compressing the internal item saem0.

nlmixr2est 4.0.2

• The loading and unloading of DLLs has been minimized in this version of nlmixr2est. This avoids loading/reloading the same DLLs and causing the CRAN mac m1 ASAN/USBAN false positive issue observed in CRAN.

• Additionally a new function nlmixr2fix(fit) has been added to nlmixr2est. It attempts to make the fit loaded from a different version of nlmixr2 compatible with nlmixr2 4.0. It also prints out the versions of nlmixr2 that were used when creating this fit. With this information you are more likely to find a way to use the fit in your current session (or in an old session). (Issue #562)

nlmixr2est 4.0.1

• Initialize lbfgsb3 error message to an empty string to address valgrind finding (as requested by CRAN).

nlmixr2est 4.0.0

• When using a model to start a new focei model, the ETAs from the last fit are used as the starting point. Now you can use foceiControl(etaMat=NA) to skip this and use eta=0 for all items.

• When using foceiControl(etaMat=fit), this will extract the ETAs from a fit for use in the next optimization.

• When using a foceiControl(etaMat=) option nlmixr2 no longer only evaluates the inner problem with the etaMat value.

• Add mceta option to "focei".

  • mceta=-1 is the default; the eta restarts at the best eta from the last step to start the inner optimization.
  • mceta=0 the eta starts at 0 to start the inner optimization.
  • mceta=1 the eta starts at either 0 or the best eta, which ever gives the lowest objective function to start the inner optimization.
  • mceta=n under the assumption of omega sample n-1 eta values and use the lowest objective function of eta sampled, last best eta and eta=0 to start the inner optimization.

• Fix Rstudio print (issue #536)

• Support rxode2’s new +var() definition in saem

• Support literal fixing of residuals (#524). All methods that support a literal fix of residuals have an option literalFixRes which defaults to TRUE. To get the behavior from older models you can use literalFixRes=FALSE

• More detailed error messages will be reported for models with errors

nlmixr2est 3.0.4

• More robust covariance calculation in focei.

• Allow hook mechanism to handle piped arguments.

• Fix for when output message from optimizing doesn’t print well (#325)

nlmixr2est 3.0.3

• Moved data check for covariates and required data items to a pre-processing step. This fixes #499. Each method that needs to have a covariate check needs to have a property covPresent. For example to apply the covariate data check to the focei method you need attr(nlmixr2Est.focei, "covPresent") <- TRUE.

• Bug fix for non-mu referenced etas when combined with mu referenced covariate values. (See #498)

• Changed option for "saem" to have literalFix=FALSE. This makes mu-referencing work better when fixing a population value.

nlmixr2est 3.0.2

• Fix bug where models where omega boundary warnings caused problems in estimation (#490)

• Created a new api for pre-processing ui, allowing adding arbitrary hooks. As written now, this includes literal fix and zero omega as well as added the new rxode2 ui processing.

• Fixed compilation to only use -I in most systems for maximum compatibility

nlmixr2est 3.0.1

New features

• Now when optimizing only a single parameter with focei-family, will change to use stats::optimize() for the outer problem (#481)

• When estimating with all fixed population parameters, do a posthoc estimation.

• Internally removed assignInMyNamespace() replacing with nlmixr2global, which fixes some edge case bugs where the nlmixr2 environment was not reset properly.

• Treated edge case where all initial parameters are zero and change scaling from scaled to unscaled (#486)

• Added mu4 referencing that will change string expressions to rxode2 numeric values. This allows derived strings to also be treated as mu expressions (#484)

Bug Fixes

• Fix focei covariance step when many omega values are fixed #482

nlmixr2est 3.0.0

• No binary linking to rxode2, lbfgsb3c and n1q1, which means that updating these will not make nlmixr2est crash without recompiling.

• New mu3 referencing will take context from the model to see if the algebraic expression can be completed from defined model variables; These variable would have to be unique.

nlmixr2est 2.2.2

Breaking changes

• Saem non-mu reference input parameters/covariates were fixed so they work correctly with fixed parameters (Issue #445)

• Focei changed back to having a lower bound for standard deviations when not specified. This means that best model fits may change. You can revert to the old settings by using foceiControl(sdLowerFact=0.0). You can also change the factors to other values than the default value, that is foceiControl(sdLowerFact=0.000001) for instance which would multiply the initial value by 0.000001 when either the lower bound isn’t specified or the lower bound is specified as zero for the error estimates related to error-based standard deviations.

• In nlmixr2, expressions are optimized. Because of that optimization, numerical rounding differences can cause different directions in optimization when fixing parameters in the model vs. fixing the parameters manually.

  This means that the fixed parameters in a model vs hard-coded fixed parameters could give different values in the final model.

  A new option literalFix was introduced which change the fixed population parameters to constants in the model while running the optimization. This makes the output of fixing within the model and fixing manually the same (which is what is likely expected). The default is for this to be turned on (ie. literalFix=TRUE). You can get back the old behavior by using the option literalFix=FALSE.

• In saem, the monte-carlo sampling occurs for all parameters including non-informative ETAs. A fix ensure that non-informative etas in saem are fixed to zero while sampling the phi values. This may change results for models with uninformative etas. To ignore the uninformative etas with saem you ca use use the prior saem handling with saemControl(handleUninformativeEtas=FALSE).

New features

• Gracefully degrade when $cov is not in the right form (see #423)

• Add support for PopED in place solving (used in babelmixr2)

• If est=foceiControl() or other nlmixr2 control with the class foceiControl infer the estimation method is focei

• Add back the warnings when estimation methods ignore the boundaries

• When using rxSolve, now respects the values from tableControl() (#465 and #297)

Bug fixes

• Will emit warnings when the return object is not a nlmixr2 fit (#453)

Other things

• Moved actual code of some matrix libraries to lotri and import them via function pointers

nlmixr2est 2.2.1

• Align with the possibility that linCmt sensitivities may not be present (like intel c++)

Bug fix

• focei cache needs to be based on the parameter order as well as the model information (#415)

nlmixr2est 2.2.0

New Features

• Algebraic mu referencing has been implemented in nlme and saem.

• New estimation method “nlm” has been added to estimate population only likelihoods using stats::nlm and possibly return a standardized nlmixr2 fit.

• New estimation method “nls” has been added to estimate population only problems. This uses minpack.lm::nlsNM by default if present, or the stats::nls

• New estimation method “optim” has been added to estimate population only likelihoods. This uses stats::optim and returns a standardized nlmixr2 fit.

• New estimation method “nlminb” has been added to estimate population only likelihoods. This uses stats::nlminb and returns a standardized nlmixr2 fit.

• New estimation methods from the minqa package: “bobyqa”, “uobyqa” and “newuoa” have been added to estimate population only likelihoods. These methods returns a standardized nlmixr2 fit.

• New estimation method “lbfgsb3c” to estimate population only likelihoods. This returns a standardized nlmixr2 fit.

• New estimation method “n1qn1” to estimate population only likelihoods. This returns a standardized nlmixr2 fit.

• Added new feature for vpcSim() where a minimum number of subjects are simulated from the model when trying to fill in ODEs that were not solved successfully. By default this is 10. This also works-around a bug when there is only one subject simulated and the data.frame has a slightly different output.

Breaking changes

• Removed fit$saemTransformedData since it isn’t actually used in saem anymore (but will break anyone’s code who is using it)

• Now the internal function .foceiPreProcessData() requires the rxode2 control rxControl() because some of the new steady state lag features need to translate the data differently based on rxControl() options.

Bug fixes

• Printing models with correlated omega values and omega values fixed to zero no longer fails (#359)

• Add back values for $parHistData (#368)

• This requires a new rxode2 which will fix multiple endpoint issues observed (#394)

• Manual back-transformed values in $parFixed are now displaying correctly and are calculated based on the confidence interval in the control instead of 95% confidence no matter what (#397)

Other changes

• An as.rxUi() method was added for fit models (#377)

nlmixr2est 2.1.8

• Version bump and a minor documentation update (same as nlmixr2est 2.1.7). This version bump is to simply allow correct binary linkage to rxode2 2.0.14. Otherwise nlmixr2 models will crash R.

nlmixr2est 2.1.7

• As requested by CRAN, remove Rvmmin

• Values in $parFixed for BSV without exponential transformation are now correctly shown (#366)

nlmixr2est 2.1.6

Breaking changes

• Since rxode2 now allows simulation with omega having diagonal zero elements, $omega and $omegaR now reflects this information including the zero omega elements in the output. On the other hand, the other eta-information and standard error information for zero etas are still excluded in $phiR, $phiSE, $eta etc.

Bug fixes

• vpcSim() works when an eta value is fixed to 0 (#341)

• augPred() now consistently uses the simulation model (instead of the inner model used for CWRES calculation).

Other changes

• Dropped dependence on orphaned package ucminf

nlmixr2est 2.1.5

• Add $fitMergeFull, $fitMergInner, $fitMergeLeft, $fitMergeRight as a complement to $dataMergeFull, $dataMergInner, $dataMergeLeft, $dataMergeRight. The fit variants prefer columns in the fit dataset instead of the original dataset. This is useful for goodness of fit plots with censoring since the DV in the fit simulates values under the ipred/residual assumption and will give more appropriate goodness of fits, otherwise these values are the limit of whatever censoring is applied

• Moved the mu reference fix for the split mu referenced model here (from babelmixr2)

nlmixr2est 2.1.4

• Breaking change, now calculate condition number based on covariance and correlation, the names have changed to be more explicit. conditionNumber changed to conditionNumberCov and a new metric conditionNumberCor has been added.

• A bug in boundary value detection prevented automatic covariance calculation with FOCEi estimation (#318)

• Fix vpcSim so that it will be a bit more robust when it is difficult to simulate.

• A bug in model piping which did not allow models to be appended to was fixed (rxode2#364)

• An internal change was made in nlmixr2.rxUi() to better support the babelmixr2 PKNCA estimation method (babelmixr2#75)

• Fixed bug where $iniUi did not return the initial ui when running non focei related methods. Also added alias of $uiIni to the same function.

• Dropped Stan headers for this package, also updated to C++17

nlmixr2est 2.1.3

• Allows $etaH and related family to be integrated into a saem fit if cwres is calculated.

• Fixed a bug where nlmixrLlikObs in the merged dataset is sometimes named llikObs, now it is always named nlmixrLlikObs

• Fixed a bug where nlmixrLlikObs shows up in merged dataset when cwres is not calculated (it was always 0), also allow cwres calculation to pick up nlmixrLlikObs in merged dataset.

• Dropped dparser dependency

nlmixr2est 2.1.2

• Fixes $etaH memory corruption so the standard errors of etas are now correct

• Removed the memory requirements for focei by neta*neta*nsub

• Fixed character based covariates so the work correctly (again) with focei. Added a test for this as well.

nlmixr2est 2.1.1

• Fixes $dataMergeInner so that observation-based log-likelihoods work with infusions. Should fix tests with ggPMX

• Fixes $etaSE and $etaRSE to work correctly when there is only 1 eta.

• Fixes npde valgrind observed on CRAN machines

nlmixr2est 2.1.0

Breaking changes

FOCEi

• Gill forward differences will not repeat now (by default), You can change back to prior behavior with foceiControl(repeatGillMax=3)

• Number of sticky recalculation is reduced to 4; to have the old behavior use foceiControl(stickyRecalcN=5)

• n2ll has been changed to ll to specify individual log-likelihoods. This was only used in simulation and was not well documented.

• Generalized log-likelihood is only supported with rxode2 2.0.8 or later.

FOCEi covariance calculation

• The S matrix calculation was made a bit more robust to errors in individual gradients. When there are errors in the individual gradient calculation, assume the gradient is the same as the overall gradient. In the tests cases, were reasonable using this adjusted S matrix. This means if some individuals do not have very much data to support a specific parameter, a S matrix calculation for the population will still be generated. When there is some patients/subject combinations that do not have sufficient data, we will add the following to the run information: S matrix had problems solving for some subject and parameters. The S matrix calculation will still fail if the percentage of parameters that are being reset is lower than foceiControl(smatPer=0.6) or whatever you specify.

• The r,s covariance matrix will now also check for unreasonably small values (controlled by foceiControl(covSmall=...)) and select a different covariance estimate method even when the “r” and “s” matrices are calculated “correctly”.

New features

• What type(s) censoring (if any) is now stored in fit$censInformation

• Standard errors of $etas can now be obtained with fit$phiSE, also available are fit$phiRSE (relative standard error), fit$phiH, (individual hessian), fit$phiC (individual covariances), fit$phiR (individual correlation matrices)

• Can also use Shi 2021 differences in addition to Gill differences. In our tests (using the same datasets as CPT) these produced worse estimates than the Gill 1983, though it is unclear why since it should be a faster more accurate method. A modified version is used in calculating the individual Hessians of numerically for the generalized likelihood approach.

• Generalized likelihood estimation is now present in nlmixr2est for focei, foce and posthoc

• nmNearPD() is a function you may use for nearest positive definite matrix. This is derived from Matrix::nearPD() but is implemented in C/C++ to be used in (possibly threaded) optimization.

• Individual Hessians can be accessed by $phiH, covariance by $phiC, eta standard errors by $phiSE and eta RSEs can be accessed by $phiRSE. There are eta aliases for these as well ($etaH, $etaC, $etaSE, and $etaRSE).

• Can now access the individual point’s contribution to the overall likelihood when merging to the original dataset. These merges can be accessed with $dataMergeFull, $dataMergeLeft, $dataMergeRight, and $dataMergeInner. The columns with the individual data column is nlmixrLlikObs.

  To calculate the total focei/foce objective function, the sum of the likelihoods still need to be adjusted by the omega/eta contribution, and the individual Hessians, and possibly the NONMEM objective function offset constant.

Censoring fixes

• Fixed bug where datasets with censoring that are not lower case cens and limit do not produce the correct table output (#180)

FOCEi updates

• Resets now scale properly when a value is simulated outside the limit
• Models with zero gradients on the first step now switch to bobyqa by default. With this, it is more important to examine the model parameters and fits for plausibility.

nlmixr2est 2.0.8

New features

• Add pd/npd as an output as well as npd/npde

SAEM bug fix

• When loading a nlmixr2 “saem” fit from another R session, nlmixr2 will no longer crash with fit$objf

NPDE/NPD fixes

• NPDE was identical to NPD even with correlated models, this was fixed (prior output was actually NPDE).

Censoring fixes

• FOCEi censoring fixes:
  • M4 method equation bug fix
  • M4 method derivative change based on equation fix
  • M2 method added missing derivative
  • Censoring already dTBS
• SAEM Censoring fixes:
  • SAEM method M4 method equation bug fix
  • Censoring limit changed to dTBS
• Censoring handling was unified

Internal changes

• Added ui$getSplitMuModel which is used in babelmixr2 and will be used in the refined stepwise covariate selection of nlmixr2extra

• Added work-around to remove _nlmixr2est_RcppExport_registerCCallable since the registering of C callable are handled manually at the moment.

nlmixr2est 2.0.7

• Use .zeros() for the matrices in armadillo in addition to relying on calloc to give zero matrices.

• Fixed one uninitialized object

• Fix for augPred so it works on population only models

• nlme no longer sets options to treat all covariates as non mu-referenced covariates, but directly calls a function that can turn on or off the mu-reference covariate selection.

• vpcSim now tries to simulate IDs that didn’t simulate correctly (with a warning)

• Export nmObjHandleControlObject

nlmixr2est 2.0.6 – new package

nlmixr2est contains the estimation functions within nlmixr2.

FOCEI family changes

• Remove lower level foceiFit function. Focei, foce, fo, foi, and posthoc now directly takes rxode2 ui objects

• New error types are supported in focei including mixing theta and etas in residual errors and different types of proportional errors

• Different types of additive and proportional errors can be used for each endpoint using + combined1() or + combined2() otherwise it takes the supplied addProp option to figure out which type of combined model is run (by default combined2())

• Focei model cache is now named focei-md5Digest.qs and uses qs compression/saving/loading.

• foceiControl() aligned between other methods.

• foceiControl(adjLik=TRUE) uses the NONMEM-style objective function throughout. foceiControl(adjLik=FALSE) uses the adjusted objective function throughout, and adjusts it back to the NONMEM objective function.

• Lag time and other between subject variability differences no longer calculate an ideal relative step size, but an absolute step size when using Gill differences (default)

• Objective function checks for infinite/NaN/NA values for the entire solving space and ensures no overflow occurs when calculating the inner hessian

SAEM changes

• mu referencing is no longer required for saem; Internally non mu-referenced values are converted to mu referenced values and the converted back when calculating the nlmixr2 object.

• nlmixr2 forced the parameter ordering to (1) population effects,

  (2) non mu-referenced between subject effects (3) omega estimates and (4) residual effects. This changes the order that nlmixr2 sees the parameters. Since this is based on a random number generator, the optimization trajectory will be different and have different results than nlmixr

• Components of omega can now be fixed.

• Residual error components can also be fixed.

• When optimizing only one residual value, nlmixr2’s saem uses nlm from R, which is more efficient than the nealder-meade method.

• Lower level saem functions (like configsaem()) are not exported because they are increasingly difficult to use and convert to something standard; a few methods (like print, summary etc) are maintained to view the lower level object and for debugging it.

• Parameter history and print-out no longer includes fixed parameters.

• The model to calculate the residuals more closely matches the model used for estimation to remove small rounding differences that may occur in the models.

• Different types of additive and proportional errors can be used for each endpoint using + combined1() or + combined2() otherwise it takes the supplied addProp option to figure out which type of combined model is run (by default combined2())

• Parameter history and printout now uses standard deviation for additive only components, matching the estimation of the components.

• rxode2 solving options are now saved in the rxControl part of the saemControl(). That is saemControl(rxControl=rxControl(...)); This fixes any conflicting option names as well as allowing alignment between the control structures in focei, nlme and saem

• saemControl() aligned between other methods.

nlme changes

• nlme has been completely rewritten to directly run from the rxode2 UI

• nlme always tries to use mu-referencing (when available)

• Internally nlme now uses parallel processing for solving so it should be faster.

• nlmixr2NlmeControl() (which will overwrite nlmeControl()) documents and adds more options to nlme. Also aligned with other methods.

• weights, fixed, random can be specified in nlmixr2NlmeControl(). If so, then the nlme object will be returned.

• returnNlme is a new option that will return the nlme object instead of the traditional nlme object.

• nlme_ode and lme_lin_cmpt are both removed.

• rxode2 solving options are now saved in the rxControl part of the saemControl(). That is nlmeControl(rxControl=rxControl(...)); This fixes any conflicting option names as well as allowing alignment between the control structures in focei, nlme and saem

nlmixr2 object change

• With saem, the nlmixr2 function now saves/compresses the phiM information. This means the gaussian and Laplacians likelihoods can be calculated when you save the nlmixr object and then restore it later.

• The nlmixr2 object compresses infrequently used and removes many unneeded objects. Even with compression, the saem objects are often a bit bigger since they include the large phiM object.

• nlmixr2 now supports non-mu referenced ETAs in the fit$parFixed and fit$parFixedDf

nlmixr2 interface change

• nlmixr2 interface changed to use rxode2 UI

• keep and drop are added to tableControl to influence the end data-frame

• $simInfo uses a quoted expression for $rx instead of a string

• $simInfo$sigma is a diagonal matrix since now the normal simulation is controlled by the variability modeled as a population value.

• nlmixr2 now allows etas that have initial omega estimates of zero to be dropped from the model (instead of issuing an error about a non-positive definite $omega matrix)

NPDE changes

• Fixed a bug where the number of simulations for a NPDE calculation are correctly passed by addNpde(fit, table=tableControl(nsim=500))

VPC changes

• vpc function rewritten and split out to vpcSim() and vpcPlot() (which is a replacement for vpc()).

• There were too many mismatches between vpc::vpc and nlmixr::vpc which caused inconsistencies in code based on load order of vpc and nlmixr. This way both coexist, and you can use the vpc simulation for other packages more easily (like ggPMX) without creating or summarizing data since ggPMX has its own methods for summarizing and creating plots.

• VPC now directly uses rxode2::rxSolve

augPred() changes

• augPred() has been written to use the new fit object.

• nlmixr2AugPred was changed to nlmixr2AugPredSolve()

• augPred uses the new interface and supports multiple endpoints. The endpoint name is now always on the plot(augPred(fit)).

getFitMethod() change

• Internally, fit estimation method is saved in fit$est, and now getFitMethod(fit) simply returns fit$est

Delete methods

• Many methods lower level utility functions have been deleted.

• nmDocx, nmLst and nmSave have been removed.

Bug fixes

• Now will reset the cache when items cannot be loaded. In the past error messages like function 'rx_0ba247452048de33b1ffb8af516714fc__calc_lhs' not provided by package 'rx_0ba247452048de33b1ffb8af516714fc_' would cause the estimation to stop. Now rxode2::rxClean() is run when this occurs.