Package 'pegs'

Pseudo-Expectation Gauss-Seidel

Description

A lightweight, dependency-free, and simplified implementation of the Pseudo-Expectation Gauss-Seidel (PEGS) algorithm. It fits the multivariate ridge regression model for genomic prediction Xavier and Habier (2022) <doi:10.1186/s12711-022-00730-w> and Xavier et al. (2025) <doi:10.1093/genetics/iyae179>, providing heritability estimates, genetic correlations, breeding values, and regression coefficient estimates for prediction. This package provides an alternative to the 'bWGR' package by Xavier et al. (2019) <doi:10.1093/bioinformatics/btz794> by using 'LAPACK' for its algebraic operations.

Author(s)

Alencar Xavier

Examples

# Load simulated data
data(test);

# Fit multivariate ridge regression
fit = pegs(Y,X);

# Heritability
print(fit$h2);

# Genetic correlations
print(fit$GC);

# Accuracy
print(mean(diag(cor(fit$hat,simu$tbv))));

---

Pseudo-Expectation Gauss-Seidel for Multi-Trait Models

Description

This function implements the Pseudo-Expectation Gauss-Seidel (PEGS) algorithm for multi-trait genomic prediction models. It supports fitting one or more random effects and can apply factor-analytic approximations to the genetic covariance structure.

Usage

pegs(Y, X, maxit = 100, logtol = -4, NNC = FALSE, 
     covbend = 1.1, covMinEv = 10e-4, XFA = -1)

Arguments

Y	A numeric matrix (n x k) of dependent variables, where n is the number of observations and k is the number of response variables. Missing values should be NA.
X	A numeric matrix (n x p) of predictors (i.e., markers) for a single random effect, or a named list of such matrices for multiple random effects.
maxit	An integer specifying the maximum number of iterations.
logtol	A numeric value for the convergence tolerance on a log10 scale. The algorithm stops when the log10 of the sum of squared changes in coefficients is less than this value.
NNC	A logical value. If TRUE, imposes a constraint of non-negative correlations on the genetic covariance matrix.
covbend	A numeric bending factor used to inflate the diagonal of the covariance matrix if it is not positive semi-definite. The bending amount is calculated as abs(min_eigenvalue) * covbend.
covMinEv	The minimum allowable eigenvalue for the genetic covariance matrix. If the smallest eigenvalue falls below this threshold, the matrix is "bent" (diagonally inflated) to make it positive semi-definite.
XFA	An integer specifying the number of factors (principal components) to use in the factor-analytic approximation of the genetic covariance matrix. XFA = -1 (Default): Fits a full-rank model (no approximation). XFA = 0: Fits a diagonal model, assuming independent genetic effects among traits. XFA > 0: Fits a reduced-rank model using the specified number of factors.

Value

A list containing the following components:

mu

Vector of intercepts for each response variable.

b

A matrix of regression coefficients (p x k) if a single X matrix is provided, or a named list of such matrices if a list is provided for X.

hat

Matrix of predicted values (n x k).

h2

Vector of heritability estimates for each trait.

GC

The genetic correlation matrix (k x k) if a single X matrix is provided, or a named list of such matrices for each random effect.

bend

The inflation factor added to the diagonal of the genetic variance matrix. This will be a single value if one random effect is fit, or a named vector if multiple effects are fit.

numit

The number of iterations until convergence.

cnv

The final convergence value.

Examples

# Load simulated data with one random effect
data(test)

# Fit a standard multi-trait model
fit <- pegs(Y, X)

# Heritability
print(fit$h2)

# Genetic correlations
print(fit$GC)

# Accuracy
print(mean(diag(cor(fit$hat, simu$tbv))))

---