mmes: *m*ixed *m*odel *e*quations *s*olver

In sommer: Solving Mixed Model Equations in R

mixed model equations solver

Description

Fits a generalized linear mixed model using model formulas, covariance structures and observation data. Special covariance structures can be specified with vsm(). Numerical fitting choices are supplied through the strictly validated control object; see mmes_control and the sommer.mmes.options vignette.

Usage

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, control=mmes_control(), .internal=NULL, ...)

Arguments

fixed	A formula specifying the response variable(s) and fixed effects, i.e: response ~ covariate By default, the fixed-effect design must have full column rank after observation filtering. Before fitting, mmes() reports all-zero columns (for example, unused levels or unobserved interaction combinations) or linearly dependent columns (for example, a parent factor together with unrestricted nested cell effects). Aliased coefficients are not automatically dropped or constrained to zero unless explicitly requested with control=mmes_control(representation=list(fixedRank="drop")). This option selects a full-rank subset using shared native sparse QR for the Henderson and direct engines. The original coefficient layout is preserved, with omitted coefficients returned as zero identification constraints and flagged in summary(fit)$betas; these zeros are not estimates of independently identifiable effects. Details, mappings and supported inference workflows are documented in mmes_control. Apply droplevels() to the retained data to remove unused factor levels. For a sparse interaction, an observed-cell factor created with interaction(..., drop=TRUE) and a formula such as response ~ 0 + cell can represent the same fitted cell means without redundant columns. Adding the parent factor again reintroduces aliasing. Dropping unused levels alone does not remove all unobserved interaction combinations. The diagnostic setup mode (control=mmes_control(output=list(setup=TRUE))) still returns the design without requiring full column rank, allowing it to be inspected. With fixedRank="drop", setup instead returns the reduced design and its fixedRank mapping.
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 a reduced-rank random covariance with specific variances fixed at zero; use fam(..., variance="homogeneous") for an estimated common specific variance; 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 control$formulation="henderson", random coefficients are represented in the eigenbasis and the diagonal relationship precision is used in the mixed model equations. With control$formulation="direct", 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 control$output$inverse=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.
control	A sommer_mmes_control object created by mmes_control. It owns numerical task, equation formulation, representation, solver, determinant, score, curvature, covariance update, proposal, convergence, output, execution and diagnostic settings. Its default is mmes_control(). Removed numerical arguments are not accepted by mmes(); configure those choices in this object. See the control help and vignette for defaults, mathematical ownership, eligibility and migration.
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.
contrasts	an optional list. See the contrasts.arg of model.matrix.default.
vcc	Optional data frame of equality/scaling constraints between covariance parameters (columns parameter, group and optionally scale), identified through the vcParams table returned by control$output$setup=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. PQL uses outer-loop controls from control$task$pql; compatibility with update and solver requests is model-dependent. For other families the reported dispersion is the Pearson estimate. The negative binomial \theta is estimated (conditional on the random effects) with control$task$pql=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.
covPar	Known covariance parameters for the solve task selected by control$task$mode="solve". Accepted forms are a compatible fitted mmes object, a named list of term parameters, or a parameter data frame. When NULL, covariance values must be fixed in the model.
.internal	Reserved for package-internal recursive fitting state. Users should not set this argument.
...	Deprecated legacy arguments. The supported arguments are getPEV, tolParInv, and verbose; they map respectively to control$output$inverse, control$execution$inverseTol, and control$diagnostics$verbose. These arguments override matching settings in control and issue a deprecation warning.

Details

Structurally unestimable REML components

For Gaussian identity-link REML fits, mmes() checks whether an independently scaled random component contributes anything outside the fixed-effect design space, using the observations retained by the common missing-data filter. A component whose entire random design lies in that space is invisible to the restricted likelihood. Such components are removed before optimization, with a warning, rather than assigned an arbitrary variance by a regularized information-matrix solve.

The conservative check handles a whole random term with a free variance scale and fixed covariance shape, or individual levels of a single dsm() or atm() shaping factor when its scale and all variance ratios are free. For example, with fixed locale means, a locale-specific independent block variance is omitted when that locale has only one observed block. A shared block variance is retained if other locales provide information about it. Correlated, multiply shaped, constrained or partly fixed structures are not automatically reduced at individual levels. The check does not diagnose general aliasing between covariance terms or weak identification near a variance boundary.

Automatic omission is not applied to ML, PQL, known-covariance solve tasks, vcc constraints, observation rotation or factor-score augmentation. For supported fits, omitted components retain their original positions in random-effect and covariance outputs, with zero covariance and zero random effects. These zeros denote an omission convention, not variances estimated to be zero. Their native-scale standard errors and Z-ratios are NA; summary() and covparams_mmes() label them "unestimable (omitted)". Only retained free parameters contribute to optimization and information criteria. Omission preserves the restricted likelihood information about retained parameters; it does not establish zero between-block variation in the population. Inference on fixed means and predictions is conditional on the omission convention.

Fixed factors with many levels

In Henderson fits the fixed-effect design X is stored as a sparse matrix. The coefficient matrix C of the mixed model equations is factorized with sparse methods at every REML iteration, whatever the control$output$inverse setting. The fill-reducing ordering of the sparse factorization eliminates the fixed-effect equations within the same factor as the random ones, which is equivalent to absorbing them. A fixed factor with thousands of levels, such as contemporary groups, does not create a dense fixed-effect block during estimation. Variance components, BLUEs and BLUPs are identical to those of an explicit absorption of the fixed effects.

control$output$inverse only controls work done after convergence:

0 (default)

No inverse of C is formed. Fixed-effect standard errors in summary() are NA.

1

Elements of the inverse are computed only within the sparsity pattern of the factor (Takahashi selected inverse), giving random-effect PEVs. The fixed part stays sparse. Fixed-effect standard errors in summary() remain NA.

2

The complete inverse of C, including the fixed block, is formed as a dense matrix. This needs about 8(p+q)^2 bytes for p fixed and q random effects (around 1.4 GB for 13,000 effects) and is the mode that reports fixed-effect standard errors in summary().

With many fixed levels, keep inverse=0 or 1. When standard errors of particular fixed-effect combinations are needed, use predict.mmes, which computes exact standard errors from the sparse C without forming Ci. postPEV can add PEVs or the full inverse after fitting.

When even the sparse REML fit is too slow, a practical alternative is to estimate the variance components with the large factor fitted as random, optionally on a well-connected subset of the data. Then refit with the factor as fixed using control=mmes_control(task=list(mode="solve")) and covPar to obtain BLUEs and BLUPs with a single factorization. Simply omitting the factor during REML is not equivalent: its variation is then attributed to the residual and to confounded random terms.

Numerical controls

The opt-in control$representation$equations="conditional" path fits Gaussian FA models using sparse latent equations with exact REML AI updates by default, or conditional-moment EM through control$update$method="em". The earlier control$update$method="conditional-fa" selector retains its EM behavior. Both optimizers support incomplete and unbalanced observations. It accepts one FA relationship term and diagonal or FA residuals. Matching fam(shared="group") tags request common trait loadings across genetic and residual sources. The default AI path is unchanged. See mmes_control and the conditional FA section of vignette("sommer.mmes.options") for eligibility and output restrictions. Conditional fits retain factorized covariance descriptors, not dense theta matrices, public C, coefficient uncertainty, or covariance-parameter uncertainty. Use covmatrix_mmes(..., se=FALSE) to explicitly reconstruct a bounded-size covariance.

mmes() keeps model specification separate from numerical policy. The control argument defaults to mmes_control(), whose modules have strict schemas: unknown field names are errors. The default configuration retains the former standard bundle: Henderson equations, REML, exact AI curvature, EM recovery, no extrapolation, 30 inner iterations, full covariance information, coefficient-matrix retention, and inherited solver/resource defaults. The mmes_control help gives every field, default, module owner and validation rule.

Ownership and eligibility

Task

control$task selects fitting or known-covariance solving, REML or ML, and outer PQL controls.

Formulation

control$formulation selects Henderson or direct observation-space equations.

Representation and solve

control$representation describes equation/covariance representation; control$solve owns linear solver, preconditioning, reuse and elimination requests.

Evaluation

control$logdet owns determinant methods; control$score owns score traces and covariance checks; control$curvature owns exact or hybrid AI and reduced-direction work.

Optimization

control$update owns EM/PXEM recovery, blending and common-scale profiling. control$proposal owns step scaling and Aitken extrapolation. Aitken and hybrid AI are mutually exclusive.

Convergence and outputs

control$convergence owns inner iteration, tolerances and stochastic accuracy; control$output owns coefficient/covariance uncertainty and setup output.

Execution and diagnostics

control$execution owns memory, hardware, threading, precision and seeds. control$diagnostics owns verbosity, package-age warnings and profiling.

These modules do not form a Cartesian product of interchangeable algorithms. Model-dependent admission remains authoritative; a requested mode may be ineligible. Consult mmes_control for prerequisites, exclusions and current behavior. No new method combinations are implied by the modular API. The internal .mmes_capabilities() catalog is descriptive and does not resolve eligibility or predict execution.

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.
remlOmissions	when present, a data frame identifying structurally unestimable random covariance components omitted from REML, with columns term, parameter and reason. It is also available in engineDiagnostics$remlOmissions.
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 effects on the original response scale.
u	the vector of random effects on the original response scale, consistent with uList.
bu	the vector of fixed and random effects together on the original response scale, formed by concatenating b and u.
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.
control	the validated modular numerical controls used for the fit. See mmes_control.
pcgRecycleState	present only when pcg$exportRecycle=TRUE. A versioned serializable Krylov basis with coefficient-coordinate metadata, for explicit import into a compatible PCG fit.
covarianceInformation	"none" when covariance information was explicitly omitted. The default retains full covariance information.
retainCoefficientMatrix	FALSE when an opt-in C-free saddle-point fit omitted the final compatibility matrix. In that mode C is NULL, and C-based post-fit methods require a refit retaining C. The default retains C.
engineDiagnostics	path-dependent numerical diagnostics. compression reports method, original observations, groups, residual classes, retained statistic bytes and setup time. pcg reports iterative work, reuse and stochastic accuracy diagnostics for REML PCG. profile contains detailed phase timings when explicitly enabled. Retained statistic bytes are not process peak memory.
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. The table includes delta-method StdError and Zratio columns when covariance information is retained, matching covparams_mmes_se(). With covarianceInformation="none" these columns are NA; covparams_mmes() still returns point estimates.
theta_se	covariance-parameter uncertainty matrix; NULL when covarianceInformation="none".
InfMat	information matrix; NULL when covarianceInformation="none".
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 control$convergence$parameterTol).
normMonitor	the matrix of the three norms calculated (see control$convergence$parameterTol).
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.

Author(s)

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

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.

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, control=mmes_control(task=list(mode="solve")),
             covPar=ans1)
head(sol1$uList[[1]])



####===========================================####
#### 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= ~ dsumm(vsm(ism(units)), by=Env),
             data=DT)
summary(ans2)

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

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






## simulated MET
nInds=300
nMarks=1000
nEnvs=3

#random population of nInds lines with nMarks markers
M <- matrix(rep(0,nInds*nMarks),nInds,nMarks)
for (i in 1:nInds) {
  M[i,] <- ifelse(runif(nMarks)<0.5,-1,1)
}
# compute the relationship matrix
MMT <- tcrossprod(M)
m <- sum(diag(MMT))/nrow(MMT)
MMT <- MMT/m
MMT <- MMT + diag(1e-05, ncol(MMT), ncol(MMT))
# get phenotypes
rg <- enhancer::simGECorMat(nEnvs, 1, mu=0.2, v=0.2)
rg <- Matrix::nearPD(rg)$mat
rg # desired genetic correlation
u <- MASS::mvrnorm(n = nMarks, rep(0, nEnvs), Sigma=rg) # marker effects
g <- crossprod(t(M),u) # genetic values
h2 <- 0.2  #heritability 
envEffect <- rnorm(nEnvs, mean = abs(rnorm(1))) # environment means
e <- MASS::mvrnorm(n = nInds, envEffect, Sigma=(1-h2)/h2*diag(diag(var(g)))) # errors
p = g + e # phenotypes
y <- as.vector(p) # phenotypes as vector
env <- sort(rep(paste0("e",1:nEnvs), nrow(g))) # environment labels
ind <- rep(paste0("i",1:nrow(g)), nEnvs ) # individuals labels
colnames(MMT) <- rownames(MMT) <- paste0("i",1:nrow(g)) # assign dimnames

MMTi <- solve(MMT + diag(1e-4,ncol(MMT),ncol(MMT)))
MMTi[lower.tri(MMTi)] <- t(MMTi)[lower.tri(MMTi)] # fill the lower triangular
MMTi <- as(as(as( MMTi,  "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(MMTi, 'inverse')=TRUE

DT <- data.frame(ind=ind, env=env, y=y)
# fit the models
mixm2 <- mmes(y ~ env,
              random=~ vsm(dsm(env), ism(ind), Gu=MMTi),
              data=DT
)
summary(mixm2)$varcomp

sommer documentation built on Oct. 25, 2026, 5:06 p.m.