Package 'sommer'

Description

Sommer is a structural univariate and multivariate linear mixed-model package for fitting models with multiple random effects and flexible covariance structures. Variance parameters are estimated by restricted maximum likelihood (REML). The package provides two complementary numerical formulations.

mmes uses Henderson's mixed model equations with an Average Information REML algorithm. The mixed-model coefficient matrix is represented as a sparse system and factorized with Eigen/Armadillo routines. Sparse LDLT factorizations are reused throughout the REML calculations. Selected inverse elements required for likelihood derivatives and prediction error variances can be obtained with Takahashi sparse-inverse recursions, so a complete dense coefficient-matrix inverse is not routinely formed.

mmer provides the marginal-covariance MNR formulation. Its REML calculations are organized around the covariance matrix of the observations and the corresponding REML projection matrix, with Newton-Raphson/Fisher-scoring and Average Information updates for variance parameters.

Both formulations support multiple random effects and structured covariance models. Through vsm, covariance structures can be composed from an arbitrary number of covariance factors, including identity, diagonal, unstructured, autoregressive, moving-average, Toeplitz, factor-analytic, reduced-rank, antedependence, correlation and spatial structures, as well as user-defined covariance functions.

The marginal formulation in mmer performs its principal REML calculations in observation space, whereas mmes performs them through the sparse mixed-model coefficient system. Their relative computational performance therefore depends on the number of observations, the number of mixed-model coefficients, sparsity, and the covariance structures in the fitted model.

The numerical algorithms are coded primarily in C++ using Armadillo and Eigen. Sommer returns REML variance-covariance estimates, BLUEs, BLUPs, residuals, fitted values, information matrices and, when requested, prediction error variance information.

Functions for genetic analysis

The package provides kernels to estimate additive (A.mat), dominance (D.mat), epistatic (E.mat), single-step (H.mat) relationship matrices for diploid and polyploid organisms. It also provides flexibility to fit other genetic models such as full and half diallel models and random regression models.

A good converter from letter code to numeric format is implemented in the function atcg1234, which supports higher ploidy levels than diploid. Additional functions for genetic analysis have been included such as build a genotypic hybrid marker matrix (build.HMM), plot of genetic maps (map.plot), creation of manhattan plots (manhattan). If you need to use pedigree you need to convert your pedigree into a relationship matrix (use the 'getA' function from the pedigreemm package).

Functions for statistical analysis and S3 methods

The vpredict function can be used to estimate standard errors for linear combinations of variance components (e.g. ratios like h2). The r2 function calculates reliability. S3 methods are available for some parameter extraction such as:

+ predict.mmes

+ fitted.mmes

+ residuals.mmes

+ summary.mmes

+ coef.mmes

+ anova.mmes

+ plot.mmes

Functions for trial analysis

Recently, spatial modeling has been added added to sommer using the two-dimensional spline (spl2Dc).

Keeping sommer updated

The sommer package is updated on CRAN every 4-months due to CRAN policies but you can find the latest source at https://github.com/covaruber/sommer. This can be easily installed typing the following in the R console:

library(devtools)

install_github("covaruber/sommer")

This is recommended if you reported a bug, was fixed and was immediately pushed to GitHub but not in CRAN until the next update.

Tutorials

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

vignette("sommer.qg")

vignette("sommer.gxe")

vignette("sommer.vs.lme4")

vignette("sommer.spatial")

Getting started

The package has been equiped with several datasets to learn how to use the sommer package (and almost to learn all sort of quantitative genetic analysis):

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

Differences of sommer >= 4.4.1 with previous versions

Since version 4.4.1, I have unified the use of the two different solving algorithms into the mmes function by just using the new argument henderson which by default is set to FALSE. Other than that the rest is the same with the addition that now the identity terms needs to be encapsulated in the ism function. In addition, now the multi-trait models need to be fitted in the long format. This are few but major changes to the way sommer models are fitted.

Differences of sommer >= 4.1.7 with previous versions

Since version 4.1.7 I have introduced the mmes-based average information function 'mmec' which is much faster when dealing with the r > c problem (more records than coefficients to estimate). This introduces its own covariance structure functons such as vsc(), usc(), dsc(), atc(), csc(). Please give it a try, although is in early phase of development.

Differences of sommer >= 3.7.0 with previous versions

Since version 3.7 I have completly redefined the specification of the variance-covariance structures to provide more flexibility to the user. This has particularly helped the residual covariance structures and the easier combination of custom random effects and overlay models. I think that although this will bring some uncomfortable situations at the beggining, in the long term this will help users to fit better models. In esence, I have abandoned the asreml formulation (not the structures available) given it's limitations to combine some of the sommer structures but all covariance structures can now be fitted using the 'vsm' functions.

Differences of sommer >= 3.0.0 with previous versions

Since version 3.0 I have decided to focus in developing the multivariate solver and for doing this I have decided to remove the M argument (for GWAS analysis) from the mmes function and move it to it's own function GWAS.

Before the mmes solver had implemented the usm(trait), diag(trait), at(trait) asreml formulation for multivariate models that allow to specify the structure of the trait in multivariate models. Therefore the MVM argument was no longer needed. After version 3.7 now the multi-trait structures can be specified in the Gt and Gtc arguments of the vsm function.

The Average Information algorithm had been removed in the past from the package because of its instability to deal with very complex models without good initial values. Now after 3.7 I have brought it back after I noticed that starting with NR the first three iterations gives enough flexibility to the AI algorithm.

Keep in mind that sommer uses direct inversion (DI) algorithm which can be very slow for datasets with many observations (big 'n'). The package is focused in problems of the type p > n (more random effect(s) levels than observations) and models with dense covariance structures. For example, for experiment with dense covariance structures with low-replication (i.e. 2000 records from 1000 individuals replicated twice with a covariance structure of 1000x1000) sommer will be faster than MME-based software. Also for genomic problems with large number of random effect levels, i.e. 300 individuals (n) with 100,000 genetic markers (p). On the other hand, for highly replicated trials with small covariance structures or n > p (i.e. 2000 records from 200 individuals replicated 10 times with covariance structure of 200x200) asreml or other MME-based algorithms will be much faster and I recommend you to use that software.

Models Enabled

General linear mixed model

Both numerical engines fit the Gaussian linear mixed-model family

$y = X\beta + Zu + e,$

with

$u \sim N(0,G), \qquad e \sim N(0,R),$

and therefore

$V = \mathrm{Var}(y) = ZGZ^\prime + R.$

For several independent random terms this becomes

$V = \sum_{k=1}^{q} Z_k G_k Z_k^\prime + R.$

Here $y$ is the response vector, $X$ and $Z$ are design matrices, $\beta$ contains fixed effects, $u$ contains random effects, and $e$ contains residual effects. The same representation covers univariate and multivariate models after the corresponding responses, design matrices and covariance structures are assembled.

Covariance models specified through vsm use a product-level variance scale and an arbitrary number of covariance-shaping factors:

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

For a random effect with known relationship covariance $A$,

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

The CovarianceFactor interface separates evaluation of each $K_j$, its derivatives, parameter reporting and trust limits from the REML solver. New covariance structures can therefore be added without changing the core mmes optimization algorithm.

Marginal MNR formulation in mmer

The MNR engine used by mmer organizes the REML calculations around the marginal covariance $V$. Define the REML projection matrix

$P = V^{-1} - V^{-1}X(X^\prime V^{-1}X)^{-1}X^\prime V^{-1}.$

Apart from constants independent of the variance parameters, the Gaussian REML log-likelihood is

$\ell_R(\theta) = -\frac{1}{2}\{\log|V|+\log|X^\prime V^{-1}X|+y^\prime P y\}.$

For variance parameter $\theta_i$, let

$V_i = \frac{\partial V}{\partial\theta_i}.$

For covariance models linear in the fitted variance parameters, the REML score has the standard form

$s_i = \frac{1}{2}\{y^\prime P V_i P y-\mathrm{tr}(P V_i)\}.$

The MNR C++ implementation forms the projected covariance derivatives $PV_i$. Its Average Information matrix is

$AI_{ij} = \frac{1}{2}y^\prime P V_i P V_j P y,$

while its Newton-Raphson/Fisher-scoring path uses trace products of the projected covariance derivatives. The resulting information system determines the variance-parameter search direction, and the implementation can blend the information calculation with an EM-type diagonal stabilization during early iterations.

This formulation is referred to as marginal because its principal covariance calculations are carried out in observation space through $V$ and $P$. This describes the numerical formulation more accurately than naming it after a particular matrix operation.

Sparse Henderson Average Information formulation in mmes

The mmes engine uses Henderson's mixed model equations. With $W=[X\;Z]$, the coefficient system is

$C = \left[ \begin{array}{cc} X^\prime R^{-1}X & X^\prime R^{-1}Z \\ Z^\prime R^{-1}X & Z^\prime R^{-1}Z+G^{-1} \end{array} \right],$

and the right-hand side is

$r = \left[ \begin{array}{c} X^\prime R^{-1}y \\ Z^\prime R^{-1}y \end{array} \right].$

The BLUE and BLUP solutions satisfy

$Cb=r, \qquad b=(\widehat{\beta}^\prime,\widehat{u}^\prime)^\prime.$

The current C++ implementation stores $C$ sparsely and uses an Eigen SimplicialLDLT factorization. Symbolic factorization information is reused when the sparsity pattern remains unchanged. Random-effect precision contributions are inserted blockwise. Structured residual covariance models use optimized diagonal or repeated-block/Kronecker paths when applicable and a generic sparse residual-factorization path otherwise.

The REML criterion is evaluated from determinant contributions for the random covariance structures, residual covariance and mixed-model coefficient matrix, together with the quadratic term $y^\prime P y$. This is a mixed-model-equation representation of the same REML objective optimized by the marginal formulation.

For covariance parameter $\eta_i$, the generic covariance engine supplies

$\Sigma_i = \frac{\partial\Sigma}{\partial\eta_i}.$

For

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

a parameter belonging to factor $K_j$ has derivative

$\frac{\partial\Sigma}{\partial\eta_{jk}} = \sigma^2 K_1\otimes\cdots\otimes \frac{\partial K_j}{\partial\eta_{jk}} \otimes\cdots\otimes K_m.$

For the working scale $\tau=\log(\sigma^2)$,

$\frac{\partial\Sigma}{\partial\tau}=\Sigma.$

The score trace terms use selected elements associated with the sparse coefficient system. These are obtained from the LDLT factors with Takahashi sparse-inverse recursions; a complete $C^{-1}$ is not required during REML iterations.

Average Information is calculated by analytically differentiating the factorized mixed-model equations. Differentiating

$Cb=r$

with respect to $\eta_i$ gives

$C b_i = r_i-C_i b,$

where

$C_i=\frac{\partial C}{\partial\eta_i}, \qquad r_i=\frac{\partial r}{\partial\eta_i}.$

The sensitivity systems therefore reuse the same sparse factorization of $C$. For covariance $\Sigma$ with precision $\Lambda=\Sigma^{-1}$ and derivative $B_i=\partial\Sigma/\partial\eta_i$,

$\Lambda_i=-\Lambda B_i\Lambda.$

This factor-differentiation formulation avoids constructing an observation-space working variate for every random covariance parameter. Residual/residual and random/residual Average Information terms use the corresponding mixed-model/Schur-complement identities.

Covariance constructors may provide analytic derivatives in C++ or R. If a constructor uses numerical differentiation, only its small covariance factor $K_j(\eta)$ is differentiated by central finite differences; the REML likelihood and downstream score and Average Information calculations are not numerically differentiated.

Sparse inverse subsets and prediction error variances

During optimization, mmes does not routinely materialize the complete $C^{-1}$. Takahashi recursions provide the selected inverse subset required for REML trace calculations. After convergence, computeCi=1 obtains the diagonal information needed for prediction error variances through the sparse inverse subset, whereas computeCi=2 requests the complete coefficient-matrix inverse.

Parameter updates and numerical safeguards

Covariance parameters are optimized in working coordinates chosen to respect their parameter spaces, such as logarithms for positive scales and hyperbolic-tangent or bounded-logit mappings for correlations. Each CovarianceFactor descriptor supplies its natural-scale reporting transformation and parameter-specific trust limit.

The information system determines a proposed REML update. In the current mmes implementation, the proposal is limited in working-coordinate space and then checked against the exact REML likelihood. If the proposed global step decreases the accepted likelihood beyond numerical tolerance, the step is geometrically reduced and reevaluated. A likelihood decrease is therefore not interpreted as convergence. Boundary handling and robust information solves are retained for difficult or weakly identified variance components.

Choosing between the two formulations

Both engines estimate the same class of linear mixed models but organize the numerical work differently. mmer uses the marginal covariance representation and is attractive when observation-space calculations are moderate in size. mmes uses the sparse Henderson coefficient system and is particularly attractive when the mixed-model equations are sparse and their factorization can be reused efficiently. The best choice depends on the dimensions, sparsity and covariance structures of the fitted model rather than on a universal rule based only on the number of records or coefficients.

See mmes, mmer, and vsm for the interfaces and available covariance structures.

Bug report and contact

If you have any questions or suggestions please post it in https://stackoverflow.com or https://stats.stackexchange.com

I'll be glad to help or answer any question. I have spent a valuable amount of time developing this package. Please cite this package in your publication. Type 'citation("sommer")' to know how to cite it.

Author(s)

Giovanny Covarrubias-Pazaran

References

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

Covarrubias-Pazaran G. 2018. Software update: Moving the R package sommer to multivariate mixed models for genome-assisted prediction. doi: https://doi.org/10.1101/354639

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

Bernardo Rex. 2010. Breeding for quantitative traits in plants. Second edition. Stemma Press. 390 pp.

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

Henderson C.R. 1975. Best Linear Unbiased Estimation and Prediction under a Selection Model. Biometrics vol. 31(2):423-447.

Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.

Lee et al. 2015. MTG2: An efficient algorithm for multivariate linear mixed model analysis based on genomic information. Cold Spring Harbor. doi: http://dx.doi.org/10.1101/027201.

Maier et al. 2015. Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder. Am J Hum Genet; 96(2):283-294.

Searle. 1993. Applying the EM algorithm to calculating ML and REML estimates of variance components. Paper invited for the 1993 American Statistical Association Meeting, San Francisco.

Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208.

Tunnicliffe W. 1989. On the use of marginal likelihood in time series model estimation. JRSS 51(1):15-27.

Examples

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

####=========================================####
#### EXAMPLES
#### Different models with sommer
####=========================================####

data(DT_example, package="enhancer")





DT <- DT_example
DT=DT[with(DT, order(Env)), ]
head(DT)

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

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

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

Description

Calculates the realized additive relationship matrix. Currently is the C++ implementation of van Raden (2008).

Usage

A.mat(X,min.MAF=0,return.imputed=FALSE)

Arguments

Details

For vanraden method: the marker matrix is centered by subtracting column means $M= X - ms$ where ms is the coumn means. Then $A=M M'/c$, where $c = \sum_k{d_k}/k$, the mean value of the diagonal values of the $M M'$ portion.

Value

If return.imputed = FALSE, the $n \times n$ additive relationship matrix is returned.

If return.imputed = TRUE, the function returns a list containing

$A

the A matrix

$X

the imputed marker matrix

References

Endelman, J.B., and J.-L. Jannink. 2012. Shrinkage estimation of the realized relationship matrix. G3:Genes, Genomes, Genetics. 2:1405-1413. doi: 10.1534/g3.112.004259

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

Examples

####=========================================####
#### random population of 200 lines with 1000 markers
####=========================================####
X <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  X[i,] <- ifelse(runif(1000)<0.5,-1,1)
}

A <- A.mat(X)

####=========================================####
#### take a look at the Genomic relationship matrix 
#### (just a small part)
####=========================================####
# colfunc <- colorRampPalette(c("steelblue4","springgreen","yellow"))
# hv <- heatmap(A[1:15,1:15], col = colfunc(100),Colv = "Rowv")
# str(hv)

Description

Identifies problem configurations (factors, levels, dimensions) where experimental solver paths provide speedups over baseline methods, fitting empirical crossover conditions.

Usage

analyze_mme_crossovers(
  benchmark,
  baseline_engine = NULL,
  candidate_engine = NULL
)

Arguments

Value

An object of class "sommer_mme_crossover".

Description

anova method for class "mmes".

Usage

## S3 method for class 'mmes'
anova(object, object2=NULL, ...)

Arguments

Value

vector of anova

Author(s)

Giovanny Covarrubias

See Also

anova, mmes, wald_mmes

Description

Antedependence covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

antem(x, order = 1L, beta = NULL, innovations = NULL, fixed = NULL)

Arguments

Details

The modified-Cholesky representation is

$Ty=e,\qquad \mathrm{Cov}(e)=D,\qquad K=T^{-1}DT^{-\mathsf T}.$

The matrix $T$ is unit lower triangular, with nonzero regression coefficients only within the requested number of preceding levels. The first innovation variance is the scale reference; the remaining innovation variances are positive ratios.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

ar1m, toeplitzm, vsm.

Examples

## Not run: 
vsm(antem(time, order=2), ism(id))

## End(Not run)

Description

Constructs the Algorithm for Proven and Young (APY) approximation for a positive-definite relationship matrix and a user-specified core set. In the genomic evaluation literature, the input is the genomic relationship matrix $G$ for genotyped animals, usually calculated from SNP marker data. The core and non-core partition is within that genotyped set; it is not a pedigree-versus-genomic partition. The non-core conditional covariance is approximated by its diagonal. The returned matrix is a sparse precision matrix suitable for the Gu argument of mmes() in a genomic BLUP model.

Usage

APY(G, core, tol = 1e-10, return.details = FALSE)

Arguments

Details

Partition $G$ into core and non-core blocks. With $B = G_{nc}G_{cc}^{-1}$ and $D = \operatorname{diag}(G_{nn} - BG_{cn})$, APY uses the precision $\left[\begin{smallmatrix} G_{cc}^{-1} + B' D^{-1} B & -B' D^{-1} \\ -D^{-1} B & D^{-1} \end{smallmatrix}\right]$. It preserves the core block and core-to-non-core relationships while replacing non-core conditional correlations with zero. Thus APY is exact when the non-core conditional covariance is diagonal; otherwise it is an approximation to the supplied relationship model.

APY was developed to provide a computationally economical approximation to the inverse genomic relationship matrix for large genomic evaluations (Misztal, 2016; Fragomeni et al., 2015). The algorithm partitions genotyped animals into core and non-core individuals. The core contains the relationships that are modeled jointly; each non-core animal is regressed on the core, with its remaining conditional variance treated as independent of other non-core animals. Thus non-core animals are not independent of the core: their core-to-non-core relationships are retained.

Pedigree relationships are a distinct input. For pedigree-only models, the pedigree relationship matrix $A$ and its sparse inverse $A^{-1}$ are normally used directly; APY is not required. In single-step genomic BLUP, pedigree and genomic information are combined, conventionally through $H^{-1} = A^{-1} + \left[\begin{smallmatrix}0&0\\0&G^{-1}-A_{22}^{-1}\end{smallmatrix}\right]$, where genotyped animals are ordered last and $A_{22}$ is their pedigree relationship block. APY approximates the genomic term $G^{-1}$ in that correction; it does not replace the pedigree term or construct $H^{-1}$. This function returns only the APY genomic precision. Do not use it as the complete single-step precision for a population that also contains non-genotyped animals.

This function constructs APY from an already available relationship matrix. It does not avoid the cost of constructing a full dense relationship matrix from markers, and it does not change the solver algorithm. Marker-based block construction and single-step matrix assembly are not implemented here.

Value

A sparse symmetric APY precision matrix with individual dimnames and attr(x, "inverse") = TRUE, ready for Gu. Metadata describing the core and conditional variances is stored in attr(x, "APY"). If return.details = TRUE, returns the diagnostic list described above.

References

Misztal I (2016). Inexpensive Computation of the Inverse of the Genomic Relationship Matrix in Populations with Small Effective Population Size. Genetics, 202(2), 401–409. doi:10.1534/genetics.115.182089.

Fragomeni BO, Lourenco DA, Tsuruta S, Masuda Y, Aguilar I, Legarra A, Lawlor TJ, Misztal I (2015). Use of genomic recursions and algorithm for proven and young animals for single-step genomic BLUP analyses–a simulation study. Journal of Animal Breeding and Genetics, 132(5), 340–345. doi:10.1111/jbg.12161.

See Also

mmes, A.mat, H.mat

Examples

## Marker-derived GBLUP for genotyped animals.
set.seed(2)
ids <- paste0("animal", seq_len(12))
markers <- matrix(sample(c(-1, 0, 1), 12 * 80, replace = TRUE), nrow = 12,
                  dimnames = list(ids, NULL))
G <- A.mat(markers)

## Blend slightly with the identity so this small illustrative G is positive
## definite, then construct APY precision for a six-animal core.
G <- 0.95 * G + 0.05 * diag(nrow(G))
dimnames(G) <- list(ids, ids)
Ginv.apy <- APY(G, core = ids[seq_len(6)])

phenotypes <- data.frame(
  id = factor(rep(ids, each = 3), levels = ids),
  y = rnorm(length(ids) * 3)
)
fit <- mmes(y ~ 1, random = ~vsm(ism(id), Gu = Ginv.apy),
            rcov = ~units, data = phenotypes, verbose = FALSE)

## In a single-step evaluation, this APY approximation supplies G^{-1} for
## the genotyped block correction; it is not itself the complete H^{-1}.

Description

First-order autoregressive covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

ar1m(x, rho = 0.30, fixed = FALSE,
     variance = c("homogeneous", "heterogeneous"), values = NULL)

Arguments

Details

The homogeneous AR(1) shape is

$K_{ij}=\rho^{|i-j|}.$

Heterogeneous mode uses

$K=D R(\rho)D,$

with the first variance fixed as the relative-scale reference. The working correlation coordinate is $\eta=\operatorname{atanh}(\rho)$, ensuring $|\rho|<1$.

Factor levels define the lag order, so numeric-looking levels must be in increasing numeric order; ar1m() warns otherwise. Inside vsm() one correlation is shared by every level of the other Kronecker factors. For a separate correlation per trial, wrap the vsm() term in dsumm.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

ar2m, ar3m, toeplitzm, vsm, dsumm.

Examples

## Not run: 
vsm(ar1m(time), ism(id))

# Separable row-by-range residual shared across trials:
# rcov = ~ vsm(dsm(trial), ar1m(range), ar1m(row), ism(units))
# Trial-specific correlations and variances:
# rcov = ~ dsumm(vsm(ar1m(range), ar1m(row), ism(units)), by=trial)

## End(Not run)

Description

Second-order autoregressive covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

ar2m(x, pacf = c(0.20, 0.10), fixed = NULL,
     variance = c("homogeneous", "heterogeneous"), values = NULL)

Arguments

Details

A stationary AR(2) covariance is parameterized through reflection coefficients / partial autocorrelations. Each PACF is mapped from an unconstrained working coordinate by tanh, then converted to stable AR coefficients by the Levinson/Durbin recursion. Heterogeneous mode uses

$K=D R D$

and estimates q-1 positive variance ratios in addition to the PACFs.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

ar1m, ar3m, toeplitzm, vsm.

Examples

## Not run: 
vsm(ar2m(time), ism(id))

## End(Not run)

Description

Third-order autoregressive covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

ar3m(x, pacf = c(0.20, 0.10, 0.05), fixed = NULL,
     variance = c("homogeneous", "heterogeneous"), values = NULL)

Arguments

Details

A stationary AR(3) covariance is parameterized through three partial autocorrelations. The PACF parameterization guarantees stationarity while allowing unconstrained working coordinates. Heterogeneous mode uses

$K=D R D$

and estimates q-1 positive variance ratios in addition to the PACFs.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

ar1m, ar2m, toeplitzm, vsm.

Examples

## Not run: 
vsm(ar3m(time), ism(id))

## End(Not run)

Description

Selected-level diagonal covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

atm(x, levs, values = NULL, fixed = NULL)

Arguments

Details

atm is a selected-level version of dsm. The design is restricted to levs; observations outside those levels receive zeros for this factor. The first selected level is the scale reference and the remaining selected levels are positive variance ratios. This representation is most naturally used for random effects; residual factors must still assign exactly one covariance-product coordinate to every observation.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

dsm, vsm, mmes.

Examples

## Not run: 
vsm(atm(environment, c("E1","E3")), ism(genotype))

## End(Not run)

Description

Creates a diagonal covariance-component selector for specified columns or levels of a covariance dimension. It is intended for the mmer covariance-model interface, typically inside vsr.

Usage

atr(x, levs)

Arguments

Details

atr() constructs a design matrix $Z$ together with a diagonal thetaC selector. Let $q$ be the number of columns of $Z$, and define

$a_j = \begin{cases} 1, & \text{if level }j\text{ is included in \code{levs}},\\ 0, & \text{otherwise}. \end{cases}$

The returned constraint matrix is

$\mathrm{thetaC}=\mathrm{diag}(a_1,\ldots,a_q).$

Thus atr() is useful when a covariance component is to be associated only with a specified subset of levels or design columns. Levels not selected in levs receive zero on the corresponding diagonal of thetaC.

This function belongs to the mmer/vsr covariance-structure interface. The related atm function belongs to the newer mmes/vsm CovarianceFactor interface and uses a different parameterization.

For a matrix input, its column names define the available levels. For factor or character input, the incidence-matrix column names define the levels. If levs is missing, all available columns are selected.

Value

A list with components:

• Z: the design/incidence representation of x.

• thetaC: a diagonal zero/one matrix selecting the levels specified by levs.

See Also

vsr, dsr, usr, csr, atm, mmer

Examples

# Select levels B and C
A <- atr(factor(c("A","B","A","C")), levs = c("B","C"))
A$Z
A$thetaC

Description

Builds a single proportion response for binomial mmes fits from counts of successes and failures (or trials). The number of trials travels with the proportion as an attribute and is used as the binomial prior weight, so the single-response interface of mmes is preserved.

Usage

binm(successes, failures = NULL, trials = NULL)

Arguments

Details

The result is successes/trials with class "binm" and attribute "trials". Subsetting keeps the trials aligned, so the vector can be stored as a data-frame column. Missing counts give a missing response, handled by the naMethodY rules of mmes. Fitting a proportion without binm treats every value as one trial, which is wrong unless all trials equal one.

Value

A numeric vector of proportions of class "binm".

See Also

mmes

Examples

set.seed(1)
d <- data.frame(x = rnorm(60), g = factor(rep(1:12, each = 5)))
d$n <- sample(5:10, 60, TRUE)
d$s <- rbinom(60, d$n, plogis(0.3 * d$x))
fit <- mmes(binm(s, trials = n) ~ x, random = ~ g, data = d,
            family = binomial(), blendWeight = 0, verbose = FALSE)
fit$b

Description

Proper conditional autoregressive spatial covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

car(x, W, rho = 0.10, fixed = FALSE)

Arguments

Details

Let $D=\mathrm{diag}(W\mathbf{1})$. The proper CAR precision and covariance are

$Q=D-\rho W,\qquad M=Q^{-1},\qquad K=M/M_{11}.$

The constructor requires positive row sums, so isolated levels are not allowed. The admissible open interval for rho is obtained from the eigenvalues of $D^{-1/2}WD^{-1/2}$ and is enforced through a bounded-logit working coordinate. This is a proper, nonsingular CAR model; an intrinsic singular CAR is not used by the current precision-based Henderson implementation. The covariance derivative is analytic.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

sar, maternm, vsm.

Examples

## Not run: 
vsm(car(location, W), ism(genotype))

## End(Not run)

Description

coef method for class "mmes".

Usage

## S3 method for class 'mmes'
coef(object, ...)

Arguments

Value

vector of coef

Author(s)

Giovanny Covarrubias

See Also

coef, mmes

Description

General positive-definite correlation structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

corgm(x, theta = NULL, fixed = NULL,
      variance = c("homogeneous", "heterogeneous"), values = NULL)

Arguments

Details

An unrestricted SPD correlation matrix is represented using a unit-diagonal lower factor $A$. The intermediate matrix $S=AA^{\mathsf T}$ is standardized to correlation scale,

$K_{ij}=S_{ij}/\sqrt{S_{ii}S_{jj}}.$

The unrestricted working parameters are the strict lower-triangular entries of $A$.

With variance="heterogeneous" the shape is $K = DRD$ with $D=\mathrm{diag}(1, \sqrt{r_2}, \ldots, \sqrt{r_q})$, where $r_l=\sigma^2_l/\sigma^2_1 = \exp(\phi_l)$ are variance ratios relative to the first level; the vsm scale is $\sigma^2_1$. The model is the same as usm but parameterized as correlations and variances, which allows individual variances to be fixed. This is what mixed-family multi-trait models (familym) need: the residual variance of a binomial or Poisson trait is fixed to one while the correlations and the variances of Gaussian traits are estimated. Reported parameters (covparams_mmes) are the correlations and the variances of each level.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

csm, usm, vsm, covmatrix_mmes, familym.

Examples

data(DT_example, package = "enhancer")
DT <- DT_example
# genetic correlations between environments with a homogeneous variance
m1 <- mmes(Yield ~ Env, random = ~ vsm(corgm(Env), ism(Name)),
           rcov = ~ units, data = DT, verbose = FALSE)
covparams_mmes(m1, 1)

# heterogeneous variances (equivalent to usm, parameterized as corr/var)
m2 <- mmes(Yield ~ Env,
           random = ~ vsm(corgm(Env, variance = "heterogeneous"), ism(Name)),
           rcov = ~ units, data = DT, verbose = FALSE)
covparams_mmes(m2, 1)
covmatrix_mmes(m2, 1)$correlation

Description

corImputation imputes missing data based on the correlation that exists between row levels.

Usage

corImputation(wide, Gu=NULL, nearest=10, roundR=FALSE)

Arguments

Value

$res

a list with the imputed matrix and the original matrix.

Author(s)

Giovanny Covarrubias-Pazaran

References

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

Examples

####################################
### imputing genotype data example
####################################
# data(DT_cpdata, package="enhancer")
# X <- GT_cpdata
# # add missing data
# v <- sample(1:length(X), 500)
# Xna <- X
# Xna[v]<- NA
# ## impute (can take some time)
# Y <- corImputation(wide=Xna, Gu=NULL, nearest=20, roundR=TRUE) 
# cm <- table(Y$imputed[v],X[v])
# ## calculate accuracy
# sum(diag(cm))/length(v)
####################################
### imputing phenotypic data example
####################################
# data(DT_h2, package="enhancer")
# X <- reshape(DT_h2[,c("Name","Env","y")], direction = "wide", idvar = "Name",
#                 timevar = "Env", v.names = "y", sep= "_")
# rownames(X) <- X$Name
# X <- as.matrix(X[,-1])
# head(X)
# # add missing data
# v <- sample(1:length(X), 50)
# Xna <- X
# Xna[v]<- NA
# ## impute
# Y <- corImputation(wide=Xna, Gu=NULL, nearest=20, roundR=TRUE)
# plot(y=Y$imputed[v],x=X[v], xlab="true",ylab="predicted")
# cor(Y$imputed[v],X[v], use = "complete.obs")

Description

covm combines two random-effect structures created with vsm into a single correlated random-effect structure. It fits an unstructured $2 \times 2$ covariance matrix between the two random effects while allowing them to have different incidence matrices. The two effects must act on the same coefficient space and use the same relationship or precision matrix. The resulting structure can be fitted with the mmes solver.

Usage

covm(ran1, ran2, thetaC = NULL, theta = NULL,
fixed = NULL, fixedSigma2 = FALSE,
labels = c("ran1", "ran2"), tol = 1e-10)

Arguments

```

diag(2) * 0.15 + matrix(0.015, 2, 2)

Details

covm is intended for models in which two distinct random effects are expected to be correlated. A common example is the joint modeling of direct and indirect genetic effects.

If the coefficient vector for the first random effect is $u_1$ and that for the second is $u_2$, covm represents their covariance structure as

$\mathrm{Var}([u_1^\prime, u_2^\prime]^\prime) = \sigma^2 K_{\mathrm{effect}} \otimes A,$

where $A$ represents the covariance structure associated with the coefficient levels, and $K_{\mathrm{effect}}$ is a $2 \times 2$ unstructured covariance shape describing the relationship between the two random effects.

Internally, covm uses the same CovarianceFactor interface as vsm. The $2 \times 2$ effect covariance is represented using a normalized-Cholesky parameterization. This guarantees a positive-definite covariance matrix during optimization.

The variance of the first random effect is represented by the overall scale parameter sigma2. The remaining variance ratio and covariance are represented through the normalized-Cholesky factor.

The two random effects may have different incidence matrices. However, they must refer to identical coefficient levels in the same order and must use the same relationship/precision matrix.

When a relationship matrix is supplied through Gu, the Henderson solver requires it to be an inverse/precision matrix with attr(Gu, "inverse") = TRUE. For example:

attr(Ai, "inverse") <- TRUE

covm(
vsm(ism(focal), Gu = Ai),
vsm(ism(neighbour), Gu = Ai)
)

The current implementation combines simple vsm random effects only. Each input must contain one covariance-product coordinate. More complex covariance-factor products should therefore not currently be supplied independently inside ran1 and ran2.

Value

A list describing a correlated random-effect structure compatible with the CovarianceFactor-v2 interface used by vsm and mmes. The main elements are:

“'

Z

A list containing the two incidence matrices, one for each random effect.

Gu

The common sparse relationship precision matrix associated with the coefficient levels. It is stored with attr(Gu, "inverse") = TRUE for use by the Henderson solver.

covStruct

A CovarianceFactor-v2 descriptor containing the overall variance scale and the normalized-Cholesky representation of the $2 \times 2$ unstructured covariance between the two random effects.

residualLocalIndex

NULL, since covm describes a random-effect structure rather than a residual covariance structure.

productDesign

NULL for the current implementation.

partitionsR

NULL for the current implementation.

covm

Logical value equal to TRUE, identifying the object as one constructed by covm.

covm_labels

The two labels used for the correlated random effects.

“'

Author(s)

Giovanny Covarrubias-Pazaran

References

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

Bijma, P. (2014). The quantitative genetics of indirect genetic effects: a selective review of modelling issues. Heredity, 112(1), 61–69.

See Also

vsm for constructing random-effect covariance structures and mmes for fitting mixed models using the Henderson solver.

Examples

data(DT_ige, package = "enhancer")
DT <- DT_ige

## Correlate two random effects with identity relationship matrices

covRes <- with(
DT,
covm(
vsm(ism(focal)),
vsm(ism(neighbour))
)
)

str(covRes)

## Custom initial variance-covariance matrix

covRes2 <- with(
DT,
covm(
vsm(ism(focal)),
vsm(ism(neighbour)),
theta = matrix(
c(0.5, 0.1,
0.1, 0.3),
2, 2
)
)
)

## Relationship precision matrix

## Ai must have matching row and column names and be marked as inverse.

##

## attr(Ai, "inverse") <- TRUE

##

## covRes3 <- with(

## DT,

## covm(

## vsm(ism(focal), Gu = Ai),

## vsm(ism(neighbour), Gu = Ai)

## )

## )

## See the DT_ige help page for a complete model-fitting example.

Description

Returns the full covariance matrix of a random or residual term of a fitted mmes model over its covariance coordinates (for example the genetic trait-by-trait covariance matrix of a multi-trait model), together with the implied correlation matrix and delta-method standard errors of every element.

Usage

covmatrix_mmes(object, term, se = TRUE, max.dim = 500L, rel_step = 1e-5)

Arguments

Details

The matrix is $\sigma^2 K_1 \otimes \cdots \otimes K_m$ built from all the covariance factors of the vsm term (the relationship matrix Gu is not included), so for vsm(usm(trait), ism(id)) it is the trait covariance matrix, and for vsm(dsm(env), usm(trait), ism(id)) the environment-by-trait matrix. Row and column names are the coordinate levels.

Standard errors use the delta method, $\mathrm{Var}(g(\hat\theta)) \approx J V J^\top$, where $V$ is object$theta_se (the covariance of the reported parameters, as used by covparams_mmes_se) and $J$ the Jacobian of each element with respect to the reported parameters. Fixed parameters contribute no variance. Diagonal standard errors of the correlation matrix are zero. Terms split by dsumm sections are not supported.

Value

A list with covariance and correlation and, if se=TRUE, covariance.se and correlation.se.

See Also

covparams_mmes, covparams_mmes_se, vpredict, stackTraits

Examples

set.seed(17)
n <- 60
DT <- data.frame(id = factor(rep(seq_len(n), each = 3)))
u <- matrix(rnorm(n * 2), n) 
DT$trait_a <- u[DT$id, 1] + rnorm(nrow(DT), sd = 0.7)
DT$trait_b <- u[DT$id, 2] + rnorm(nrow(DT), sd = 1)
DTL <- stackTraits(DT, traits = c("trait_a", "trait_b"), keep = "id")
fit <- mmes(value ~ trait,
            random = ~ vsm(usm(trait), ism(id)),
            rcov = ~ vsm(usm(trait), ism(record)),
            data = DTL, verbose = FALSE)
G <- covmatrix_mmes(fit, 1)       # genetic
G$covariance
G$correlation
G$correlation.se
R <- covmatrix_mmes(fit, 2)       # residual
R$correlation

Description

Returns covariance parameters in the native, model-facing scale defined by each CovarianceFactor descriptor, rather than the unconstrained working coordinates or normalized variance ratios used internally by the optimizer.

Usage

covparams_mmes(object, term = NULL)

Arguments

Details

Every CovarianceFactor descriptor carries a native reporting callback. Built-in structures use it to report meaningful quantities such as level-specific variances for dsm(), variances and covariances for usm(), and variance-correlation parameters for correlation structures. Structural zeros and internal coordinates such as log standard deviations, normalized Cholesky coefficients, and variance ratios are not returned where the covariance model provides a more direct native representation.

All built-in CovarianceFactor families provide native reporters, including identity, diagonal, unstructured, AR(1–3), compound-symmetry, moving-average, general-correlation, factor-analytic, antedependence, reduced-rank, Matern, Toeplitz, SAR, and CAR structures. Heterogeneous AR and compound-symmetry models report their dependence parameter together with level-specific variances rather than an overall variance plus variance ratios.

Both fixed and estimated covariance parameters are included. For an arbitrary Kronecker product, the product-level variance scale is absorbed by the first covariance-shaping factor. Subsequent factors are necessarily reported on their normalized relative scale because the product has only one identifiable overall scale.

User-defined ownm() structures can supply a native_report callback. If none is supplied, their overall variance and naturally transformed parameter values are returned.

Use covparams_mmes_se() when native-scale delta-method standard errors are also required.

Value

A data frame with columns term, factor, section, parameter, and estimate. section identifies the dsumm section a parameter belongs to and is NA otherwise.

See Also

covparams_mmes_se, mmes, vsm, ownm

Examples

## Not run: 
fit <- mmes(Yield ~ Env,
            random = ~ vsm(dsm(Env), ism(Name)),
            rcov = ~ units, data = DT_example)
covparams_mmes(fit)
fit$covParNative

## End(Not run)

Description

Returns descriptor-defined native covariance parameters with delta-method standard errors computed from the fitted covariance-parameter uncertainty matrix.

Usage

covparams_mmes_se(object, term = NULL, rel_step = 1e-6)

Arguments

Details

The fitted theta_se matrix describes uncertainty in the reported covPar coordinates. For each covariance structure, covparams_mmes_se numerically evaluates the Jacobian of its descriptor's native reporting callback and applies the multivariate delta method,

$V_{native} = J V_{covPar} J^{\mathsf T}.$

Consequently, coupled transformations use the complete covariance information. For example, the standard error of a dsm() level variance accounts for uncertainty and covariance in both the overall variance and its variance ratio. Likewise, usm() standard errors account for all relevant normalized Cholesky parameters.

The transformation is available for every built-in CovarianceFactor family supported by covparams_mmes(); user-defined ownm() structures use their descriptor-provided native_report callback.

Parameters marked fixed in the covariance descriptor contribute no independent sampling uncertainty: their rows and columns in theta_se are excluded before transformation. A native quantity involving both fixed and free parameters can still have a nonzero standard error through its free parameters. Both fixed and estimated native quantities are returned.

The same table is stored automatically in object$covParNativeSE by mmes().

Value

A data frame containing term, factor, parameter, estimate, StdError, and Zratio. A zero standard error is reported for a fully fixed native quantity, with Zratio=NA.

See Also

covparams_mmes, mmes, vsm

Examples

## Not run: 
fit <- mmes(Yield ~ Env,
            random = ~ vsm(dsm(Env), ism(Name)),
            rcov = ~ units, data = DT_example)
covparams_mmes_se(fit)
fit$covParNativeSE

## End(Not run)

Description

Compound-symmetry covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

csm(x, rho = 0.10, fixed = FALSE,
    variance = c("homogeneous", "heterogeneous"), values = NULL)

Arguments

Details

With variance = "homogeneous", csm uses unit diagonal and a common off-diagonal correlation,

$K_{ii}=1, \quad K_{ij}=\rho\ (i\ne j).$

With variance = "heterogeneous", it uses

$K=D C(\rho)D,$

where the first variance is the reference and the remaining q-1 variances are estimated as positive ratios. Positive definiteness requires $-1/(q-1)<\rho<1$. The product-level variance remains owned by vsm.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

corgm, dsm, vsm.

Examples

## Not run: 
vsm(csm(environment), ism(genotype))

## End(Not run)

Description

Combines a design/incidence representation with a user-supplied covariance-component constraint matrix for use by the mmer covariance-model interface, typically inside vsr.

Usage

csr(x, mm)

Arguments

Details

csr() provides the general constraint-matrix interface underlying custom covariance patterns for mmer/vsr. Unlike dsr and usr, it does not construct a predefined diagonal or unstructured pattern. Instead, the supplied matrix mm determines the covariance-component structure through thetaC.

Conceptually, thetaC identifies which elements of the covariance matrix are associated with estimated variance-covariance components and which constraints are shared among elements. The interpretation of the entries therefore follows the thetaC convention used by vsr.

When x is a matrix it is used directly as $Z$. Otherwise, factors and character vectors are converted to incidence matrices, while numeric non-factor inputs are treated as a single design column.

The function assigns the column names of $Z$ to both dimensions of mm; therefore mm should be conformable with the number of columns of the resulting design matrix.

Value

A list with components:

• Z: the design/incidence matrix associated with the covariance dimension.

• thetaC: the user-supplied covariance-component constraint matrix, with dimensions named according to Z.

See Also

vsr, dsr, usr, atr, mmer

Examples

# A custom two-level covariance-component pattern
x <- factor(c("A","B","A","B"))
mm <- matrix(c(1, 2,
               2, 1), 2, 2, byrow = TRUE)
C <- csr(x, mm)
C$Z
C$thetaC

Description

C++ implementation of the dominance matrix. Calculates the realized dominance relationship matrix. Can help to increase the prediction accuracy when 2 conditions are met; 1) The trait has intermediate to high heritability, 2) The population contains a big number of individuals that are half or full sibs (HS & FS).

Usage

D.mat(X,nishio=TRUE,min.MAF=0,return.imputed=FALSE)

Arguments

Details

Missing values are replaced by the column (marker) mean before the matrix is computed. Markers with minor allele frequency less than or equal to min.MAF (and therefore all monomorphic markers) are removed. For each marker $k$, $p_k$ is the frequency of the allele counted by the code 1, estimated as the column mean of $(X+1)/2$, and $q_k = 1-p_k$.

For the nishio method (Nishio and Satoh, 2014; Vitezica et al., 2013): genotypes are coded as $w_{ik} = -2q_k^2$ for $x_{ik}=1$, $2p_kq_k$ for $x_{ik}=0$ and $-2p_k^2$ for $x_{ik}=-1$, computed as $w_{ik} = -x_{ik}^2 + (p_k-q_k)x_{ik} + 2p_kq_k$. Fractional values are passed through the same formula; note that a mean-imputed value $x_{ik} = p_k-q_k$ receives the heterozygote code $2p_kq_k$. Then $D = W W'/c$, where $c = \sum_k (2p_kq_k)^2$. Under Hardy-Weinberg equilibrium this coding is orthogonal to the additive coding.

For the Su method (Su et al., 2012): heterozygosity is coded as $h_{ik} = 1-|x_{ik}|$ and centered as $m_{ik} = h_{ik} - 2p_kq_k$. Then $D = M M'/c$, where $c = \sum_k 2p_kq_k(1-2p_kq_k)$.

Versions of sommer earlier than 4.4.9 used a different formula for nishio=TRUE: $h_{ik}$ centered by its column mean and divided by $\sum_k (2p_kq_k)^2$. That was not the Nishio and Satoh (2014) method, so results from those versions will differ.

Value

If return.imputed = FALSE, the $n \times n$ dominance relationship matrix is returned.

If return.imputed = TRUE, the function returns a list containing

$D

the D matrix

$X

the imputed marker matrix

References

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

Nishio M and Satoh M. 2014. Including Dominance Effects in the Genomic BLUP Method for Genomic Evaluation. Plos One 9(1), doi:10.1371/journal.pone.0085792

Vitezica ZG, Varona L, Legarra A. 2013. On the Additive and Dominant Variance and Covariance of Individuals Within the Genomic Selection Scope. Genetics 195(4): 1223-1230. doi:10.1534/genetics.113.155176

Su G, Christensen OF, Ostersen T, Henryon M, Lund MS. 2012. Estimating Additive and Non-Additive Genetic Variances and Predicting Genetic Merits Using Genome-Wide Dense Single Nucleotide Polymorphism Markers. PLoS ONE 7(9): e45293. doi:10.1371/journal.pone.0045293

Examples

####=========================================####
#### EXAMPLE 1
####=========================================####
####random population of 200 lines with 1000 markers
X <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  X[i,] <- sample(c(-1,0,0,1), size=1000, replace=TRUE)
}

D <- D.mat(X)

Description

This function takes a matrix that is in data frame format and transforms it into a matrix. Other packages that allows you to obtain an additive relationship matrix from a pedigree is the 'pedigreemm' package.

Usage

dfToMatrix(x, row="Row",column="Column",
             value="Ainverse", returnInverse=FALSE, 
             bend=1e-6)

Arguments

Value

Author(s)

Giovanny Covarrubias-Pazaran

References

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

Examples

library(Matrix)
m <- matrix(1:9,3,3)
m <- tcrossprod(m)

mdf <- as.data.frame(as.table(m))
mdf

dfToMatrix(mdf, row = "Var1", column = "Var2", 
            value = "Freq",returnInverse=FALSE )

Description

Diagonal heterogeneous covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

dsm(x, values = NULL, fixed = NULL, theta = NULL)

Arguments

Details

For $q$ levels, dsm uses the first variance as the internal scale reference and estimates $q-1$ positive variance ratios,

$K=\mathrm{diag}(1,r_2,\ldots,r_q).$

The ratios are optimized on log scales and reported on their natural positive scale.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

vsm, mmes.

Examples

## Not run: 
vsm(dsm(environment), ism(genotype))

## End(Not run)

Description

Creates a diagonal variance-component structure for use by the mmer covariance-model interface, typically inside vsr. Each column or level represented by x receives its own diagonal variance component, while all off-diagonal covariance components are constrained to zero.

Usage

dsr(x)

Arguments

Details

dsr() is part of the covariance-structure interface used by mmer/vsr; it is distinct from dsm, which belongs to the newer mmes/vsm CovarianceFactor interface.

The function returns a design matrix $Z$ and a covariance-parameter constraint matrix thetaC. If $q$ columns or levels are represented by x, then

$\mathrm{thetaC}=I_q.$

Consequently, the associated covariance matrix has the form

$\Sigma=\mathrm{diag}(\sigma_1^2,\ldots,\sigma_q^2),$

so that the $q$ variances are estimated separately and all covariances are fixed to zero.

For a factor or character vector, missing observations are retained through the construction of the incidence matrix. If only one non-missing level is present, a one-column incidence matrix is constructed explicitly.

For a numeric non-factor input, x is treated as a single design column rather than expanded into factor levels.

Value

A list with components:

• Z: the design/incidence matrix associated with the covariance dimension.

• thetaC: a diagonal matrix defining the diagonal covariance-component pattern used by vsr.

See Also

vsr, usr, csr, atr, dsm, mmer

Examples

# Diagonal covariance structure across levels
D <- dsr(factor(c("A","B","A","C")))
D$Z
D$thetaC

Description

dsumm replicates a residual vsm structure independently for every level (section) of a grouping variable, such as a trial, location, or experiment. Each section receives its own variance and its own copy of every covariance parameter of the inner structure.

Usage

dsumm(x, by, levels = NULL)

Arguments

Details

Direct sum versus Kronecker product. A vsm() term is a Kronecker product with a single overall variance. Including a diagonal section factor,

$\mathrm{vsm}(\mathrm{dsm}(t), \mathrm{ar1m}(r), \mathrm{ar1m}(c), \mathrm{ism}(units)):\quad R = \sigma^2\, D_t \otimes K_r(\rho_r) \otimes K_c(\rho_c),$

gives each section its own variance, but a single correlation shared by all sections, because each factor of a Kronecker product appears only once. In contrast, dsumm builds the block-diagonal direct sum

$R = \bigoplus_{t} \sigma^2_t\, K_r(\rho_{r,t}) \otimes K_c(\rho_{c,t}),$

so every section has its own variance and its own correlations. Sections are mutually independent, so their parameters are estimated only from their own observations; with section-specific fixed effects the joint fit reproduces separate per-section fits.

Choosing between the two. Use the vsm() form when trials are small or similar and a common spatial correlation is adequate (fewer parameters). Use dsumm when correlations are expected to differ between trials. With $T$ sections and $p$ inner covariance parameters, dsumm estimates $T(p+1)$ residual parameters.

Unlisted levels. Levels omitted from levels are modelled as $\sigma^2_t I$. Their observations are retained even when the variables of the inner structure (for example row or range) are missing, so trials without spatial coordinates can be analysed together with spatially designed trials.

Computation. Each section is embedded in the rectangle spanned by its observed levels of the inner factors. Covariance applications use the small factor inverses, and missing cells (e.g. plots with a missing response) are absorbed exactly by a Schur complement, so the dense observation-level covariance matrix is never formed. This requires sections that are close to complete rectangles and no W weights.

Reporting. covparams_mmes returns one row per parameter and section, with the section in the section column. The variance rows belong to the by factor; correlation rows are labelled by their inner factor, e.g. ar1m(range).

Value

A residual covariance structure, in the same form as a residual vsm() term, for use in the rcov argument of mmes.

Limitations

dsumm is currently available only for residual (rcov) structures, without W weights, and cannot be nested.

See Also

vsm, ar1m, dsm, mmes, covparams_mmes.

Examples

## Not run: 

mix <- mmes(yield ~ trial,
            rcov = ~ dsumm(vsm(ar1m(range), ar1m(row), ism(units)), by = trial),
            data = dat)
covparams_mmes(mix)

# Spatial structure only for two trials; remaining trials are iid.
mix2 <- mmes(yield ~ trial,
             rcov = ~ dsumm(vsm(ar1m(range), ar1m(row), ism(units)),
                            by = trial, levels = c("T1", "T2")),
             data = dat)

# Site-specific unstructured covariance among environments.
mix3 <- mmes(y ~ site,
             rcov = ~ dsumm(vsm(usm(env), ism(units)), by = site),
             data = dat2)

## End(Not run)

Description

Calculates the realized epistatic relationship matrix of second order (additive x additive, additive x dominance, or dominance x dominance) using hadamard products with the C++ Armadillo library.

Usage

E.mat(X,nishio=TRUE,type="A#A",min.MAF=0.02)

Arguments

Details

it is computed as the Hadamard product of the epistatic relationship matrix; E=A#A, E=A#D, E=D#D.

Value

The epistatic relationship matrix is returned.

References

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

Endelman, J.B., and J.-L. Jannink. 2012. Shrinkage estimation of the realized relationship matrix. G3:Genes, Genomes, Genetics. 2:1405-1413. doi: 10.1534/g3.112.004259

Nishio M and Satoh M. 2014. Including Dominance Effects in the Genomic BLUP Method for Genomic Evaluation. Plos One 9(1), doi:10.1371/journal.pone.0085792

Su G, Christensen OF, Ostersen T, Henryon M, Lund MS. 2012. Estimating Additive and Non-Additive Genetic Variances and Predicting Genetic Merits Using Genome-Wide Dense Single Nucleotide Polymorphism Markers. PLoS ONE 7(9): e45293. doi:10.1371/journal.pone.0045293

Examples

####=========================================####
####random population of 200 lines with 1000 markers
####=========================================####
X <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  X[i,] <- sample(c(-1,0,0,1), size=1000, replace=TRUE)
}

E <- E.mat(X, type="A#A") 
# if heterozygote markers are present can be used "A#D" or "D#D"

Description

Factor-analytic covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

fam(x, k = 1L, loadings = NULL, specific = NULL, fixed = NULL,
    variance = c("heterogeneous", "homogeneous"))

Arguments

Details

The factor-analytic shape is based on

$M=\Lambda\Lambda^{\mathsf T}+\Psi,$

where $\Psi$ is diagonal positive, followed by $K=M/M_{11}$. For rotational identification, the leading $k\times k$ loading block is lower triangular and its diagonal loadings are positive. The first specific variance is an internal reference; remaining specific variances are positive ratios. This removes the otherwise redundant internal scale because vsm supplies the overall $\sigma^2$.

With variance="heterogeneous", fixed has length $qk-k(k-1)/2+q-1$: triangular loadings followed by q-1 specific-variance ratios. With variance="homogeneous", all specific ratios are one, $\Psi=I_q$ internally, and fixed has length $qk-k(k-1)/2$. The common specific variance on the fitted scale is estimated, not fixed at one, and is reported as common_specific_variance. This homogeneous model reproduces the covariance family previously fitted by rrm(). For a genuine reduced-rank model with all specific variances fixed at zero, use rrm.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

rrm, usm, vsm.

Examples

## Not run: 
vsm(fam(environment, k=2), ism(genotype))
vsm(fam(environment, k=2, variance="homogeneous"), ism(genotype))

## End(Not run)

Description

Assigns a stats::family() object to each trait of a long-format multi-trait mmes model, e.g. a binary disease score together with a Gaussian yield, so that their genetic and residual correlations can be estimated jointly by penalized quasi-likelihood (PQL).

Usage

familym(by, ...)

Arguments

Details

The object is passed as the family argument of mmes. Each PQL iteration builds the working response and working weights of every row with the family of its trait, and fits a weighted Gaussian multi-trait model.

The residual variance of binomial, Poisson and negative binomial traits is not a free parameter (their dispersion is one), while that of Gaussian and other families is estimated. The residual structure must therefore allow individual trait variances to be fixed: use rcov = ~ vsm(corgm(trait, variance = "heterogeneous"), ism(record)) (trait correlations and heterogeneous variances, see corgm) or rcov = ~ vsm(dsm(trait), ism(units)) (uncorrelated residuals). The first level of the trait factor must be a fixed-dispersion trait, because its variance is the vsm scale; reorder the levels otherwise. mmes then fixes the residual variance of every fixed-dispersion trait to one and estimates the rest.

As for single-family PQL fits, the likelihood, AIC and BIC are NA. The fitted object has a named dispersion vector (one per trait: one for fixed-dispersion traits, the REML residual variance for the others), and fitted() and residuals() apply the inverse link and variance function of each trait to its rows. Random-effect variances are on the link scale of each trait.

Value

An object of class "sommer_familym".

See Also

mmes, corgm, stackTraits, covmatrix_mmes

Examples

set.seed(12)
ng <- 100; nr <- 4
G <- matrix(c(1, 0.4, 0.4, 0.8), 2)
u <- matrix(rnorm(ng * 2), ng) %*% chol(G)
d <- data.frame(id = factor(rep(paste0("g", 1:ng), each = nr)))
d$disease <- rbinom(nrow(d), 1, plogis(-0.5 + u[as.integer(d$id), 1]))
d$yield <- 10 + u[as.integer(d$id), 2] + rnorm(nrow(d))
L <- stackTraits(d, traits = c("disease", "yield"))

fit <- mmes(value ~ trait,
            random = ~ vsm(usm(trait), ism(id)),
            rcov = ~ vsm(corgm(trait, variance = "heterogeneous"), ism(record)),
            family = familym("trait", disease = binomial(), yield = gaussian()),
            data = L, verbose = FALSE)
fit$dispersion
covmatrix_mmes(fit, 1)$correlation   # genetic correlation (link scale)
head(fitted(fit))                    # probabilities for disease rows

Description

fitted method for class "mmes".

Usage

## S3 method for class 'mmes'
fitted(object, type=c("response", "link"), ...)

Arguments

Value

For Gaussian fits, fitted values of the form y.hat = Xb + Zu. For non-Gaussian PQL fits, response-scale conditional means by default, or Xb + Zu with type="link".

Author(s)

Giovanny Covarrubias

See Also

fitted, mmes

Examples

# data(DT_cpdata, package="enhancer")
# DT <- DT_cpdata
# GT <- GT_cpdata
# MP <- MP_cpdata
# #### create the variance-covariance matrix
# A <- A.mat(GT) # additive relationship matrix
# #### look at the data and fit the model
# head(DT)
# mix1 <- mmes(Yield~1,
#               random=~vsm(ism(id),Gu=A)
#                       + Rowf + Colf + spl2Dc(Row,Col),
#               rcov=~units,
#               data=DT)
# 
# ff=fitted(mix1)
# 
# colfunc <- colorRampPalette(c("steelblue4","springgreen","yellow"))
# lattice::wireframe(`u:Row.fitted`~Row*Col, data=ff$dataWithFitted,  
#           aspect=c(61/87,0.4), drape=TRUE,# col.regions = colfunc,
#           light.source=c(10,0,10))
# lattice::levelplot(`u:Row.fitted`~Row*Col, data=ff$dataWithFitted, col.regions = colfunc)

Description

fixm creates a square matrix with 3's in the diagnals and off-diagonals to quickly specify a fixed constraint in the Gtc argument of the vsm function.

Usage

fixm(x, reps=NULL)

Arguments

Value

$res

a matrix or a list of matrices with the constraints to be provided in the Gtc argument of the vsm function.

Author(s)

Giovanny Covarrubias-Pazaran

References

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

Examples

fixm(4)
fixm(4,2)

Description

Instantiates reproducible synthetic mixed model datasets and formulas spanning factor dimensions, covariance structures, relationship densities, and tasks.

Usage

generate_mme_workload(name = NULL,
                      n_ind = 60L,
                      n_levels = 5L,
                      structure = c("ar1", "ar2", "us", "diag", "fa",
                                    "rr", "spatial", "sar", "car"),
                      relationship = c("none", "dense", "sparse"),
                      missing_rate = 0.0,
                      reml = TRUE,
                      computeCi = 0L,
                      seed = 42L)

Arguments

Value

A list of class "mme_workload" containing data, formulas, REML and computeCi settings, and workload metadata.

See Also

mmes, mmes_control

Examples

# Generate a small 1D autoregressive workload
wl <- generate_mme_workload(
  name = "demo_ar1",
  n_ind = 20,
  n_levels = 4,
  structure = "ar1"
)

Description

THIS FUNCTION IS DEPRECATED. Fits a multivariate/univariate linear mixed model GWAS by likelihood methods (REML), see the Details section below. It uses the mmer function and its core coded in C++ using the Armadillo library to optimize dense matrix operations common in the derect-inversion algorithms. After the model fit extracts the inverse of the phenotypic variance matrix to perform the association test for the "p" markers. Please check the Details section (Model enabled) if you have any issue with making the function run.

The sommer package is updated on CRAN every 3-months due to CRAN policies but you can find the latest source at https://github.com/covaruber/sommer . This can be easily installed typing the following in the R console:

library(devtools)

install_github("covaruber/sommer")

This is recommended since bugs fixes will be immediately available in the GitHub source. 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")

or visit https://covaruber.github.io

Usage

GWAS(fixed, random, rcov, data, weights, W,
    nIters=20, tolParConvLL = 1e-03, tolParInv = 1e-06, 
    init=NULL, constraints=NULL,method="NR", 
    getPEV=TRUE,naMethodX="exclude",
    naMethodY="exclude",returnParam=FALSE, 
    dateWarning=TRUE,date.warning=TRUE,verbose=FALSE,
    stepWeight=NULL, emWeight=NULL,
    M=NULL, gTerm=NULL, n.PC = 0, min.MAF = 0.05, 
    P3D = TRUE)

Arguments

Details

Citation

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

Models Enabled

For details about the models enabled and more information about the covariance structures please check the help page of the package (sommer). In general the GWAS model implemented in sommer to obtain marker effect is a generalized linear model of the form:

b = (X'V-X)X'V-y

with X = ZMi

where: b is the marker effect (dimensions 1 x mt) y is the response variable (univariate or multivariate) (dimensions 1 x nt) V- is the inverse of the phenotypic variance matrix (dimensions nt x nt) Z is the incidence matrix for the random effect selected (gTerm argument) to perform the GWAS (dimensions nt x ut) Mi is the ith column of the marker matrix (M argument) (dimensions u x m)

for t traits, n observations, m markers and u levels of the random effect. Depending if P3D is TRUE or FALSE the V- matrix will be calculated once and used for all marker tests (P3D=TRUE) or estimated through REML for each marker (P3D=FALSE).

vignette('sommer.start')

Bug report and contact

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

If you have any bug report please go to https://github.com/covaruber/sommer or send me an email to address it asap.

Value

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

Author(s)

Giovanny Covarrubias-Pazaran

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

Covarrubias-Pazaran G. 2018. Software update: Moving the R package sommer to multivariate mixed models for genome-assisted prediction. doi: https://doi.org/10.1101/354639

Bernardo Rex. 2010. Breeding for quantitative traits in plants. Second edition. Stemma Press. 390 pp.

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

Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.

Lee, D.-J., Durban, M., and Eilers, P.H.C. (2013). Efficient two-dimensional smoothing with P-spline ANOVA mixed models and nested bases. Computational Statistics and Data Analysis, 61, 22 - 37.

Lee et al. 2015. MTG2: An efficient algorithm for multivariate linear mixed model analysis based on genomic information. Cold Spring Harbor. doi: http://dx.doi.org/10.1101/027201.

Maier et al. 2015. Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder. Am J Hum Genet; 96(2):283-294.

Rodriguez-Alvarez, Maria Xose, et al. Correcting for spatial heterogeneity in plant breeding experiments with P-splines. Spatial Statistics 23 (2018): 52-71.

Searle. 1993. Applying the EM algorithm to calculating ML and REML estimates of variance components. Paper invited for the 1993 American Statistical Association Meeting, San Francisco.

Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208.

Tunnicliffe W. 1989. On the use of marginal likelihood in time series model estimation. JRSS 51(1):15-27.

Zhang et al. 2010. Mixed linear model approach adapted for genome-wide association studies. Nat. Genet. 42:355-360.

Examples

####=========================================####
#### For CRAN time limitations most lines in the 
#### examples are silenced with one '#' mark, 
#### remove them and run the examples using
#### command + shift + C |OR| control + shift + C
####=========================================####
#####========================================####
##### potato example
#####========================================####
# 
# data(DT_polyploid, package="enhancer")
# DT <- DT_polyploid
# GT <- GT_polyploid
# MP <- MP_polyploid
# ####=========================================####
# ####### convert markers to numeric format
# ####=========================================####
# numo <- atcg1234(data=GT, ploidy=4);
# numo$M[1:5,1:5];
# numo$ref.allele[,1:5]
# 
# ###=========================================####
# ###### plants with both genotypes and phenotypes
# ###=========================================####
# common <- intersect(DT$Name,rownames(numo$M))
# 
# ###=========================================####
# ### get the markers and phenotypes for such inds
# ###=========================================####
# marks <- numo$M[common,]; marks[1:5,1:5]
# DT2 <- DT[match(common,DT$Name),];
# DT2 <- as.data.frame(DT2)
# DT2[1:5,]
# 
# ###=========================================####
# ###### Additive relationship matrix, specify ploidy
# ###=========================================####
# A <- A.mat(marks)
# ###=========================================####
# ### run it as GWAS model
# ###=========================================####
# ans2 <- GWAS(tuber_shape~1,
#              random=~vsr(Name,Gu=A),
#              rcov=~units,
#              gTerm = "u:Name",
#              M=marks, data=DT2)
# plot(ans2$scores[,1])

Description

Given a matrix A and a matrix G returns a H matrix with the C++ Armadillo library.

Usage

H.mat(A, G, tau = 1, omega = 1, tolparinv=1e-6)

Arguments

Details

See references

Value

H Matrix with the relationship between the individuals based on pedigree and corrected by molecular information

References

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

Martini, J. W., Schrauf, M. F., Garcia-Baccino, C. A., Pimentel, E. C., Munilla, S., Rogberg-Munoz, A., ... & Simianer, H. (2018). The effect of the H-1 scaling factors tau and omega on the structure of H in the single-step procedure. Genetics Selection Evolution, 50(1), 16.

Examples

####=========================================####
####random population of 200 lines with 1000 markers
####=========================================####
M <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  M[i,] <- sample(c(-1,0,0,1), size=1000, replace=TRUE)
}
rownames(M) <- 1:nrow(M)
v <- sample(1:nrow(M),100)
M2 <- M[v,]

A <- A.mat(M) # assume this is a pedigree-based matrix for the sake of example
G <- A.mat(M2)

H <- H.mat(A,G)
# colfunc <- colorRampPalette(c("steelblue4","springgreen","yellow"))
# hv <- heatmap(H[1:15,1:15], col = colfunc(100),Colv = "Rowv")

Description

Identity covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

ism(x)

Arguments

Details

ism returns $K=I$. It has no covariance-shape parameters; the only unknown covariance scale is the $\sigma^2$ supplied by vsm. It is commonly used as the final term of vsm to indicate an independent main effect.

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

vsm, mmes.

Examples

## Not run: 
vsm(ism(id))
vsm(dsm(environment), ism(genotype))

## End(Not run)

Description

Reconstructs the normalized loadings matrix and specific variances of a fam() or rrm() covariance-shaping factor fitted inside vsm(), on the same scale as the fitted covariance matrix returned in object$theta.

Usage

loadings_mmes(object, term = NULL, varianceScale=TRUE, rotation=TRUE)

Arguments

Details

fam()/rrm() remove their otherwise redundant internal scale so that vsm() owns the single overall variance; the loadings and specific variances reported by loadings_mmes are re-normalized so that

$\Sigma = \sigma^2(\Lambda\Lambda^{\mathsf T}+\Psi)$

reproduces object$theta[[term]] exactly when varianceScale=FALSE, where $\sigma^2$ is object$covPar[[term]][1]. With varianceScale=TRUE, the returned loadings crossproduct plus the diagonal of the returned specific variances directly reproduces that covariance.

Value

A list containing:

loadings

a levels-by-factor matrix $\Lambda$.

specific

a named vector of specific variances $\Psi$ (identically zero for every level in a rrm() term).

sigma2

the fitted overall variance scale.

model

"fa" or "rr".

term

the resolved term name.

See Also

scores_mmes, fam, rrm, vsm, mmes.

Examples

## Not run: 
mix <- mmes(BLUEs ~ trial,
            random = ~ vsm(fam(trial, 2), ism(genotype)),
            rcov = ~ units, data = dt)
loadings_mmes(mix)

## End(Not run)

Description

Moving-average covariance structures of order one or two. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

mam(x, order = 1L, theta = NULL, fixed = NULL)
ma1m(x, theta = 0.15, fixed = FALSE)
ma2m(x, theta = c(0.15, 0.05), fixed = NULL)

Arguments

Details

Using the conventional polynomial $e_t+\theta_1e_{t-1}+\theta_2e_{t-2}$, the autocovariance at lag $h$ is proportional to $\sum_{j=0}^{p-h}\theta_j\theta_{j+h}$ with $\theta_0=1$. The covariance is normalized to correlation scale and is zero beyond lag $p$. Invertibility of the MA polynomial is not required merely to define this covariance matrix. ma1m and ma2m are convenience wrappers.