Changelog

Source: NEWS.md

nlmixr2est (development version)

• In the mu-referenced estimation family (mfocei/ifocei, mfoce/ifoce, mfocep/ifocep, magq/iagq, mlaplace/ilaplace), the regression-updated mu thetas (population and covariate coefficients) now appear as standard columns, in natural theta order, in the live iteration print and in fit$parHist/parHistData; gradient rows record NaN for them (shown as blank console cells) and the previous appended | mu| row was removed.

• SAEM no longer errors when the between-subject-variability covariance (Omega) collapses to a non-positive-definite matrix mid-run. It is now projected to the nearest positive-definite matrix (with a one-time user-visible warning recorded in fit$runInfo) instead of stopping with inv_sympd(): matrix is singular or not positive definite.

• Fixed the foceiControl(fast=TRUE) analytic outer gradient for FOCEI models whose residual variance depends on the prediction (prop(), add()+prop(), combined1, pow(), add()+pow()). The (f,R) determinant chain rule aliased d(dfr)/df onto pffR = d(dff)/dR. Those coincide only when the determinant coefficients are second partials of a potential, which holds for the exact censored Laplace determinant but not for the Gauss-Newton expected information used on a normal observation: there (dff,dfr,drr) = (1/R, 0, 0.5/R^2), so d(dfr)/df = 0 while d(dff)/dR = -1/R^2. The aliasing injected a spurious -a(s)/R^2 * (a_l*aR_m + aR_l*a_m) term into dHt/ddir. It is proportional to aR = dR/ddir, so additive error, lnorm, transforms on an additive endpoint, and every FOCE variant (frozen variance) were unaffected – only FOCEI with a prediction-dependent variance was wrong, by 20% to 600% depending on the error model. The bad gradient made foceif/mufoceif stall with false convergence (8) short of the focei optimum. Verified exact against Richardson-extrapolated central differences across add, prop, add+prop, combined1, pow, add+pow, lnorm, boxCox and yeoJohnson.

• The analytic-covariance augmented model now relies on rxode2::rxOptExpr()’s own chunked (and rxControl(cores=)-parallel) optimization instead of a hand-rolled line chunker. Every rxOptExpr() call site (focei, nlm, nls, saem, nlme, rpem and the analytic covariance) now passes the fit’s rxControl(cores=) through as parallel=, so the chunks are optimized with the same thread setting the solves use. Requires an rxode2 with chunked rxOptExpr() (current rxode2 main).

• Added vaeCovariates(), which returns the subject-level covariates that est = "vae" would explore during automated covariate selection (name, type and centering value), using the same discovery rules as the fit; warn=FALSE silences the time-varying-covariate exclusion warning. Exported so other packages can inspect the search space without running an estimation.

• Added an automatic differentiation variational inference method (est = "advi", Kucukelbir et al. 2017) with adviControl(). The variational gradient is obtained from the FOCEi forward sensitivities (the inner per-subject eta gradient and the outer theta-sensitivity score) rather than automatic differentiation, and the whole optimization (the adaptive step-size-scale search plus the stochastic-gradient-ascent loop) runs in a single C++ call. It supports mean-field and block full-rank variational families (adviFamily), both a point-estimate (variational-EM, output comparable to FOCEi/SAEM) and a full-Bayes (pointEstimate=FALSE, with a full-rank variational posterior over the population parameters) mode, mu-referenced and non-mu structural thetas plus residual error, multiple endpoints and BLQ censoring (via the reused inner likelihood), the paper’s adaptive step-size with a step-size-scale search (adaptEta), and a warm-resume API (adviControl(resume=)). The reparameterization noise is drawn from a counter-based stream keyed by the global iteration index, so a shorter run is a bit-for-bit prefix of a longer one and results are independent of the thread count. The run prints the standard nlmixr2 iteration table (like saem/vae, adviControl(print=)) with the step-size search and the main run shown as labeled stages, and the walk is saved as standard parHistData (fit$parHist). The inner and theta-sensitivity models are set up once and reused by the output step (no model re-build at finalize). The per-subject ELBO and gradient are computed in parallel over subjects (rxControl(cores=)); a serial id-ordered reduction keeps the result identical for any thread count.

• Added quasi-random (Sobol) importance sampling to the imp/impmap estimation methods (QRPEM, Leary & Dunlavey PAGE 2012): impmapControl(qr = TRUE) draws the E-step samples from a low-discrepancy Sobol sequence mapped through the inverse normal CDF, so the E-step integrals converge at O(1/N) instead of O(1/sqrt(N)) – more accurate at the same isample and a smoother objective trace. qrShift controls the Cranley-Patterson randomization of the point set and qrRefresh whether the shift is redrawn each iteration (FALSE makes each EM iteration a deterministic map). The Monte-Carlo covariance (impCov = TRUE) uses the quasi-random points too. Results stay reproducible and thread-count independent.

• Added SIR (sampling-importance-resampling) acceleration for the imp/impmap M-step: impmapControl(sir = TRUE, sirSample = ) runs the non-mu / residual-error Newton update on an equal-weight systematic resample of the importance samples, cutting the theta-sensitivity solves from isample to sirSample (default max(25, ceiling(isample/10))) per subject per iteration.

• Added est = "qrpem" with qrpemControl(), sugar for the impmap importance-sampling EM with qr = TRUE and sir = TRUE.

• Fixed impmapControl(impSeed = ) being ignored: the FOCEI solve setup reset the ambient seed before the importance-sampling E-step, so every impSeed produced the identical draw. The E-step / covariance / QR-shift / SIR streams now derive from impSeed directly, so a different impSeed gives a different (but reproducible and thread-count-independent) importance sample.

• The est = "vae" training loop (variational-autoencoder NLME) now runs entirely in C++ (vaeTrainCpp_): the burn-in / KL-anneal / EM / smoothing schedule, the LSTM encoder forward/backward, the FOCEi inner-likelihood evaluation, the closed-form M-step (including BICc-ELBO covariate selection), Adam, and the iteration/parameter-history walk. Each gradient step now re-parameterizes the inner problem in place (updateTheta) instead of re-running the full FOCEi setup (rxSymInvCholCreate + foceiSetup_) it did every step in R, so VAE fits are substantially faster.

• The analytic FOCEI/FOCE outer gradient and observed-information covariance (foceiControl(fast = TRUE) / covType = "analytic") now build a smaller augmented sensitivity model by reusing eta sensitivities for mu-referenced covariate coefficients. A covariate coefficient b enters through an eta’s parameter (cl = exp(tcl + eta.cl + b*cov)), so for a covariate that is constant within each subject every sensitivity of b equals the linked eta’s scaled by the covariate value (df/db = cov*df/deta). Those covariate directions are therefore emitted as algebraic scaled copies instead of integrating their own state-sensitivity ODEs, shrinking the compiled model ~25-58% and cutting its (super-linear) build time ~1.3x-4x on larger covariate models – the result is bit-identical. Detection reuses rxode2’s own mu-reference classification (muRefCovariateDataFrame + mu2RefCovariateReplaceDataFrame), covers bare and algebraic (log(WT/70), WT-70) covariates, and falls back to the full symbolic build for anything the eta-scaling identity does not cover (time-varying covariates, a shared/reused eta). The analytic covariance also now restricts itself to Gaussian endpoints (t/cauchy/count/ordinal likelihoods and multiple estimated transform lambdas fall back to the finite-difference covariance).

• Extended that reuse to a covariate on an eta-less parameter (v = exp(tv + b*cov) with no eta.v): the coefficient now reuses the structural theta’s own sensitivity direction (df/db = cov*df/dtv, since tv and b enter the mu identically) instead of building its own state-sensitivity ODEs. Previously such a model errored inside the direction map and silently fell back to the finite-difference covariance; it now stays analytic for both the covariance and the fast outer gradient (matching finite differences), and holds the integrated-direction count fixed regardless of how many covariates sit on eta-less parameters (~1.3x fewer for one, ~2.6x for three). A structural-theta occurrence guard (shared with the eta guard) keeps the reuse exact.

• Fixed the analytic (fast=TRUE) outer gradient and covariance being wrong for a model where a mu-referenced parameter’s random effect is shared across parameters (e.g. eta.cl used in both cl and v). The mu-referenced theta reused that eta’s state sensitivity, but df/dtheta (one parameter) differs from df/deta (all the parameters the eta appears in), so the gradient/covariance for that theta was incorrect. Such a theta now gets its own true-sensitivity direction (the eta keeps its own), so the gradient and covariance stay analytic and are correct.

• Added a fast-SAEM (f-SAEM, Karimi, Lavielle and Moulines 2020) simulation step: saemControl() gains fast/fastKernel/fastCov/fastIter/fastLik options and a new est = "fsaem" method (sugar for saemControl(fast = TRUE)). When enabled, the early SAEM iterations replace the random-walk Metropolis simulation with an independent Metropolis-Hastings kernel whose Gaussian proposal is centered at each subject’s conditional MAP (reusing the FOCEi inner likelihood), which reaches the MLE in fewer iterations; later iterations degrade to the standard kernels. Supports continuous single-endpoint models (additive, proportional or combined error), including non-time-varying mu-referenced covariates (absorbed into the per-subject prior mean of the inner); models outside this envelope (e.g. time-varying covariates, mixtures) transparently run standard SAEM.

• saem and fsaem now fit general log-likelihood endpoints (ll(name) ~ <expr>, e.g. Weibull/exponential time-to-event), in the style of saemix likelihood models: the model returns the per-observation log-likelihood, the simulation step uses it directly as the observation loss, and the population parameters and between-subject variances are estimated by the standard SAEM M-step (no separate residual error). Fixed-effect-only parameters of such a model are refined by a direct L-BFGS-B optimization of the observation likelihood (saemControl() gains lbfgsLmm/lbfgsFactr/lbfgsPgtol/ lbfgsMaxIter, derived from sigdig like foceiControl()), respecting the parameter bounds from the model.

• The fsaem IMH kernel now draws its proposals from rxode2’s threefry engine (seeded from the fit’s seed plus the chain/subject index), so results are reproducible regardless of the number of threads. Proposals for a bounded log-likelihood parameter are re-drawn up to saemControl(nRetry = 10) times and then clamped to the violated boundary. The whole saem/fsaem fit now runs inside rxode2::rxWithSeed() so a session’s first fit is seeded and fits never contaminate each other’s RNG state.

• New estimation methods est = "impmap" and est = "imp": importance-sampling expectation-maximization in the style of NONMEM’s METHOD=IMP, with the E-step proposal centered at each subject’s MAP mode (impmap) or at the running conditional mean (imp, NONMEM MAPITER=0). Options are set with impmapControl() / impControl(), including the number of importance samples (isample), maximum EM iterations (nIter), proposal scale adaptation (gamma, iscaleMin, iscaleMax, iaccept), a NONMEM-style windowed convergence test (ctol, nConvWindow), a reproducible thread-independent sampling seed (impSeed), and an experimental Monte-Carlo covariance (impCov=TRUE). Mu-referenced models (including covariates), mixture models, parameter bounds, and fix()ed thetas/omegas are supported; the reported objective is a FOCEi evaluation at the EM estimates so fits are comparable with the focei family.

• The “initial ETAs were nudged” warning is now only raised when an ETA actually stayed at zero and a nudge was performed, not merely when the nudge check ran (which previously produced a spurious warning on well-behaved fits).

• For a fully mu-referenced model (every eta mu-referenced, so the initial etas are all zero) a user-specified non-default foceiControl(mceta=) now falls back to the default (mceta=-2) with a warning, since the mceta starting-point search has nothing to explore.

• The reported fit timing (fit$time) again attributes the symengine model build and rxode2 compilation to setup (focei family) / configure (saem) instead of leaking it into the other bucket. The nlm-family methods (nlm/nlminb/bobyqa/newuoa/uobyqa/n1qn1/lbfgsb3c/optim/nls) now time their preprocessing and EBE model build as setup and the model build + optimizer as optimize, instead of leaving nearly all of it in other.

• Fixed nlmControl() listing eventSens/sensMethod twice, which made a control round-trip (do.call(nlmControl, ...)) error with “formal argument … matched by multiple actual arguments”.

• fast=TRUE (and the *f methods) with a linCmt() model now downgrades to fast=FALSE up front with a message instead of re-attempting the symengine augmented-model build (and silently falling back to finite differences) on every outer-gradient call; the analytic covariance likewise reports the linCmt() fallback instead of failing silently.

• Fixed a segfault in the analytic covariance step for models out of analytic scope (e.g. models using tad()/podo()): building the augmented sensitivity model frees the fit’s global solve, and the finite-difference sandwich fallback then solved against freed memory. The freed solve is now restored before the fallback runs.

• Iteration printing now labels the estimation phase on the back-transformed (X) row: est = "vae" shows Burn in/KL anneal/EM/Smooth in the objective column (with a key legend in the header), and est = "saem" tags the row SA: X / EM: X for the burn-in and EM phases.

• Fixed muModel = "lin"/"irls" erroring with “undefined columns selected” when a model has two or more mu-referenced covariates that are expressions (e.g. log(WT/70)) rather than bare data columns (#711).

• est = "vae" now assembles its fit with nlmixr2CreateOutputFromUi() and computes its covariance directly in the focei covariance step; vaeControl(covMethod=) now takes the focei choices ("analytic" (default), "r,s", "r", "s", "") instead of "linear".

• Fixed the analytic (fast=TRUE) outer gradient never being used during live optimization: the gradient hook read omega/theta/etaObf from the fit environment, which are only written at finalize, so every call silently fell back to finite differences. The C++ hook now refreshes that state each gradient call; the cached augmented-model metadata (foceiModel$outerMeta) carries all fields the batched solve needs.

• The mu-referenced families (mfoceif/ifoceif/…) now consume the analytic gradient on the profiled (mu-reduced) parameter set; a gradient/parameter-set size mismatch stops (never a silent FD fallback) and FD is only used when the sensitivity system fails to solve (with a one-time warning).

• The mu-referenced FOCEI families (mfocei/ifocei/…) now profile plain (covariate-free) mu-referenced population thetas out of the outer optimizer via the in-C++ regression as well (intercept-only groups), so outer gradients – numeric or analytic – are only calculated for the non-mu-referenced parameters (residual errors, omegas, non-mu thetas); user-fixed mu-referenced thetas remain outer-optimized. The regress/re-optimize cycle defaults were tightened (muModelTol 1e-3 -> 1e-5, muModelMaxCycles 10 -> 20) so the profiled fits converge to the same optimum as the base methods (usually faster, since the outer optimizer sees a well-converged profiled objective).

• Renamed the mu-referenced FOCEI-family estimation methods (all introduced in this development version): the irls* methods are now i* (ifocei, ifoce, ifocep, iagq, ilaplace and fast variants ifoceif/ifocef/ifocepf) and the mu* methods are now m* (mfocei, mfoce, mfocep, magq, mlaplace and mfoceif/mfocef/mfocepf), with matching control functions (ifoceiControl(), mfoceiControl(), …). The old names are removed (they never shipped in a release).

• The analytic (fast=TRUE) outer gradient and covMethod="analytic" now handle models whose residual-error parameters are all fixed. Such a fit’s outer problem is omega-only (structural thetas mu-profiled, residuals fixed), and omegas do not enter the ODE, so no model re-solve per parameter is needed; previously the analytic path fell back to finite differences (re-solving the ODE for every omega step). The general (f,R) assembler now runs with no free residual direction, reading the (fixed) variance from the solved rx_r_.

• Bounded mu-referenced parameters (population thetas and covariate coefficients) are now also profiled by the mu/irls regression in the clamped m*/i* family (mfocei/ifocei and variants, muModel != "none"), with the update clamped to the bounds (box-constrained least squares, foceiControl(muModelClampRetries=)); parameters that were clamped during the fit are reported once as a fit note. Every other (non-clamped) method keeps the previous style: a bounded population theta rejects its whole mu group (with a warning) and a bounded covariate coefficient stays an ordinary outer-optimized parameter.

• Fixed the mu-family regression’s handling of a user-fixed covariate coefficient: its (fixed) contribution was added into the regression target but never taken back out, so each regression pass shifted the group’s linear predictor by the fixed term (biasing the population theta and inflating the eta variance). The fixed contribution is now left out of the regression entirely.

• fast=TRUE now defaults the outer optimizer to lbfgsb3c (FD methods keep nlminb); an explicit outerOpt is honored.

• The iteration print and $parHistData track analytic gradients as their own type (A/"Analytic Gradient"), and the fit header reports the gradient actually used and the mu-model variant, e.g. (outer: lbfgsb3c; grad: analytic; mu: irls).

• est = "vae" iteration print now shows back-transformed (X) parameters (e.g. exp() thetas) like focei/saem, drops the redundant unscaled (U) block (vae never scales), and the development-time verbose print was removed.

• Mixture models: parallel inner loops now solve mixture components strictly one at a time (components share the physical subjects’ data structures), parallelizing over subjects within each component.

• Fixed a segmentation fault in est = "vae": the parallel inner-likelihood loop now caps its thread count at the rxode2 solve’s op->cores (1 for models rxode2 flags not thread safe, e.g. linCmt() gradients), matching the other FOCEi parallel loops.

• SAEM no longer errors with No data with ID: <id> for a subject that has a dose but no usable observation (e.g. all of its DV values are missing, which rxode2::etTrans() converts to EVID==2 records). The shared preprocessor (.foceiPreProcessData()) now drops any subject without an observation – emitting the same IDs without observations dropped message rxode2 already uses for dose-only subjects – and the shared table builder re-inserts those subjects’ rows with a population PRED (solved at eta = 0) and NA individual columns (IPRED, etas, residuals). Every estimation method now handles observation-less subjects the same way: the subject is reported in $runInfo and appears in the output with a population prediction and NA individual values, and the SAEM-specific guard in .configsaem() is removed (#687).

• Censored (M2/M3/M4/BLQ) observations now use the exact censored second derivative for the FOCEI inner Laplace Hessian (foceiControl(censOption = "laplace"), the new default); censOption = "gauss" keeps the historic uncensored Gauss-Newton curvature.

• Added foceiControl(censOption = ...) to choose the censored (M2/M3/M4/BLQ) second-derivative treatment. "gauss" (the default) keeps the historic uncensored Gauss-Newton curvature, matching common PMx tools; "laplace" uses the exact censored second derivative for the inner Laplace Hessian and analytic covariance. Non-censored fits are unchanged. The option is accepted by saemControl/nlmControl for a uniform interface but is inert there (SAEM has no Laplace inner Hessian; NLM uses a finite-difference Hessian that already reflects censoring), so their censoring text stays plain while FOCEI/FOCE note the treatment used (e.g. "M3 censoring (gauss)").

• The mu-referenced FOCEI families (mfocei/ifocei/mfoce/mfocep/…) now compute their covariance on the full corresponding focei/foce/focep model, at the mu fit’s converged theta and eta with the inner problem frozen (as if the full model had produced that point), instead of the mu->phi reduced model used during estimation. This fixes incorrect standard errors on the mu-referenced/covariate (“linear”) parameters (the sandwich covMethod="r,s" was the most affected) and honors the requested covMethod.

• covMethod="analytic" now falls back to the finite-difference sandwich ("r,s") when the analytic covariance is unavailable for a model (e.g. linCmt()), instead of the R-matrix ("r") alone.

• The analytic fast gradient and covMethod="analytic" now cover matExp() / indLin() (matrix-exponential / inductive-linearization) models for FOCEI, FOCE, and foce+ (and their mu/irls families), matching the equivalent ODE fit; an indLin() forcing state no longer misorders compartments in the augmented model.

• The analytic fast outer gradient now covers foce+ (focep), matching the finite-difference gradient (its live-conditional-R kernel was already used by the analytic covariance).

• The analytic fast gradient now covers censored M2/M3/M4 observations for FOCEI (both censOption values) and FOCE (default censOption="gauss"); the FOCE EBE is re-solved with the exact censored score at the frozen variance. The reported censoring text notes the second-derivative treatment used (e.g. "M3 censoring (laplace)").

• covMethod="analytic" now covers censored M2/M3/M4 observations for FOCEI and FOCE with the default censOption="gauss" (censored score partials with the Gauss-Newton determinant, and a censoring-aware FOCE EBE re-solve); censOption="laplace" still uses the FD covariance.

• Fixed the sign of the M2 (interval) upper-tail term in the censored inner gradient, which had shifted M2 EBEs and objective values.

• Added foceiControl(fast = TRUE) to compute the FOCEI outer (population) gradient analytically from Almquist (2015) sensitivity equations instead of by finite differences; out-of-scope models fall back to the finite-difference gradient. Under fast, the outer optimizer defaults to "lbfgsb3c"; pairing it with a derivative-free outerOpt reverts to fast = FALSE.

• The analytic fast outer gradient now solves every subject’s augmented sensitivity model in one OpenMP-threaded rxode2 population solve (for both the FOCEI and FOCE families) instead of one R solve per subject, which previously dominated the gradient cost.

• The analytic fast outer gradient now supports a both-sides transform with an estimated boxCox/yeoJohnson lambda: lambda enters as a prediction-sensitivity direction plus the DV-transform residual chain and the transform Jacobian.

• covMethod="analytic" now supports an estimated boxCox/yeoJohnson lambda for FOCEI, FOCE, and foce+ (the observed-information carries the DV-transform 2nd-order chain).

• Fixed the post-fit foceiCovAnalytic()/getVarCov() recompute falling back to the finite-difference covariance for a general (f,R) variance (multi-endpoint / non-add-prop / estimated-lambda transform): it now routes to the same general (f,R) assembly as the live fit, reproducing the installed analytic covariance instead of returning NULL.

• Fixed the analytic FOCE/foce+ covariance silently falling back to the finite-difference Hessian: the per-subject eta=0 solve list dropped its NULL foce+ slots (list[[i]] <- NULL shrinks the list), so foce+ (and some FOCE) covariances went out of bounds.

• Fixed the general (f,R) analytic covariance (multi-endpoint / non-add-prop variance / both-sides transform) reporting covMethod="r" instead of "analytic".

• Added the *f convenience estimation methods – focef, focepf, foceif, mfocef, mfocepf, mfoceif, ifocef, ifocepf, ifoceif – each equivalent to its base method with foceiControl(fast = TRUE) as the default (the analytic outer gradient + Eq-48 warm-start).

• The analytic gradient/covariance augmented model now carries modeled dosing parameters (f(), lag()/alag(), rate(), dur()) and their second-order dose-based (“jump”) sensitivities via rxode2’s eventSens = "jump"; previously such models produced an incorrect (dose-unscaled) augmented solve.

• The default outer optimizer (foceiControl(outerOpt=)) is now "nlminb" for the finite-difference methods (and "lbfgsb3c" when fast = TRUE).

• foceiControl(mceta=) now defaults to -2: when the analytic gradient supplies the EBE sensitivity (fast = TRUE), the next inner-problem starting eta is the Almquist (2015) Eq-48 extrapolation eta* + (d eta*/d theta)(theta_new - theta_old), accepted only within the standardized-eta reset bound (else the last eta, or 0 when that is also out of bound). mceta = -1 jumps between the extrapolation and eta=0. Both fall back to keeping the last eta when no analytic sensitivity is available (fast = FALSE), so the default is behavior-preserving there.

• The analytic covariance (covMethod = "analytic") now supports foce = "foce+" (the live conditional residual variance), including the focep method; previously it fell back to the finite-difference covariance.

• The FOCEI covType control was removed; the analytic-vs-finite-difference R-matrix choice is now part of foceiControl(covMethod = c("analytic", "r,s", "r", "s", "")), with "analytic" (the exact observed-information R-matrix) the default.

• foceiControl(covSolveTol=) now also tightens the ODE solves behind the finite-difference covariance methods (previously it only applied to the analytic augmented-sensitivity solves).

• The SAEM iteration history (parHistData) now records the off-diagonal Omega covariances of declared blocks (cov.<eta>.<eta>), alongside the existing diagonal variances and residual parameters.

• SAEM now reports a full theta + residual + Omega covariance by default (saemControl(covFull = TRUE)): the linearized-FIM (covMethod = "linFim") variance block (Omega variances/covariances and residual error parameters) is added to the structural-theta block, Omega rows named by the random effect (om.eta.cl / cov.eta.cl.eta.v), and the residual SEs are surfaced in the parameter table. covFull = FALSE keeps the legacy structural-theta-only covariance.

• Added saemControl(covMethod = "sa"), a stochastic-approximation Fisher Information covariance (Kuhn & Lavielle 2005, as used by Monolix): a dedicated post-estimation phase (nSaCov iterations) holds the parameters at the converged estimate and Monte-Carlo averages the Louis observed-information into a converged full theta + Omega + residual covariance.

• Fixed covMethod = "fim": the SAEM Fisher information omitted the deterministic mu-block complete Hessian, leaving its fixed-effect block indefinite (fim produced NaN standard errors). The mu Fisher information is now added, so fim inverts to a valid positive definite full theta + Omega + residual covariance.

• covMethod = "fim"/"sa" now report off-diagonal Omega covariances and proportional/combined residual standard errors by splicing them from the linearized-FIM variance block (the analytic simulation FIM keeps the structural-theta block).

• The analytic covariance now reports fit$covMethod as "analytic" (instead of "r"), and its Omega variance/covariance rows are named by the random effect – om.eta.cl / cov.eta.cl.eta.v – rather than by the mu-referenced theta (om.tcl).

• Residual (error-model) parameters are now part of the focei-family covariance: their standard errors are estimated alongside the structural thetas (skipCov no longer skips them; only fixed, IOV, and mlogit-scale mixture-probability thetas are skipped).

• foceiCovAnalytic() now caches its result on the fit and installs the covariance as fit$cov, so repeated calls and getVarCov() reuse it instead of recomputing the augmented sensitivity solve every time.

• Added the focep, mfocep, and ifocep estimation methods – the foce, mfoce, and ifoce methods with foce = "foce+" forced (the live conditional residual variance R).

• Added foceiControl(foce = c("nonmem", "foce+")) to choose how FOCE evaluates the residual variance R: "nonmem" (new default) freezes R at the eta = 0 population prediction to match NONMEM FOCE, while "foce+" keeps the live conditional R (the FOCE behavior in nlmixr2est 6.0.1 and earlier), which can be slightly more accurate but does not match NONMEM and falls back to the finite-difference covariance under covType = "analytic".

• Added foceiControl(covType = "analytic"), an exact analytic observed-information covariance for FOCEI and FOCE fits. The R-matrix is assembled in closed form – a data term from the analytic 2nd-order sensitivities (rxode2 .rxSens) and a log-determinant term whose 3rd-order tensor is recovered by Shi (2021) central differences of those sensitivities (keeping the augmented ODE at O(ndir^2)) – matching NONMEM $COV MATRIX=R. It is computed while the optimizer is live and covers additive, proportional, and combined error, and mu-referenced, covariate, and other non-mu-referenced structural parameters (and non-mu-referenced etas) as well as SD-scale inter-occasion variability; any fit outside its scope (FO, nAGQ > 1, censoring, DV-transformed error, bounded-parameter transforms, a structural theta shared by two etas, or a non-SD iovXform – and a pure-proportional variance that vanishes at a near-zero prediction) emits a message and falls back to the finite-difference Hessian. FOCE (interaction off) uses the general total-derivative Hessian and re-solves the FOCE empirical-Bayes estimates to the FOCE inner stationarity before assembly (FOCE reduces to the FOCEI result for additive error). It defaults to covMethod = "r" (the observed-information inverse); an explicit covMethod = "r,s" or "s" is honored, with the analytic R feeding the native finite-difference sandwich / S-matrix. covFull chooses the "r" fit$cov shape: the structural-theta block (default, matching the finite-difference shape) or the full theta + residual sigma + Omega matrix, which adds the Omega/residual SEs the "r" matrix does not provide (via rxode2’s rxOmegaVarCovDeriv). The augmented-solve tolerance is derived from sigdig (override with covSolveTol). The default remains covType = "fd".

• covFull = TRUE now also works with covType = "fd": the full theta + residual sigma + Omega covariance is computed by central finite differences of the objective over the same parameter set the analytic engine uses, at the finite-difference covariance seam. The Omega block is perturbed on the variance-covariance scale directly (natural scale, no Cholesky Jacobian). It reuses the same infrastructure as the theta-only foceiCalcR: the per-parameter step is chosen by the Gill-Murray-Saunders-Wright (1983) optimal finite-difference-interval routine (gill83), and the Hessian uses the 5-point diagonal / 4-point off-diagonal stencils, so the multi-random-effect Hessian stays positive-definite. It matches the analytic covariance (and NONMEM $COV MATRIX=R) to finite-difference tolerance, handles IOV (the occasion-variance SD is an ordinary theta), and gives users the residual/Omega SEs the theta-only FD cov omits.

• Fixed the FOCE (interaction = FALSE) objective function and empirical-Bayes estimates. The residual variance R entering the FOCE inner likelihood must be evaluated at the eta = 0 population prediction and held constant (the truncated Sheiner-Beal inner gradient drops dR/deta). The previous code froze R symbolically, which only removed explicit eta symbols: for ODE models a live model state stayed in R (so R remained eta-dependent), and for linCmt() models the freeze injected a second linCmt/linCmtB call that corrupted the inner sensitivities and stalled the eta optimization. R is now supplied at eta = 0 from the inner model itself, so FOCE with an ODE model agrees with the closed-form (linCmt) result and both match the NONMEM FOCE reference (Wang 2007). The eta = 0 R depends only on theta, so it is cached across the inner iterations (recomputed once per parameter update), keeping FOCE close to the FOCEI per-iteration cost.

• Fixed the FOCEI covMethod = "r" / "s" / "r,s" standard errors, which were inflated by a constant factor (√2 for "r", 2 for "s") because the R- and S-matrix covariances used 2*R^{-1} / 4*S^{-1} instead of R^{-1} / S^{-1}; they now match NONMEM $COV (#666).

• Added sensMethod to nlmControl()/foceiControl(); the nlm-family methods can now compute ODE parameter sensitivities with the in-engine discrete adjoint ("adjoint") using the matching s-method, or pick it automatically ("auto") when estimated thetas exceed ODE states, matching the forward result

• Matrix-exponential / inductive-linearization models (matExp() with optional indLin() forcing) now estimate with the focei family (focei/foce/foi/posthoc), the nlm family (nlm/nlminb/…), and SAEM, matching the equivalent ODE model. A hand-written matExp() model that used indLin() previously registered the forcing state as compartment 1, reversing it relative to the ODE and misplacing default (compartment-1) dosing; the generated cmt()/d/dt() declarations now order compartments source-first from the k_<from>_<to> graph. SAEM also now materializes the implied d/dt() for these models instead of erroring

• Added foceiControl(warm=c("calc", "save")); "calc" (new default) warm-starts each n1qn1 inner optimization from the eta Hessian calculated at the starting eta and the current theta (including ll()/dnorm() models with finite-difference Hessians); the Hessian is always recalculated rather than reused from an earlier outer evaluation, since theta moves between evaluations and a saved Hessian would be stale. "save" keeps the prior behavior

• Computing NPDE for a fit with a degenerate simulated covariance (e.g. a residual SD estimated near zero) no longer aborts the R session; the affected subject’s NPDE is set to NA instead

• matExp()/indLin() models now estimate with the focei family, the nlm family, and SAEM, matching the equivalent ODE model; compartments are ordered source-first from the k_<from>_<to> graph so default dosing is placed correctly

• Added sensMethod to the nlm-family controls (nlmControl(), nlminbControl(), optimControl(), n1qn1Control(), lbfgsb3cControl()) and to foceiControl() (focei/foce inner ETA sensitivities); ODE parameter sensitivities can be computed with the forward sensitivities ("forward") or the in-engine discrete adjoint ("adjoint"). When left at "default", the method is taken from the global option getOption("nlmixr2est.adjoint") (default "forward"), so the package-wide policy can be set in one place

• fo/foi now force forward sensitivities (adjoint does not apply to the eta=0 linearization), and the adjoint base-method restore in the focei family is a strict no-op for forward fits, fixing fo/foi tables/residuals

• Fix a fit aborting with initial 'omega' matrix inverse is non-positive definite when a degenerate fit (e.g. SAEM collapsing an uninformative random-effect variance to 0) leaves a singular omega; the sym-inv-chol setup now nearPD-corrects it so the residual/table diagnostics still run

• Fix SAEM erroring with No data with ID: <id> for a dosed subject with no usable observation; such subjects are now dropped before estimation and re-inserted into the output with a population PRED and NA individual columns, like FOCEi (#687)

• Internal consolidation of data preparation and the nlm-family control/fit functions across estimation methods; no change to any fit result

• Fix cov2cor error when omega has exactly one nonzero diagonal

• Fix SAEM linearized-FIM covariance (covMethod = "linFim") erroring when exactly one covariate-model parameter is estimated

• Internal consolidation of data preparation across estimation methods (no change to any fit result). The shared preprocessor .foceiPreProcessData() already fed every method; this removes the duplication layered on top of it: two never-called data-setup functions (.nlminbFitDataSetup, .nlsFitDataSetup) were deleted; the column-name normalization and the time-varying-covariate detection were each extracted into a single shared helper (.nmUpcaseNonCov, .nlmixrTimeVaryingCovariates); and the nine nlm-family *FamilyControl/*FamilyFit functions (nlm, nlminb, bobyqa, newuoa, uobyqa, n1qn1, lbfgsb3c, optim, nls) were collapsed onto two generics (.nlmFamilyControlGeneric, .nlmFamilyFitGeneric). SAEM’s internal event-table dv column drop in .configsaem() is now by-name with a layout assertion instead of a positional index.

• Fix FOCEi aborting R with Cube::slice(): index out of bounds when mceta >= 1 and maxInnerIterations == 0

• Fix Windows heap-corruption segfault for gradient/pooled estimator fits at more than one core

• Fix the SAEM linearized-FIM covariance (covMethod = "linFim") erroring (or falling back) when exactly one covariate-model parameter is estimated, due to a vector-collapse transpose bug in calc.COV().

• Defensively use drop = FALSE when subsetting the omega covariance matrix for the correlation (cov2cor) calculation, so an omega with exactly one nonzero diagonal element does not collapse to a scalar and trigger a “‘V’ is not a square numeric matrix” error.

• Fix cov2cor error when omega has exactly one nonzero diagonal by subsetting with drop = FALSE

• SAEM covariance no longer errors with Error in rxode2::rxInv(.tmp) : Not a matrix. for models with a single population parameter (which arises e.g. with M2/M3/M4 censoring). The covariance fallback inverts a subset of the FIM (.saem$Ha[1:.nth, 1:.nth]); for a single parameter that subset dropped from a 1x1 matrix to a scalar, which rxode2::rxInv() rejects. The subset now keeps drop = FALSE.

• The parallel test suite (Config/testthat/parallel) now gives each worker its own rxode2 model-compile directory and sizes the worker pool to the host (capped at 2 on CRAN). Previously all workers compiled models into one shared cache directory and raced, producing spurious “error building model” failures and the 6h CI timeouts.

• The parallel test suite (Config/testthat/parallel) is now robust on

  101. Three independent problems were fixed:
  • Each worker gets its own rxode2 model-compile directory; workers previously raced on one shared cache directory, producing spurious “error building model” failures and 6h timeouts.
  • The worker pool is sized to the host (detectCores()/2, 1 on CRAN) and can be pinned with the NLMIXR2_TESTTHAT_CPUS environment variable (needed inside containers, where detectCores() reports the host rather than the cpuset allotment).

• The test suite is now robust on CI. Three independent problems were fixed:

  • Each test process gets its own rxode2 model-compile directory; with Config/testthat/parallel the workers previously raced on one shared cache directory, producing spurious “error building model” failures and 6h timeouts.
  • The suite runs a single testthat worker by default (overridable with the NLMIXR2_TESTTHAT_CPUS environment variable) with within-solve threads kept at 2, so workers x threads no longer saturates every CPU – leaving a core for the CI runner’s heartbeat agent. Saturating all cores was starving the agent and killing the devel/oldrel jobs with exit 143 (“the runner has received a shutdown signal”).
  • BLAS/OpenMP thread pools are now capped in the CI environment (OPENBLAS_NUM_THREADS etc. in the workflow) instead of via Sys.setenv() in tests/testthat.R. OpenBLAS is loaded at process startup and reads that variable only then, so the in-R setting was always too late and OpenBLAS ran one thread per core.

• The test suite runs a single testthat worker on CI and on CRAN (so it does not oversubscribe a core-limited runner) and parallel (Config/testthat/parallel) elsewhere; rxode2’s within-solve threads are capped to 2 only on CRAN and left to rxode2’s own management otherwise.

• Fix SAEM covariance error (rxInv(.tmp): Not a matrix) for models with a single population parameter

• Test suite uses a single testthat worker on CI/CRAN and parallel elsewhere; rxode2’s within-solve threads capped to 2 only on CRAN

• fit$time now reports every estimation stage consistently

• foceiControl() now defaults to outerOpt = "lbfgsb3c" (previously "nlminb") and sigdig = 4 (previously 3). rxUiDeparse() of a foceiControl() correctly omits outerOpt when it is left at this default.

• Aggregated ODE-solve warnings flushed during focei/saem estimation (e.g. [lsoda -- internal t + h = t ...]: N warning(s) for subject(s): ...) now report the real subject id instead of Unknown. Estimation passes a declassed data frame to rxode2, so rxode2’s subject-id factor table was empty during the fit; the fit’s idLvl is now pushed into rxode2 via rxSetIdLvlFactors() right after estimation setup. Degrades gracefully (the warning falls back to an internal #N index) when run against an rxode2 that predates the rxSetIdLvlFactors symbol.

• 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).

• Added new mu-referenced FOCEI-family estimation methods: mfocei, ifocei (FOCEI); mfoce, ifoce (FOCE); magq, iagq (adaptive Gauss-Hermite quadrature); mlaplace, ilaplace (Laplace). For any theta/eta that participates in a mu-ref covariate relationship (e.g. cl <- exp(tcl + eta.cl + allo.cl*logWT)), the covariate-coefficient theta(s) are excluded from the outer gradient optimizer entirely and instead re-derived directly by a closed-form OLS (mu* methods) or curvature-reweighted IRLS (irls* methods) regression of each subject’s back-calculated conditional value on the covariate(s); the regression residual becomes that subject’s eta. This runs natively in C++ inside the same inner-optimization pass every outer iteration already uses – no restart loop, no repeated model setup/compilation, no R-level round trip – and converges to comparable objective function values and parameter estimates as the corresponding standard method, typically in less time. New foceiControl() options control the mechanism: muModel ("none"/"lin"/"irls", default "none"), muRefCovAlg, muModelTol, muModelMaxCycles. Ordinary focei/foce/agq/ laplace/fo/foi default to muModel="none" and are completely unaffected; fo/foi are intentionally not given mu*/irls* counterparts (their first-order approximation has no per-subject conditional estimate to regress against).

• foceiControl() now defaults to outerOpt = "lbfgsb3c" and sigdig = 4

• Added mu-referenced FOCEI-family estimation methods: mfocei/ ifocei, mfoce/ifoce, magq/iagq, mlaplace/ilaplace, with new foceiControl() options muModel, muRefCovAlg, muModelTol, muModelMaxCycles

• Errors during estimation are now collected and reported together instead of only the last one (new collectErr argument to .collectWarn()).

• FOCEi no longer aborts R with Cube::slice(): index out of bounds when mceta >= 1 is combined with a pure evaluation that runs with maxInnerIterations == 0 (the covariance step, or nlmixr2extra::linearize()). The Monte-Carlo ETA-sample cube is filled only when maxInnerIterations > 0, so the inner-loop read indexed an empty cube and the resulting exception, thrown inside the OpenMP parallel region, was uncaught and aborted R. The read is now guarded (id < n_slices) and skipped when the cube holds no slice for the subject. instead of only the last one

• Fix Windows heap-corruption segfault building (focei, foce, fo, laplace, agq, bobyqa, nlm, optim, nls, nlminb, lbfgsb3c, n1qn1, newuoa, uobyqa) fits at more than one core. On Windows each package statically links its own OpenMP runtime, so when the parallel inner loop called rxode2’s solver across threads rxode2 saw every worker as thread 0 and collapsed its per-thread solve buffers onto a single slot, racing and corrupting the heap. The inner loop now hands rxode2 the real thread id via setRxThreadId() from rxode2 api (requires the matching rxode2).

• Fixed covMethod = "r" and covMethod = "s" standard errors, which were inflated by constant factors (sqrt(2) and 2, respectively). With the objective on the -2*logLik scale the R matrix is the observed information and the S matrix is the score cross-product, so the covariances are R^-1 and S^-1; they were returned as 2*R^-1 and 4*S^-1. The default sandwich method "r,s" (R^-1 S R^-1) was already correct and is unchanged. covR/covS are fixed at source so the sandwich-selection heuristic also compares consistently-scaled covariances (#666).

• The iteration-time progress output emitted by every estimator (focei, saem, bobyqa, nlm, optim, nls, nlminb, lbfgsb3c, n1qn1, newuoa, uobyqa) now flows through a single shared printer (scaleApplyIterPrintControl/scalePrintFun in src/scale.h). Each estimator’s iteration trace has the same #/U/X row layout, column wrapping, ANSI handling, periodic header re-emit cadence, and per-iteration user-interrupt check.

• Fix issue 641: FOCEI now updates additive mu-referenced population parameters with large-magnitude initial estimates

• Iteration-time progress output for all estimators now flows through a shared printer; new iterPrintControl() bundles the every, ncol, headerEvery, useColor, simple options

• focei (and the foce/fo/foi/posthoc family) again shows the Function Val. objective-function column in its iteration trace. The column had silently disappeared when the shared printer gained its showOfv flag, because focei sets up its scaling struct by hand and the flag defaulted to off — dropping the objective from every outer optimizer (including foceiControl(outerOpt = "bobyqa")) and misaligning the gradient (G/F/C/M) rows, whose method label lives in that same column slot.

• Restore the Function Val. objective column for focei/foce/fo/foi/ posthoc, which had dropped out once the shared printer gained showOfv

• Added sensMethod to the nlm-family controls (nlmControl(), nlminbControl(), optimControl(), n1qn1Control(), lbfgsb3cControl()) and to foceiControl() (focei/foce inner ETA sensitivities); ODE parameter sensitivities can be computed with the forward sensitivities ("forward") or the in-engine discrete adjoint ("adjoint"). When left at "default", the method is taken from the global option getOption("nlmixr2est.adjoint") (default "forward"), so the package-wide policy can be set in one place

• fo/foi now force forward sensitivities (adjoint does not apply to the eta=0 linearization), and the adjoint base-method restore in the focei family is a strict no-op for forward fits, fixing fo/foi tables/residuals

• Fix a fit aborting with initial 'omega' matrix inverse is non-positive definite when a degenerate fit (e.g. SAEM collapsing an uninformative random-effect variance to 0) leaves a singular omega; the sym-inv-chol setup now nearPD-corrects it so the residual/table diagnostics still run

• Fix SAEM erroring with No data with ID: <id> for a dosed subject with no usable observation; such subjects are now dropped before estimation and re-inserted into the output with a population PRED and NA individual columns, like FOCEi (#687)

• Internal consolidation of data preparation and the nlm-family control/fit functions across estimation methods; no change to any fit result

• Fix cov2cor error when omega has exactly one nonzero diagonal

• Fix SAEM linearized-FIM covariance (covMethod = "linFim") erroring when exactly one covariate-model parameter is estimated

• Fix FOCEi aborting R with Cube::slice(): index out of bounds when mceta >= 1 and maxInnerIterations == 0

• Fix Windows heap-corruption segfault for gradient/pooled estimator fits at more than one core

• Fix SAEM covariance error (rxInv(.tmp): Not a matrix) for models with a single population parameter

• Test suite uses a single testthat worker on CI/CRAN and parallel elsewhere; rxode2’s within-solve threads capped to 2 only on CRAN

• Restore the Function Val. objective column for focei/foce/fo/foi/ posthoc, which had dropped out once the shared printer gained showOfv

• Periodic header re-emits now repeat only the column labels, not the full Key: legend

• fit$time now reports every estimation stage consistently

• foceiControl() now defaults to outerOpt = "lbfgsb3c" and sigdig = 4

• Added mu-referenced FOCEI-family estimation methods: mfocei/ ifocei, mfoce/ifoce, magq/iagq, mlaplace/ilaplace, with new foceiControl() options muModel, muRefCovAlg, muModelTol, muModelMaxCycles

• Errors during estimation are now collected and reported together instead of only the last one

• Fix issue 641: FOCEI now updates additive mu-referenced population parameters with large-magnitude initial estimates

• Iteration-time progress output for all estimators now flows through a shared printer; new iterPrintControl() bundles the every, ncol, headerEvery, useColor, simple options

• 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

• 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

• Fix heap-buffer overflow and wrong back-transform in SAEM lambda (Box-Cox) residual-error models

• Fix heap-buffer overflow and wrong back-transform in SAEM lambda (Box-Cox) residual-error models

• Guard against null pointer arithmetic in inner.cpp

• Use OpenMP threading for S matrix calculation

• Use OpenMP threading wile calculating NPDEs

• Bounded the adaptive finite-difference step size (Shi et al. 2021) used for the FOCEi eta gradient to a reasonable region. When the finite differences carried no detectable curvature the step-size search could grow the step without bound and probe a parameter far outside the region where the local model holds, producing a degenerate ODE solve that corrupted the shared solver state and collapsed the eta finite-difference sensitivity on later inner iterations (the eta could get stuck near 0). The step is now clamped both above and below.

• saemControl(covMethod = "") no longer errors. "" is a documented choice that skips the covariance step, but saemControl() rejected it because match.arg() selects choices with pmatch(), and pmatch("") matches nothing; it is now handled explicitly.

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

CRAN release: 2026-06-03

• 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

CRAN release: 2026-05-31

• 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

CRAN release: 2025-10-09

• 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

CRAN release: 2025-08-29

• 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

CRAN release: 2025-07-24

• 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

CRAN release: 2025-07-19

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

nlmixr2est 4.0.0

CRAN release: 2025-07-15

• 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

CRAN release: 2025-02-18

• 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

CRAN release: 2025-01-18

• 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

CRAN release: 2024-11-23

• 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

CRAN release: 2024-10-22

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

CRAN release: 2024-09-18

• 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

CRAN release: 2024-05-28

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

CRAN release: 2024-01-31

• 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

CRAN release: 2023-12-12

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

CRAN release: 2023-10-08

• 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

CRAN release: 2023-09-18

• As requested by CRAN, remove Rvmmin

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

nlmixr2est 2.1.6

CRAN release: 2023-05-25

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

CRAN release: 2023-04-22

• 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

CRAN release: 2023-04-02

• 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

CRAN release: 2022-11-10

• 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

CRAN release: 2022-11-02

• 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

CRAN release: 2022-10-22

• 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

CRAN release: 2022-10-19

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

CRAN release: 2022-06-22

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

CRAN release: 2022-05-23

• 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

CRAN release: 2022-05-12

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.