Using covariance structures with sommer

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

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,
  control = mmes_control(diagnostics = list(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,
  control = mmes_control(diagnostics = list(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,
  control = mmes_control(diagnostics = list(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,
  control = mmes_control(diagnostics = list(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.

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,
  control = mmes_control(output = list(setup = 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.

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.