---
title: "Using covariance structures with sommer"
author: "Giovanny Covarrubias-Pazaran"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using covariance structures with sommer}
  %\VignetteEngine{knitr::knitr}
  %\VignetteEncoding{UTF-8}
---

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:

```{r, eval=FALSE}
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

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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:

```{r, eval=FALSE}
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.

## 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.

```{r, eval=FALSE}
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`:

```{r, eval=FALSE}
fit_environment_shape <- function(shape){
  mmes(Yield ~ Env,
  random = ~ vsm(shape, ism(DT$Name)),
  rcov = ~ vsm(ism(units)),
  data = DT,
  control = mmes_control(diagnostics = list(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 factor-analytic covariance with heterogeneous specific variances by default, or an estimated common specific variance with `variance="homogeneous"`. `rrm()` fits a genuine reduced-rank covariance with all specific variances fixed at zero.
* `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) Matrix representation and number of parameters

Every structure below is a $q \times q$ covariance shape $K$ among the $q$ levels of the factor (for example, environments). The term covariance is $\Sigma = \sigma^2 K$, so each `vsm()` term estimates **one** additional overall scale $\sigma^2$. The shape $K$ is scaled so that its reference entry is fixed ($K_{11}=1$, or the first specific variance, or the first innovation variance, depending on the structure). This removes the redundancy between $\sigma^2$ and $K$, so *the covariance factor parameters plus $\sigma^2$ equal the number of free parameters of the usual textbook covariance*. For instance, `usm()` has $q(q+1)/2 - 1$ factor parameters plus $\sigma^2$, which equals $q(q+1)/2$ in total.

The matrices are written for $q=3$ with $\sigma^2$ omitted. The counts apply to a general $q$ and are the defaults; entries fixed through `fixed=` are not estimated.

**Identity, `ism()`** (no factor parameters; 1 parameter in total, $\sigma^2$)

$$
K=\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}
$$

**Diagonal, `dsm()`** ($q-1$ variance ratios; $q$ in total). `atm()` is the same for the $q_s$ selected levels ($q_s-1$ ratios).

$$
K=\begin{pmatrix}1&0&0\\0&r_2&0\\0&0&r_3\end{pmatrix}
$$

**Compound symmetry, `csm()`** (1 correlation $\rho$; 2 in total)

$$
K=\begin{pmatrix}1&\rho&\rho\\\rho&1&\rho\\\rho&\rho&1\end{pmatrix}
$$

**Heterogeneous compound symmetry, `csm(variance="heterogeneous")`** ($1+(q-1)=q$ parameters; $q+1$ in total). Entries are $K_{ij}=\sqrt{r_ir_j}\,\rho$ for $i\neq j$, and $K_{ii}=r_i$ with $r_1=1$:

$$
K=D\,C(\rho)\,D=\begin{pmatrix}1&\rho\sqrt{r_2}&\rho\sqrt{r_3}\\\rho\sqrt{r_2}&r_2&\rho\sqrt{r_2r_3}\\\rho\sqrt{r_3}&\rho\sqrt{r_2r_3}&r_3\end{pmatrix}
$$

**Autoregressive, `ar1m()`, `ar2m()`, `ar3m()`** (order $p$ partial autocorrelations: 1, 2 or 3 parameters; $p+1$ in total; heterogeneous variance adds $q-1$ ratios). For AR(1):

$$
K=\begin{pmatrix}1&\rho&\rho^2\\\rho&1&\rho\\\rho^2&\rho&1\end{pmatrix}
$$

For AR($p$) the correlations $\rho_h$ at lag $h$ follow the Yule-Walker recursion from the partial autocorrelations, and $K_{ij}=\rho_{|i-j|}$.

**Moving average, `mam(order=m)`** ($m$ parameters; $m+1$ in total). Correlation is zero beyond lag $m$. For MA(1):

$$
K=\begin{pmatrix}1&\rho_1&0\\\rho_1&1&\rho_1\\0&\rho_1&1\end{pmatrix}
$$

**Toeplitz, `toeplitzm()`** ($q-1$ lag correlations, parameterized through partial autocorrelations; $q$ in total)

$$
K=\begin{pmatrix}1&\rho_1&\rho_2\\\rho_1&1&\rho_1\\\rho_2&\rho_1&1\end{pmatrix}
$$

**Unstructured, `usm()`** ($q(q+1)/2-1$ parameters; $q(q+1)/2$ in total, with $K_{11}=1$)

$$
K=\begin{pmatrix}1&k_{12}&k_{13}\\k_{12}&k_{22}&k_{23}\\k_{13}&k_{23}&k_{33}\end{pmatrix}
$$

**General correlation, `corgm()`** ($q(q-1)/2$ correlations; $q(q-1)/2+1$ in total; `variance="heterogeneous"` adds $q-1$ ratios)

$$
K=\begin{pmatrix}1&\rho_{12}&\rho_{13}\\\rho_{12}&1&\rho_{23}\\\rho_{13}&\rho_{23}&1\end{pmatrix}
$$

**Factor analytic, `fam(k=k)`**

$$
K=\frac{\Lambda\Lambda^\top+\Psi}{(\Lambda\Lambda^\top+\Psi)_{11}},
\qquad \Lambda\in\mathbb{R}^{q\times k},\quad \Psi=\mathrm{diag}(\psi_1,\ldots,\psi_q)
$$

The leading $k\times k$ block of $\Lambda$ is lower triangular for rotational identification, and the first specific variance is the reference. The default heterogeneous mode estimates $kq-k(k-1)/2$ loadings plus $q-1$ specific-variance ratios; with $\sigma^2$ this is $kq-k(k-1)/2+q$ in total. With `variance="homogeneous"` only the $kq-k(k-1)/2$ loadings are estimated ($\Psi=\psi I$).

**Reduced rank, `rrm(k=k)`** ($K=\Lambda\Lambda^\top$, with $\Lambda_{11}=1$, so $kq-k(k-1)/2-1$ parameters; $kq-k(k-1)/2$ in total). It is singular for $k<q$ because $\Psi=0$.

**Antedependence, `antem(order=p)`** ($K=T^{-1}DT^{-\top}$ with $T$ unit lower triangular holding the regression coefficients in the first $p$ subdiagonals and $D=\mathrm{diag}(1,d_2,\ldots,d_q)$). For $q=3$, $p=1$:

$$
T=\begin{pmatrix}1&0&0\\-\beta_{21}&1&0\\0&-\beta_{32}&1\end{pmatrix}
$$

It estimates $\sum_{i=2}^{q}\min(i-1,p)$ coefficients plus $q-1$ innovation ratios (for $q=3$, $p=1$: $2+2=4$ parameters; 5 in total).

**Spatial and user-defined structures.** `maternm()` has range and smoothness ($K_{ij}=\frac{2^{1-\nu}}{\Gamma(\nu)}(d_{ij}/\text{range})^\nu \mathcal{K}_\nu(d_{ij}/\text{range})$; up to 2 parameters, 3 in total); `sar()` and `car()` have one spatial parameter $\rho$ that controls the strength of dependence through the weights matrix $W$ (2 in total); `ownm(K=)` with a known matrix has no factor parameters (1 in total, $\sigma^2$), whereas a user function adds `length(par)` parameters.

| Structure | Factor parameters | Total with $\sigma^2$ ($q=3$) | Total ($q=5$) |
|---|---|---|---|
| `ism()` | 0 | 1 | 1 |
| `dsm()` | $q-1$ | 3 | 5 |
| `csm()` | 1 | 2 | 2 |
| `csm()` heterogeneous | $q$ | 4 | 6 |
| `ar1m()` | 1 | 2 | 2 |
| `ar1m()` heterogeneous | $q$ | 4 | 6 |
| `ar2m()` | 2 | 3 | 3 |
| `mam(order=1)` | 1 | 2 | 2 |
| `toeplitzm()` | $q-1$ | 3 | 5 |
| `corgm()` | $q(q-1)/2$ | 4 | 11 |
| `fam(k=1)` | $2q-1$ | 6 | 10 |
| `rrm(k=1)` | $q-1$ | 3 | 5 |
| `antem(order=1)` | $2(q-1)$ | 5 | 9 |
| `usm()` | $q(q+1)/2-1$ | 6 | 15 |

Use the table to judge whether the data can support a structure: the number of free parameters should be small relative to the number of observations per environment pair. Because $\sigma^2$ is shared, models differing only in these structures can be compared with `logLik` and information criteria, and with likelihood-ratio tests when they are nested (for example `csm` within heterogeneous `csm`, and heterogeneous `csm` and `fam` within `usm`).

### 3) 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.

```{r, eval=FALSE}
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.

### 4) 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:

```{r, eval=FALSE}
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)`.

### 5) 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:

```{r, eval=FALSE}
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`.

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

```{r, eval=FALSE}
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`:

```{r, eval=FALSE}
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.