Package {sommer}

Type:	Package
Title:	Solving Mixed Model Equations in R
Version:	4.4.87
Date:	2026-09-30
Maintainer:	Giovanny Covarrubias-Pazaran <cova_ruber@live.com.mx>
Description:	Structural multivariate-univariate linear mixed model solver for estimation of multiple random effects with unknown variance-covariance structures (e.g., unstructured, diagonal, autoregressive, factor analytic, reduced rank, among others) and known covariance among levels of random effects (e.g., pedigree and genomic relationship matrices) (Covarrubias-Pazaran, 2016 <doi:10.1371/journal.pone.0156744>; Maier et al., 2015 <doi:10.1016/j.ajhg.2014.12.006>; Jensen et al., 1997). REML estimates can be obtained using the Direct-Inversion Newton-Raphson and Direct-Inversion Average Information algorithms for the problems r x r (r being the number of records) or using the Henderson-based average information algorithm for the problem c x c (c being the number of coefficients to estimate).
Depends:	R (≥ 3.5.0), Matrix (≥ 1.6-2), methods, stats, MASS, crayon, enhancer
LazyLoad:	yes
License:	GPL-3
Imports:	Rcpp (≥ 0.12.19)
BugReports:	https://github.com/covaruber/sommer/issues
URL:	https://github.com/covaruber/sommer
LinkingTo:	Rcpp, RcppArmadillo, RcppProgress, RcppEigen, Matrix
SystemRequirements:	METIS graph partitioning library (optional; enables nested-dissection sparse ordering for large relationship matrices, auto-detected at configure time, falls back to AMD ordering when absent)
Suggests:	rmarkdown, knitr, plyr, parallel, orthopolynom, RSpectra, lattice, emmeans (≥ 1.4), estimability, testthat (≥ 3.0.0)
VignetteBuilder:	knitr
Config/testthat/edition:	3
RoxygenNote:	7.3.2
NeedsCompilation:	yes
Packaged:	2026-10-04 04:23:15 UTC; giovannycovarrubias
Author:	Giovanny Covarrubias-Pazaran [aut, cre]
Repository:	CRAN
Date/Publication:	2026-10-04 03:20:02 UTC

Additive relationship matrix

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)

Construct an APY Relationship Precision Matrix

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

Dominance relationship matrix

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)

Epistatic relationship matrix

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"

Genome wide association study analysis

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])

Combined relationship matrix H

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")

Multivariate Newton-Raphson algorithm

Description

Multivariate Newton-Raphson algorithm used behind mmes when the henderson argument is set to FALSE. Algorithm made available for users that want to avoid the user-friendly interface and provide their matrices directly.

Usage

MNR(Y, # response
    X, Gx, # fixed effects
    Z, K, # random effects
    R,  # residual effects
    Ge, GeI, # initial values and constraints
    W, isInvW, # weights matrix
    iters, tolpar, tolparinv, # other params
    ai, pev,
    verbose, retscaled,
    stepweight, emweight
    )

Arguments

Details

This is the Rcpp-coded Direct-Inversion LMM REML algorithm used behind the mmer function.

Value

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

References

Mrode, R. A. (2014). Linear models for the prediction of animal breeding values. Cabi.

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

See Also

mmes – the core function of the package

Examples

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


# response matrix
y <- cbind(imputev(DT$Yield), imputev(DT$Firmness))
# fixed effect incidence matrix
X <- model.matrix(~Rowf, data=DT)
# random effect incidence matrix
Z <- Matrix::sparse.model.matrix(~id - 1, data=DT)
colnames(Z) <- gsub("id","",colnames(Z))
# covariance for random effect
GT <- GT_cpdata
A <- A.mat(GT) # additive relationship matrix
A <- A[colnames(Z),colnames(Z)] # make sure of the order
A <- A + diag(1e-4,nrow(A), nrow(A))
# residual effect incidence matrix (dimensions equal nrow(y))
R1 <- Matrix::Diagonal(n=nrow(y)) 
# weights matrix (dimensions equal nrow(y))
W <- diag(nrow(y)) 
# some other parameters
maxIter=3
stepWeight <- rep(0.9, maxIter)
stepWeight[1:2] <- c(0.5, 0.7)
emWeights <- rep(0,maxIter)
# model fit
res <- MNR( Y=y, # multi-trait response
            X=list(X), Gx=list(diag(2)), # fixed effects
            Z=list(Z), K=list(A), # random effects
            R=list(R1), # residual effects
            Ge=list( (diag(2)*.3)+.15 , diag(2)*0.75 ), # inital vc 
            GeI=list(unsm2(2),diag(2)), # vc constraints
            W=W, isInvW=TRUE, # weights for records
            iters=maxIter,
            tolpar=1e-4, tolparinv=1e-6, # tolerance
            ai=FALSE, pev=FALSE, # algorithm specifics
            verbose=TRUE, 
            retscaled=FALSE, 
            stepweight=stepWeight, # second derivatives weights
            emweight=emWeights # em update weights
)
res$sigma
res$Beta
res$U[[1]]

anova form a GLMM fitted with mmes

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

Antedependence covariance structure

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)

First-order autoregressive covariance structure

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)

Second-order autoregressive covariance structure

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)

Third-order autoregressive covariance structure

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)

Selected-level diagonal covariance structure

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)

Selected-level diagonal covariance structure for mmer and vsr

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

Binomial proportion response with its number of trials

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(), verbose = FALSE)
fit$b

Proper conditional autoregressive spatial covariance structure

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)

coef form a GLMM fitted with mmes

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

Imputing a matrix using correlations

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")

General positive-definite correlation structure

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

Covariance Between Two Random Effects

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.

Covariance and correlation matrix of a fitted term with standard errors

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

Extract native-scale covariance parameters from an mmes fit

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)

Native-scale covariance parameters and standard errors

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)

Compound-symmetry covariance structure

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)

User-specified covariance-component structure for mmer and vsr

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

data frame to matrix

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 )

Diagonal heterogeneous covariance structure

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)

Diagonal covariance structure for mmer and vsr

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

Direct sum of residual covariance structures across sections

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)

Factor-analytic covariance structure

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)

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.

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

## End(Not run)

Different distribution families for the traits of a multi-trait model

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

fitted form a LMM fitted with mmes

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)

fixed indication matrix

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)

Identity covariance structure

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)

Extract factor-analytic loadings from a fitted mmes model

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, where \sigma^2 is object$covPar[[term]][1].

Value

A list containing:

loadings

a levels-by-factor matrix \Lambda.

specific

a named vector of specific variances \Psi (identically 1 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)

Moving-average covariance structures of order one or two

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.

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(ma1m(time), ism(id))
vsm(ma2m(time), ism(id))

## End(Not run)

Matern spatial covariance structure

Description

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

Usage

maternm(x, range = NULL, nu = 0.5, fixed = c(FALSE, FALSE),
        distance = NULL, anisotropy = c("none", "geometric"),
        angle = pi/2, ratio = 1.5)

Arguments

Details

For distance d, the correlation is

K(d)=\frac{2^{1-\nu}}{\Gamma(\nu)}z^{\nu}K_{\nu}(z),\qquad z=\frac{\sqrt{2\nu}d}{\phi},

where \phi is range. Both range and smoothness are positive and are optimized on log scales. The implementation uses the generic R evaluator and central factor-level numerical derivatives. Repeated coordinate rows are mapped to the same covariance 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

metricm, sar, car, vsm.

Examples

## Not run: 
vsm(maternm(cbind(xcoord,ycoord)), ism(genotype))
vsm(maternm(time), ism(id))

## End(Not run)

Marker effects recovered from a fitted relationship-based term

Description

Back-transforms the BLUPs of a random term fitted with a genomic relationship matrix (GBLUP) into BLUPs of marker effects, optionally with standard errors and back-solved GWAS-type tests. This is the post-fit equivalent of ASReml's mef argument and does not require refitting the model.

Usage

meffects_mmes(object, term, M, method = c("A.mat", "custom"),
              min.MAF = 0, Z = NULL, scale = NULL, blend = 0, Gu = NULL,
              se = FALSE, chunk = 1000)

Arguments

Details

If the fitted relationship matrix is G=(1-\lambda)sZZ^\top+\lambda I, the BLUP of the marker effects is

\hat a=(1-\lambda)s\,Z^\top G^{-1}\hat u,

using the same G^{-1} as the fit, so Z\hat a=\hat u-\lambda G^{-1}\hat u (exactly \hat u when \lambda=0). For structured terms such as vsm(dsm(env), ism(id), Gu=Gi) effects are returned per coordinate (environment, trait, ...).

With se=TRUE, following Gualdron-Duarte et al. (2014),

\mathrm{Var}(\hat a)=H^\top(\sigma^2_c G - C^{uu})H,\qquad H=(1-\lambda)sG^{-1}Z,

where \sigma^2_c is the variance of the coordinate and C^{uu} the prediction error covariance of the term; z=\hat a/\mathrm{SE} gives back-solved GWAS tests. Only diagonal elements are formed: the C^{uu} part uses solves with the stored coefficient matrix in blocks of chunk markers. The cost grows as O(mn^2) for m markers and n individuals.

The rebuild of G^{-1} requires the data used in the original call; supply Gu otherwise. Markers removed by the filters are returned as NA.

Value

A data frame with columns marker, coordinate (structured terms only), effect and, if se=TRUE, se, z and p.value. Attributes scale and markersKept.

References

Gualdron-Duarte JL, Cantet RJC, Bates RO, Ernst CW, Raney NE, Steibel JP (2014). Rapid screening for phenotype-genotype associations by linear transformations of genomic evaluations. BMC Bioinformatics 15: 246.

See Also

A.mat, mmes, vsm

Examples

set.seed(7)
M <- matrix(sample(c(-1, 0, 1), 40 * 80, TRUE), 40)
rownames(M) <- paste0("g", 1:40); colnames(M) <- paste0("snp", 1:80)
G <- 0.95 * A.mat(M) + 0.05 * diag(40)
Gi <- solve(G); attr(Gi, "inverse") <- TRUE
d <- data.frame(id = factor(rep(rownames(M), each = 2), levels = rownames(M)))
d$y <- rnorm(80)
fit <- mmes(y ~ 1, random = ~ vsm(ism(id), Gu = Gi), data = d, verbose = FALSE)
head(meffects_mmes(fit, 1, M, blend = 0.05, se = TRUE))

Distance-based (metric) spatial correlation structures

Description

Exponential, power, gaussian, spherical and circular correlation structures computed from distances between locations, isotropic or anisotropic. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

metricm(x, model = c("exponential", "power", "gaussian", "spherical", "circular"),
        range = NULL, rho = NULL, anisotropy = c("none", "product", "geometric"),
        angle = pi/2, ratio = 1.5, metric = c("euclidean", "manhattan"),
        distance = NULL, fixed = NULL)

Arguments

Details

With h=d/\phi:

Product anisotropy (exponential, power and gaussian) corresponds to ASReml aexp/agau; geometric anisotropy corresponds to ASReml's anisotropic Matern and is also available in maternm. The power model is the exponential model with \phi=-1/\log\rho, and the exponential model equals maternm(nu=0.5).

Ranges and ratios are optimized on log scales, correlations and angles through bounded logits, with analytic derivatives. The structure is dense over the q locations, so it is intended for up to a few thousand locations. On a complete row-by-column grid, product anisotropy is a Kronecker product of two 1-D structures; vsm(metricm(row, "power"), metricm(col, "power"), ism(units)) fits the same model while keeping the compact grid engine. The gaussian model is nearly singular for close locations; keep a residual (nugget) term. With geometric anisotropy the angle is not identifiable when the ratio is close to one, and the ratio can drift to extreme values when the data are isotropic.

Value

A list with the incidence matrix Z and a CovarianceFactor descriptor covFactor, for use inside vsm.

See Also

maternm, ar1m, sar, car, vsm.

Examples

set.seed(1)
co <- cbind(x = runif(30, 0, 10), y = runif(30, 0, 10))
d <- data.frame(x = rep(co[, 1], each = 2), y = rep(co[, 2], each = 2), one = factor(1))
d$yield <- rnorm(60)
fit <- mmes(yield ~ 1, random = ~ vsm(metricm(cbind(x, y), model = "power"), ism(one)),
            data = d, verbose = FALSE)
covparams_mmes(fit)

mixed model equations for r records

Description

THIS FUNCTION HAS BEEN DEPRECATD PLEASE EXPLORE mmes. The mmer function uses the Direct-Inversion Newton-Raphson or Average Information coded in C++ using the Armadillo library to optimize dense matrix operations common in genomic selection models. These algorithms are intended to be used for problems of the type c > r (more coefficients to estimate than records in the dataset) and/or dense matrices. For problems with sparse data, or problems of the type r > c (more records in the dataset than coefficients to estimate), the MME-based algorithm in the mmes function is faster and we recommend to shift to use that function.

Usage

mmer(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=TRUE, reshapeOutput=TRUE, stepWeight=NULL,
     emWeight=NULL, contrasts=NULL)

Arguments

Details

The use of this function requires a good understanding of mixed models. Please review the 'sommer.quick.start' vignette and pay attention to details like format of your random and fixed variables (e.g. character and factor variables have different properties when returning BLUEs or BLUPs, please see the 'sommer.changes.and.faqs' vignette).

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")

Citation

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

Special variance structures

vsr(atr(x,levels),y)

can be used to specify heterogeneous variance for the "y" covariate at specific levels of the covariate "x", e.g., random=~vsr(at(Location,c("A","B")),ID) fits a variance component for ID at levels A and B of the covariate Location.

vsr(dsr(x),y)

can be used to specify a diagonal covariance structure for the "y" covariate for all levels of the covariate "x", e.g., random=~vsr(dsr(Location),ID) fits a variance component for ID at all levels of the covariate Location.

vsr(usr(x),y)

can be used to specify an unstructured covariance structure for the "y" covariate for all levels of the covariate "x", e.g., random=~vsr(usr(Location),ID) fits variance and covariance components for ID at all levels of the covariate Location.

vsr(overlay(...,rlist=NULL,prefix=NULL))

can be used to specify overlay of design matrices between consecutive random effects specified, e.g., random=~vsr(overlay(male,female)) overlays (overlaps) the incidence matrices for the male and female random effects to obtain a single variance component for both effects. The 'rlist' argument is a list with each element being a numeric value that multiplies the incidence matrix to be overlayed. See overlay for details.Can be combined with vsr().

vsr(leg(x,n),y)

can be used to fit a random regression model using a numerical variable x that marks the trayectory for the random effect y. The leg function can be combined with the special functions dsr, usr at and csr. For example random=~vsr(leg(x,1),y) or random=~vsr(usr(leg(x,1)),y).

vsr(x,Gtc=fcm(v))

can be used to constrain fixed effects in the multi-response mixed models. This is a vector that specifies if the fixed effect is to be estimated for such trait. For example fixed=cbind(response.i, response.j)~vsr(Rowf, Gtc=fcm(c(1,0))) means that the fixed effect Rowf should only be estimated for the first response and the second should only have the intercept.

gvsr(x,y)

can be used to fit variance and covariance parameters between two or more random effects. For example, indirect genetic effect models.

spl2Da(x.coord, y.coord, at.var, at.levels))

can be used to fit a 2-dimensional spline (e.g., spatial modeling) using coordinates x.coord and y.coord (in numeric class) assuming a single variance component. The 2D spline can be fitted at specific levels using the at.var and at.levels arguments. For example random=~spl2Da(x.coord=Row,y.coord=Range,at.var=FIELD).

spl2Db(x.coord, y.coord, at.var, at.levels))

can be used to fit a 2-dimensional spline (e.g., spatial modeling) using coordinates x.coord and y.coord (in numeric class) assuming multiple variance components. The 2D spline can be fitted at specific levels using the at.var and at.levels arguments. For example random=~spl2Db(x.coord=Row,y.coord=Range,at.var=FIELD).

S3 methods

S3 methods are available for some parameter extraction such as fitted.mmer, residuals.mmer, summary.mmer, randef, coef.mmer, anova.mmer, plot.mmer, and predict.mmer to obtain adjusted means. In addition, the vpredict function (replacement of the pin function) can be used to estimate standard errors for linear combinations of variance components (e.g., ratios like h2).

Additional Functions

Additional functions for genetic analysis have been included such as relationship matrix building (A.mat, D.mat, E.mat, H.mat), build a genotypic hybrid marker matrix (build.HMM), plot of genetic maps (map.plot), and manhattan plots (manhattan). If you need to build a pedigree-based relationship matrix use the getA function from the pedigreemm package.

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, just make sure you have read the vignettes carefully before sending your question.

Example Datasets

The package has been equiped with several datasets to learn how to use the sommer package:

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

* DT_ige dataset to show how to fit indirect genetic effect models.

Models Enabled

For details about the models enabled and more information about the covariance structures please check the help page of the package (sommer).

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

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

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

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

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

####===========================================####
#### Univariate heterogeneous variance models  ####
####===========================================####

## Compound simmetry (CS) + Diagonal (DIAG) model
ans2 <- mmer(Yield~Env,
             random= ~Name + vsr(dsr(Env),Name),
             rcov= ~ vsr(dsr(Env),units),
             data=DT)
summary(ans2)

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

ans3 <- mmer(Yield~Env,
             random=~ vsr(usr(Env),Name),
             rcov=~vsr(dsr(Env),units), 
             data=DT)
summary(ans3)



####==========================================####
#### Multivariate homogeneous variance models ####
####==========================================####

## Multivariate Compound simmetry (CS) model
DT$EnvName <- paste(DT$Env,DT$Name)
ans4 <- mmer(cbind(Yield, Weight) ~ Env,
              random= ~ vsr(Name, Gtc = unsm(2)) + vsr(EnvName,Gtc = unsm(2)),
              rcov= ~ vsr(units, Gtc = unsm(2)),
              data=DT)
summary(ans4)

####=============================================####
#### Multivariate heterogeneous variance models  ####
####=============================================####

## Multivariate Compound simmetry (CS) + Diagonal (DIAG) model
ans5 <- mmer(cbind(Yield, Weight) ~ Env,
              random= ~ vsr(Name, Gtc = unsm(2)) + vsr(dsr(Env),Name, Gtc = unsm(2)),
              rcov= ~ vsr(dsr(Env),units, Gtc = unsm(2)),
              data=DT)
summary(ans5)

####===========================================####
#### Multivariate unstructured variance models ####
####===========================================####

ans6 <- mmer(cbind(Yield, Weight) ~ Env,
              random= ~ vsr(usr(Env),Name, Gtc = unsm(2)),
              rcov= ~ vsr(dsr(Env),units, Gtc = unsm(2)),
              data=DT)
summary(ans6)

####=========================================####
####=========================================####
#### EXAMPLE SET 2
#### 2 variance components
#### one random effect with variance covariance structure
####=========================================####
####=========================================####

data("DT_cpdata", package="enhancer")
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
head(DT)
GT[1:4,1:4]
#### create the variance-covariance matrix
A <- A.mat(GT)
#### look at the data and fit the model
mix1 <- mmer(Yield~1,
             random=~vsr(id, Gu=A) + Rowf,
             rcov=~units,
             data=DT)
summary(mix1)$varcomp

#### multi trait example
mix2 <- mmer(cbind(Yield,color)~1,
              random = ~ vsr(id, Gu=A, Gtc = unsm(2)) + # unstructured at trait level
                            vsr(Rowf, Gtc=diag(2)) + # diagonal structure at trait level
                                vsr(Colf, Gtc=diag(2)), # diagonal structure at trait level
              rcov = ~ vsr(units, Gtc = unsm(2)), # unstructured at trait level
              data=DT)
summary(mix2)

mixed model equations solver

Description

Fits linear mixed models by restricted maximum likelihood (REML) or maximum likelihood (ML), and generalized linear mixed models by penalized quasi-likelihood (PQL). Random and residual covariance structures are specified with vsm and covariance constructors. Gaussian fits can use Henderson mixed-model equations or direct inversion in observation space; solver selects LDLT, CHOLMOD or PCG only for the Henderson engine. Alternatively, solveOnly=TRUE solves the mixed-model equations at known covariance parameters without estimating them. Numerical engines are implemented in C++ using Armadillo and Eigen.

Usage

mmes(fixed, random, rcov, data, W, weights=NULL,
     nIters=30, tolParConvLL=1e-04,
     tolParConvNorm=1e-04, tolParInv=1e-06,
     naMethodX="exclude", naMethodY="exclude",
     naMethodRandom="exclude", naMethodR="exclude",
     returnParam=FALSE, dateWarning=TRUE,
     verbose=TRUE, stepWeight=NULL, emWeight=NULL,
     contrasts=NULL, getPEV=TRUE, henderson=TRUE,
    computeCi=0, solver="auto", pcgTol=1.0e-8,
    pcgMaxIters=0, pcgTraceProbes=8,
    pcgLanczosSteps=20, REML=TRUE, vcc=NULL,
    family=stats::gaussian(), pqlControl=list(),
    .pqlInner=FALSE, .pqlFixedDispersion=FALSE,
    .pqlWorkingPrecision=NULL, .pqlBaseW=NULL,
    .pqlBaseFactor=NULL, acceleration="none", .pqlStart=NULL,
     factorScoreAugmentation="none",
     .factorScoreParameters=NULL,
    pcgPreconditioner="diagonal", pcgNystromRank=32L,
    solveOnly=FALSE, covPar=NULL)

Arguments

Details

Choosing a fitting mode

The default, henderson=TRUE, estimates covariance parameters with mixed-model equations in coefficient space. Supply relationship precisions as Gu, with attr(Gu, "inverse")=TRUE. Set henderson=FALSE to estimate parameters with the direct observation-space engine; Gu is still supplied as a relationship precision. The direct engine does not use the Henderson solver or PCG controls. Non-Gaussian families use an outer PQL iteration around Gaussian working fits, not exact GLMM maximum likelihood.

For known-parameter prediction, use solveOnly=TRUE and covPar. This mode always uses its own matrix-free PCG solver, irrespective of solver and henderson; it does not estimate covariance parameters or return likelihoods or PEVs. returnParam=TRUE, by contrast, prepares the design and covariance descriptors without fitting or solving the model.

Starting values and convergence

When vsm() variance scales are omitted, initial scales are chosen from fixed-effect residual variation when possible. Explicit sigma2 values are retained; fixedSigma2=TRUE prevents estimation of that scale. Covariance constructors control their own starting shape parameters and which of those parameters are fixed.

Henderson fits stop when either the absolute change in log-likelihood is below tolParConvLL or the parameter convergence norm is below tolParConvNorm. Reaching nIters alone is not convergence; inspect convergence, llik and normMonitor. A nearly unchanged likelihood need not imply equally precise covariance estimates. When comparing equivalent parameterizations, tighten both convergence tolerances and allow enough iterations. PCG trace and log-determinant approximations introduce an additional accuracy tradeoff controlled by the PCG arguments.

Evaluation environments and observation filtering

The current interface follows R-style model evaluation. Variables referenced by fixed, random, and rcov may be supplied as columns of data or resolved from the formula/calling environment. This makes it possible, for example, to keep relationship matrices, external grouping variables, or spatial coordinates outside the analysis data frame.

Missingness is handled through a common observation-selection mechanism. The response, fixed effects, random effects, and residual covariance specification are evaluated against the original observation set and their validity is combined into a single inclusion mask. The same retained rows are then used for the response vector, fixed-effect design matrix, random-effect design matrices, optional weight matrix, and residual covariance indexing. This avoids independently subsetting different parts of the mixed model.

The fitted object contains obsInfo, an observation-level audit table describing the filtering process. It records the original row index, validity of the response, fixed, random, and residual components, whether the observation was retained, and, when applicable, the reason for exclusion. This can be useful for diagnosing model specifications involving missing values.

Estimated marginal means

When the suggested emmeans package is installed, univariate mmes fits can be supplied directly to emmeans::emmeans(). The resulting means are population-level fixed-effect marginal means with covariance obtained from the mixed-model equations. Inference uses asymptotic degrees of freedom (df = Inf). Multivariate-response fits are not currently supported by this interface.

Random and residual formulas are evaluated as R language objects rather than by splitting their textual representation at plus signs. Consequently, arithmetic or nested expressions inside model terms are not interpreted as separate top-level random or residual terms.

The use of this function requires a good understanding of mixed models. Please review the covariance-structure vignette and pay attention to details like format of your random and fixed variables (e.g. character and factor variables have different properties when returning BLUEs or BLUPs).

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

vignette("sommer.covariance.structures", package="sommer")

vignette("sommer.genetic.evaluation", package="sommer")

vignette("sommer.qg", package="sommer")

vignette("sommer.gxe", package="sommer")

Citation

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

Special variance structures

vsm(atm(x,levels),ism(y))

can be used to specify heterogeneous variance for the "y" covariate at specific levels of the covariate "x", e.g., random=~vsm(atm(Location,c("A","B")),ism(ID)) fits a variance component for ID at levels A and B of the covariate Location.

vsm(dsm(x),ism(y))

can be used to specify a diagonal covariance structure for the "y" covariate for all levels of the covariate "x", e.g., random=~vsm(dsm(Location),ism(ID)) fits a variance component for ID at all levels of the covariate Location.

vsm(usm(x),ism(y))

can be used to specify an unstructured covariance structure for the "y" covariate for all levels of the covariate "x", e.g., random=~vsm(usm(Location),ism(ID)) fits variance and covariance components for ID at all levels of the covariate Location.

vsm(rrm(Location, 2), ism(ID))

fits the reduced-rank covariance model with two factors across locations for ID. Use fam(Location, 2) instead to estimate location-specific variances as well as factor loadings.

vsm(ism(overlay(...,rlist=NULL,prefix=NULL)))

can be used to specify overlay of design matrices between consecutive random effects specified, e.g., random=~vsm(ism(overlay(male,female))) overlays (overlaps) the incidence matrices for the male and female random effects to obtain a single variance component for both effects. The 'rlist' argument is a list with each element being a numeric value that multiplies the incidence matrix to be overlayed. See overlay for details.Can be combined with vsm().

vsm(ism(redmm(x,M,nPC)))

can be used to create a reduced model matrix of an effect (x) assumed to be a linear function of some feature matrix (M), e.g., random=~vsm(ism(redmm(x,M))) creates an incidence matrix from a very large set of features (M) that belong to the levels of x to create a reduced model matrix. See redmm for details.Can be combined with vsm().

vsm(leg(x,n),ism(y))

can be used to fit a random regression model using a numerical variable x that marks the trayectory for the random effect y. The leg function can be combined with the special functions dsm, usm at and csm. For example random=~vsm(leg(x,1),ism(y)) or random=~vsm(usm(leg(x,1)),ism(y)).

spl2Dc(x.coord, y.coord, at.var, at.levels)

can be used to fit a 2-dimensional spline (e.g., spatial modeling) using coordinates x.coord and y.coord (in numeric class) assuming multiple variance components. The 2D spline can be fitted at specific levels using the at.var and at.levels arguments. For example random=~spl2Dc(x.coord=Row,y.coord=Range,at.var=FIELD).

Covariance between random effects

covm( vsm(ism(ran1)), vsm(ism(ran2)) )

can be used to specify covariance between two different random effects, e.g., random=~covm( vsm(ism(x1)), vsm(ism(x2)) ) where two random effects in their own vsm() structure are encapsulated. Only applies for simple random effects.

S3 methods

S3 methods are available for some parameter extraction such as fitted.mmes, residuals.mmes, summary.mmes, randef, coef.mmes, anova.mmes, plot.mmes, and predict.mmes to obtain adjusted means. In addition, the vpredict function (replacement of the pin function) can be used to estimate standard errors for linear combinations of variance components (e.g., ratios like h2). The r2 function calculates reliability.

Additional Functions

Additional functions for genetic analysis have been included such as relationship matrix building (A.mat, D.mat, E.mat, H.mat), build a genotypic hybrid marker matrix (build.HMM), plot of genetic maps (map.plot), and manhattan plots (manhattan). If you need to build a pedigree-based relationship matrix use the getA function from the pedigreemm package.

Bug report and contact

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

If you have any bug report please go to https://github.com/covaruber/sommer or send me an email to address it asap, just make sure you have read the vignettes carefully before sending your question. If you need professional consulting support visit https://covaruber.github.io/sommer/

Example Datasets

The package has been equiped with several datasets to learn how to use the sommer package:

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

* DT_ige dataset to show how to fit indirect genetic effect models.

Models Enabled

For details about the models enabled and more information about the covariance structures please check the help page of the package (sommer).

Value

A Gaussian or PQL fit returns an S3 object of class mmes. The main components are listed below; some are engine-specific or depend on computeCi. PQL fits additionally return family, dispersion, deviance and pqlMonitor; their likelihood and information criteria are not GLMM likelihood quantities. With returnParam=TRUE, a setup list is returned instead, including design matrices, covariance descriptors, vcParams and optimizer controls. With solveOnly=TRUE, the return value has class mmesSolve; see the solveOnly argument for its components and limitations.

Author(s)

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

References

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

Jensen, J., Mantysaari, E. A., Madsen, P., and Thompson, R. (1997). Residual maximum likelihood estimation of (co) variance components in multivariate mixed linear models using average information. Journal of the Indian Society of Agricultural Statistics, 49, 215-236.

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

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

Examples

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

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

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

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

## Known variances: solve the MME only (e.g. genetic evaluation)
sol1 <- mmes(Yield~Env,
             random= ~ Name + Env:Name,
             rcov= ~ units,
             data=DT, solveOnly=TRUE, covPar=ans1)
head(sol1$uList[[1]])



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

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

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

User-defined covariance structure

Description

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

Usage

ownm(x, K = NULL, fun = NULL, par = numeric(), fixed = NULL,
  dfun = NULL, par_names = NULL, native_report = NULL)

Arguments

Details

There are two modes. With K, ownm creates a fixed covariance factor. With fun, the matrix is evaluated at every optimizer point. The returned matrix is symmetrized and normalized by its first diagonal element so that vsm retains the unique overall variance scale. If dfun is supplied, its derivative is normalized analytically using the quotient rule; otherwise a central numerical derivative of the normalized factor is used. This is the general extension mechanism for covariance structures that do not require a native C++ evaluator. When native_report is omitted, the overall variance and naturally transformed par values are reported.

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, covparams_mmes, rrm, maternm, sar, car.

Examples

## Not run: 
K <- matrix(c(1,.3,.3,1),2,2)
vsm(ownm(group, K=K), ism(id))

cfun <- function(p) {
  r <- tanh(p[1]); matrix(c(1,r,r,1),2,2)
}
native <- function(scale, par, factor, absorb_scale=TRUE) {
  c(variance=scale, correlation=tanh(par[1]))
}
vsm(ownm(group, fun=cfun, par=0, native_report=native), ism(id))

## End(Not run)

plot form a LMM plot with mmes

Description

plot method for class "mmes".

Usage

## S3 method for class 'mmes'
plot(x,stnd=TRUE, ...)

Arguments

Value

vector of plot

Author(s)

Giovanny Covarrubias covarrubiasp@wisc.edu

See Also

plot, mmes

Examples

data(DT_yatesoats, package="enhancer")
DT <- DT_yatesoats
head(DT)
m3 <- mmes(fixed=Y ~ V + N + V:N,
           random = ~ B + B:MP,
           rcov=~units,
           data = DT)

plot the change of VC across iterations

Description

plot for monitoring.

Usage

pmonitor(object, ...)

Arguments

Value

vector of plot

Author(s)

Giovanny Covarrubias

See Also

plot, mmes

Examples

data(DT_yatesoats, package="enhancer")
DT <- DT_yatesoats
head(DT)
m3 <- mmes(fixed=Y ~ V + N + V:N,
           random = ~ B + B:MP,
           rcov=~units,
           data = DT)
pmonitor(m3)

Post-fit prediction error variances and inverse coefficient matrix

Description

Computes prediction error variances (PEVs) or the complete inverse of the mixed model coefficient matrix after a model has been fitted with mmes. For PEV-only calculations, the function uses the Takahashi sparse inverse subset approach, avoiding computation of the complete inverse of the coefficient matrix.

Usage

postPEV(object, mode = 1L)

Arguments

Details

The function is intended for post-processing models fitted with mmes, allowing prediction error variances or the complete inverse of the mixed model coefficient matrix to be calculated only when they are needed.

When mode = 1, the sparse mixed model coefficient matrix stored in the fitted object is factorized and the Takahashi equations are used to obtain the selected elements of its inverse required for the diagonal prediction error variances. The complete inverse is not formed.

When mode = 2, the complete inverse of the mixed model coefficient matrix is computed. The prediction error variances are then obtained from the corresponding diagonal elements of this inverse.

When mode = 0, no matrix factorization or inverse calculation is performed and previously stored inverse-derived results are cleared.

predict.mmes does not require any particular mode: it computes exact standard errors for arbitrary linear combinations of fixed and random effects directly from the stored coefficient matrix C, without using Ci at all (see predict.mmes for details). postPEV() is only needed when the diagonal prediction error variances (uPevList) or the complete inverse (Ci) are wanted for their own sake.

For this post-fit calculation to be available, the mmes model must have been fitted with a version of the Henderson mixed model solver that stores the final sparse coefficient matrix and its scaling factor in the fitted model object. Models can therefore be initially fitted with computeCi = 0 and the PEVs or complete inverse calculated afterwards only if required.

Value

An object of class "mmes" containing the original fitted model together with the requested post-fit inverse information.

With mode = 1, uPevList contains the prediction error variances for the random effects while Ci remains empty.

With mode = 2, Ci contains the complete inverse of the mixed model coefficient matrix and uPevList contains the corresponding prediction error variances for the random effects.

With mode = 0, inverse-derived outputs are cleared.

The returned object also contains CiComputed and CiMode, indicating whether the complete inverse was computed and which calculation mode was used, respectively.

References

Takahashi, K., Fagan, J., and Chin, M. S. (1973). Formation of a sparse bus impedance matrix and its application to short circuit study. Proceedings of the 8th PICA Conference, Minneapolis, Minnesota.

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

####=========================================####
#### Fit a mixed model without computing PEVs
#### or the complete inverse during estimation
####=========================================####

# mod <- mmes(
#   fixed = Yield ~ 1,
#   random = ~ vsm(ism(Genotype)),
#   rcov = ~ units,
#   data = DT,
#   computeCi = 0
# )

####=========================================####
#### Compute PEVs afterwards using the
#### Takahashi sparse inverse subset
####=========================================####

# mod <- postPEV(mod, mode = 1)
# mod$uPevList

####=========================================####
#### Compute the complete inverse afterwards
####=========================================####

# mod <- postPEV(mod, mode = 2)
# mod$Ci
# mod$uPevList

Predict form of a LMM fitted with mmes

Description

predict method for class "mmes".

Usage

## S3 method for class 'mmes'
predict(object, Dtable=NULL, D, levels=NULL, sed=FALSE,
        pairwise=FALSE, adjust="none", df=Inf, ...)

Arguments

Details

This function allows to produce predictions specifying those variables that define the margins of the hypertable to be predicted (argument D). Predictions are obtained for each combination of values of the specified variables that is present in the data set used to fit the model. See vignettes for more details.

Standard errors are exact for any computeCi setting (0, 1, or 2) and do not require object$Ci: Var(D %*% bu) = D C.inv() D.t() is instead obtained by solving C x = d for each row d of D against the stored coefficient matrix object$C, which mmes() always returns regardless of computeCi. This is cheaper than forming the complete inverse (computeCi=2) and, unlike the Takahashi selected-inverse subset (computeCi=1), is exact for arbitrary linear combinations, not just diagonal PEVs.

For predicted values the pertinent design matrices X and Z together with BLUEs (b) and BLUPs (u) are multiplied and added together.

predicted.value equal Xb + Zu.1 + ... + Zu.n

For computing standard errors for predictions the parts of the coefficient matrix:

C11 equal (X.t() V.inv() X).inv()

C12 equal 0 - [(X.t() V.inv() X).inv() X.t() V.inv() G Z]

C22 equal PEV equal G - [Z.t() G[V.inv() - (V.inv() X X.t() V.inv() X V.inv() X)]G Z.t()]

In practive C equals ( W.t() V.inv() W ).inv()

when both fixed and random effects are present in the inclusion set. If only fixed and random effects are included, only the respective terms from the SE for fixed or random effects are calculated.

Value

Author(s)

Giovanny Covarrubias-Pazaran

References

Welham, S., Cullis, B., Gogel, B., Gilmour, A., and Thompson, R. (2004). Prediction in linear mixed models. Australian and New Zealand Journal of Statistics, 46, 325 - 347.

See Also

predict, mmes

Examples

data(DT_yatesoats, package="enhancer")
DT <- DT_yatesoats
m3 <- mmes(fixed=Y ~ V + N + V:N ,
           random = ~ B + B:MP,
           rcov=~units,
           data = DT)

m3 <- postPEV(m3, mode = 2)
 
#############################
## predict means for nitrogen
#############################
Dt <- m3$Dtable; Dt
# first fixed effect just average
Dt[1,"average"] = TRUE
# second fixed effect include
Dt[2,"include"] = TRUE
# third fixed effect include and average
Dt[3,"include"] = TRUE
Dt[3,"average"] = TRUE
Dt

pp=predict(object=m3, Dtable=Dt, D="N")
pp$pvals

#############################
## predict means for variety
#############################

Dt <- m3$Dtable; Dt
# first fixed effect include
Dt[1,"include"] = TRUE
# second fixed effect just average
Dt[2,"average"] = TRUE
# third fixed effect include and average
Dt[3,"include"] = TRUE
Dt[3,"average"] = TRUE
Dt

pp=predict(object=m3, Dtable=Dt, D="V")
pp$pvals

#############################
## predict means for nitrogen:variety
#############################
# prediction matrix D based on (equivalent to classify in asreml)
Dt <- m3$Dtable; Dt
# first fixed effect include and average
Dt[1,"include"] = TRUE
Dt[1,"average"] = TRUE
# second fixed effect include and average
Dt[2,"include"] = TRUE
Dt[2,"average"] = TRUE
# third fixed effect include and average
Dt[3,"include"] = TRUE
Dt[3,"average"] = TRUE
Dt

pp=predict(object=m3, Dtable=Dt, D="N:V")
pp$pvals

Reliability

Description

Calculates the reliability of BLUPs in a sommer model.

Usage

r2(object, object2=NULL)

Arguments

Details

The reliability method calculated is the classical animal model: R2=(G-PEV)/G

Value

result

a list with as many elements as random effects fitted containing reliabilities for individual BLUPs.

References

Mrode, R. A. (2014). Linear models for the prediction of animal breeding values. Cabi.

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

See Also

mmes – the core function of the package

Examples

####=========================================####
#### Example population
####=========================================####
data(DT_example, package="enhancer")
DT <- DT_example
head(DT)
ans1 <- mmes(Yield~Env, 
             random= ~ Name + Env:Name,
             rcov= ~ units,
             data=DT)
ans1 <- postPEV(ans1, mode = 1)
rel=r2(ans1)

extracting random effects

Description

This function is extracts the random effects from a mixed model fitted by mmer.

Usage

randef(object)

Arguments

Value

$randef

a list structure with the random effects or BLUPs.

Examples

# randef(model)

Residuals form a GLMM fitted with mmes

Description

residuals method for class "mmes".

Usage

## S3 method for class 'mmes'
residuals(object,
                          type=c("response", "deviance", "pearson", "working"), ...)

Arguments

Value

For Gaussian fits, residuals of the form e = y - Xb - Zu. For non-Gaussian PQL fits, residuals on the requested scale.

Author(s)

Giovanny Covarrubias

See Also

residuals, mmes

Back-transform rotated mmes random effects

Description

Back-transforms the selected relationship random effect from the solver eigenbasis to the original relationship-level basis.

Usage

rotate_back_mmes(object)

Arguments

Details

If Gu=U\Lambda U' and a denotes the solver-basis effect, the reported effect is g=Ua. Fitted mmes objects already expose this original-basis effect in uList; this function provides the explicit transformation from the retained engine output.

Value

A matrix of original-basis random effects, with relationship levels in rows and covariance-product coordinates in columns.

See Also

mmes, vsm

Reduced-rank covariance structure

Description

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

Usage

rrm(x, k = 1L, loadings = NULL, fixed = NULL)

Arguments

Details

The reduced-rank structure uses

M=\Lambda\Lambda^{\mathsf T}+I_q,\qquad K=M/M_{11}.

The rank-k term captures the dominant covariance pattern and the identity term provides a common isotropic remainder, keeping the covariance strictly positive definite for precision-based Henderson calculations. The leading k\times k loading block is lower triangular for rotational identification, and its diagonal loadings are positive. rrm is compiled through the generic CovarianceFactor callback interface rather than requiring model-specific solver code.

Value

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

See Also

fam, ownm, vsm.

Examples

## Not run: 
vsm(rrm(environment, k=2), ism(genotype))

## End(Not run)

Simultaneous autoregressive spatial covariance structure

Description

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

Usage

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

Arguments

Details

The simultaneous autoregressive covariance is constructed from

B=I-\rho W,\qquad M=B^{-1}B^{-\mathsf T},\qquad K=M/M_{11}.

For a general real weights matrix, rho is restricted to the conservative interval (-1/r(W),1/r(W)), where r(W) is the spectral radius. A bounded-logit working coordinate enforces the interval. The first covariance derivative is supplied analytically to the generic CovarianceFactor engine.

Value

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

See Also

car, maternm, vsm.

Examples

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

## End(Not run)

Predict factor-analytic scores from a fitted mmes model

Description

Predicts latent factor scores for each level of the main random-effect term crossed with a fam() or rrm() covariance-shaping factor, combining the fitted loadings (see loadings_mmes), the fitted covariance matrix, and the BLUPs already stored in object$uList.

Usage

scores_mmes(object, term = NULL, method = c("regression", "bartlett"), 
            varianceScale=TRUE, rotation=TRUE)

Arguments

Details

Let U be the levels-by-levels BLUP matrix in object$uList[[term]] and \Lambda the loadings from loadings_mmes. The regression (Thomson) scores are

F = U\Sigma^{-1}\Lambda,

and the Bartlett scores are

F = U\Psi^{-1}\Lambda(\Lambda^{\mathsf T}\Psi^{-1}\Lambda)^{-1}.

Because U already contains shrunken (BLUP) effects, the resulting scores are similarly regularized.

Value

A matrix with one row per level of the main random-effect term and one column per latent factor.

See Also

loadings_mmes, fam, rrm, vsm, mmes.

Examples

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

## End(Not run)

Solving Mixed Model Equations in R

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)

Two-dimensional penalised tensor-product of marginal B-Spline basis.

Description

Auxiliary function used for modelling the spatial or environmental effect as a two-dimensional penalised tensor-product (isotropic approach) based on Lee et al. (2013) and Rodriguez-Alvarez et al. (2018). This is a modified wrapper of some portions of the SpATS package to build a single incidence matrix containing all the columns from tensor products of the x and y coordinates and it fits such matrix as a single random effect. Then the heterogeneous covariances structure capabilities of sommer can be used to enhance the model fit. You may be interested in reading and citing not only sommec but also Wageningen publications if using this 2D spline methodology.

Usage

spl2Dc(x.coord,y.coord,at.var=NULL,at.levels=NULL, type="PSANOVA", 
      nsegments = c(10,10), penaltyord = c(2,2), degree = c(3,3), 
      nestorder = c(1,1), thetaC=NULL, theta=NULL, sp=FALSE)

Arguments

Details

The following documentation is taken from the SpATS package. Please refer to this package and associated publications if you are interested in going deeper on this technique:

Within the P-spline framework, anisotropic low-rank tensor-product smoothers have become the general approach for modelling multidimensional surfaces (Eilers and Marx 2003; Wood 2006). In the original SpATS package, was proposed to model the spatial or environmental effect by means of the tensor-product of B-splines basis functions. In other words, was proposed to model the spatial trend as a smooth bivariate surface jointly defined over the the spatial coordinates. Accordingly, the current function has been designed to allow the user to specify the spatial coordinates that the spatial trend is a function of. There is no restriction about how the spatial coordinates shall be specified: these can be the longitude and latitude of the position of the plot on the field or the column and row numbers. The only restriction is that the variables defining the spatial coordinates should be numeric (in contrast to factors).

As far as estimation is concerned, we have used in this package the equivalence between P-splines and linear mixed models (Currie and Durban, 2002). Under this approach, the smoothing parameters are expressed as the ratio between variance components. Moreover, the smooth components are decomposed in two parts: one which is not penalised (and treated as fixed) and one with is penalised (and treated as random). For the two-dimensional case, the mixed model representation leads also to a very interesting decomposition of the penalised part of the bivariate surface in three different components (Lee and Durban, 2011): (a) a component that contains the smooth main effect (smooth trend) along one of the covariates that the surface is a function of (as, e.g, the x-spatial coordinate or column position of the plot in the field), (b) a component that contains the smooth main effect (smooth trend) along the other covariate (i.e., the y-spatial coordinate or row position); and (c) a smooth interaction component (sum of the linear-by-smooth interaction components and the smooth-by-smooth interaction component).

The original implementation of SpATS assumes two different smoothing parameters, i.e., one for each covariate in the smooth component. Accordingly, the same smoothing parameters are used for both, the main effects and the smooth interaction. However, this approach can be extended to deal with the ANOVA-type decomposition presented in Lee and Durban (2011). In their approach, four different smoothing parameters are considered for the smooth surface, that are in concordance with the aforementioned decomposition: (a) two smoothing parameter, one for each of the main effects; and (b) two smoothing parameter for the smooth interaction component.

It should be noted that, the computational burden associated with the estimation of the two-dimensional tensor-product smoother might be prohibitive if the dimension of the marginal bases is large. In these cases, Lee et al. (2013) propose to reduce the computational cost by using nested bases. The idea is to reduce the dimension of the marginal bases (and therefore the associated number of parameters to be estimated), but only for the smooth-by-smooth interaction component. As pointed out by the authors, this simplification can be justified by the fact that the main effects would in fact explain most of the structure (or spatial trend) presented in the data, and so a less rich representation of the smooth-by-smooth interaction component could be needed. In order to ensure that the reduced bivariate surface is in fact nested to the model including only the main effects, Lee et al. (2013) show that the number of segments used for the nested basis should be a divisor of the number of segments used in the original basis (nsegments argument). In the present function, the divisor of the number of segments is specified through the argument nestorder. For a more detailed review on this topic, see Lee (2010) and Lee et al. (2013). The "PSANOVA" approach represents an alternative method. In this case, the smooth bivariate surface (or spatial trend) is decomposed in five different components each of them depending on a single smoothing parameter (see Lee et al., 2013).

—————–

As mentioned at the beginning, the piece of documentation stated above was taken completely from the SpATS package in order to provide a deeper explanation. In practice, sommec uses some pieces of code from SpATS to build the design matrix containing all the columns from tensor products of the x and y coordinates and it fits such matrix as a single random effect. As a result the same variance component is assumed for the linear, linear by linear, linear by spline, and spline by spline interactions. This results in a less flexible approach than the one proposed by Rodriguez-Alvarez et al. (2018) but still makes a pretty good job to model the spatial variation. Use under your own risk.

References

Rodriguez-Alvarez, M.X, Boer, M.P., van Eeuwijk, F.A., and Eilers, P.H.C. (2018). SpATS: Spatial Analysis of Field Trials with Splines. R package version 1.0-9. https://CRAN.R-project.org/package=SpATS.

Rodriguez-Alvarez, M.X., et al. (2015) Fast smoothng parameter separaton n multdmensonal generalzed P-splnes: the SAP algorthm. Statistics and Computing 25.5: 941-957.

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.

Gilmour, A.R., Cullis, B.R., and Verbyla, A.P. (1997). Accounting for Natural and Extraneous Variation in the Analysis of Field Experiments. Journal of Agricultural, Biological, and Environmental Statistics, 2, 269 - 293.

See Also

mmes

Examples

## ============================ ##
## example to use spl2Dc() 
## ============================ ## 
data(DT_cpdata, package="enhancer")
# DT <- DT_cpdata
# GT <- GT_cpdata
# MP <- MP_cpdata
# A <- A.mat(GT)
## ============================ ##
## mimic 2 fields
## ============================ ## 
# aa <- DT; bb <- DT
# aa$FIELD <- "A";bb$FIELD <- "B"
# set.seed(1234)
# aa$Yield <- aa$Yield + rnorm(length(aa$Yield),0,4)
# DT2 <- rbind(aa,bb)
# head(DT2)
# mix <- mmes(Yield~1, henderson = F,
#             random=~
#               vsm(dsm(FIELD),ism(Rowf)) +
#               vsm(dsm(FIELD),ism(Colf)) +
#                 spl2Dc(Row,Col,at.var=FIELD),
#             rcov=~vsm(dsm(FIELD),ism(units)),
#             data=DT2)
# 
# # extract spatial effects
# blup <- mix$uList$`spl2Dc(Row, Col, at.var = FIELD`
# head(blup) # 2 fields
# # recreate the incidence matrices
# xx=with(DT2, spl2Dc(Row,Col,at.var=FIELD))
# # get fitted values Zu for spatial effects and add them to the dataset
# field1 <- xx$Z$`A:all` %*% blup[,1]
# field2 <- xx$Z$`B:all` %*% blup[,2]
# DT2$spat <- field1+field2
# # plots the spatial effects
# lattice::levelplot(spat~Row*Col|FIELD, data=DT2)

Get Tensor Product Spline Mixed Model Incidence Matrices

Description

spl2Dmats gets Tensor-Product P-Spline Mixed Model Incidence Matrices for use with sommer and its main function mmes. We thank Sue Welham for making the TPSbits package available to the community. If you're using this function for your research please cite her TPSbits package :) this is mostly a wrapper of her tpsmmb function to enable the use in sommer.

Usage

spl2Dmats(
  x.coord.name,
  y.coord.name,
  data,
  at.name,
  at.levels, 
  nsegments=NULL,
  minbound=NULL,
  maxbound=NULL,
  degree = c(3, 3),
  penaltyord = c(2,2), 
  nestorder = c(1,1),
  method = "Lee"
)

Arguments

Value

List of length 7 elements:

1. data = the input data frame augmented with structures required to fit tensor product splines in asreml-R. This data frame can be used to fit the TPS model.

   Added columns:

   • TP.col, TP.row = column and row coordinates

   • TP.CxR = combined index for use with smooth x smooth term

   • TP.C.n for n=1:(diff.c) = X parts of column spline for use in random model (where diff.c is the order of column differencing)

   • TP.R.n for n=1:(diff.r) = X parts of row spline for use in random model (where diff.r is the order of row differencing)

   • TP.CR.n for n=1:((diff.c*diff.r)) = interaction between the two X parts for use in fixed model. The first variate is a constant term which should be omitted from the model when the constant (1) is present. If all elements are included in the model then the constant term should be omitted, eg. y ~ -1 + TP.CR.1 + TP.CR.2 + TP.CR.3 + TP.CR.4 + other terms...

   • when asreml="grp" or "sepgrp", the spline basis functions are also added into the data frame. Column numbers for each term are given in the grp list structure.

2. fR = Xr1:Zc

3. fC = Xr2:Zc

4. fR.C = Zr:Xc1

5. R.fC = Zr:Xc2

6. fR.fC = Zc:Zr

7. all = Xr1:Zc | Xr2:Zc | Zr:Xc1 | Zr:Xc2 | Zc:Zr

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

M <- spl2Dmats(x.coord.name = "Col", y.coord.name = "Row", data=DT, nseg =c(14,21))
head(M$data)
# m1g <- mmes(Yield~1+TP.CR.2+TP.CR.3+TP.CR.4,
#             random=~Rowf+Colf+vsm(ism(M$fC))+vsm(ism(M$fR))+
#               vsm(ism(M$fC.R))+vsm(ism(M$C.fR))+vsm(ism(M$fC.fR))+
#               vsm(ism(id),Gu=A),
#             data=M$data, tolpar = 1e-6,
#             iters=30)
# 
# summary(m1g)$varcomp

Stack several trait columns into the long format used by multi-trait models

Description

Converts a wide data frame (one column per trait) into the long format used by multi-trait mmes models: one row per record and trait, a factor identifying the trait, a single response column and a record key that pairs the traits measured on the same experimental unit.

Usage

stackTraits(data, traits, keep = NULL, trait = "trait", value = "value",
            record = "record")

Arguments

Details

Rows are stacked trait by trait, so row i of trait t is row (t-1)n+i of the output. The record key is taken from the row names of data when these are unique and not the default 1..n; otherwise it is "r1", "r2", .... Missing trait values are kept as NA (they are handled by the naMethodY rules of mmes, so incomplete records are allowed).

Using the record key in the residual, rcov = ~ vsm(usm(trait), ism(record)), pairs the residuals of the traits of each unit explicitly. With ism(units) the pairing is implied by the order of the rows within each trait instead.

Value

A data frame with the keep columns and the record, trait and value columns, with length(traits) * nrow(data) rows.

See Also

mmes, covmatrix_mmes, familym

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")
head(DTL)
table(DTL$trait)

# bivariate model with the residual paired through the record key
fit <- mmes(value ~ trait,
            random = ~ vsm(usm(trait), ism(id)),
            rcov = ~ vsm(usm(trait), ism(record)),
            data = DTL, verbose = FALSE)
covmatrix_mmes(fit, 1)$correlation   # genetic correlation

Covariance structure spanning several random terms

Description

Combines two or more random terms built with vsm into one random structure with an estimated covariance between the terms. Typical uses are direct-maternal animal models, correlated permanent-environment effects and indirect genetic effects.

Usage

strm(..., cov = usm, Gu = NULL, labels = NULL, sigma2 = NULL,
     fixedSigma2 = FALSE)

Arguments

Details

For terms u_1,\dots,u_T with incidence matrices Z_1,\dots,Z_T the structure is

\mathrm{Var}(u_1,\dots,u_T)=\sigma^2\,K_{terms}\otimes K_{inner}\otimes A,

where K_{terms} is given by cov, K_{inner} is the Kronecker product of the covariance factors used inside the terms (e.g. dsm(env)), and A is the relationship matrix. The requirements are those of a direct product: all terms share the coefficient levels (their levels are aligned on the union, or on the levels of Gu), the same Gu, and the same inner covariance constructors with the same levels; the parameters of the inner factors are shared. The sigma2 of individual terms is ignored. BLUPs are returned in one uList element with a column per term (and inner level). covm(ran1, ran2) is the two-term special case with cov = usm.

The terms do not need to be adjacent in the model formula and their order is checked by level names rather than assumed.

Value

A random structure for the random formula of mmes.

See Also

covm, vsm, usm, mmes

Examples

set.seed(3)
ids <- paste0("i", 1:30)
A <- diag(30) + diag(0.2, 30); dimnames(A) <- list(ids, ids)
Ai <- solve(A); attr(Ai, "inverse") <- TRUE
d <- data.frame(id = factor(sample(ids, 200, TRUE), levels = ids),
                dam = factor(sample(ids[1:15], 200, TRUE), levels = ids))
d$y <- rnorm(200)
fit <- mmes(y ~ 1, random = ~ strm(dir = vsm(ism(id)), mat = vsm(ism(dam)), Gu = Ai),
            data = d, verbose = FALSE)
covparams_mmes(fit, 1)

summary form a GLMM fitted with mmer

Description

summary method for class "mmer".

Usage

## S3 method for class 'mmer'
summary(object, ...)

Arguments

Value

vector of summary

Author(s)

Giovanny Covarrubias-Pazaran

See Also

summary, mmer

summary form a GLMM fitted with mmes

Description

summary method for class "mmes".

Usage

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

Arguments

Value

vector of summary

Author(s)

Giovanny Covarrubias-Pazaran

See Also

summary, mmes

General positive-definite Toeplitz correlation structure

Description

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

Usage

toeplitzm(x, pacf = NULL, fixed = NULL)

Arguments

Details

A full q by q Toeplitz correlation matrix is parameterized by q-1 reflection coefficients. Each PACF is mapped from an unconstrained working coordinate with tanh; the reflection coefficients are converted to a stable autoregressive representation and then to the implied Toeplitz correlation sequence. This parameterization preserves positive definiteness while allowing every lag correlation to vary indirectly.

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, ar3m, mam, vsm.

Examples

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

## End(Not run)

Get Tensor Product Spline Mixed Model Incidence Matrices

Description

tpsmmbwrapper is a wrapper of tpsmmb function from the TPSbits package to avoid version dependencies but if you're using this function for your research please cite the TPSbits package. This function is internally used by the spl2Dmatrices function to get Tensor-Product P-Spline Mixed Model Bits (design matrices) for use with sommer.

Usage

tpsmmbwrapper(
  columncoordinates,
  rowcoordinates,
  data,
  nsegments=NULL,
  minbound=NULL,
  maxbound=NULL,
  degree = c(3, 3),
  penaltyord = c(2, 2),
  nestorder = c(1, 1),
  asreml = "mbf",
  eigenvalues = "include",
  method = "Lee",
  stub = NULL
)

Arguments

Value

List of length 7, 8 or 9 (according to the asreml and eigenvalues parameter settings).

1. data = the input data frame augmented with structures required to fit tensor product splines in asreml-R. This data frame can be used to fit the TPS model.

   Added columns:

   • TP.col, TP.row = column and row coordinates

   • TP.CxR = combined index for use with smooth x smooth term

   • TP.C.n for n=1:(diff.c) = X parts of column spline for use in random model (where diff.c is the order of column differencing)

   • TP.R.n for n=1:(diff.r) = X parts of row spline for use in random model (where diff.r is the order of row differencing)

   • TP.CR.n for n=1:((diff.c*diff.r)) = interaction between the two X parts for use in fixed model. The first variate is a constant term which should be omitted from the model when the constant (1) is present. If all elements are included in the model then the constant term should be omitted, eg. y ~ -1 + TP.CR.1 + TP.CR.2 + TP.CR.3 + TP.CR.4 + other terms...

   • when asreml="grp" or "sepgrp", the spline basis functions are also added into the data frame. Column numbers for each term are given in the grp list structure.

2. mbflist = list that can be used in call to asreml (so long as Z matrix data frames extracted with right names, eg BcZ<stub>.df)

3. BcZ.df = mbf data frame mapping onto smooth part of column spline, last column (labelled TP.col) gives column index

4. BrZ.df = mbf data frame mapping onto smooth part of row spline, last column (labelled TP.row) gives row index

5. BcrZ.df = mbf data frame mapping onto smooth x smooth term, last column (labelled TP.CxR) maps onto col x row combined index

6. dim = list structure, holding dimension values relating to the model:

   1. "diff.c" = order of differencing used in column dimension

   2. "nbc" = number of random basis functions in column dimension

   3. "nbcn" = number of nested random basis functions in column dimension used in smooth x smooth term

   4. "diff.r" = order of differencing used in column dimension

   5. "nbr" = number of random basis functions in column dimension

   6. "nbrn" = number of nested random basis functions in column dimension used in smooth x smooth term

7. trace = list of trace values for ZGZ' for the random TPspline terms, where Z is the design matrix and G is the known diagonal variance matrix derived from eigenvalues. This can be used to rescale the spline design matrix (or equivalently variance components).

8. grp = list structure, only added for settings asreml="grp", asreml="sepgrp" or asreml="own". For asreml="grp", provides column indexes for each of the 5 random components of the 2D splines. For asreml="sepgrp", provides column indexes for each of the X and Z component matrices for the 1D splines, plus the composite smooth x smooth interaction term. For asreml="own", provides column indexes for the composite random model. Dimensions of the components can be derived from the values in the dim item. The Z terms are scaled by the associated eigenvalues when eigenvalues="include", but not when eigenvalues="omit".

9. eigen = list structure, only added for option setting eigenvalues="omit". Holds the diagonal elements of the inverse variance matrix for the terms Xc:Zr (called diagr), Zc:Xr (called diagc) and Zc:Zr (called diagcr).

unstructured indication matrix

Description

unsm creates a square matrix with ones in the diagonals and 2's in the off-diagonals to quickly specify an unstructured constraint in the Gtc argument of the vsm function.

Usage

  unsm(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

unsm(3)
unsm(3,2)

Unstructured positive-definite covariance structure

Description

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

Usage

usm(x, theta = NULL, fixed = NULL)

Arguments

Details

The covariance shape is parameterized as

K=LL^{\mathsf T},

with lower-triangular L and L_{11}=1. Diagonal elements after the first are positive and represented on log scales; lower off-diagonal elements are unrestricted. This guarantees positive definiteness without repeatedly repairing an unconstrained covariance matrix.

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(usm(environment), ism(genotype))

## End(Not run)

Unstructured covariance structure for mmer and vsr

Description

Creates an unstructured variance-covariance component pattern for use by the mmer covariance-model interface, typically inside vsr.

Usage

usr(x)

Arguments

Details

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

After constructing the design matrix Z, the function calls unsm(q), where q is the number of columns of Z, to create the thetaC pattern for an unstructured covariance matrix. Thus the associated covariance model permits a separate variance for every level and a separate covariance for every pair of levels:

\Sigma = \begin{bmatrix} \sigma_1^2 & \sigma_{12} & \cdots & \sigma_{1q}\\ \sigma_{12} & \sigma_2^2 & \cdots & \sigma_{2q}\\ \vdots & \vdots & \ddots & \vdots\\ \sigma_{1q} & \sigma_{2q} & \cdots & \sigma_q^2 \end{bmatrix}.

The exact component coding is stored in the matrix returned by unsm(q) and passed to the vsr machinery through thetaC.

For a factor or character vector, missing observations are retained through the construction of the incidence matrix. For a numeric non-factor input, x is treated as a single design column.

Value

A list with components:

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

• thetaC: the unstructured covariance-component pattern produced by unsm().

See Also

vsr, dsr, csr, atr, usm, mmer

Examples

# Unstructured covariance pattern across three levels
U <- usr(factor(c("A","B","A","C")))
U$Z
U$thetaC

Equality and scaling constraints between variance parameters

Description

The vcc argument of mmes constrains groups of covariance parameters to be equal or to keep fixed ratios. Parameters are identified through the vcParams table.

Details

Parameters are estimated on transformed scales that guarantee valid covariance matrices (log for variances, atanh or bounded logits for correlations). A constraint is enforced exactly when it is affine on that scale, which covers equalities of any parameter and ratios of positive parameters (a log-scale offset). Relations such as \sigma^2_3=\sigma^2_1+\sigma^2_2 are not affine on the log scale and are not supported.

Internally the free working parameters are written as w=o+T\phi. Each average-information update is projected onto this set in the information metric, which equals the constrained Newton step, and then repaired exactly; both the Henderson and the direct engines support constraints. The covariance of the estimates is T(T^\top\mathcal{I}T)^{-1}T^\top, so constrained parameters share their standard errors.

The vcParams table

mmes(..., returnParam = TRUE)$vcParams (before fitting) and fit$vcParams (after fitting) list every covariance parameter with columns index, term, parameter, kind, start (natural scale), free, structure and position (and group when constraints were used). kind is "scale" for the variance scale (sigma2) owned by each vsm() term, or the transform used to estimate the parameter: "exp" (positive quantities such as ranges and variance ratios), "identity" (e.g. Cholesky or factor-analytic loadings), "tanh" or "bounded_logit" (correlation-type parameters).

Specifying constraints

vcc is a data frame with columns

parameter

the index of the parameter, or its name "term:parameter" (e.g. "vsm(ism(B)):sigma2"); the parameter name alone is accepted when unique.

group

parameters with the same group value are constrained together.

scale

optional, default 1: each member equals scale times the first member of its group (natural scale). Scaling is allowed for "scale", "exp" and "identity" parameters; correlation-type parameters can only be equated.

All members of a group must have the same kind: variance scales can be constrained with other variance scales, and other parameters with parameters estimated through the same transform. If any member of a group is fixed, the whole group is fixed at the value implied by that member. Starting values of the members are made consistent with the first (or fixed) member.

See Also

mmes, vsm, covparams_mmes

Examples

data(DT_yatesoats, package = "enhancer")
DT <- DT_yatesoats
p <- mmes(Y ~ V * N, random = ~ B + B:MP, rcov = ~ units, data = DT,
          returnParam = TRUE)
p$vcParams
# block and whole-plot variances equal
m1 <- mmes(Y ~ V * N, random = ~ B + B:MP, rcov = ~ units, data = DT,
           vcc = data.frame(parameter = c(1, 2), group = 1), verbose = FALSE)
# block variance twice the whole-plot variance
m2 <- mmes(Y ~ V * N, random = ~ B + B:MP, rcov = ~ units, data = DT,
           vcc = data.frame(parameter = c("vsm(ism(B:MP)):sigma2", "vsm(ism(B)):sigma2"),
                            group = 1, scale = c(1, 2)), verbose = FALSE)
unlist(m2$covPar)

vpredict form of a LMM fitted with mmes

Description

vpredict method for class "mmes".

Post-analysis procedure to calculate linear combinations of variance components. Its intended use is when the variance components are either simple variances or are variances and covariances in an unstructured matrix. The functions covered are linear combinations of the variance components (for example, phenotypic variance), a ratio of two components (for example, heritabilities) and the correlation based on three components (for example, genetic correlation).

The calculations are based on the estimated variance parameters and their variance matrix as represented by the inverse of the Fisher or Average information matrix. Note that this matrix has zero values for fixed variance parameters including those near the parameter space boundary.

The transform is specified with a formula. On the left side of the formula is a name for the transformation. On the right side of the formula is a transformation specified with shortcut names like 'V1', 'V2', etc. The easiest way to identify these shortcut names is to use 'summary(object)$varcomp'. The rows of this object can referred to with shortcuts 'V1', 'V2', etc. See the example below.

Usage

vpredict(object, transform)
## S3 method for class 'mmes'
vpredict(object, transform)

Arguments

Details

The delta method (e.g., Lynch and Walsh 1998, Appendix 1; Ver Hoef 2012) uses a Taylor series expansion to approximate the moments of a function of parameters. Here, a second-order Taylor series expansion is implemented to approximate the standard error for a function of (co)variance parameters. Partial first derivatives of the function are calculated by algorithmic differentiation with deriv.

Though vpredict can calculate standard errors for non-linear functions of (co)variance parameters from a fitted mmes model, it is limited to non-linear functions constructed by mathematical operations such as the arithmetic operators +, -, *, / and ^, and single-variable functions such as exp and log. See deriv for more information.

Value

Author(s)

Giovanny Covarrubias

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

Lynch, M. and B. Walsh 1998. Genetics and Analysis of Quantitative Traits. Sinauer Associates, Inc., Sunderland, MA, USA.

Ver Hoef, J.M. 2012. Who invented the delta method? The American Statistician 66:124-127. DOI: 10.1080/00031305.2012.687494

See Also

vpredict, mmes

Examples

####=========================================####
####=========================================####
#### EXAMPLE 1
#### simple example with univariate models
####=========================================####
####=========================================####
# data(DT_cpdata, package="enhancer")
# DT <- DT_cpdata
# GT <- GT_cpdata
# MP <- MP_cpdata
# #### create the variance-covariance matrix 
# A <- A.mat(GT)
# #### look at the data and fit the model
# head(DT)
# mix1 <- mmes(Yield~1,
#               random=~vsm(ism(id),Gu=A), 
#               data=DT)
# summary(mix1)$varcomp
# #### run the vpredict function
# vpredict(mix1, h2 ~ V1 / ( V1 + V2 ) )

variance structure specification

Description

vs DEPRECATED NOW. was the main function to build the variance-covariance structure for the random effects to be fitted in the mmer solver.

Usage

  vs(..., Gu=NULL, Gti=NULL, Gtc=NULL, reorderGu=TRUE, buildGu=TRUE)

Arguments

Details

When providing initial values in the Gti argument the user has to provide scaled variance component values. The user can provide values from a previous model by accessing the sigma_scaled output from an mmer model or if an specific value is desired the user can obtain the scaled value as:

m = x/var(y)

where x is the desired initial value and y is the response variable. You can find an example in the DT_cpdata dataset.

Value

$res

a list with all neccesary elements (incidence matrices, known var-cov structures, unknown covariance structures to be estimated and constraints) to be used in the mmer solver.

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

See Also

The core function of the package: mmer

Examples

data(DT_example, package="enhancer")
DT <- DT_example
A <- A_example

## ============================ ##
## example to without structure
## ============================ ##

mix <- mmer(Yield~Env, 
            random= ~ vs(Name),
            rcov=~ vs(units),
            data=DT)

Variance-structure model with arbitrary Kronecker covariance factors

Description

vsm constructs random-effect or residual covariance structures for mmes. Covariance-shaping terms are represented by CovarianceFactor v2 descriptors and may be combined in an arbitrary-depth Kronecker product. A single product-level variance, \sigma^2, is owned by vsm, which avoids scale confounding among the individual covariance factors.

Usage

vsm(..., Gu = NULL, sigma2 = NULL, fixedSigma2 = FALSE,
  rotation = FALSE, isFixed = FALSE, verbose = TRUE)

Arguments

Details

Let K_1,\ldots,K_m be the dimensionless covariance shapes supplied by the covariance constructors. The covariance represented by one vsm term is

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

For a random effect with known relationship precision Gu, the factor product determines the covariance among the crossed covariance coordinates and Gu determines the relationship structure among the levels of the final main effect.

The final argument in ... is the main-effect incidence term. Every term before it is a covariance-shaping factor. Thus vsm(dsm(location), ar1m(row), ism(genotype), Gu=Ainv) represents a location-by-row covariance shape crossed with genotype. Kronecker ordering is left to right: earlier factors are the slow index and later factors are the fast index.

All optimizer coordinates are stored on unconstrained working scales. Each CovarianceFactor descriptor also carries its evaluator, derivative strategy, reporting transformation, and trust-region caps. Built-in high-use structures may use native C++ covariance primitives, while newer or user-defined structures can use generic R callbacks without requiring changes to ai_mme_sp2.

For residual covariance models, the covariance-shaping factors must define exactly one product coordinate for each observation. mmes combines that local coordinate with an independently constructed residual block index, so a large observation-level identity design is not materialized.

When rotation=TRUE, Gu is decomposed as

Gu=U\Lambda U',

where \Lambda is diagonal. The aligned original precision, diagonal precision, and eigenvectors are retained in the returned descriptor. The mmes front end validates the observation layout and applies the engine-specific rotation after observation filtering. Residual covariance must satisfy R_{(i,a),(l,b)} = c_{ab}\delta_{il}, where i,l index relationship levels and a,b index rotation blocks, so that U'RU = R; otherwise mmes stops with an error.

Value

A list containing the random-effect design blocks in Z, the precision matrix in Gu, the flattened Kronecker descriptor in covStruct, the product design, and, for residual models, the local covariance-coordinate index. When rotation is requested, GuRot contains the diagonal precision and rotation contains its eigen metadata. The covStruct object contains descriptor version 2 CovarianceFactor objects.

See Also

mmes, ism, dsm, atm, usm, csm, ar1m, ar2m, ar3m, mam, corgm, fam, antem, rrm, maternm, toeplitzm, sar, car, and ownm. See dsumm for section-specific (direct-sum) residual structures.

Examples

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

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


str(with(DT, vsm(dsm(Env),ism(Name))))
str(with(DT, vsm(csm(Env),ism(Name))))
str(with(DT, vsm(toeplitzm(Env),ism(Name))))

## Gu in 3-column (row, column, value) format with level names in rowNames
Ainv <- data.frame(Row=c(1,2,2,3), Column=c(1,1,2,3),
                   Ainverse=c(2,-1,2,1))
attr(Ainv, "rowNames") <- c("a","b","c")
id <- factor(c("a","b","c","a"))
str(vsm(ism(id), Gu=Ainv)$Gu)

variance structure specification

Description

vsr DEPRECATED NOW. was the main function to build the variance-covariance structure for the random effects to be fitted in the mmer solver.

Usage

  vsr(..., Gu=NULL, Gti=NULL, Gtc=NULL, reorderGu=TRUE, buildGu=TRUE)

Arguments

Details

When providing initial values in the Gti argument the user has to provide scaled variance component values. The user can provide values from a previous model by accessing the sigma_scaled output from an mmer model or if an specific value is desired the user can obtain the scaled value as:

m = x/var(y)

where x is the desired initial value and y is the response variable. You can find an example in the DT_cpdata dataset.

Value

$res

a list with all neccesary elements (incidence matrices, known var-cov structures, unknown covariance structures to be estimated and constraints) to be used in the mmer solver.

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

Examples

data(DT_example, package="enhancer")
DT <- DT_example
A <- A_example

## ============================ ##
## example to without structure
## ============================ ##

mix <- mmer(Yield~Env, 
            random= ~ vsr(Name),
            rcov=~ vsr(units),
            data=DT)

Wald tests for fixed effects with optional small-sample denominator df

Description

Pseudo analysis of variance for the fixed terms of an mmes fit, based on incremental or conditional Wald statistics, with chi-square, residual, Satterthwaite or Kenward-Roger reference distributions.

Usage

wald_mmes(object, ssType = c("incremental", "conditional"),
          denDF = c("none", "residual", "satterthwaite", "kr"), terms = NULL)

Arguments

Details

Let \Phi be the covariance of the fixed effects and M=\Phi^{-1}. For a given order of the terms, the upper Cholesky factor U of M defines z=U b; the incremental Wald statistic of a term is the sum of the squared elements of z belonging to it, i.e. a test of L b=0 with L the corresponding rows of U. Conditional tests use the order [terms not containing the tested term, tested term, terms containing it].

For Satterthwaite and Kenward-Roger df the model inputs are rebuilt from the stored call and X^\top V^{-1}X is re-evaluated at perturbed covariance parameters (on the reported scale of covPar); its numerical first and second derivatives give the quantities P_i and the Kenward-Roger adjustment. The covariance of the covariance parameters is theta_se, the inverse average information. The second-derivative term of V enters through the Hessian of X^\top V^{-1}X; it vanishes for variance-component models, where the method is the exact Kenward-Roger approximation. Because the average (rather than expected) information is used, df can differ slightly from software using the expected information. The re-evaluation requires the data and objects used in the original call to be available, REML estimation, and a Gaussian fit.

anova(object) with a single model returns wald_mmes(object, ...).

Value

A data frame of class "wald.mmes" with columns Df, denDF, Wald, F.value and p.value, one row per term.

References

Kenward MG and Roger JH (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics 53: 983-997.

See Also

mmes, predict.mmes, anova.mmes

Examples

data(DT_yatesoats, package = "enhancer")
m <- mmes(Y ~ V * N, random = ~ B + B:MP, rcov = ~ units,
          data = DT_yatesoats, verbose = FALSE)
wald_mmes(m)
wald_mmes(m, ssType = "conditional", denDF = "kr")