--- title: "Modular numerical controls for sommer::mmes" author: "sommer development team" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{sommer.mmes.options} %\VignetteEngine{knitr::knitr} %\VignetteEncoding{UTF-8} --- # Purpose For a clickable decision tree with searchable controls, defaults, examples, and compatibility checks, open the [interactive mmes_control guide](https://covaruber.github.io/sommer/mmes-options.html). `mmes()` separates statistical model specification from numerical policy. Formulas, data, weights, family, constraints, and known covariance parameters remain `mmes()` arguments. Numerical choices belong in one validated object: ```r fit <- mmes( fixed = trait ~ environment, random = ~ genotype, rcov = ~ units, data = trial_data, control = mmes_control() ) ``` This vignette documents the modular API, defaults, mathematical ownership, eligibility constraints, and migration from the removed numerical formals. It is explanation only: code is displayed, not run. Building this vignette does not fit models, benchmark methods, or change hardware settings. # Public interfaces The current public numerical interface is: ```r mmes(fixed, random, rcov, data, W, weights = NULL, naMethodX = "exclude", naMethodY = "exclude", naMethodRandom = "exclude", naMethodR = "exclude", contrasts = NULL, family = stats::gaussian(), vcc = NULL, covPar = NULL, parameters = list(), curves = list(), start = NULL, method = mmes_method(), control = mmes_control(), familyParameters = list()) mmes_method(objective = "reml", integration = "auto", approximation = "auto", task = "fit") mmes_control( formulation = "henderson", optimizer = list(), mode = list(), working = list(), likelihood = list(), population = list(), inference = list(), representation = list(), solve = list(), logdet = list(), score = list(), curvature = list(), update = list(), proposal = list(), convergence = list(), output = list(), execution = list(), diagnostics = list() ) ``` `fixed`, `random`, `rcov`, `data`, and `W` remain model inputs. `weights` declares independent blocks in a supplied `W`; it does not create weights. `family`, `vcc`, `covPar`, missing-data policies, and `contrasts` retain their statistical or data-preparation roles. `parameters` declares unknown quantities with `parameter()`; `curves` declares named native curves referenced by `curve_term("name")` in the formula. `start` supplies natural-scale model starts; `familyParameters` configures fixed/estimated family nuisance parameters. These are model inputs, not numerical controls. `mmes_method()` declares objective, integration, approximation, and fit/solve task. Its defaults preserve Gaussian REML and historical PQL for ordinary non-Gaussian mixed models, with REML Gaussian working fits. This is **not** an exact non-Gaussian REML likelihood. Explicit native Laplace requires `mmes_method(objective="ml", integration="laplace")` and an eligible model. Non-Gaussian mixed exact integration, PQL with explicit exact/Laplace integration, and nonlinear PQL are rejected rather than silently substituted. `mmes_control()` validates top-level and nested module names. Unknown fields, duplicate names, and malformed named lists error. Removed options are not accepted through `...`, aliases, or a legacy `engine` argument. Controls record requests; they do not fit a model or establish model-specific eligibility. Fitting validation remains authoritative. # Traditional nonlinear fixed means `curves=list(curve=nonlinear_curve("time", "exponential"))`, referenced as `curve_term("curve")` in the formula, adds the fixed mean `asymptote * exp(-rate * time)` to ordinary formula contributions. `"logistic"` and `"saturation"` provide logistic growth and Michaelis-Menten saturation, respectively. These are traditional parametric curves. ```r fit <- mmes( fixed = y ~ 0 + curve_term("curve"), random = ~ group, data = dat, curves = list(curve = nonlinear_curve("time", "exponential")), start = list(parameters = c(asymptote = 4, rate = 0.6)), method = mmes_method(objective = "ml"), control = mmes_control( optimizer = list(algorithm = "nlminb", gradient = "analytic"), convergence = list(maxIters = 500L, likelihoodTol = 1e-8, parameterTol = 1e-6), output = list(inverse = 0L) ) ) coef(fit) predict(fit, newdata = dat) predict(fit, newdata = data.frame(time = c(1, 2)), re.form = NA) ``` The curve parameters are fixed; random terms remain additive Gaussian effects. Gaussian integration is exact. The residual variance and random covariance parameters are estimated jointly with the fixed mean. Positive curve parameters use log coordinates; starts and coefficients use natural coordinates. Automatic starts are deterministic, with no RNG use or automatic input scaling. Nonlinear input missingness follows `naMethodX` and shares the ordinary observation mask. The exact Gaussian path requires univariate Gaussian identity ML, Henderson, at least one random term, and independent homogeneous residuals. Unsupported controls error rather than selecting a different model. REML, direct/PQL fitting, structured residuals, and weights are not supported on this route. Fixed-only curves use the reference specification route. Multiple independent or jointly correlated random curve parameters are supported as described below. All existing native response families are also available, including ordinal, nominal, zero-inflated and `familym()` mixtures. These use complete analytic Laplace ML. The traditional curve contributes to the **link predictor**; it is shared across traits and every nonreference nominal logit. Family nuisance choices use the `mmes(familyParameters=)` argument; outer/mode/quadrature choices use `optimizer`, `mode`, and `population`. For example, `family=nbinom2(theta=3)` with `method=mmes_method(objective="ml"), familyParameters=list(theta=3, estimate=FALSE)` fixes the NB2 shape. A zero-component regression uses `family=zipoisson(zero=~treatment)`. For these families `type="response"` returns response means or category probabilities and `type="link"` returns predictors. Categorical mixtures return trait-named lists. Population prediction integrates additive and curve deviations jointly, including covariance between nominal logits; it is not inverse-link at zero random effects. Accuracy and workspace failures are explicit. `simulate(re.form="conditional")` draws observations at fitted modes; `simulate(re.form="population")` draws fresh prior random effects before evaluating the nonlinear predictor. A supplied `seed` is local to simulation. Neither operation computes uncertainty. Fitting computes no SE, PEV, parameter covariance, or outer observed Hessian. Uncertainty remains strictly post-fit: ```{r, eval=FALSE} post <- post_nonlinear(fit, mode=2) vcov(post) # Natural fixed/curve covariance vcov(post, type="conditional") # Conditional mode covariance vcov(post, type="random") # Empirical Bayes approximate PEV confint(post) # Working-normal, positive curves transformed predict(post, newdata=dat, interval="confidence", cov.fit=TRUE) predict(post, newdata=dat, re.form="population", interval="prediction") ``` `postPEV()` dispatches to the same nonlinear implementation. Conditional mode covariance and estimated-parameter approximate PEV are distinct; neither is an exact unconditional frequentist PEV. Numerical choices use `inference` in `mmes_control()`. Prediction intervals are normal approximations for scalar response-scale outputs, while categorical probabilities support confidence intervals. `summary()` displays saved uncertainty without computing it; `summary(inference=TRUE)` requests it. `AIC()`/`BIC()` count all free coordinates. Model comparisons require regular nested ordinary fixed terms with the same nonlinear architecture; automatic nonlinear or variance-boundary tests are rejected. New conditional predictions must use represented random levels and the fitted covariance basis. For these additive Gaussian models, the population mean equals the zero-random-effect mean. See `?nonlinear_curve` for curve equations, parameter names, restrictions, and identification diagnostics. Nonlinear scalability can be inspected without changing the model: ```{r, eval=FALSE} ctrl <- mmes_control(diagnostics=list(profile="force")) fit <- mmes(y~0+curve_term("curve"), random=~group, data=dat, curves=list(curve=nonlinear_curve("time")), method=mmes_method(objective="ml"), control=ctrl) profile <- profile_nonlinear(fit, replications=3L) profile$table fit$executionPlan$dimensions fit$executionPlan$memory ``` Preparation, outer optimization, mode work, initialization and convergence diagnostics are recorded separately. Kernel times are inclusive, overlapping, and not additive. Memory reports retained R object sizes and the declared admission policy, not native peak RSS. The nonlinear sparse workspace guard bounds factor fill by connected random components before Hessian assembly: independent groups retain linear admission, whereas a large connected relationship graph may be rejected conservatively. Prediction and post-fit inference keep their separate memory checks. The reproducible installed benchmark is `system.file("benchmarks", "nonlinear-qualification.R", package="sommer")`; it varies observations and groups across exact Gaussian and Gaussian/Poisson random-curve paths. Timings are descriptive, never portable performance guarantees. ## Random curve parameters (Laplace ML) ```r model <- nonlinear_curve("time", "exponential", random = list(rate = ~ vsm(ism(subject), sigma2 = 0.04))) fit <- mmes(y ~ 0 + curve_term("curve"), data = dat, curves = list(curve = model), method = mmes_method(objective = "ml"), start = list(parameters = c(asymptote = 4, rate = 0.6)), control = mmes_control( mode = list(maxit = 100L, tol = 1e-8, stepTol = 1e-8, warmStart = "off", starts = "deterministic"), population = list(order = 20L, maxOrder = 160L, tol = 1e-7), convergence = list(maxIters = 500L), diagnostics = list(verbose = FALSE))) ``` Here `log(rate_j) = log(rate) + b_j`. Any one parameter of the three traditional curves can vary; positive parameters use log coordinates, and logistic midpoint uses identity coordinates. An optional ordinary `random=~subject` adds independent additive effects. Supported covariance descriptors and known relationship precision retain their existing conventions. One unit-incidence curve coordinate per row is required. Exact zero random variance is unsupported; with additive effects retained, remove its declaration for the zero-curve-variance model. The mode solver uses exact conditional curvature, safeguards Newton steps, and compares deterministic starts. Laplace scores include mode motion and third curve derivatives. No uncertainty or PEV calculations are introduced. Mode diagnostics report failures and competing candidates; they do not certify a global optimum. Conditional predictions use fitted curve modes (`re.form=NULL`); `re.form=NA` uses shared, zero-deviation parameters. `re.form="population"` integrates against the fitted prior variance with increasing-order Gaussian quadrature, and errors if the requested agreement is not reached. This is a population mean, not an uncertainty calculation. New-data conditional/population predictions currently require represented groups and the fitted covariance basis. The fitted object retains `groupParameters` and `curveRandomMode`. ## Structured and reduced-rank random covariance The nonlinear paths have dedicated numerical checks for `ism()`, `dsm()`, `usm()`, `ar1m()`, `csm()`, `fam()` (heterogeneous or homogeneous specifics), and representative Kronecker combinations. `fam()` requires positive specific variances. Known positive-definite relationship precision keeps the usual `Gu` convention. Other constructors are eligible only when their evaluator and derivative contracts match; eligibility alone is not comprehensive validation. Pure `rrm()` factors are supported in both ordinary additive effects and random curve parameters: ```r model <- nonlinear_curve("time", "exponential", random = list(rate = ~ vsm(rrm(environment, k = 1), ism(subject)))) ``` Reduced-rank effects use `b = L v`, with a nonsingular Gaussian latent prior. The existing normalized loading convention is retained; there is no singular inverse, added ridge, or artificial specific variance. Loading-dependent designs and their analytic derivatives enter the exact Gaussian or full Laplace score. The declared curve incidence is still one unit per observation; weighted latent loadings are internal, not a new random-slope API. `randomMode`, `curveRandomMode` and `groupParameters` use restored original coordinates. `latentRandomMode` and `randomCoordinates` retain the latent point mode and dimensions. Predictions rebuild the saved basis; population means use the implied marginal prior variance. Neither latent modes nor these population integrals are PEV or parameter uncertainty. Existing post-fit inference guards remain in force. No new numerical controls are needed for this parameterization. ## Multiple random curve parameters Existing named formulas and all `vsm()` arguments keep their syntax and meaning. Separate named formulas specify independent Gaussian parameter blocks: ```r model <- nonlinear_curve("time", random=list( asymptote=~vsm(ism(subject), sigma2=0.04), rate=~vsm(ism(subject), sigma2=0.01))) ``` For correlated deviations, use the separate `joint` argument: ```r model <- nonlinear_curve("time", joint=list( parameters=c("asymptote","rate"), random=~vsm(usm(.parameter), ism(subject), sigma2=0.04))) ``` The internal `.parameter` axis follows `joint$parameters` in that exact order. It is a covariance-shaping factor, not a user data column or the final grouping term. `usm()` models cross-parameter covariance; `dsm()`, `fam()`, `rrm()` and other eligible native descriptors retain their existing semantics. Known relationships remain `Gu` precisions. No response replication or user-data reshaping is needed. A joint block can coexist with independent formulas for other parameters, but the parameter names cannot overlap. Ordinary additive effects remain separate independent `mmes(random=~...)` terms; their cross-covariance with curve parameters is not supported. Reports qualify covariance terms by `curve[parameter]` or the joint parameter names, avoiding ambiguous names when grouping formulas coincide. `vcc` retains its existing indexed/named constraint semantics. Positive parameters vary on their log coordinates; logistic midpoint varies on its identity coordinate. All cross-parameter second and third derivatives enter exact conditional curvature and complete analytic Laplace scores. Population predictions integrate the *joint prior*, rather than integrating each parameter separately. Quadrature cost grows with its effective rank; existing `population` and `execution$denseMemoryMB` controls apply and failure is explicit. Conditional prediction, generative simulation and serialization retain the joint basis. PEV, group-parameter standard errors and prediction uncertainty remain explicitly post-fit. `vcov(fit, type="group")` explicitly requests natural-scale group-parameter covariance, including cross-parameter terms; it is not computed during fitting. ## Curve and parameter declarations These are distinct random-effect locations: * `mmes(random=~group)` declares ordinary additive response-predictor effects. * `parameters=list(rate=parameter(fixed=~treatment, random=~group, link="log"))` declares a regression for an unknown parameter used in a nonlinear formula. * A native `nonlinear_curve(random=list(rate=~group))` declares curve-coordinate deviations; `joint` supplies an admitted correlated block. `optimizer$gradient="auto"` preserves native analytic scores and selects the reported numerical method for reference expression models. An explicit `"analytic"` request is rejected for reference expressions rather than silently replaced by numerical derivatives. Native curves reject explicit `"central"`. Native curves use the public formula bridge: ```r fit <- mmes(y ~ 0 + curve_term("component"), random = ~ group, data = dat, curves = list(component = model), start = list(parameters = natural_starts), method = mmes_method(objective = "ml"), control = mmes_control(mode = list(maxit = 100L, tol = 1e-8))) ``` Declared `parameter()` regression models have their own supported family, expression, and numerical workspace limits (30 declarations, 60 mean coefficients, 512 modes for nonexact integration). The dense reference route supports nuisance `familyParameters`, checked inference controls, simulation, and bounded joint predictor quadrature. It rejects RR covariance, weights, `W`, `vcc`, custom contrasts and nondefault representation/solver/update controls. Ordinary Gaussian models retain `representation=list(fixedRank="drop")`; this is not a license to drop unidentified nonlinear weights. Consult the model specification vignette rather than assuming every native curve layout is interchangeable. Unsupported combinations error; they do not enable a different estimator. # Ownership and defaults The default control retains the established numerical bundle: Henderson formulation, REML fitting, exact AI curvature, EM recovery, no extrapolation, 30 inner iterations, full covariance information, coefficient-matrix retention, and inherited solver/resource defaults. These defaults preserve the prior numerical bundle; they do not authorize new combinations. | Module | Mathematical ownership | Fields and defaults | |:--|:--|:--| | `formulation` | Equation domain | `"henderson"`; `"henderson"` or `"direct"`. | | `method` (not a control module) | Objective, integration, approximation and task | `mmes_method(objective="reml", integration="auto", approximation="auto", task="fit")`. Ordinary non-Gaussian mixed auto fits preserve PQL/REML working fits. Laplace requires ML. | | `optimizer` | Outer optimizer policy | `maxit=NULL`, `gradient="auto"`, `algorithm="nlminb"`, `scoreTol=NULL`; gradient choices are `"auto"`, `"analytic"`, `"central"`. Auto preserves native analytic defaults and selects reported numerical derivatives for reference expressions; unsupported explicit requests error. `scoreTol` sets the nonlinear specification outer-score threshold; NULL uses 0.001. Gaussian AI/update modules retain their role. Do not duplicate conflicting `convergence$maxIters`. | | `mode` | Conditional random-effect mode | `maxit=100L`, `tol=1e-8`, `stepTol=1e-8`, `warmStart="auto"`, `starts="deterministic"`. Auto preserves native GLMM auto reuse and deterministic cold nonlinear modes; nonlinear explicit `"on"` errors. Non-default `stepTol` is supported only by nonlinear modes; native GLMM modes use `tol`. | | `working` | Outer PQL working fits | `maxit=20L`, `tol=1e-5`, `estimateTheta=FALSE`, `warmStart=TRUE`. `convergence` controls inner Gaussian fits. | | `likelihood` | Native likelihood algebra | `backend="auto"`, `trace="panels"`, `blocks="off"`; backend selects random/observation algebra, not the Gaussian solver. Conditional-mode reuse belongs only to `mode$warmStart`. | | `population` | Nonlinear population quadrature | `order=20L`, `maxOrder=160L`, `tol=1e-7`; failure to meet accuracy/workspace limits errors. | | `inference` | Strictly post-fit nonlinear uncertainty | Native inference tolerances and conditioning safeguards; see the generated catalog for constructor defaults. | | `representation` | Equation/covariance representation | `compression="off"`, `factorScores="none"`, `fixedRank="error"`, `equations="normal"`, `sparseFactors="auto"`, `gridAugment="auto"`, `matrixFree="auto"`; modes can inherit global options, and `equations` defaults to `"saddle-point"` if that mode is globally forced. | | `solve` | Linear solves, preconditioning, reuse, elimination | `method="auto"`; PCG fields below; `blockSchur`, `blockBanded`, `factorSchurMultiFactor`, and `dFirstML` default to `"auto"`. | | `logdet` | Determinant approximation | `scaling="none"`, `lanczosSteps=20L`. | | `score` | Covariance derivatives and trace estimation | `traceProbes=8L`, `estimator="hutchinson"`, `rank=16L`, `covarianceChecks="auto"`, `traceBatch="off"`. | | `curvature` | AI metric and reduced-direction assembly | `method="exact-ai"`, `assembly="sensitivity"`, `reduced="off"`. | | `update` | Covariance-coordinate update and recovery | `method="ai"`, `recovery="em"`, `blendWeight=NULL`, `commonScale="off"`. | | `proposal` | Proposal scaling and extrapolation | `stepWeight=NULL`, `extrapolation="none"`. | | `convergence` | Inner iterations and stochastic accuracy | `maxIters=30L`, `likelihoodTol=1e-4`, `parameterTol=1e-4`, `accuracy="fixed"`, `maxTraceProbes=NULL`, `maxLanczosSteps=NULL`, `logDetTol=1e-3`, `scoreTol=1e-3`. | | `output` | Coefficient and covariance uncertainty; setup output | `inverse=0L`, `covariance="full"`, `coefficientMatrix=TRUE`, `setup=FALSE`, `selectedInverse="auto"`, `panelSize=64L`. | | `execution` | Memory, hardware, threads, precision, seed | `denseMemoryMB=1000`, `memoryPolicy="speed"`, `gpu=FALSE`, `gpuMinSize=256L`, `blasThreads=NULL`, `inverseTol=1e-6`, `seed=0L`. | | `diagnostics` | Messages, package-age warning, profiling | `verbose=TRUE`, `dateWarning=TRUE`, `profile="off"`. | The PCG fields under `solve` default to `tol=1e-8`, `maxIters=0L`, `preconditioner="diagonal"`, `nystromRank=32L`, `warmStart=FALSE`, `recycleRank=0L`, `recycleMethod="coarse"`, `recycleState=NULL`, `exportRecycle=FALSE`, `rhsMethod="independent"`, `blockSize=16L`, and `residualRefresh=0L`. `maxIters=0L` selects the solver's automatic limit. PCG-specific requests matter only when an eligible PCG path uses them. `output$inverse` controls final coefficient-inverse work and fit-time PEVs. Mode `0L` computes neither. Mode `1L` computes random-effect PEVs from the selected inverse without forming full `Ci`. Mode `2L` forms full `Ci` and computes PEVs. To defer this work, retain the default mode `0L` and `output$coefficientMatrix=TRUE`, then call `postPEV(fit, mode=1L)` for PEVs only or `postPEV(fit, mode=2L)` for PEVs and full `Ci`. Native Laplace GLMMs require `output$inverse=0L`; setting it to `1L` or `2L` is an error because all GLMM PEV work is post-fit. Use `postPEV()` after fitting. Laplace fixed-effect and covariance-parameter standard errors are likewise calculated on demand by `vcov()` or `summary()`. For Laplace prediction uncertainty, `predict(..., se.fit=TRUE)` uses bounded solve panels; `cov.fit=TRUE` requests the full covariance. Both use `execution$denseMemoryMB` and `output$panelSize`, and the default simulation memory budget also uses `execution$denseMemoryMB`. Options may supply defaults for experimental modes, memory budget/policy, and GPU settings. For example, a non-`NULL` `sommer.mme.denseMemoryMB` option replaces the `1000` MB default. Explicit control values override option-derived mode defaults. Use them for reproducible intent. Hardware and library availability can still affect execution. Constructors and valid values: * `solve$method`: `"auto"`, `"ldlt"`, `"cholmod"`, or `"pcg"`. * `representation$compression`: `"off"`, `"groups"`, or `"classes"`. * `representation$factorScores`: `"none"`, `"fixed-shape"`, or `"profile"`. * `representation$fixedRank`: `"error"` (default) or opt-in `"drop"`. * `representation$equations`: `"normal"`, `"saddle-point"` or experimental `"conditional"`. * `logdet$scaling`: `"none"` or `"diagonal"`. * `score$estimator`: `"hutchinson"` or `"control-variate"`. * `curvature$method`: `"exact-ai"` or `"hybrid-ai"`. * `curvature$assembly`: `"sensitivity"` (default) or `"augmented-working"`. * `update$method`: `"ai"`, conditional `"em"`, or the earlier `"conditional-fa"` EM selector. * `update$recovery`: `"em"` or `"pxem"`. * `proposal$extrapolation`: `"none"` or `"aitken"`. * `convergence$accuracy`: `"fixed"` or `"adaptive"`. * `output$inverse`: `0`, `1`, or `2`; `output$covariance`: `"full"` or `"none"`. * `execution$memoryPolicy`: `"speed"` or `"limit"`. * Experimental mode fields accept `"auto"`, `"off"`, `"force"`, or logical `TRUE`/`FALSE`. A forced request still obeys mathematical, structural, and memory checks. The complete field reference and mapping from old names are in `?mmes_control`. # Decision graph and constraints ## Coefficient-space augmented working-variate AI ```{r, eval=FALSE} control <- mmes_control( solve = list(method = "cholmod"), curvature = list(assembly = "augmented-working"), diagnostics = list(profile = "force") ) fit <- mmes(y ~ environment, random = ~ id, data = dat, control = control) fit$engineDiagnostics$aiAssembly fit$executionPlan$admission$aiAssembly ``` Both AI assemblies process directions together through the existing joint coefficient solver. The augmented route uses a prior-shifted border and its Schur complement, not a second factorization of enlarged equations. Random working directions remain in coefficient space: no observation-by-random- parameter matrix or observation covariance is formed. Structured residual derivatives retain their existing diagonal, block, grid or sparse operators. Compression retains within-group residual variation. The option changes assembly, not the likelihood, score, optimizer or coefficient/inverse output coordinates. It supports deterministic Gaussian Henderson REML/ML, including admitted specialized coefficient solve paths, weights, positive-specific FA, tied/reduced directions, EM/PXEM recovery, hybrid exact refreshes and common-scale profiling. Existing restrictions on combinations of those methods continue to apply. PCG, PQL/Laplace, conditional FA, factor-score augmentation and zero-specific latent reduced rank are explicitly rejected. The option is preserved in stored controls and refits. Existing residual-engine limitations remain: weighted incomplete serial blocks can fail selected-inverse demand closure with either AI assembly. Choosing augmented assembly does not bypass that trace preparation. Workspace admission uses `execution$denseMemoryMB`. Multi-column tiles, where required, preserve every off-diagonal AI block. No approximate diagonal AI or silent sensitivity substitution is used. Numerical products are rebuilt for each coefficient/covariance evaluation; structural maps and scratch allocations are reused where valid. Neither formulation is universally faster or more stable. Use matched backends, starts, tolerances and output settings when comparing them; the benchmark runner reports separate frozen-parameter AI parity as well as complete-fit timings. ## Rank-deficient fixed effects By default, `mmes()` stops with diagnostics when the fixed design contains zero columns or linear dependencies. For Gaussian identity-link fits, both the Henderson and direct engines can instead use a shared sparse C++ QR basis selected once after observation filtering: ```{r, eval=FALSE} fit <- mmes( response ~ locale + cell, data = observations, control = mmes_control(representation = list(fixedRank = "drop")) ) fit$fixedRank[c("rank", "retained", "omitted")] summary(fit)$betas ``` The original coefficient layout is retained. Omitted coefficients are constrained to zero and explicitly marked aliased, with no inferential standard errors or tests. Retained coefficients are conditional on those constraints. This preserves the fitted-mean space but need not reproduce ASReml's particular coefficient constraints. Different fixed bases can shift REML likelihoods by a constant. `fit$fixedRank$relation` is a sparse map from the retained design to the original columns. Direct-engine `VarBeta` stays in retained coordinates; it is not expanded to a potentially huge original-column covariance. Henderson prediction contrasts are checked against the original model's estimability relations. Non-estimable rows are rejected. Wald tests, likelihood-ratio comparisons, Satterthwaite/Kenward-Roger inference, and `emmeans` require an explicitly full-rank refit. Existing direct-engine prediction limitations are unchanged. Character classify predictions for rank-reduced fits stop before a Cartesian level grid is constructed. To explicitly request observed combinations: ```{r, eval=FALSE} cells <- predict(fit, D = "locale:cell", observed = TRUE) cells$observed ``` This path averages the fixed predictor over retained records within each observed group and excludes random contributions by default. It is not an equally weighted balanced-grid marginal mean. Random contributions can be requested explicitly through `Dtable`, without changing its fixed rows: ```{r, eval=FALSE} rules <- fit$Dtable rules$include[rules$type == "random"] <- TRUE conditional <- predict(fit, D = "locale:cell", observed = TRUE, Dtable = rules) ``` For random terms, include alone uses empirical within-group design averages; include plus average uses equally weighted co-occurring indicator levels; average alone uses equally weighted selected fitted coefficient levels. Random level restrictions are given through entries of `levels` named by the fitted random term. The returned `Dtable` contains only term-level rules; selections are returned separately in `levels` and can be reused with those rules. A `levels` entry named exactly as `D` may select and order observed labels. Unobserved requests produce an error. Numeric classify variables and custom fixed-effect rules require explicit estimable contrasts instead. See `?predict.mmes` for random-level filtering and weighting details. Observed combinations can still produce a large dense covariance. `maxMemoryMB` defaults to the fitted dense memory budget (1000 MB otherwise) and checks estimated prediction workspace before constructing grouped contrasts or solving. Restrict `levels` or use smaller explicit contrast batches if necessary. `PEV = FALSE` does not remove the mandatory `vcov`. This option is not supported with PQL, solve mode, conditional FA, factor scores, observation rotation or compression; these combinations fail explicitly rather than silently ignore the rank policy. ## Conditional FA: reusable posterior moments with AI or EM The default ordinary AI workflow is unchanged. The opt-in `representation$equations="conditional"` selects sparse latent FA equations using the existing formula preparation, observation map, covariance descriptors, model methods, controls, and execution-plan reporting. The optimizer is a separate block: `update$method="ai"` (the default) uses exact REML scores and average-information curvature on these equations; `update$method="em"` uses conditional-moment EM. The earlier `update$method="conditional-fa"` API retains its EM behavior even without an explicit representation request. ```r long <- stackTraits(wide, traits = trait_names, keep = "id") Ainv <- solve(A) attr(Ainv, "inverse") <- TRUE fit <- mmes( value ~ trait, random = ~ vsm(fam(trait, k = 2), ism(id), Gu = Ainv), rcov = ~ vsm(fam(trait, k = 1), ism(record)), data = long, control = mmes_control( representation = list(equations = "conditional"), update = list(method = "ai"), convergence = list(maxIters = 300), output = list(covariance = "none", coefficientMatrix = FALSE, inverse = 0, panelSize = 32) ) ) fit$convergence fit$conditionalFA$likelihood fit$conditionalFA$updates # accepted AI steps and any EM recovery fit$conditionalFA$ai # final score norm, curvature rank and workspace fit$executionPlan$admission$conditionalFA covmatrix_mmes(fit, term = 1, se = FALSE) ``` Missing responses are excluded by the usual observation map, not imputed. Repeated individuals and incomplete residual pairing blocks are supported; `record` must identify the actual experimental unit. The genetic incidence must select exactly one individual and one trait per observed value. The fixed design may contain covariates and interactions but must have full column rank. Each trait must retain observations. Known positive-definite relationship **precision** matrices retain the usual `Gu` contract. Each covariance source must be entirely free or entirely fixed. FA specifics can be heterogeneous or homogeneous; diagonal and homogeneous residuals are also supported. ### Mathematical blocks Write the latent model as $$ y=Xb+Z_g(L_g f_g+s_g)+L_e f_e+\epsilon,\quad \operatorname{Var}(s_g)=\Psi_g\otimes A,\quad \operatorname{Var}(\epsilon)=\operatorname{diag}(\psi_{e,t(i)}). $$ Residual factors are indexed by pairing key, genetic factors by individual. The prepared sparse design is denoted by $D$ and its random-effect precision by $Q$ (with zeros on fixed-effect coordinates). At every iteration: 1. Assemble $H=D^\mathsf{T}R_0^{-1}D+Q$ and factor it numerically. 2. Compute conditional means and the entries of $H^{-1}$ needed by the complete-data quadratics, including fixed-effect uncertainty. 3. Select the update block: exact REML score/AI contractions, or loading regressions and expected quadratic forms for EM. Do **not** treat predicted factors as observed data. 4. Evaluate the observed restricted log likelihood and reject a material decrease. EM convergence requires likelihood and parameter stabilization; AI requires a small diagonally scaled REML score. On the original response scale, omitting the Gaussian normalizing constant, the restricted objective uses $$ -2\ell_R=\log|R_0|-\log|Q_{\mathrm{random}}|+\log|H| +(y-D\hat a)^\mathsf{T}R_0^{-1}(y-D\hat a) +\hat a^\mathsf{T}Q\hat a . $$ Integrating the fixed coefficients with a flat measure gives the restricted objective; their conditional second moments are essential to the EM update. Specific variances have a machine-precision positive floor to keep the conditional distribution well-defined during EM. AI uses log-specific coordinates and, for shared factors, logit allocations. ### Exact score and AI blocks For covariance coordinate $\theta_i$, let $V_i=\partial V/\partial\theta_i$, $P=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}$, and $z=Py$. The score and average information are $$ s_i=\tfrac12(z^\mathsf{T}V_i z-\operatorname{tr}(PV_i)),\qquad \mathcal A_{ij}=\tfrac12(V_i z)^\mathsf{T}P(V_j z). $$ The implementation differentiates the sparse latent design, residual diagonal and random precision. Selected posterior moments give the exact score traces; relationship-factor solves apply random prior covariances without forming a dense relationship inverse. The already factored $H$ projects derivative right-hand sides in bounded panels. Neither $V$ nor $P$ is constructed. Shared loading derivatives affect both covariance sources in the same coordinate. Independent loading matrices retain rotationally redundant coordinates. A diagonally scaled, rank-revealing curvature eigensolve drops numerically null directions rather than inverting a singular information matrix. Steps have a bounded coordinate displacement and a backtracking observed-likelihood Armijo test. If no AI proposal is accepted, an EM recovery step uses the same posterior moments; `conditionalFA$updates` records `"ai"` or `"em-recovery"` explicitly. This changes optimization, not the REML objective. `conditionalFA$ai$scoreNorm` is the diagonally scaled score checked against `convergence$parameterTol`; a small displacement alone cannot certify AI convergence. Near-zero conditional noise can make the latent equations ill-conditioned even when the marginal covariance remains well-conditioned. Conditional AI therefore has explicit numerical bounds: writing $\tau=\max(10^{-8},\texttt{execution\$inverseTol})$, free physical specifics must be at least $\tau\operatorname{var}(y)$ and shared allocations lie in $[\tau,1-\tau]$ (`inverseTol` must be below 0.5). Initial values must satisfy these bounds. An active-set solve excludes outward directions at a bound; convergence uses the projected KKT score, not an unconstrained zero-gradient claim. The actual `varianceFloor`, `allocationFloor`, and `activeBounds` are reported in `conditionalFA$ai`. These physical-coordinate constraints are not identical to the ordinary engine's descriptor-coordinate bounds. EM recovery stays within the AI feasible set. Explicit conditional EM retains its earlier machine-precision floor, so optima at numerical boundaries need not match exactly between optimizers. `conditionalFA$rawLikelihood` retains this original-scale objective. For compatibility with the ordinary solver's standardized-response reporting, `llik` and `conditionalFA$likelihood` add $(n-\operatorname{ncol}(X))\log(\operatorname{sd}(y))$. This constant does not affect updates or convergence. `AIC` and `BIC` retain the existing solver's fixed-design-column penalty convention; the covariance-model parameter count is recorded separately rather than silently changing those established fields. ### Reuse and memory policy The native conditional-moment block reuses symbolic LDLT analysis only while the sparse pattern is identical. Numeric factors are recomputed at every covariance update. Fixed sparse quadratic contractions are compiled once; identical factor incidence/missingness patterns share their Gram contractions. Selected inverse entries are obtained in bounded panels. No dense observation covariance or full posterior inverse is constructed. EM does not assemble an AI matrix; conditional AI does assemble a dense parameter-by-parameter matrix and observation-by-parameter derivative RHS. `output$panelSize` controls inverse/projection-panel width; `execution$denseMemoryMB` admits their dense workspace and, for AI, the retained RHS/curvature/eigensolve workspace (not a process-memory cap). Oversized AI requests fail explicitly rather than switching optimizers. Sparse factor fill, latent genetic specifics, and the number of exact inverse solves can still be expensive, especially with dense relationships. This first engine prioritizes exactness and bounded dense intermediates; it is **not a claim of MegaLMM-like thousand-trait runtime**. EM can converge slowly near boundaries. Compare converged objectives and predictions, not equal iteration counts. ### Shared-loading covariance model For common trait axes, use matching tags: ```r random = ~ vsm(fam(trait, k = 2, shared = "axes"), ism(id), Gu = Ainv) rcov = ~ vsm(fam(trait, k = 2, shared = "axes"), ism(record)) ``` This changes the covariance family to $G=L D_g L^\mathsf{T}+\Psi_g$ and $R=L D_e L^\mathsf{T}+\Psi_e$. Both sources must be free with equal rank. Diagonal factor variances are updated from conditional quadratics and normalized so $D_g+D_e=I$; the compensating scale is absorbed into $L$. Starting loading columns must be positive scalar multiples between sources. Defaults satisfy this. The common basis is in `fit$conditionalFA$state$loading`, with allocations `weightG` and `weightE`; sign and ordering are not uniquely identified. `loadings_mmes()` instead reports triangular marginal loadings for each source. Other engines reject shared tags rather than ignoring the restriction. ### Output and compatibility Fits remain `mmes` objects. `fitted()`, `residuals()`, `coef()`, covariance parameter reports, and explicit `covmatrix_mmes(..., se=FALSE)` remain available. `theta` has named empty entries rather than dense trait matrices. The covariance descriptors and conditional state are the authoritative factorized covariance representation. Covariance-parameter SEs, retained public `C`, coefficient PEVs, and uncertainty-dependent methods are not available on this path; they fail explicitly. No implicit handoff to another estimator is made. For shared models, the usual marginal parameter reports are redundant representations; the independent parameter count is reported in `fit$conditionalFA$parameterCount`. The public `llik` retains its usual one-row history shape. `monitor` contains the final marginal covariance parameters only, not every intermediate EM state; the likelihood history is retained in `conditionalFA$likelihood`. Weights, non-Gaussian families, ML, multiple random covariance sources, rotation, user `vcc`, random-regression/overlapping incidence, `rrm()`, spatial residuals, and combinations with AI acceleration/augmentation are not admitted. Such models continue to use the unchanged ordinary engine unless the conditional engine is explicitly requested. ### Validation and measured limits The posterior kernel is checked against dense means, inverse entries and log determinants, including symbolic reuse/invalidation. End-to-end tests compare incomplete-record REML, GLS and genetic BLUPs with dense reference equations, check shared and independent EM monotonicity, and verify convergence against a separately optimized one-dimensional restricted likelihood with a fixed genetic covariance. Unequal genetic/residual ranks and homogeneous specifics have separate coverage. Conditional AI scores and curvature are also checked against independently finite-differenced marginal covariance matrices, including rank-two shared loadings and homogeneous specifics. Convergence is compared with a separately optimized restricted likelihood; numerical-bound KKT convergence, EM recovery and legacy-EM parity have separate tests. Baseline-versus-candidate checks used the enhancer split-plot, compound-symmetry and diagonal GxE, multivariate, random-regression, field row/column, and genomic relationship models. Covariance parameters, BLUPs, likelihood histories and convergence flags matched exactly for all seven. This is regression evidence, not an independent certification of a model optimum. The earlier conditional EM selector and the explicit conditional EM control also produce identical histories, parameters and predictions. Measured R allocations changed by less than 0.03%; elapsed times were noisy on the shared machine, so unchanged ordinary paths are not claimed as a speedup. A separate single-thread, synthetic score benchmark with 2,000 traits, rank 5 and 200 individuals compared the new regression score routine with `U %*% solve(Sigma, L)`: maximum absolute difference was below $5\,10^{-14}$, median elapsed time was 0.003 versus 0.146 seconds (seven repetitions), and profiled R allocations were 3.84 versus 30.67 MiB. This measures **score calculation only**, not model fitting or process peak memory. Conditional EM is not uniformly competitive with AI. In a four-trait, 24-individual, four-replicate incomplete-data experiment with a near-boundary genetic specific variance, AI converged in about 5 seconds; conditional EM had not met both convergence tolerances after 1,500 updates and about 123 seconds. The engine reports nonconvergence explicitly. These measurements motivate future acceleration/reduced-latent-system work; they do not establish thousand-trait fitting scalability. On the same near-boundary fixture, conditional AI with active numerical bounds converged in 16 updates rather than the 1,500 unfinished EM updates. Conditional AI took about 1.5 seconds versus 0.6 seconds for ordinary AI. These shared-machine timings do not establish a general speed advantage; the numerical bounds and stopping rules must also be considered when comparing optima. ## Ordinary fitting paths The sequence below shows dependencies, not a sequence of arbitrary method switches. Several paths dispatch early, and later choices remain gated by the model, task, solver, and output request. ```text Model formulas, covariance descriptors, data and weights | +-- method$task = "solve" | +-- Known covariance parameters; no covariance estimation | +-- Gaussian identity and restricted solve/output requests | +-- method$task = "fit" | +-- Non-Gaussian family | +-- Outer PQL uses working; method defaults preserve REML working fits | +-- Inner fits remain subject to Gaussian-path restrictions | +-- Gaussian identity / inner Gaussian fit | +-- formulation = "direct" | +-- Observation-space fitting | +-- No Henderson solver selection | +-- formulation = "henderson" +-- solve$method = "pcg" | +-- REML estimation; iterative solves and stochastic | determinant/score evaluation +-- solve$method = "auto", "ldlt", or "cholmod" +-- Deterministic sparse-factorization paths +-- Exact determinant/score paths when eligible Representation, evaluation, update, proposal, convergence, output, execution, and diagnostics requests apply within those branches only when their model-dependent admission checks pass. ``` Model-independent combinations are checked by one compatibility table when `mmes_control()` builds the object and again in `mmes()`. A supplied setting that conflicts with another, or that the selected path would ignore, stops the call. The error names every offending setting and suggests a fix. Restating a default is always accepted. Rules that depend on the data, family, `vcc` or rotation are checked at fit time. See `?mmes_control` (section Compatibility) for the rule list. ```r mmes_control(solve = list(method = "ldlt"), score = list(traceProbes = 32L)) #> Error: Incompatible mmes_control() settings: #> * PCG preconditioning, recycling, stochastic log-determinant and #> trace-estimation settings only apply when solve$method='pcg' ... #> Offending: score$traceProbes. Remove them or set solve$method='pcg'. ``` Examples of coupled prerequisites: * `formulation="direct"` uses observation-space equations. Henderson settings such as `solve$method`, `representation`, `curvature` or PCG estimators have no effect there and stop the call when supplied. * For covariance estimation, PCG is a Henderson REML bundle. It combines iterative linear solves with stochastic determinant and score work, and supports `output$inverse` 0 or 2, not selected-inverse mode 1. The known- covariance solve task has its own restricted PCG path and requires `output$inverse=0`. Solve tolerance does not bound stochastic estimator error. * `proposal$extrapolation="aitken"` and `curvature$method="hybrid-ai"` cannot be requested together. * `representation$compression` requires an eligible Gaussian Henderson model, independent residuals, fixed diagonal weights, and matching design rows. PQL, known-covariance solving, rotation, correlated residuals, saddle-point equations, factor-score augmentation, and latent reduced rank are excluded. * `representation$equations="saddle-point"` requires an eligible deterministic Gaussian Henderson path. It is not compatible with PCG, PQL, known- covariance solving, acceleration, rotation, compression, factor-score augmentation, or latent reduced rank. * `representation$factorScores` requires a compatible deterministic Henderson REML model and eligible covariance structure. It excludes known-covariance solving, PCG, user covariance constraints, rotation, compression, positive EM blending, and latent reduced rank. * `representation$sparseFactors` applies only to eligible supported structured factors, including selected AR, ANTE, SAR, and CAR structures. It does not introduce new covariance models. On complete grid residual blocks it also keeps banded AR/ANTE precisions sparse in the coefficient matrix. * `representation$gridAugment` completes grid residual blocks such as `dsumm(vsm(ar1m(range), ar1m(row), ism(units)), by=trial)` by adding missing cells as pseudo-observations with dummy fixed effects (exact for REML). It requires a Gaussian identity REML Henderson fit with the `ldlt` or `cholmod` solver and no weights, rotation, factor scores, saddle-point equations or known-covariance solving; `"force"` reports why a model is ineligible. `representation$matrixFree` avoids eligible tensor-prior assembly for PCG; it does not make every model operation matrix-free. * Native Laplace `likelihood` execution controls are distinct from Gaussian solver controls: `mode$warmStart="auto"` and `optimizer$gradient="auto"` are defaults (auto retains native analytic scores); `mode$warmStart="off"` restores zero starts and `optimizer$gradient="central"` selects reference central differences. `likelihood$trace="panels"` and `likelihood$blocks="off"` remain defaults; exact `trace="selected"` and `blocks="on"` are opt-in random-backend paths. Failed warm modes retry zero. ZIP auto does not reuse modes; explicit on cold-verifies the historical root. Missing inverse-subset entries use exact solve panels. Work/choices persist through refits and saved models, and all public SE/PEV remain post-fit. Non-default settings are rejected for PQL. * `solve$blockSchur`, `blockBanded`, `factorSchurMultiFactor`, and `dFirstML` request specialized elimination only for eligible coefficient graphs, likelihoods, backends, and workspaces. Forcing a request does not bypass structural or memory checks. * `score$traceBatch` requires an eligible exact LDLT score path and panel workspace. `score$covarianceChecks` can use available certificates, but no certificate bypasses positive-definiteness validation. * `curvature$reduced` needs an eligible reducing constraint map and supported update path. `update$commonScale` needs free product scales, a full-rank fixed design, and scale-preserving constraints. Both exclude several accelerated, PCG, PQL, and special-representation paths. * `output$covariance="none"` is restricted to supported Gaussian Henderson fits; it does not disable curvature needed by optimization. Setting For Henderson, `output$coefficientMatrix=FALSE` requires the eligible saddle-point path and `output$inverse=0`. For direct fits, `coefficientMatrix=TRUE` retains sparse prepared inputs for post-fit covariance, not Henderson equations; `FALSE` omits those inputs and disables post-fit uncertainty methods. These output requests do not substitute for one another. `output$selectedInverse` matters only for eligible Henderson selected-inverse output (`output$inverse=1`); `output$panelSize` controls full-inverse panels (`output$inverse=2`), and diagonal `postPEV()` panels for direct fits. * GPU requests require a CUDA-enabled Henderson build and usable CUDA runtime and device. Direct inversion and the known-covariance solve path remain CPU paths. * `output$inverse` controls coefficient-matrix uncertainty. `output$covariance` controls covariance-parameter information. `output$coefficientMatrix` controls retained coefficient-matrix output. These are separate requests with separate eligibility conditions. * `"auto"` records a request for policy-based selection, not a promise that a specialized method will be selected. `"force"` cannot bypass mathematical or structural restrictions. For the descriptive capability inventory and conceptual decision graph, inspect the internal catalog: ```r capabilities <- sommer:::.mmes_public_catalog() capabilities$controls ``` This inventory does not inspect a fitted model, resolve eligibility, or predict the selected execution path. It complements, but does not replace, control validation and fit-time diagnostics. # Migration from removed numerical formals Move numerical policy into `control`, statistical choices into `method`, and model starts/family nuisance choices into their `mmes()` arguments. Removed `task` control paths and `nonlinear=` are not compatibility aliases. | Former argument or path | Current setting | |:--|:--| | `henderson=TRUE/FALSE` | `formulation="henderson"/"direct"` | | `solver` | `solve$method` | | `nIters` | `convergence$maxIters` | | `tolParConvLL` | `convergence$likelihoodTol` | | `tolParConvNorm` | `convergence$parameterTol` | | `tolParInv` | `execution$inverseTol` | | `REML=TRUE/FALSE`, `task$objective` | `method=mmes_method(objective="reml"/"ml")` | | `solveOnly=TRUE/FALSE`, `task$mode` | `method=mmes_method(task="solve"/"fit")` | | `pqlControl`, `task$pql` | `control=mmes_control(working=list(...))` | | `task$glmm$method="laplace"` | `method=mmes_method(objective="ml", integration="laplace")` | | `task$glmm$familyParameters` | `mmes(familyParameters=list(...))` | | `task$nonlinear$start` | `mmes(start=list(parameters=...))` | | `nonlinear=curve` | `curves=list(curve=curve)` and formula `y~0+curve_term("curve")` | | `task$nonlinear$mode/population/inference` | `control` modules `mode`, `population`, `inference` | | `computeCi` | `output$inverse` | | `getPEV` | No separate replacement; Gaussian fits use `output$inverse` for fit-time PEVs. Native Laplace requires `output$inverse=0` and uses `postPEV()` after fitting. | | `returnParam` | `output$setup` | | `blend`, `blendWeight` | `update$recovery`, `update$blendWeight` | | `stepWeight` | `proposal$stepWeight` | | `acceleration="aitken"` | `proposal$extrapolation="aitken"` | | `acceleration="hybrid-ai"` | `curvature$method="hybrid-ai"` | | `compression` | `representation$compression` | | `factorScoreAugmentation` | `representation$factorScores` | | `experimental$sparseFactors`, `matrixFreePcg`, `saddlePoint` | `representation$sparseFactors`, `representation$matrixFree`, `representation$equations` | | `experimental$blockSchur`, `blockBanded`, `factorSchurMultiFactor`, `dFirstML` | Same-named `solve` fields | | `experimental$pdCertificates`, `traceBatch` | `score$covarianceChecks`, `score$traceBatch` | | `experimental$reducedAI`, `commonScale` | `curvature$reduced`, `update$commonScale` | | `experimental$cholmodFinalSelected`, `panelCi2` | `output$selectedInverse`, `output$panelSize` | | `experimental$profile` | `diagnostics$profile` | | `covarianceInformation` | `output$covariance` | | `retainCoefficientMatrix` | `output$coefficientMatrix` | | PCG `tol`, `maxIters`, preconditioner, reuse | `solve` fields | | PCG determinant scaling and Lanczos depth | `logdet$scaling`, `logdet$lanczosSteps` | | PCG trace probes, estimator, rank | `score$traceProbes`, `score$estimator`, `score$rank` | | PCG covariance checks and trace batching | `score$covarianceChecks`, `score$traceBatch` | | PCG accuracy budgets | `convergence` fields | | PCG seed | `execution$seed` | | Experimental elimination requests | Corresponding `solve` fields | | Experimental sparse factors, matrix-free, saddle point | Corresponding `representation` fields | | Experimental trace checks/batching | `score$covarianceChecks`, `score$traceBatch` | | Experimental reduced AI | `curvature$reduced` | | Experimental common scale | `update$commonScale` | | Experimental inverse selection/panel size | `output$selectedInverse`, `output$panelSize` | | Experimental profiling | `diagnostics$profile` | | `denseMemoryMB`, `memoryPolicy`, `gpu`, `gpuMinSize`, `blasThreads` | Same-named `execution` fields | | `verbose`, `dateWarning` | Same-named `diagnostics` fields | | `recordDiagnostics`, `verboseDiagnostics` | Removed flags; use `diagnostics$verbose` for fit output and `diagnostics$profile` for profiling | | `engine` | Pass the complete object as `control` | | `...` legacy/pass-through options | No replacement; set documented module fields only | For example, former settings spread across top-level formals and engine lists now have one canonical owner: ```r control <- mmes_control( formulation = "henderson", solve = list(method = "pcg", tol = 1e-8, maxIters = 0L), score = list(traceProbes = 16L), convergence = list(accuracy = "adaptive", maxIters = 40L), output = list(inverse = 0L, covariance = "full") ) method <- mmes_method(objective = "reml") ``` The statistical model is unchanged by choosing a numerical representation or solver only when that method is mathematically equivalent and eligible. Do not interpret the module structure as permission to combine unsupported methods.