mixed model equations solver

Fits a generalized linear mixed model coded in C++ using the Armadillo and Eigen libraries to optimize matrix operations. Special covariance structures can be specified with the vsm() function. Two engines are available; henderson formulation and direct reml (henderson argument) and three solvers in each; ldlt, cholmod and pcg (solver argument) to adjust to the modeling situation.

Usage

mmes(fixed, random, rcov, data, W, weights=NULL,
     nIters=30, tolParConvLL=1e-04,
     tolParConvNorm=1e-04, tolParInv=1e-06,
     naMethodX="exclude", naMethodY="exclude",
     naMethodRandom="exclude", naMethodR="exclude",
     returnParam=FALSE, dateWarning=TRUE,
     verbose=TRUE, stepWeight=NULL, emWeight=NULL,
     contrasts=NULL, getPEV=TRUE, henderson=TRUE,
    computeCi=0, solver="auto", pcgTol=1.0e-8,
    pcgMaxIters=0, pcgTraceProbes=8,
    pcgLanczosSteps=20, REML=TRUE, vcc=NULL,
    family=stats::gaussian(), pqlControl=list(),
    .pqlInner=FALSE, .pqlFixedDispersion=FALSE,
    .pqlWorkingPrecision=NULL, .pqlBaseW=NULL,
    .pqlBaseFactor=NULL, acceleration="none", .pqlStart=NULL,
    factorScoreAugmentation="none", .factorScoreParameters=NULL,
    pcgPreconditioner="diagonal", pcgNystromRank=32L,
    solveOnly=FALSE, covPar=NULL, blasThreads=NULL)

Arguments

fixed

A formula specifying the response variable(s) and fixed effects, i.e:

response ~ covariate

random

A formula specifying the random effects, e.g., random = ~ genotype + year. A simple random term that is not wrapped in vsm is treated as an identity-covariance random effect.

The vsm function is the main interface for specifying structured covariance models for random effects. Its current covariance model is

$$G = \sigma^2 (K_1 \otimes K_2 \otimes \cdots \otimes K_m) \otimes A,$$

where $\sigma^2$ is the single variance scale owned by vsm(), $K_1,\ldots,K_m$ are covariance-shaping factors supplied by covariance constructors, and $A$ is the known covariance relationship for the final random effect. In the Henderson implementation Gu supplies the inverse of $A$; it must have row and column names matching the levels of the final random effect and attr(Gu,"inverse")=TRUE.

Within vsm(...), the last constructor supplies the incidence matrix for the main random effect, while all preceding constructors define covariance-shaping factors. The number of covariance factors is not limited. For example,

random = ~ vsm(dsm(Location), ism(Name))

fits heterogeneous variances across Location for the random effect Name, while

random = ~ vsm(dsm(Environment), ar1m(Row), ar1m(Column), ism(Name), Gu=Ainv)

specifies an arbitrary Kronecker product of covariance structures before the relationship matrix for Name.

Covariance constructors currently available include:

ism for identity covariance;

dsm and atm for diagonal and selected-level diagonal covariance structures;

usm for an unstructured covariance matrix;

csm and corgm for correlation structures;

ar1m, ar2m, and ar3m for autoregressive covariance structures;

mam (including ma1m and ma2m) for moving-average covariance structures;

toeplitzm for a general positive-definite Toeplitz correlation structure;

fam for factor-analytic covariance;

rrm for the reduced-rank covariance approximation;

antem for antedependence covariance;

maternm for Matern spatial covariance;

metricm for exponential, power, gaussian, spherical and circular distance-based correlations (isotropic or anisotropic);

sar and car for simultaneous and conditional autoregressive spatial covariance models; and

ownm for a fixed or user-defined covariance function.

Design-matrix utilities such as overlay, spl2Dc, leg, and redmm can still be combined with vsm() when appropriate.

See vsm and the individual covariance-constructor help pages for details on parameterizations and examples.

A single relationship term may request Lee–van der Werf rotation with vsm(..., Gu=Gu, rotation=TRUE). With henderson=TRUE, random coefficients are represented in the eigenbasis and the diagonal relationship precision is used in the mixed model equations. With henderson=FALSE, the equivalent transformed observation-covariance parameterization is used. Public random effects, fitted values, and residuals are returned in the original basis. The current implementation requires a complete balanced Gaussian layout, one rotated relationship term, and computeCi=0 during fitting. Residual structures are allowed when they are rotation invariant: residuals may be structured across rotation blocks (e.g., rcov=~vsm(usm(trait), ism(units)) or rcov=~vsm(dsm(Env), ism(units))), but must be independent and identically structured across the levels of the rotated relationship term. Residual structures that correlate observations of different levels (e.g., ar1m(), maternm(), sar(), or car() over plots), residual variances that differ among levels within a rotation block, and W matrices that are not diagonal with a constant weight within each rotation block are rejected, as are non-Gaussian PQL fits.

rcov

A formula specifying the residual covariance structure. The default is rcov = ~ units, which fits a homogeneous residual variance.

Structured residual covariance is specified with the same vsm and CovarianceFactor constructors used for random effects. The residual term should be expressed as a single vsm() term; arbitrary covariance complexity is represented by placing multiple covariance factors inside that term rather than by summing multiple residual terms.

The final term is normally ism(units), where units is the observation-level incidence generated internally by mmes. For example,

rcov = ~ vsm(dsm(Environment), ism(units))

fits heterogeneous residual variances across environments, while

rcov = ~ vsm(ar1m(Row), ar1m(Column), ism(units))

fits a separable row-by-column AR1 residual covariance, with one correlation per direction shared by every other Kronecker factor. For section-specific variances and correlations, use

rcov = ~ dsumm(vsm(ar1m(Row), ar1m(Column), ism(units)), by = Trial)

(see dsumm). More generally, structures such as usm, csm, toeplitzm, maternm, sar, car, and ownm can be used when their incidence layout is appropriate for the residual model.

Residual covariance factors must identify exactly one covariance-product coordinate for each observation. mmes() constructs the residual block and local-coordinate indexing internally, so the data do not need to be manually sorted by the variables defining the residual structure.

Residual pairing key (multi-trait models). For a long-format multi-trait model the residual must know which records belong to the same experimental unit. With ism(units) records are paired implicitly by their order of appearance within each trait; a warning is issued when repeated residual coordinates must be paired by row order. An explicit key is given by replacing units with a column identifying the unit, e.g.

rcov = ~ vsm(usm(trait), ism(record))

where record is the column created by stackTraits. Records sharing a key value form one residual block (missing traits are simply absent from it); a trait may appear only once per key value.

As with random effects, vsm() owns one residual variance scale $\sigma^2$; the preceding covariance constructors provide dimensionless covariance shapes that are combined through an arbitrary Kronecker product.

See vsm and the individual covariance-constructor help pages for details.

data

Optional data frame containing variables used by the model. Model expressions are evaluated using normal R scoping rules: variables are first sought in data and may also be obtained from the environments associated with the model formulas or from the calling environment. Consequently, auxiliary objects such as relationship matrices, coordinate vectors, or grouping variables do not have to be copied into data. If data is omitted, variables are evaluated from the formula/calling environment.

Observation-level variables used by the response, fixed effects, random effects, and residual covariance specification must nevertheless have compatible lengths. Internally, mmes() constructs a common observation mask and applies it consistently to all observation-level model components.

W

Weights matrix (e.g., when covariance among plots exist). Internally W is squared and inverted as Wsi = solve(chol(W)), then the residual matrix is calculated as R = Wsi*O*Wsi.t(), where * is the matrix product, and O is the original residual matrix. Weights can be combined with any residual structure, including Kronecker (multi-trait) residuals.

weights

Optional one-sided formula describing independent row blocks of the supplied W matrix, for example weights = ~trial. Rows with different formula groups must have zero cross-group entries in W. This is structural metadata only: it does not generate numeric weights. Declared blocks are processed independently by the Henderson weighting path; without this formula, W is handled without a known grouping.

nIters

Maximum number of iterations allowed.

tolParConvLL

Convergence criteria based in the change of log-likelihood between iteration i and i-1.

tolParConvNorm

When using the Henderson method this argument is the convergence criteria based in the norm proposed by Jensen, Madsen and Thompson (1997):

e1 = || InfMatInv.diag()/sqrt(N) * dLu ||

where InfMatInv.diag() is the diagonal of the inverse of the information matrix, N is the total number of variance components, and dLu is the vector of first derivatives.

tolParInv

Tolerance parameter for matrix inverse used when singularities are encountered in the estimation procedure. By default the value is 1e-06. This parameter should be fairly small because the it is used to bend matrices like the information matrix in the henderson algorithm or the coefficient matrix when it is not positive-definite.

naMethodX

Missing-data policy for variables entering the fixed-effect design matrix. The default, "exclude", marks observations with missing fixed-effect information for removal from the common observation set. "include" retains the historical sommer behavior for supported cases by imputing missing fixed-effect covariates. This option should be used carefully, because imputation changes the fitted design matrix.

naMethodY

Missing-data policy for the response variable(s). The default, "exclude", removes observations with missing response information from the common observation set. The historical "include" and "include2" behaviors are retained where supported: "include" imputes missing responses, whereas "include2" is intended for multivariate-response situations in which an observation may be retained when at least one response is observed.

naMethodRandom

Missing-data policy for observation-level variables used to construct random-effect terms. The default is "exclude". Under this policy, an observation for which a required random-effect classification or covariate cannot be evaluated is excluded from the common observation set. Random terms are evaluated before the final observation subset is applied, so the same inclusion mask is used consistently for the response, fixed-effect design, and all random-effect design matrices.

naMethodR

Missing-data policy for observation-level variables used to define the residual covariance structure. The default is "exclude". Observations lacking information required to assign their residual covariance-product coordinate are excluded from the common observation set. The residual block and local-coordinate layout are then constructed on the retained observations.

returnParam

A TRUE/FALSE value to indicate if the program should return the parameters to be used for fitting the model instead of fitting the model.

dateWarning

A TRUE/FALSE value to indicate if the program should warn you when is time to update the sommer package.

verbose

A TRUE/FALSE value to indicate if the program should return the progress of the iterative algorithm.

stepWeight

A vector of values (of length equal to the number of iterations) indicating the weight used to multiply the update (delta) for variance components at each iteration. If NULL the 1st iteration will be multiplied by 0.5, the 2nd by 0.7, and the rest by 0.9. This argument can help to avoid that variance components go outside the parameter space in the initial iterations which happens very often with the AI method but it can be detected by looking at the behavior of the likelihood. In that case you may want to give a smaller weight.

emWeight

A vector of values (of length equal to the number of iterations) indicating with values between 0 and 1 the weight assigned to the EM information matrix; the values 1 - emWeight are applied to the AI information matrix to produce the joint information matrix used in the variance-component update. This creates an EM-heavy warm start early in the fit and gradually transitions to an AI-dominated update as the optimizer approaches convergence. The default schedule is a logarithmic decay from 1 to 0.03 over the first min(nIters, 20) iterations, followed by a 0.03 floor. Values outside [0, 1] are rejected.

contrasts

an optional list. See the contrasts.arg of model.matrix.default.

getPEV

A logical value indicating whether prediction error variance results should be organized and returned when they are requested through computeCi.

henderson

A logical value indicating which REML engine to use. TRUE (the default) uses the Henderson mixed-model-equations algorithm, efficient when there are many more records than coefficients to estimate. FALSE uses a direct-inversion engine that inverts the phenotypic covariance matrix in observation space instead; it is less efficient than henderson=TRUE when there are many more records than coefficients, but can be more efficient when there are many more coefficients to estimate than records available (e.g., marker/SNP-BLUP-style models). Both engines share the same vsm() covariance-structure parameterization. The direct-inversion engine currently supports computeCi values 0 and 2 only (no Takahashi selected-inverse mode), has no solver/PCG/CHOLMOD options, and requires a single response column (use the long-format vsm(usm(trait), ...) convention for multi-trait models, as with the Henderson engine).

computeCi

An integer indicating whether post-fit prediction error variance information and/or the inverse of the mixed model coefficient matrix should be computed when the Henderson algorithm is used (henderson=TRUE). The available options are:

0: no additional inverse-related computation is performed after convergence. The full inverse of the coefficient matrix (Ci) is not formed and uPevList is not computed. This is the fastest option and is the default.

1: prediction error variances are computed using the Takahashi sparse inverse subset algorithm. This method uses the sparse LDLT factorization of the coefficient matrix and recursively evaluates only those elements of the inverse that belong to the sparsity pattern induced by the factorization. In particular, the diagonal elements required for prediction error variances are obtained without constructing the complete inverse matrix. Therefore, uPevList is returned while Ci is not fully materialized.

With solver="cholmod" (see below), the supernodal factorization used during REML iterations has no Takahashi selected-inverse equivalent. In that case computeCi=1 automatically performs one additional sparse LDLT factorization of the converged coefficient matrix after REML convergence, purely to extract the selected inverse subset. This is a single one-time cost, not repeated every iteration, so the iterative speed benefit of solver="cholmod" is retained.

2: the complete inverse of the coefficient matrix is computed. The prediction error variances in uPevList are then extracted from the appropriate diagonal elements of Ci. This option is the most computationally and memory intensive, especially for large mixed model equation systems. With solver="cholmod" the full inverse is obtained directly from the supernodal factor (no extra LDLT refactorization is needed for this mode).

For large models, computeCi=0 is recommended during model fitting. If prediction error variances or the full inverse are needed afterwards, they can be obtained with postPEV without refitting the model.

solver

Linear-system solver used by the Henderson implementation. "auto" (default) picks a solver automatically based on the density of the random-effect relationship matrices (Gu) actually supplied: if any random effect uses a dense Gu (e.g. a genomic or marker-based relationship matrix, which is typically close to fully dense), "cholmod" is selected; otherwise (e.g. sparse pedigree-based relationship matrices, or no Gu at all) "ldlt" is selected. Passing an explicit value ("ldlt", "pcg", or "cholmod") always overrides the automatic choice. "ldlt" uses the sparse direct simplicial LDLT factorization path. "pcg" selects the iterative preconditioned conjugate-gradient path. For REML fits with computeCi=0, random-effect precisions produced by vsm() are applied matrix-free using their Kronecker factors, avoiding explicit assembly of $K^{-1} \otimes A^{-1}$ during optimization. The complete coefficient matrix is materialized once after convergence for compatibility with prediction and post-fit methods. Other PCG configurations transparently use the assembled sparse path. "cholmod" uses the supernodal (blocked, BLAS-3) sparse Cholesky factorization bundled with R's Matrix package (SuiteSparse CHOLMOD), accessed through Matrix's public C API - no additional system library needs to be installed, since Matrix is already a dependency of sommer. Supernodal factorization can be faster than the simplicial LDLT path for large coefficient matrices (e.g., models built on large pedigree or genomic relationship matrices), because it reorganizes the sparse factorization into dense blocks that are handled with threaded BLAS-3 kernels rather than column-by-column updates; the speed-up depends on the BLAS linked to the running R installation. solver="cholmod" currently supports all computeCi values (see above for how computeCi=1 is handled for this solver). The covariance-model specification is independent of this choice: all solvers receive the same generic CovarianceFactor descriptors produced by vsm().

pcgTol

Convergence tolerance used by the PCG linear solver. The default is 1e-8. This argument is used when solver="pcg". Smaller values request more accurate iterative solves but may increase computation time.

pcgMaxIters

Maximum number of PCG iterations for an individual linear solve. A value of 0 lets the C++ solver choose its internal iteration limit from the problem dimension. This argument is used only by the PCG path.

pcgTraceProbes

Number of stochastic probe vectors used by the PCG-based machinery for trace quantities required during REML calculations. The default is 8. Increasing this value can reduce stochastic approximation variability at additional computational cost. This argument is relevant to solver="pcg".

pcgLanczosSteps

Number of Lanczos steps used by the PCG-based stochastic log-determinant calculations. The default is 20. Larger values can improve the approximation at additional computational cost. This argument is relevant to solver="pcg".

REML

Logical, TRUE (default) fits variance/covariance parameters by restricted maximum likelihood (REML), the classical choice for BLUP-style variance component estimation. Setting REML=FALSE fits by maximum likelihood (ML) instead: the log-likelihood, score equations, and reported monitor/log-lik values drop the REML "restriction" term that accounts for the degrees of freedom used to estimate the fixed effects. ML estimates of variance components are known to be biased downward relative to REML, but unlike REML log-likelihoods, ML log-likelihoods are directly comparable (via a likelihood-ratio test) across models that differ in their fixed-effects specification - REML likelihoods are only comparable across models that share the same fixed effects. REML=FALSE currently requires solver="ldlt" or solver="cholmod" (solver= "auto" already resolves to one of these). solver="pcg" is not yet supported with REML=FALSE. Prediction error variances (uPevList/computeCi) are still computed from the REML-style inverse regardless of REML.

vcc

Optional data frame of equality/scaling constraints between covariance parameters (columns parameter, group and optionally scale), identified through the vcParams table returned by returnParam=TRUE and stored in the fit. See vcc.

family

A stats::family() object. The default Gaussian identity family uses the ordinary linear mixed-model implementation. Other families are fitted by penalized quasi-likelihood (PQL): each outer iteration forms an IRLS working response and weights, then fits the resulting weighted Gaussian mixed model with the Henderson solver. The initial implementation supports one numeric response. Any family with a valid link can be used, e.g. binomial() (logit, probit, cloglog), poisson(), Gamma(), inverse.gaussian(), gaussian(link="log"), quasipoisson() and MASS::negative.binomial(theta). Binomial counts out of a known number of trials are supplied as a single proportion response built with binm, whose trials are used as prior weights. Binomial, Poisson and negative binomial working residual dispersions are fixed to one; for other families the reported dispersion is the Pearson estimate. The negative binomial $\theta$ is estimated (conditional on the random effects) with pqlControl=list(estimateTheta=TRUE) and returned in theta.nb and theta.nb.se. Because the working-model log-likelihood is not a GLMM likelihood, llik, AIC and BIC are NA for PQL fits and anova() likelihood ratio tests are refused; use deviance and Wald tests instead. Formula offsets are included on the link scale. A supplied symmetric positive definite W is combined with the IRLS precision after observation filtering, and may be sparse and non-diagonal.

Long-format multi-trait data with a different family per trait are fitted by supplying family = familym("trait", Disease = binomial(), Yield = gaussian()). See familym for the required residual structure; in that case dispersion is a named vector with one value per trait.

pqlControl

A named list controlling non-Gaussian PQL fits. Supported entries are maxit (outer-iteration limit, default 20), tol (relative deviance convergence tolerance, default 1e-5) and estimateTheta (estimate the negative binomial $\theta$, default FALSE) and warmStart (reuse covariance parameters from the preceding working fit when the observation and model layouts agree, default TRUE). The pqlMonitor also reports warmStarted and innerIterations and symbolicAnalyses. With the LDLT solver, compatible warm fits reuse symbolic factorizations and selected-inverse topology across PQL iterations while recomputing numeric factors. Other solvers reuse their symbolic factorizations within each working fit.

.pqlInner

Internal-use logical indicating that mmes() is being called for a Gaussian working-model fit within the PQL iteration. Users should not set this argument directly.

.pqlFixedDispersion

Internal-use logical indicating that the working residual dispersion is fixed, as required for binomial and Poisson PQL models. Users should not set this argument directly.

.pqlWorkingPrecision

Internal-use vector containing the current IRLS working precisions after observation filtering. Users should not set this argument directly.

.pqlBaseW

Internal-use base observation precision matrix supplied by the outer PQL fit before combination with the IRLS working precisions. Users should not set this argument directly.

.pqlBaseFactor

Internal-use cached factor of .pqlBaseW, reused across PQL iterations. Users should not set this argument directly.

acceleration

Optional optimizer acceleration: "none" (default) or "aitken". Aitken proposals are capped, applied every third accepted iteration during EM blending, and checked by the existing likelihood line search. Requires a deterministic Henderson solver; not available for PCG or direct inversion. The Henderson result's engineDiagnostics reports coefficient evaluations, line-search halvings, accelerated proposals, symbolic-analysis and ML selected-inverse topology counters, and active block-Schur partition dimensions, grouped weight-block counts and assembled PCG batch solve/iteration counts, plus peak sparse fill in the irregular residual fallback.

.pqlStart

Internal-use covariance warm-start state for successive PQL working fits.

factorScoreAugmentation

Latent-factor parameterization for FA/RR and compound-symmetry covariance factors with nonnegative correlation. The default "none" uses the marginal covariance MME. "fixed-shape" uses the augmented MME with fixed loading and specific-variance parameters; it is currently limited to Henderson REML, computeCi=0, and one eligible low-rank-plus-diagonal shaping factor per random term. It replaces each $q d$ marginal coefficient block with $(k+q)d$ factor and specific-effect coefficients constrained to share one variance scale. "profile" estimates free eligible shape parameters with a bounded L-BFGS-B outer profile optimizer. It starts from the values in the constructors (so it does not first factor the expensive marginal FA/RR system); at each candidate shape it rebuilds the augmented design and refits shared random and residual scales by REML. The outer iteration budget is dimension-adjusted, with at least five and at most twelve outer iterations; nIters controls each inner AI fit. Inner covariance starts and symbolic factorizations are reused when the augmented layout matches; numeric factors are recomputed. The outer objective uses finite-difference gradients, is experimental, can require many augmented REML fits, and can converge to a local optimum; inspect factorScoreProfile. Shape-parameter standard errors are returned as NA; the inner AI information matrix is not their uncertainty estimate.

Both augmented modes currently require Gaussian identity Henderson REML, computeCi=0, one eligible low-rank-plus-diagonal shaping factor per random term, no rotation, and no user vcc. With solver="auto", dense Gu selects CHOLMOD; its exact block-Schur path eliminates independent specific-effect blocks while retaining factor effects in the border. The returned uList, bu, fitted values, and prediction contrasts use the original effects; C remains augmented and is tagged with CRepresentation="factor-score-augmented".

.factorScoreParameters

Internal-use FA/RR shape parameters on the working scale, indexed by random term. The profile optimizer uses these to fix shapes in its inner fits; users should not set this argument directly.

pcgPreconditioner

PCG preconditioner: "diagonal" (default) or "nystrom". The Nyström option builds a deterministic landmark approximation to assembled $C$ and applies its Woodbury inverse as an SPD preconditioner; PCG continues to multiply by the exact $C$ and verifies the true residual. Nyström currently uses assembled-C PCG, so it disables the matrix-free shortcut.

pcgNystromRank

Positive number of deterministic landmarks for pcgPreconditioner="nystrom"; the default is 32 and the value is capped at the coefficient dimension.

solveOnly

If TRUE, no variance components are estimated: the mixed-model equations are solved once at known covariance parameters (genetic evaluation). The coefficient matrix is never assembled; a matrix-free preconditioned conjugate gradient iterates on the data with block-Jacobi preconditioning (one block for the fixed effects when there are at most 1000 of them, and one covariance-dimension block per random-effect level). pcgTol is the relative-residual tolerance and pcgMaxIters the iteration limit (0 = automatic). Requires a Gaussian identity model, computeCi=0, no rotation, no vcc, a diagonal W if any, and a block-diagonal residual covariance whose blocks have at most 10000 records. The result is an object of class mmesSolve with b, u, bu, uList, fitted, residuals, theta (covariance matrices used), vcParams (the natural-scale parameters used, reusable as covPar) and pcg diagnostics. No PEV or likelihood is computed.

covPar

Known covariance parameters for solveOnly=TRUE. One of: a fitted mmes object with the same random and residual terms (its final estimates are used exactly); a list with one element per term, named by term label or in formula order (random terms, then residual), each holding either the natural-scale parameter vector (as in fit$covPar) or, for ism(), dsm() and usm() terms, the full covariance matrix (e.g. $G_0$ and $R_0$); or a data frame with columns term, parameter and value (as vcParams of an mmesSolve result, or the vcParams table of returnParam=TRUE with a value column). When NULL, every parameter must already be fixed in the formula (e.g. vsm(ism(id), Gu=Ainv, sigma2=0.3, fixedSigma2=TRUE)).

blasThreads

Optional BLAS thread limit for the duration of this call: NULL (default) preserves the existing setting, a positive integer requests that count, and "auto" uses allocation-aware CPU availability from the optional parallelly package. Runtime control requires the optional flexiblas package for FlexiBLAS, or RhpcBLASctl for a recognized supported BLAS. An unsupported or rejected explicit request is an error. The previous setting is restored on return or error. This is a process-wide BLAS setting, not an OpenMP limit or multi-node allocation. OpenMP workers may also invoke BLAS: avoid oversubscription by choosing appropriate OMP_NUM_THREADS before launching R. More BLAS threads need not improve performance. Verbose output distinguishes configured limits from measured CPU use. The result's threadingDiagnostics records the backend, configured BLAS threads (NA when unknown), allocation-aware CPUs, environment, numerical path, and CPU/elapsed times. Average CPUs means aggregate process CPU seconds divided by elapsed seconds, not a sampled active-thread count. Henderson fits additionally record numeric factorization measurements, excluding symbolic analysis; block-Schur measurements include block extraction and Schur assembly. For PQL/profile fits these native measurements describe the returned inner fit, while the fit-level measurement covers the outer call.

Details

Iteration output

With verbose=TRUE and henderson=TRUE, the iteration table reports:

iteration: The iteration number.
LogLik: The log likelihood at the current accepted parameter values (restricted log likelihood when REML=TRUE).
wall: The local clock time, displayed as HH:MM:SS.
cpu(sec): Cumulative elapsed seconds since the native fitting routine started; despite the column name, this is not process CPU time.
restrained: The number of variance parameters restrained at a boundary during the iteration.
EMweight: The weight of the EM information in the blended AI/EM update; 1 uses EM information and 0 uses AI information.
notPD: The one-based covariance structure indices whose proposed updates failed the positive-definiteness check and required local fallback. Multiple indices are comma-separated (e.g., 2,3); - means none failed. Indices follow the random terms, with the residual structure last.
trustAlpha: The multiplier applied to the joint proposed step in working parameter coordinates after local safeguards. A value of 1 means no trust scaling was needed; smaller values shorten the step. This is distinct from any subsequent likelihood line-search halving.
pivot: For solver="ldlt" only, the smallest diagonal pivot in the LDLT factorization of the mixed-model coefficient matrix, a numerical diagnostic of positive definiteness.

The safeguard columns describe the proposed update, whereas LogLik describes the current accepted fit. Floating-point values are displayed to three decimal places (scientific notation is used when needed to fit the column). This formatting does not round stored estimates; very small positive values may therefore appear as 0.000.

Evaluation environments and observation filtering

The current interface follows R-style model evaluation. Variables referenced by fixed, random, and rcov may be supplied as columns of data or resolved from the formula/calling environment. This makes it possible, for example, to keep relationship matrices, external grouping variables, or spatial coordinates outside the analysis data frame.

Missingness is handled through a common observation-selection mechanism. The response, fixed effects, random effects, and residual covariance specification are evaluated against the original observation set and their validity is combined into a single inclusion mask. The same retained rows are then used for the response vector, fixed-effect design matrix, random-effect design matrices, optional weight matrix, and residual covariance indexing. This avoids independently subsetting different parts of the mixed model.

The fitted object contains obsInfo, an observation-level audit table describing the filtering process. It records the original row index, validity of the response, fixed, random, and residual components, whether the observation was retained, and, when applicable, the reason for exclusion. This can be useful for diagnosing model specifications involving missing values.

Estimated marginal means

When the suggested emmeans package is installed, univariate mmes fits can be supplied directly to emmeans::emmeans(). The resulting means are population-level fixed-effect marginal means with covariance obtained from the mixed-model equations. Inference uses asymptotic degrees of freedom (df = Inf). Multivariate-response fits are not currently supported by this interface.

Random and residual formulas are evaluated as R language objects rather than by splitting their textual representation at plus signs. Consequently, arithmetic or nested expressions inside model terms are not interpreted as separate top-level random or residual terms.

The use of this function requires a good understanding of mixed models. Please review the 'sommer.quick.start' vignette and pay attention to details like format of your random and fixed variables (e.g. character and factor variables have different properties when returning BLUEs or BLUPs).

For tutorials on how to perform different analysis with sommer please look at the vignettes by typing in the terminal:

vignette("v1.sommer.quick.start")

vignette("v2.sommer.changes.and.faqs")

vignette("v3.sommer.qg")

vignette("v4.sommer.gxe")

Citation

Type citation("sommer") to know how to cite the sommer package in your publications.

Special variance structures

vsm(atm(x,levels),ism(y))

can be used to specify heterogeneous variance for the "y" covariate at specific levels of the covariate "x", e.g., random=~vsm(at(Location,c("A","B")),ism(ID)) fits a variance component for ID at levels A and B of the covariate Location.

vsm(dsm(x),ism(y))

can be used to specify a diagonal covariance structure for the "y" covariate for all levels of the covariate "x", e.g., random=~vsm(dsm(Location),ism(ID)) fits a variance component for ID at all levels of the covariate Location.

vsm(usm(x),ism(y))

can be used to specify an unstructured covariance structure for the "y" covariate for all levels of the covariate "x", e.g., random=~vsm(usm(Location),ism(ID)) fits variance and covariance components for ID at all levels of the covariate Location.

vsm(usm(rrmat(x,y,z,nPC)),ism(y))

can be used to specify an unstructured covariance structure for the "y" effect for all levels of the covariate "x", and a response variable "z", e.g., random=~vsm(rrmat(Location,ID,response, nPC=2),ism(ID)) fits a reduced-rank factor analytic covariance for ID at 2 principal components of the covariate Location.

vsm(ism(overlay(...,rlist=NULL,prefix=NULL)))

can be used to specify overlay of design matrices between consecutive random effects specified, e.g., random=~vsm(ism(overlay(male,female))) overlays (overlaps) the incidence matrices for the male and female random effects to obtain a single variance component for both effects. The `rlist` argument is a list with each element being a numeric value that multiplies the incidence matrix to be overlayed. See overlay for details.Can be combined with vsm().

vsm(ism(redmm(x,M,nPC)))

can be used to create a reduced model matrix of an effect (x) assumed to be a linear function of some feature matrix (M), e.g., random=~vsm(ism(redmm(x,M))) creates an incidence matrix from a very large set of features (M) that belong to the levels of x to create a reduced model matrix. See redmm for details.Can be combined with vsm().

vsm(leg(x,n),ism(y))

can be used to fit a random regression model using a numerical variable x that marks the trayectory for the random effect y. The leg function can be combined with the special functions dsm, usm at and csm. For example random=~vsm(leg(x,1),ism(y)) or random=~vsm(usm(leg(x,1)),ism(y)).

spl2Dc(x.coord, y.coord, at.var, at.levels))

can be used to fit a 2-dimensional spline (e.g., spatial modeling) using coordinates x.coord and y.coord (in numeric class) assuming multiple variance components. The 2D spline can be fitted at specific levels using the at.var and at.levels arguments. For example random=~spl2Dc(x.coord=Row,y.coord=Range,at.var=FIELD).

Covariance between random effects

covm( vsm(ism(ran1)), vsm(ism(ran2)) )

can be used to specify covariance between two different random effects, e.g., random=~covm( vsm(ism(x1)), vsm(ism(x2)) ) where two random effects in their own vsm() structure are encapsulated. Only applies for simple random effects.

S3 methods

S3 methods are available for some parameter extraction such as fitted.mmes, residuals.mmes, summary.mmes, randef, coef.mmes, anova.mmes, plot.mmes, and predict.mmes to obtain adjusted means. In addition, the vpredict function (replacement of the pin function) can be used to estimate standard errors for linear combinations of variance components (e.g., ratios like h2). The r2 function calculates reliability.

Additional Functions

Additional functions for genetic analysis have been included such as relationship matrix building (A.mat, D.mat, E.mat, H.mat), build a genotypic hybrid marker matrix (build.HMM), plot of genetic maps (map.plot), and manhattan plots (manhattan). If you need to build a pedigree-based relationship matrix use the getA function from the pedigreemm package.

Bug report and contact

If you have any technical questions or suggestions please post it in https://stackoverflow.com or https://stats.stackexchange.com for the community to help you.

If you have any bug report please go to https://github.com/covaruber/sommer or send me an email to address it asap, just make sure you have read the vignettes carefully before sending your question. If you need professional consulting support visit https://covaruber.github.io/sommer/

Example Datasets

The package has been equiped with several datasets to learn how to use the sommer package:

* DT_halfdiallel, DT_fulldiallel and DT_mohring datasets have examples to fit half and full diallel designs.

* DT_h2 to calculate heritability

* DT_cornhybrids and DT_technow datasets to perform genomic prediction in hybrid single crosses

* DT_wheat dataset to do genomic prediction in single crosses in species displaying only additive effects.

* DT_cpdata dataset to fit genomic prediction models within a biparental population coming from 2 highly heterozygous parents including additive, dominance and epistatic effects.

* DT_polyploid to fit genomic prediction and GWAS analysis in polyploids.

* DT_gryphon data contains an example of an animal model including pedigree information.

* DT_btdata dataset contains an animal (birds) model.

* DT_legendre simulated dataset for random regression model.

* DT_sleepstudy dataset to know how to translate lme4 models to sommer models.

* DT_ige dataset to show how to fit indirect genetic effect models.

Models Enabled

For details about the models enabled and more information about the covariance structures please check the help page of the package (sommer).

Value

If all parameters are correctly indicated the program will return a list with the following information:

data: the dataset used in model fitting after application of the common observation mask.
obsInfo: an observation-level audit table describing the common missing-data filter. It includes the original row index, validity indicators for response, fixed, random, and residual model components, the final inclusion indicator, and the reason for exclusion when applicable.
Dtable: the table to be used for the predict function to help the program recognize the factors available.
llik: the vector of log-likelihoods across iterations
b: the vector of fixed effect.
u: the vector of random effect.
bu: the vector of fixed and random effects together.
Ci: the inverse of the coefficient matrix.
Ci_11: the inverse of the coefficient matrix pertaining to the fixed effects.
theta: a list of estimated variance covariance matrices. Each element of the list corresponds to the different random and residual components
covPar: the fitted covariance parameters in their descriptor-defined reported coordinates.
pcgMatrixFree: logical indicating whether the matrix-free Kronecker PCG path was used during optimization.
residualCovarianceCompact: logical indicating that a complete, unweighted, non-diagonal Kronecker residual covariance was evaluated directly from its factors without materializing the full covariance matrix. In this case the residual element of theta is a $0 \\times 0$ matrix; use covparams_mmes() for fitted residual covariance parameters.
covParNative: a data frame containing all fixed and estimated covariance parameters in the native model-facing scale defined by each CovarianceFactor. For example, dsm() entries are level-specific variances rather than variance ratios. This is the same result returned by covparams_mmes().
covParNativeSE: the native-scale covariance-parameter table augmented with delta-method standard errors and Z ratios. This is the same result returned by covparams_mmes_se().
theta_se: inverse of the information matrix.
InfMat: information matrix.
monitor: The values of the variance-covariance components across iterations during the REML estimation.
AIC: Akaike information criterion
BIC: Bayesian information criterion
convergence: a TRUE/FALSE statement indicating if the model converged.
partitions: a list where each element contains a matrix indicating where each random effect starts and ends.
partitionsX: a list where each element contains a matrix indicating where each fixed effect starts and ends.
percDelta: the matrix of percentage change in deltas (see tolParConvNorm argument).
normMonitor: the matrix of the three norms calculated (see tolParConvNorm argument).
toBoundary: the matrix of variance components that were forced to the boundary across iterations.
Cchol: the Cholesky decomposition of the coefficient matrix.
y: the response vector.
W: the column binded matrix W = [X Z]
uList: a list containing the BLUPs in data frame format where rows are levels of the random effects and column the different factors at which the random effect is fitted. This is specially useful for diagonal and unstructured models.
uPevList: a list containing the BLUPs in data frame format where rows are levels of the random effects and column the different factors at which the random effect is fitted. This is specially useful for diagonal and unstructured models.
args: the fixed, random and residual formulas from the mmes model.
constraints: The vector of constraints.

References

Covarrubias-Pazaran G. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 2016, 11(6): doi:10.1371/journal.pone.0156744

Jensen, J., Mantysaari, E. A., Madsen, P., and Thompson, R. (1997). Residual maximum likelihood estimation of (co) variance components in multivariate mixed linear models using average information. Journal of the Indian Society of Agricultural Statistics, 49, 215-236.

Sanderson, C., & Curtin, R. (2025). Armadillo: An Efficient Framework for Numerical Linear Algebra. arXiv preprint arXiv:2502.03000.

Gilmour et al. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51(4):1440-1450.

Author

Coded by Giovanny Covarrubias-Pazaran with contributions of Christelle Fernandez Camacho, Johan Aparicio-Arce, and Claudio Flavio 5.1 jaja (laughing in Spanish).

Examples


####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples
####=========================================####

data(DT_example, package="enhancer")
DT <- DT_example
head(DT)

####=========================================####
#### Univariate homogeneous variance models  ####
####=========================================####

## Compound simmetry (CS) model
ans1 <- mmes(Yield~Env,
             random= ~ Name + Env:Name,
             rcov= ~ units,
             data=DT)
summary(ans1)

## Known variances: solve the MME only (e.g. genetic evaluation)
sol1 <- mmes(Yield~Env,
             random= ~ Name + Env:Name,
             rcov= ~ units,
             data=DT, solveOnly=TRUE, covPar=ans1)
head(sol1$uList[[1]])

# \donttest{

####===========================================####
#### Univariate heterogeneous variance models  ####
####===========================================####
DT=DT[with(DT, order(Env)), ]
## Compound simmetry (CS) + Diagonal (DIAG) model
ans2 <- mmes(Yield~Env,
             random= ~Name + vsm(dsm(Env),ism(Name)),
             rcov= ~ vsm(dsm(Env),ism(units)),
             data=DT)
summary(ans2)

####===========================================####
####  Univariate unstructured variance models  ####
####===========================================####

ans3 <- mmes(Yield~Env,
             random=~ vsm(usm(Env),ism(Name)),
             rcov=~vsm(dsm(Env),ism(units)),
             data=DT)
summary(ans3)


# }