Using covariance structures with sommer

The sommer package fits mixed models with structured covariance models for random effects and residuals. In mmes(), covariance structures are supplied through vsm(). This vignette explains the relation between vsm(), the covariance constructors, and the internal CovarianceFactor descriptor used by the Henderson mixed-model solver.

SECTION 1: Covariance structures in vsm()

1. The product-level variance scale
2. Covariance-shaping factors
3. Homogeneous and heterogeneous marginal variances

SECTION 2: Fitting structured models

1. Compound symmetry
2. Autoregressive covariance
3. Fixing covariance parameters

SECTION 3: CovarianceFactor descriptors

1. Complete covariance-structure catalog
2. Combining two simple effects with covm()

SECTION 4: CovarianceFactor descriptors

1. Descriptor fields
2. Working and reported parameters

SECTION 1: Covariance structures in vsm()

1) The product-level variance scale

Each vsm() term owns one overall variance parameter, sigma2. The covariance constructors inside vsm() define dimensionless covariance shapes. If the supplied factors are \(K_1,\ldots,K_m\), the covariance represented by one vsm() term is

$$ \Sigma = \sigma^2(K_1 \otimes K_2 \otimes \cdots \otimes K_m). $$

This convention avoids confounding several absolute variance scales in the same Kronecker product. The sigma2 argument supplies a starting value, and fixedSigma2 = TRUE fixes that overall scale.

2) Covariance-shaping factors

The final argument to vsm() is the main-effect incidence term. All preceding arguments are covariance-shaping factors. For example, the following structure has an environment covariance factor, an AR(1) row covariance factor, and an identity covariance among genotypes:

vsm(
  dsm(Env),
  ar1m(Row),
  ism(Name)
)

The factors are combined from left to right. Earlier factors index the slow dimension of the Kronecker product and later factors index the fast dimension. A known relationship precision matrix can be supplied through Gu for the levels of the final main-effect term.

Common covariance constructors include ism() for an identity shape, dsm() for diagonal relative variances, csm() for compound symmetry, ar1m() through ar3m() for autoregressive correlations, usm() for an unstructured covariance, and corgm() for a general correlation matrix.

3) Homogeneous and heterogeneous marginal variances

Correlation structures support a common variance argument. The default is "homogeneous". For compound symmetry with correlation matrix \(C(\rho)\), the two modes are

$$ K = C(\rho) $$

and

$$ K = D C(\rho) D, \qquad D = \mathrm{diag}(1,\sqrt{r_2},\ldots,\sqrt{r_q}). $$

The heterogeneous mode estimates \(q-1\) positive variance ratios. The first level is the reference; its absolute variance is represented by vsm(..., sigma2 = ...). The same convention applies to ar1m(), ar2m(), and ar3m().

SECTION 2: Fitting structured models

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

1) Compound symmetry

The following model fits a homogeneous compound-symmetry covariance among environments for genotype effects:

fit_cs <- mmes(
  Yield ~ Env,
  random = ~ vsm(csm(Env), ism(Name)),
  rcov = ~ vsm(ism(units)),
  data = DT,
  verbose = FALSE
)

To estimate environment-specific variance ratios while retaining the same common correlation, select heterogeneous variance mode. values supplies positive starting marginal variances:

n_env <- nlevels(factor(DT$Env))

fit_csh <- mmes(
  Yield ~ Env,
  random = ~ vsm(
    csm(Env, variance="heterogeneous", values=rep(1, n_env)),
    ism(Name)
  ),
  rcov = ~ vsm(ism(units)),
  data = DT,
  verbose = FALSE
)

The homogeneous and heterogeneous models have one overall random-effect variance from vsm(). The heterogeneous model additionally estimates one common correlation and \(q-1\) environment variance ratios.

2) Autoregressive covariance

For ordered levels, use an autoregressive correlation. The following example fits an AR(1) covariance with heterogeneous marginal variances:

DT$EnvOrder <- factor(DT$Env, levels=unique(DT$Env), ordered=TRUE)
n_env <- nlevels(DT$EnvOrder)

fit_ar1h <- mmes(
  Yield ~ Env,
  random = ~ vsm(
    ar1m(EnvOrder, rho=0.30, variance="heterogeneous",
         values=rep(1, n_env)),
    ism(Name)
  ),
  rcov = ~ vsm(ism(units)),
  data = DT,
  verbose = FALSE
)

For AR(2) and AR(3), use ar2m() or ar3m() and provide starting partial autocorrelations through pacf. Their heterogeneous mode uses the same variance and values arguments.

3) Fixing covariance parameters

The fixed argument belongs to the covariance constructor and fixes its factor parameters. In a heterogeneous compound-symmetry model, its order is the common correlation followed by the \(q-1\) variance ratios. The following keeps all environment variance ratios at one while estimating the common correlation:

n_env <- nlevels(factor(DT$Env))

environment_cs <- csm(
  DT$Env,
  variance="heterogeneous",
  values=rep(1, n_env),
  fixed=c(FALSE, rep(TRUE, n_env - 1L))
)

To fix a residual variance to one, fix the product-level vsm() scale rather than a covariance-factor parameter:

fit_fixed_residual <- mmes(
  Yield ~ Env,
  random = ~ vsm(ism(Name)),
  rcov = ~ vsm(ism(units), sigma2=1, fixedSigma2=TRUE),
  data = DT,
  verbose = FALSE
)

For a heterogeneous residual structure, fixedSigma2=TRUE fixes the reference residual variance. Fixing the associated variance-ratio entries to TRUE fixes the other marginal residual variances relative to that reference.

SECTION 3: Complete covariance-structure catalog

The following catalog uses the same DT_example data throughout. It creates an ordered environment factor for time-series structures, a one-dimensional environment coordinate for the Matérn structure, and a named chain adjacency matrix for SAR and CAR structures. The code is shown with eval=FALSE because fitting every candidate model is illustrative rather than a recommended model-selection workflow.

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

env_levels <- unique(as.character(DT$Env))
DT$EnvOrder <- factor(DT$Env, levels=env_levels, ordered=TRUE)
DT$EnvCoordinate <- as.numeric(DT$EnvOrder)
n_env <- length(env_levels)

# DT_example has three environments. AR(3) needs at least four ordered
# levels, so this demonstration-only partition is used for the AR(3) rows
# below. Replace it with a scientific time, distance, or ordered factor.
DT$CatalogOrder <- factor(rep(seq_len(4L), length.out=nrow(DT)), ordered=TRUE)
n_catalog <- nlevels(DT$CatalogOrder)

# Named first-neighbour adjacency among ordered environments.
W_env <- matrix(0, n_env, n_env,
                dimnames=list(env_levels, env_levels))
W_env[cbind(seq_len(n_env - 1L), 2:n_env)] <- 1
W_env[cbind(2:n_env, seq_len(n_env - 1L))] <- 1

# A known positive-definite covariance shape for ownm().
K_env <- 0.40 ^ abs(outer(seq_len(n_env), seq_len(n_env), "-"))

# Every object below has the same observation layout and can be used as
# a covariance factor in vsm(shape, ism(Name)). Most use EnvOrder; AR(3)
# uses CatalogOrder because EnvOrder has too few levels for that structure.
environment_shapes <- list(
  identity = ism(DT$EnvOrder),
  diagonal = dsm(DT$EnvOrder),
  selected_diagonal = atm(DT$EnvOrder, levs=env_levels[1:3]),
  compound_symmetry = csm(DT$EnvOrder),
  compound_symmetry_heterogeneous = csm(
    DT$EnvOrder, variance="heterogeneous", values=rep(1, n_env)
  ),
  ar1 = ar1m(DT$EnvOrder),
  ar1_heterogeneous = ar1m(
    DT$EnvOrder, variance="heterogeneous", values=rep(1, n_env)
  ),
  ar2 = ar2m(DT$EnvOrder),
  ar2_heterogeneous = ar2m(
    DT$EnvOrder, variance="heterogeneous", values=rep(1, n_env)
  ),
  ar3 = ar3m(DT$CatalogOrder),
  ar3_heterogeneous = ar3m(
    DT$CatalogOrder, variance="heterogeneous", values=rep(1, n_catalog)
  ),
  ma1 = mam(DT$EnvOrder, order=1L),
  ma2 = mam(DT$EnvOrder, order=2L),
  unstructured = usm(DT$EnvOrder),
  general_correlation = corgm(DT$EnvOrder),
  factor_analytic = fam(DT$EnvOrder, k=1L),
  antedependence = antem(DT$EnvOrder, order=1L),
  user_defined = ownm(DT$EnvOrder, K=K_env),
  reduced_rank = rrm(DT$EnvOrder, k=1L),
  matern = maternm(DT$EnvCoordinate),
  toeplitz = toeplitzm(DT$EnvOrder),
  sar = sar(DT$EnvOrder, W=W_env),
  car = car(DT$EnvOrder, W=W_env)
)

mam(order=1L) and mam(order=2L) are the canonical moving-average interfaces; ma1m() and ma2m() are convenience wrappers. atm() is a selected-level diagonal structure: observations outside levs have zero incidence in that covariance factor, so it is appropriate only when that selected-level interpretation is intended. The AR(3) entries use CatalogOrder only because DT_example has three environment levels; an AR(3) analysis requires at least four scientifically meaningful ordered levels.

The common fitting pattern is identical for all entries in environment_shapes:

fit_environment_shape <- function(shape){
  mmes(
    Yield ~ Env,
    random = ~ vsm(shape, ism(DT$Name)),
    rcov = ~ vsm(ism(units)),
    data = DT,
    verbose = FALSE
  )
}

fit_identity <- fit_environment_shape(environment_shapes$identity)
fit_ar1 <- fit_environment_shape(environment_shapes$ar1)
fit_matern <- fit_environment_shape(environment_shapes$matern)
fit_car <- fit_environment_shape(environment_shapes$car)

The constructors differ in their assumptions and parameters:

• ism() has no factor parameters and specifies independent, equal-variance levels.
• dsm() estimates positive relative variances; atm() does the same for a chosen subset of levels.
• csm() estimates a common correlation, with optional heterogeneous variance ratios through variance="heterogeneous".
• ar1m(), ar2m(), and ar3m() use stationary ordered-level correlations. AR(2) and AR(3) use partial autocorrelations, and all AR functions offer homogeneous or heterogeneous marginal variances.
• mam() specifies an MA(1) or MA(2) correlation with zero correlation beyond its order.
• usm() estimates a fully unstructured positive-definite covariance, whereas corgm() estimates a general correlation matrix and leaves marginal scale to vsm().
• fam() fits a reduced factor-analytic covariance plus specific variances; rrm() fits a reduced-rank covariance with an identity remainder.
• antem() uses a modified-Cholesky antedependence covariance with ordered levels.
• ownm() accepts either a known positive-definite matrix through K or a user-supplied covariance function.
• maternm() uses numeric spatial coordinates and estimates range and smoothness.
• toeplitzm() estimates a general stationary Toeplitz correlation using partial autocorrelations.
• sar() and car() use a named spatial weights matrix. SAR accepts a general square weights matrix; CAR requires a symmetric, non-negative, zero-diagonal adjacency matrix with no isolated levels.

The choice should be driven by the scientific design: use ordered structures only when the level order is meaningful, spatial structures only with defensible coordinates or adjacency, and flexible structures such as usm() or corgm() only when the data support their larger number of parameters.

2) Combining two simple effects with covm()

covm() is not a covariance-shaping factor. It combines two simple vsm() random-effect structures that use the same main-effect levels and the same relationship precision matrix. The combined effect has a normalized \(2 \times 2\) unstructured covariance between the two effects and one product-level scale.

effect_1 <- vsm(ism(DT$Name))
effect_2 <- vsm(ism(DT$Name))
joint_effect <- covm(effect_1, effect_2, labels=c("effect_1", "effect_2"))

Use covm() when two effects must be correlated. For crossed covariance factors such as environments by genotypes, use multiple factors directly inside vsm() instead.

3) Covariance across several terms with strm()

strm() generalizes covm() to any number of random terms and any covariance constructor over the terms (ASReml str()). All terms must share the coefficient levels, the relationship matrix, and the same inner covariance factors; their variance is \(\sigma^2 K_{terms}\otimes K_{inner}\otimes A\). A direct-maternal-permanent environment animal model is:

fit <- mmes(y ~ 1,
            random = ~ strm(dir = vsm(ism(id)), mat = vsm(ism(dam)),
                            pe = vsm(ism(pe)), cov = usm, Gu = Ainv),
            data = animals)
covparams_mmes(fit, 1)  # variances and covariances among dir, mat and pe

cov can be any constructor applied to the term index, e.g. cov = dsm (independent terms), cov = corgm, or cov = function(x) fam(x, k = 1). Inner factors are shared, e.g. strm(vsm(dsm(env), ism(id)), vsm(dsm(env), ism(dam)), Gu = Ainv).

4) Equality and scaling constraints with vcc

Covariance parameters can be constrained to be equal or to keep fixed ratios (ASReml vcc). Parameters are identified in the vcParams table:

p <- mmes(Y ~ V * N, random = ~ B + B:MP, rcov = ~ units, data = DT_yatesoats,
          returnParam = TRUE)
p$vcParams
fit <- mmes(Y ~ V * N, random = ~ B + B:MP, rcov = ~ units, data = DT_yatesoats,
            vcc = data.frame(parameter = c("vsm(ism(B:MP)):sigma2", "vsm(ism(B)):sigma2"),
                             group = 1, scale = c(1, 2)))

Here the block variance is twice the whole-plot variance. Variance scales can only be grouped with other variance scales, correlation-type parameters can only be equated, and a fixed member fixes its whole group. See ?vcc for the rules and the limitations with respect to ASReml’s vcm.

SECTION 4: CovarianceFactor descriptors

1) Descriptor fields

Every covariance constructor returns a list containing an incidence matrix Z and a compiled covFactor descriptor. vsm() combines those descriptors, adds log(sigma2) as the first optimizer parameter, and passes the resulting structure to mmes().

structure <- csm(DT$Env, variance="heterogeneous")
str(structure$covFactor)

The principal CovarianceFactor fields are:

• dim and levels: covariance dimension and the matching level order.
• par: starting values on an unconstrained working scale.
• free: logical indicators specifying which entries of par are estimated.
• par_names: readable names for the covariance parameters.
• evaluator: the function or native operation that evaluates the covariance shape.
• derivative: analytic or numerical derivative information used by the AI REML algorithm.
• report: transformations used to report natural-scale parameters.
• trust_cap: maximum proposal sizes for optimizer coordinates.
• structurally_diagonal: whether the covariance shape is always diagonal.

These fields are created and validated by the covariance constructors. They are useful for inspecting a model, but users should normally select a constructor and its documented arguments rather than edit a descriptor manually.

2) Working and reported parameters

The optimizer works on unconstrained coordinates. For example, a correlation in \((-1,1)\) is represented using atanh(rho), positive variance ratios use logarithms, and compound-symmetry correlations use a bounded-logit transformation to remain in their positive-definite interval. The report component stores the inverse transformation so model summaries can display natural-scale correlations and variance ratios.

The resulting vsm() object stores the flattened covariance descriptor in covStruct:

random_structure <- vsm(csm(DT$Env, variance="heterogeneous"), ism(DT$Name))
random_structure$covStruct$par_names
random_structure$covStruct$free

This separation between a single product-level scale and normalized covariance factors makes homogeneous and heterogeneous structures identifiable, composable, and usable in the same mmes() interface.

Literature

Covarrubias-Pazaran G. 2016. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 11(6):1-15.

Gilmour AR, Thompson R, and Cullis BR. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51:1440-1450.