Package 'selection.index'

Description

Implements the Base Index where coefficients are set equal to economic weights. This is a simple, non-optimized approach that serves as a baseline comparison.

Unlike the Smith-Hazel index which requires matrix inversion, the Base Index is computationally trivial and robust when covariance estimates are unreliable.

Usage

base_index(
  pmat,
  gmat,
  wmat,
  wcol = 1,
  selection_intensity = 2.063,
  compare_to_lpsi = TRUE,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation:

Index coefficients: $b = w$

The Base Index is appropriate when: - Covariance estimates are unreliable - Computational simplicity is required - A baseline for comparison is needed

Value

List with:

• summary - Data frame with coefficients and metrics

• b - Vector of Base Index coefficients (equal to w)

• w - Named vector of economic weights

• Delta_G - Named vector of expected genetic gains per trait

• lpsi_comparison - Optional comparison with Smith-Hazel LPSI

Examples

## Not run: 
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
weights <- c(10, 8, 6, 4, 2, 1, 1)

result <- base_index(pmat, gmat, weights, compare_to_lpsi = TRUE)
print(result)

## End(Not run)

Description

Implements the Combined Linear Genomic Selection Index where selection combines both phenotypic observations and Genomic Estimated Breeding Values (GEBVs). This is used for selecting candidates with both phenotype and genotype data (e.g., in a training population).

Usage

clgsi(
  phen_mat = NULL,
  gebv_mat = NULL,
  pmat,
  gmat,
  wmat,
  wcol = 1,
  P_y = NULL,
  P_g = NULL,
  P_yg = NULL,
  reliability = NULL,
  selection_intensity = 2.063,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation:

The CLGSI combines phenotypic and genomic information:

$I = \mathbf{b}_y' \mathbf{y} + \mathbf{b}_g' \hat{\mathbf{g}}$

Coefficients are obtained by solving the partitioned system:

$\begin{bmatrix} \mathbf{b}_y \\ \mathbf{b}_g \end{bmatrix} = \begin{bmatrix} \mathbf{P} & \mathbf{P}_{y\hat{g}} \\ \mathbf{P}_{y\hat{g}}' & \mathbf{P}_{\hat{g}} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{G} \\ \mathbf{C}_{\hat{g}g} \end{bmatrix} \mathbf{w}$

Where: - $\mathbf{P}$ = Var(phenotypes) - $\mathbf{P}_{\hat{g}}$ = Var(GEBVs) - $\mathbf{P}_{y\hat{g}}$ = Cov(phenotypes, GEBVs) - $\mathbf{G}$ = Genotypic variance-covariance - $\mathbf{C}_{\hat{g}g}$ = Cov(GEBV, true BV)

Value

List with components:

• b_y - Coefficients for phenotypes

• b_g - Coefficients for GEBVs

• b_combined - Full coefficient vector [b_y; b_g]

• P_combined - Combined variance matrix

• Delta_H - Expected genetic advance per trait

• GA - Overall genetic advance

• PRE - Percent relative efficiency

• hI2 - Index heritability

• rHI - Index accuracy

• summary - Data frame with all metrics

References

Cerón-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing.

Examples

## Not run: 
# Generate example data
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate phenotypes and GEBVs
set.seed(123)
n_genotypes <- 100
n_traits <- ncol(gmat)

phen_mat <- matrix(rnorm(n_genotypes * n_traits, mean = 15, sd = 3),
  nrow = n_genotypes, ncol = n_traits
)
gebv_mat <- matrix(rnorm(n_genotypes * n_traits, mean = 10, sd = 2),
  nrow = n_genotypes, ncol = n_traits
)
colnames(phen_mat) <- colnames(gmat)
colnames(gebv_mat) <- colnames(gmat)

# Economic weights
weights <- c(10, 5, 3, 3, 5, 8, 4)

# Calculate CLGSI
result <- clgsi(phen_mat, gebv_mat, pmat, gmat, weights, reliability = 0.7)
print(result$summary)

## End(Not run)

Description

Generic mathematical operations optimized with C++/Eigen. No design-specific logic - purely mathematical primitives that can be orchestrated by R code to implement any experimental design.

This architecture allows: - Easy addition of new experimental designs (in R only) - C++ speed for heavy computation - Single source of truth (design_stats.R) - Better maintainability and testability

Description

Implements the CPPG-LGSI which combines phenotypic and genomic information while achieving predetermined proportional gains between traits. This is the most general constrained genomic index.

Usage

cppg_lgsi(
  T_C = NULL,
  Psi_C = NULL,
  d,
  phen_mat = NULL,
  gebv_mat = NULL,
  pmat = NULL,
  gmat = NULL,
  wmat = NULL,
  wcol = 1,
  U = NULL,
  reliability = NULL,
  k_I = 2.063,
  L_I = 1,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation (Chapter 6, Section 6.4):

Coefficient vector: $beta_CP = beta_CR + theta_CP * delta_CP$

Where beta_CR is from CRLGSI and:

$theta_CP = (beta_C' * Phi_C * (Phi_C' * T_C^{-1} * Phi_C)^{-1} * d_C) / (d_C' * (Phi_C' * T_C^{-1} * Phi_C)^{-1} * d_C)$

Selection response: $R_CP = (k_I / L_I) * sqrt(beta_CP' * T_C * beta_CP)$

Value

List with:

• summary - Data frame with coefficients and metrics

• b - Vector of CPPG-LGSI coefficients ($\beta_{CP}$)

• b_y - Coefficients for phenotypes

• b_g - Coefficients for GEBVs

• E - Expected genetic gains per trait

• theta_CP - Proportionality constant

• gain_ratios - Ratios of achieved to desired gains

Examples

## Not run: 
# Simulate data
set.seed(123)
n_genotypes <- 100
n_traits <- 5

phen_mat <- matrix(rnorm(n_genotypes * n_traits, 15, 3), n_genotypes, n_traits)
gebv_mat <- matrix(rnorm(n_genotypes * n_traits, 10, 2), n_genotypes, n_traits)

gmat <- cov(phen_mat) * 0.6
pmat <- cov(phen_mat)

# Desired proportional gains
d <- c(2, 1, 1, 0.5, 0)

w <- c(10, 8, 6, 4, 2)

result <- cppg_lgsi(
  phen_mat = phen_mat, gebv_mat = gebv_mat,
  pmat = pmat, gmat = gmat, d = d, wmat = w,
  reliability = 0.7
)
print(result$summary)
print(result$gain_ratios)

## End(Not run)

Description

Implements the CRLGSI which combines phenotypic and genomic information with restrictions on genetic gains. This extends CLGSI to include constraints.

Usage

crlgsi(
  T_C = NULL,
  Psi_C = NULL,
  phen_mat = NULL,
  gebv_mat = NULL,
  pmat = NULL,
  gmat = NULL,
  wmat,
  wcol = 1,
  restricted_traits = NULL,
  U = NULL,
  reliability = NULL,
  k_I = 2.063,
  L_I = 1,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation (Chapter 6, Section 6.3):

The CRLGSI combines phenotypic and genomic data with restrictions.

Coefficient vector: $beta_CR = K_C * beta_C$

Where K_C incorporates the restriction matrix.

Selection response: $R_CR = (k_I / L_I) * sqrt(beta_CR' * T_C * beta_CR)$

Expected gains: $E_CR = (k_I / L_I) * (Psi_C * beta_CR) / sqrt(beta_CR' * T_C * beta_CR)$

Value

List with:

• summary - Data frame with coefficients and metrics

• b - Vector of CRLGSI coefficients ($\beta_{CR}$)

• b_y - Coefficients for phenotypes

• b_g - Coefficients for GEBVs

• E - Expected genetic gains per trait

• R - Overall selection response

Examples

## Not run: 
# Simulate data
set.seed(123)
n_genotypes <- 100
n_traits <- 5

phen_mat <- matrix(rnorm(n_genotypes * n_traits, 15, 3), n_genotypes, n_traits)
gebv_mat <- matrix(rnorm(n_genotypes * n_traits, 10, 2), n_genotypes, n_traits)

gmat <- cov(phen_mat) * 0.6 # Genotypic component
pmat <- cov(phen_mat)

w <- c(10, 8, 6, 4, 2)

# Restrict traits 2 and 4
result <- crlgsi(
  phen_mat = phen_mat, gebv_mat = gebv_mat,
  pmat = pmat, gmat = gmat, wmat = w,
  restricted_traits = c(2, 4), reliability = 0.7
)
print(result$summary)

## End(Not run)

Description

Implements the Pesek & Baker (1969) Desired Gains Index where breeders specify target genetic gains instead of economic weights. This enhanced version includes calculation of implied economic weights and feasibility checking.

Usage

dg_lpsi(
  pmat,
  gmat,
  d,
  return_implied_weights = TRUE,
  check_feasibility = TRUE,
  selection_intensity = 2.063
)

Arguments

Details

Mathematical Formulation:

1. Index coefficients: $\mathbf{b} = \mathbf{G}^{-1}\mathbf{d}$

2. Expected response: $\Delta \mathbf{G} = (i/\sigma_I) \mathbf{G}\mathbf{b}$

CRITICAL: Scale Invariance Property

The achieved gains $\Delta\mathbf{G}$ are determined by selection intensity (i), genetic variance (G), and phenotypic variance (P), NOT by scaling $\mathbf{b}$. If you multiply $\mathbf{b}$ by constant c, $\sigma_I$ also scales by c, causing complete cancellation in $\Delta\mathbf{G} = (i/(c\sigma_I))\mathbf{G}(c\mathbf{b}) = (i/\sigma_I)\mathbf{G}\mathbf{b}$.

What DG-LPSI Actually Achieves:

- Proportional gains matching the RATIOS in d (not absolute magnitudes) - Achieved magnitude depends on biological/genetic constraints - Use feasibility checking to verify if desired gains are realistic

3. Implied economic weights (Section 1.4 of Chapter 4):

$\hat{\mathbf{w}} = \mathbf{G}^{-1} \mathbf{P} \mathbf{b}$

The implied weights represent the economic values that would have been needed in a Smith-Hazel index to achieve the desired gain PROPORTIONS. Large implied weights indicate traits that are "expensive" to improve (low heritability or unfavorable correlations), while small weights indicate traits that are "cheap" to improve.

Feasibility Checking:

The function estimates maximum possible gains as approximately 3.0 * sqrt(G_ii) (assuming very intense selection with i ~ 3.0) and warns if desired gains exceed 80% of these theoretical maxima.

Value

List with:

• summary - Data frame with coefficients, metrics, and implied weights

• b - Vector of selection index coefficients

• Delta_G - Named vector of achieved genetic gains per trait

• desired_gains - Named vector of desired gains (input d)

• gain_errors - Difference between desired and achieved gains

• implied_weights - Economic weights that would achieve these gains in Smith-Hazel LPSI

• implied_weights_normalized - Normalized implied weights (max absolute = 1)

• feasibility - Data frame with feasibility analysis per trait

• hI2 - Index heritability

• rHI - Index accuracy

References

Pesek, J., & Baker, R. J. (1969). Desired improvement in relation to selection indices. Canadian Journal of Plant Science, 49(6), 803-804.

Examples

# Load data
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Specify desired gains (e.g., increase each trait by 1 unit)
desired_gains <- rep(1, ncol(gmat))

# Calculate Desired Gains Index with all enhancements
result <- dg_lpsi(pmat, gmat, d = desired_gains)

# View summary
print(result$summary)

# Extract implied weights to understand relative "cost" of gains
print(result$implied_weights_normalized)

# Check feasibility
print(result$feasibility)

Description

Implements the ESIM by maximising the squared accuracy $\rho_{HI}^2$ through the generalised eigenproblem of the multi-trait heritability matrix $\mathbf{P}^{-1}\mathbf{C}$.

Unlike the Smith-Hazel LPSI, **no economic weights are required**. The net genetic merit vector $\mathbf{w}_E$ is instead implied by the solution.

Usage

esim(pmat, gmat, selection_intensity = 2.063, n_indices = 1L)

Arguments

Details

Eigenproblem (Section 7.1):

$(\mathbf{P}^{-1}\mathbf{C} - \lambda_E^2 \mathbf{I})\mathbf{b}_E = 0$

The solution $\lambda_E^2$ (largest eigenvalue) equals the maximum achievable index heritability $h^2_{I_E}$.

Key metrics:

$R_E = k_I \sqrt{\mathbf{b}_E^{\prime}\mathbf{P}\mathbf{b}_E}$

$\mathbf{E}_E = k_I \frac{\mathbf{C}\mathbf{b}_E}{\sqrt{\mathbf{b}_E^{\prime}\mathbf{P}\mathbf{b}_E}}$

Implied economic weights:

$\mathbf{w}_E = \frac{\sqrt{\lambda_E^2}}{\mathbf{b}_E^{\prime}\mathbf{P}\mathbf{b}_E} \mathbf{C}^{-1}\mathbf{P}\mathbf{b}_E$

Uses cpp_symmetric_solve and cpp_quadratic_form_sym from math_primitives.cpp for efficient matrix operations, and R's eigen() for the eigendecomposition.

Value

Object of class "esim", a list with:

summary

Data frame with b coefficients, hI2, rHI, sigma_I, Delta_G, and lambda2 for each index requested.

b

Named numeric vector of optimal ESIM coefficients (1st index).

Delta_G

Named numeric vector of expected genetic gains per trait.

sigma_I

Standard deviation of the index $\sigma_I$.

hI2

Index heritability $h^2_{I_E}$ (= leading eigenvalue).

rHI

Accuracy $r_{HI_E}$.

lambda2

Leading eigenvalue (maximised index heritability).

implied_w

Implied economic weights $\mathbf{w}_E$.

all_eigenvalues

All eigenvalues of $\mathbf{P}^{-1}\mathbf{C}$.

selection_intensity

Selection intensity used.

References

Ceron-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Section 7.1.

Examples

## Not run: 
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

result <- esim(pmat, gmat)
print(result)
summary(result)

## End(Not run)

Description

Estimates and imputes missing values in randomized complete block design (RCBD), Latin square design (LSD), or split plot design (SPD) experimental data.

Uses one of six methods: REML, Yates, Healy, Regression, Mean, or Bartlett.

Usage

estimate_missing_values(
  data,
  genotypes,
  replications,
  columns = NULL,
  main_plots = NULL,
  design = c("RCBD", "LSD", "SPD"),
  method = c("REML", "Yates", "Healy", "Regression", "Mean", "Bartlett"),
  tolerance = 1e-06
)

Arguments

Details

The function handles missing values by iteratively estimating them based on the experimental design structure:

**RCBD:** 2-way blocking (genotypes × blocks) **LSD:** 3-way blocking (genotypes × rows × columns) **SPD:** Nested structure (blocks > main plots > sub-plots)

Method Availability:

• RCBD: All methods (REML, Yates, Healy, Regression, Mean, Bartlett)

• LSD: All methods (REML, Yates, Healy, Regression, Mean, Bartlett)

• SPD: Mean only (other methods fall back to Mean)

Method Selection Guide:

• Use REML for complex missing patterns or when precision is critical

• Use Yates for balanced designs with few missing values

• Use Healy when multiple values missing from same treatment/block

• Use Regression for fast, deterministic estimation

• Use Mean for quick estimation when precision is less critical

• Use Bartlett when traits are highly correlated

The function uses the centralized design_stats engine for all ANOVA computations, ensuring consistency with gen_varcov(), phen_varcov(), and mean_performance().

Value

Matrix of the same dimensions as data with all missing values replaced by estimates.

References

Yates, F. (1933). The analysis of replicated experiments when the field results are incomplete. Empire Journal of Experimental Agriculture, 1, 129-142.

Healy, M. J. R., & Westmacott, M. (1956). Missing values in experiments analysed on automatic computers. Applied Statistics, 5(3), 203-206.

Bartlett, M. S. (1937). Some examples of statistical methods of research in agriculture and applied biology. Supplement to the Journal of the Royal Statistical Society, 4(2), 137-183.

Examples

# RCBD example with missing values
data(seldata)
test_data <- seldata[, 3:5]
test_data[c(1, 10, 25), 1] <- NA
test_data[c(5, 15), 2] <- NA

# Impute using Yates method
imputed <- estimate_missing_values(test_data, seldata$treat, seldata$rep, method = "Yates")

# Check that no NA remain
anyNA(imputed) # Should be FALSE

## Not run: 
# Latin Square Design example
# lsd_data should have genotypes, rows, and columns
imputed_lsd <- estimate_missing_values(
  data = lsd_data[, 3:7],
  genotypes = lsd_data$treat,
  replications = lsd_data$row,
  columns = lsd_data$col,
  design = "LSD",
  method = "REML"
)

# Split Plot Design example
# spd_data should have sub-plots, blocks, and main plots
imputed_spd <- estimate_missing_values(
  data = spd_data[, 3:7],
  genotypes = spd_data$subplot,
  replications = spd_data$block,
  main_plots = spd_data$mainplot,
  design = "SPD",
  method = "Mean"
)

## End(Not run)

Description

Genetic Advance for PRE

Usage

gen_advance(phen_mat, gen_mat, weight_mat)

Arguments

Value

Genetic advance of character or character combinations

Examples

gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
gen_advance(phen_mat = pmat[1, 1], gen_mat = gmat[1, 1], weight_mat = weight[1, 2])

Description

Genotypic Variance-Covariance Analysis

Usage

gen_varcov(
  data,
  genotypes,
  replication,
  columns = NULL,
  main_plots = NULL,
  design_type = c("RCBD", "LSD", "SPD"),
  method = c("REML", "Yates", "Healy", "Regression", "Mean", "Bartlett")
)

Arguments

Value

A Genotypic Variance-Covariance Matrix

Examples

# RCBD example
gen_varcov(data = seldata[, 3:9], genotypes = seldata$treat, replication = seldata$rep)

# Latin Square Design example (requires columns parameter)
# gen_varcov(data=lsd_data[,3:7], genotypes=lsd_data$treat,
#           replication=lsd_data$row, columns=lsd_data$col, design_type="LSD")

# Split Plot Design example (requires main_plots parameter)
# gen_varcov(data=spd_data[,3:7], genotypes=spd_data$subplot,
#           replication=spd_data$block, main_plots=spd_data$mainplot, design_type="SPD")

Description

Computes the genetic-genomic covariance matrix (A) as defined in Chapter 8 (Equation 8.12) for GESIM and related genomic eigen selection indices.

Structure: A = [[C, C_g-gamma], [C_gamma-g, $\Gamma$]] (2t x 2t, square symmetric)

where: - C = Var(g) = true genotypic variance-covariance (t x t) - $\Gamma$ = Var($\gamma$) = genomic variance-covariance (t x t) - C_g-gamma = Cov(g, $\gamma$) = covariance between true BVs and GEBVs (t x t) - C_gamma-g = Cov($\gamma$, g) = transpose of C_g-gamma (t x t)

Usage

genetic_genomic_varcov(
  gmat,
  Gamma = NULL,
  reliability = NULL,
  C_gebv_g = NULL,
  square = TRUE
)

Arguments

Details

The genetic-genomic matrix relates selection on phenotypes + GEBVs to expected genetic gains.

**For GESIM (Chapter 8):** Requires the full (2t × 2t) square matrix for the eigenproblem: ($\Phi$^(-1) A - $\lambda$I)b = 0

**For LMSI/CLGSI (Chapter 4):** Can use the rectangular (2t × t) form in the equation: b = P^(-1) G w, where G is (2t × t).

When reliability is provided: - $C_{\gamma g}$ = diag($\sqrt{r^2}$)

When reliability is NULL: - $C_{\gamma g}$ = Gamma (assumes unbiased GEBVs, perfect prediction)

Value

Genetic-genomic covariance matrix: - If square = TRUE: (2t × 2t) symmetric matrix for GESIM/eigen indices - If square = FALSE: (2t × t) rectangular matrix for LMSI where t is the number of traits

References

Cerón-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Chapters 4 & 8.

Examples

## Not run: 
# Generate example data
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate genomic covariance
Gamma <- gmat * 0.8

# For GESIM: Get square (2t × 2t) matrix
A_square <- genetic_genomic_varcov(gmat, Gamma, reliability = 0.7)
print(dim(A_square)) # Should be 14 x 14 (2t × 2t)

# For LMSI: Get rectangular (2t × t) matrix
A_rect <- genetic_genomic_varcov(gmat, Gamma, reliability = 0.7, square = FALSE)
print(dim(A_rect)) # Should be 14 x 7 (2t × t)

## End(Not run)

Description

Computes genomic variance-covariance matrix ($\Gamma$ or Gamma) from a matrix of Genomic Estimated Breeding Values (GEBVs).

$\gamma$ (gamma) represents GEBV vectors obtained from genomic prediction models (e.g., GBLUP, rrBLUP, Genomic BLUP). This function computes Var($\gamma$) = $\Gamma$.

Usage

genomic_varcov(gebv_mat, method = "pearson", use = "complete.obs")

Arguments

Details

The genomic variance-covariance matrix $\Gamma$ captures genetic variation as predicted by molecular markers. It is computed as:

where $\gamma_i$ is the GEBV vector for genotype i and $\mu_{\gamma}$ is the mean GEBV vector.

**Missing Value Handling:** - "complete.obs": Uses only complete observations (recommended) - "pairwise.complete.obs": Uses pairwise-complete observations (may not be PSD) - "everything": Fails if any NA present

When using pairwise deletion, the resulting matrix may not be positive semi-definite (PSD), which can cause numerical issues in selection indices.

**Applications:** In selection index theory: - Used in LGSI (Linear Genomic Selection Index) - Component of $\Phi$ (phenomic-genomic covariance) - Component of A (genetic-genomic covariance)

Value

Symmetric genomic variance-covariance matrix (n_traits x n_traits)

References

Cerón-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Chapters 4 & 8.

Examples

## Not run: 
# Simulate GEBVs
set.seed(123)
n_genotypes <- 100
n_traits <- 5
gebv_mat <- matrix(rnorm(n_genotypes * n_traits),
  nrow = n_genotypes, ncol = n_traits
)
colnames(gebv_mat) <- paste0("Trait", 1:n_traits)

# Compute genomic variance-covariance
Gamma <- genomic_varcov(gebv_mat)
print(Gamma)

## End(Not run)

Description

Implements the GESIM by maximising the squared accuracy through the generalised eigenproblem combining phenotypic data with GEBVs (Genomic Estimated Breeding Values). No economic weights are required.

Usage

gesim(pmat, gmat, Gamma, selection_intensity = 2.063)

Arguments

Details

Eigenproblem (Section 8.2):

$(\mathbf{\Phi}^{-1}\mathbf{A} - \lambda_G^2 \mathbf{I}_{2t})\boldsymbol{\beta}_G = 0$

where:

$\mathbf{\Phi} = \begin{bmatrix} \mathbf{P} & \mathbf{\Gamma} \\ \mathbf{\Gamma} & \mathbf{\Gamma} \end{bmatrix}$

$\mathbf{A} = \begin{bmatrix} \mathbf{C} & \mathbf{\Gamma} \\ \mathbf{\Gamma} & \mathbf{\Gamma} \end{bmatrix}$

Implied economic weights:

$\mathbf{w}_G = \mathbf{A}^{-1}\mathbf{\Phi}\boldsymbol{\beta}$

Selection response:

$R_G = k_I \sqrt{\boldsymbol{\beta}_G^{\prime}\mathbf{\Phi}\boldsymbol{\beta}_G}$

Expected genetic gain per trait:

$\mathbf{E}_G = k_I \frac{\mathbf{A}\boldsymbol{\beta}_G}{\sqrt{\boldsymbol{\beta}_G^{\prime}\mathbf{\Phi}\boldsymbol{\beta}_G}}$

Value

Object of class "gesim", a list with:

summary

Data frame with coefficients and metrics.

b_y

Coefficients for phenotypic data.

b_gamma

Coefficients for GEBVs.

b_combined

Combined coefficient vector.

E_G

Expected genetic gains per trait.

sigma_I

Standard deviation of the index.

hI2

Index heritability (= leading eigenvalue).

rHI

Accuracy $r_{HI}$.

R_G

Selection response.

lambda2

Leading eigenvalue.

implied_w

Implied economic weights.

selection_intensity

Selection intensity used.

References

Ceron-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Section 8.2.

Examples

## Not run: 
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate GEBV covariance (in practice, compute from genomic predictions)
Gamma <- gmat * 0.8 # Assume 80% GEBV-phenotype covariance

result <- gesim(pmat, gmat, Gamma)
print(result)

## End(Not run)

Description

Implements the GW-ESIM by incorporating genome-wide marker effects directly into the eigen selection index framework. Uses N marker scores alongside phenotypic data.

Usage

gw_esim(pmat, gmat, G_M, M, selection_intensity = 2.063)

Arguments

Details

Eigenproblem (Section 8.3):

$(\mathbf{Q}^{-1}\mathbf{X} - \lambda_W^2 \mathbf{I}_{(t+N)})\boldsymbol{\beta}_W = 0$

where:

$\mathbf{Q} = \begin{bmatrix} \mathbf{P} & \mathbf{G}_M \\ \mathbf{G}_M^{\prime} & \mathbf{M} \end{bmatrix}$

$\mathbf{X} = \begin{bmatrix} \mathbf{C} & \mathbf{G}_M \\ \mathbf{G}_M^{\prime} & \mathbf{M} \end{bmatrix}$

Selection response:

$R_W = k_I \sqrt{\boldsymbol{\beta}_W^{\prime}\mathbf{Q}\boldsymbol{\beta}_W}$

Expected genetic gain per trait:

$\mathbf{E}_W = k_I \frac{\mathbf{X}\boldsymbol{\beta}_W}{\sqrt{\boldsymbol{\beta}_W^{\prime}\mathbf{Q}\boldsymbol{\beta}_W}}$

Value

Object of class "gw_esim", a list with:

summary

Data frame with key metrics.

b_y

Coefficients for phenotypic data.

b_m

Coefficients for marker scores.

b_combined

Combined coefficient vector.

E_W

Expected genetic gains per trait.

sigma_I

Standard deviation of the index.

hI2

Index heritability (= leading eigenvalue).

rHI

Accuracy.

R_W

Selection response.

lambda2

Leading eigenvalue.

selection_intensity

Selection intensity used.

References

Ceron-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Section 8.3.

Examples

## Not run: 
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate marker data
N_markers <- 100
n_traits <- nrow(gmat)
G_M <- matrix(rnorm(n_traits * N_markers, sd = 0.5), n_traits, N_markers)
M <- diag(N_markers) + matrix(rnorm(N_markers^2, sd = 0.1), N_markers, N_markers)
M <- (M + t(M)) / 2 # Make symmetric

result <- gw_esim(pmat, gmat, G_M, M)
print(result)

## End(Not run)

Description

Implements the GW-LMSI which uses all available genome-wide markers directly as predictors in the selection index. Unlike LMSI which uses aggregated marker scores per trait, GW-LMSI treats each individual marker as a separate predictor.

Usage

gw_lmsi(
  marker_mat,
  trait_mat = NULL,
  gmat,
  P_GW = NULL,
  G_GW = NULL,
  wmat,
  wcol = 1,
  lambda = 0,
  selection_intensity = 2.063,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation:

The GW-LMSI maximizes the correlation between the index $I_{GW} = \mathbf{b}_{GW}^{\prime}\mathbf{m}$ and the breeding objective $H = \mathbf{w}^{\prime}\mathbf{g}$.

Marker covariance matrix:

$\mathbf{P}_{GW} = Var(\mathbf{m})$

Covariance between markers and traits:

$\mathbf{G}_{GW} = Cov(\mathbf{m}, \mathbf{g})$

Index coefficients (with optional Ridge regularization):

$\mathbf{b}_{GW} = (\mathbf{P}_{GW} + \lambda \mathbf{I})^{-1} \mathbf{G}_{GW} \mathbf{w}$

Accuracy:

$\rho_{HI} = \sqrt{\frac{\mathbf{b}_{GW}^{\prime} \mathbf{G}_{GW} \mathbf{w}}{\mathbf{w}^{\prime} \mathbf{G} \mathbf{w}}}$

Selection response:

$R_{GW} = k \sqrt{\mathbf{b}_{GW}^{\prime} \mathbf{P}_{GW} \mathbf{b}_{GW}}$

Expected genetic gain per trait:

$\mathbf{E}_{GW} = k \frac{\mathbf{G}_{GW}^{\prime} \mathbf{b}_{GW}}{\sigma_{I_{GW}}}$

Note on Singularity Detection and Regularization: The function automatically detects problematic cases:

1. **High-dimensional case**: When n_markers > n_genotypes, P_GW is mathematically singular (rank-deficient). The function issues a warning and suggests an appropriate lambda value.

2. **Ill-conditioned case**: When P_GW has a high condition number (> 1e10), indicating numerical instability.

3. **Numerical singularity**: When P_GW has eigenvalues near zero.

Ridge regularization adds $\lambda$I to P_GW, ensuring positive definiteness. Recommended lambda values are 0.01-0.1 times the average diagonal element of P_GW. Users can also set lambda = 0 to force generalized inverse (less stable but sometimes needed).

Value

List of class "gw_lmsi" with components:

b

Index coefficients for markers (n_markers).

P_GW

Marker covariance matrix (n_markers x n_markers).

G_GW

Covariance between markers and traits (n_markers x n_traits).

rHI

Index accuracy (correlation between index and breeding objective).

sigma_I

Standard deviation of the index.

R

Selection response (k * sigma_I).

Delta_H

Expected genetic gain per trait (vector of length n_traits).

GA

Overall genetic advance in breeding objective.

PRE

Percent relative efficiency (if GAY provided).

hI2

Index heritability.

lambda

Ridge regularization parameter used.

n_markers

Number of markers.

high_dimensional

Logical indicating if n_markers > n_genotypes.

condition_number

Condition number of P_GW (if computed).

summary

Data frame with metrics.

References

Lande, R., & Thompson, R. (1990). Efficiency of marker-assisted selection in the improvement of quantitative traits. Genetics, 124(3), 743-756.

Cerón-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Chapter 4.

Examples

## Not run: 
# Load data
data(seldata)
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate marker data
set.seed(123)
n_genotypes <- 100
n_markers <- 200
n_traits <- ncol(gmat)

# Marker matrix (coded as 0, 1, 2)
marker_mat <- matrix(sample(0:2, n_genotypes * n_markers, replace = TRUE),
  nrow = n_genotypes, ncol = n_markers
)

# Trait matrix
trait_mat <- matrix(rnorm(n_genotypes * n_traits, mean = 15, sd = 3),
  nrow = n_genotypes, ncol = n_traits
)

# Economic weights
weights <- c(10, 5, 3, 3, 5, 8, 4)

# Calculate GW-LMSI with Ridge regularization
result <- gw_lmsi(marker_mat, trait_mat, gmat,
  wmat = weights, lambda = 0.01
)
print(result$summary)

## End(Not run)

Description

Converts genetic distance (in Morgans) to recombination fraction using Haldane's mapping function. This function models the relationship between genetic distance and the probability of recombination between loci.

Usage

haldane_mapping(distance)

Arguments

Details

Mathematical Formula (Chapter 10, Section 10.1):

The relationship between recombination fraction (r) and genetic distance (d):

$r = \frac{1}{2}(1 - e^{-2d})$

Where: - d = Genetic distance in Morgans - r = Recombination fraction (probability of recombination per meiosis)

This function assumes no crossover interference beyond that implied by the mapping function itself.

Value

Recombination fraction (r) ranging from 0 to 0.5. - r = 0 indicates complete linkage (no recombination) - r = 0.5 indicates independent assortment (unlinked loci)

Examples

# Zero distance means complete linkage (no recombination)
haldane_mapping(0) # Returns 0

# 1 Morgan distance
haldane_mapping(1) # Returns ~0.43

# Large distance approaches 0.5 (independent assortment)
haldane_mapping(10) # Returns ~0.5

# Vector of distances
distances <- c(0, 0.1, 0.5, 1.0, 2.0)
haldane_mapping(distances)

Description

Converts recombination fraction back to genetic distance (in Morgans). This is the inverse of Haldane's mapping function.

Usage

inverse_haldane_mapping(recombination_fraction)

Arguments

Details

Mathematical Formula:

Solving Haldane's equation for d:

$d = -\frac{1}{2} \ln(1 - 2r)$

Value

Genetic distance in Morgans.

Examples

# Convert recombination fraction to distance
inverse_haldane_mapping(0.25) # Returns ~0.347 Morgans
inverse_haldane_mapping(0.5) # Returns Inf (unlinked)

Description

Implements the Linear Genomic Selection Index where selection is based solely on Genomic Estimated Breeding Values (GEBVs). This is used for selecting candidates that have been genotyped but not phenotyped (e.g., in a testing population).

Usage

lgsi(
  gebv_mat,
  gmat,
  wmat,
  wcol = 1,
  reliability = NULL,
  selection_intensity = 2.063,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation:

The LGSI maximizes the correlation between the index I = b' * gebv and

Index coefficients: $\mathbf{b} = \mathbf{P}_{\hat{g}}^{-1} \mathbf{C}_{\hat{g}g} \mathbf{w}$

Where: - $\mathbf{P}_{\hat{g}}$ = Var(gebv) - variance-covariance of GEBVs - $\mathbf{C}_{\hat{g}g}$ = Cov(gebv, g) - covariance between GEBVs and true breeding values

If reliability (r) is known: $\mathbf{C}_{\hat{g}g} = \text{diag}(r) \mathbf{P}_{\hat{g}}$

Expected response: $\Delta \mathbf{H} = \frac{i}{\sigma_I} \mathbf{C}_{\hat{g}g} \mathbf{b}$

Value

List with components:

• b - Index coefficients

• P_gebv - GEBV variance-covariance matrix

• reliability - Reliability values used

• Delta_H - Expected genetic advance per trait

• GA - Overall genetic advance in the index

• PRE - Percent relative efficiency (if GAY provided)

• hI2 - Index heritability

• rHI - Index accuracy

• sigma_I - Standard deviation of the index

• summary - Data frame with all metrics

References

Cerón-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing.

Examples

## Not run: 
# Generate example data
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate GEBVs (in practice, these come from genomic prediction)
set.seed(123)
n_genotypes <- 100
n_traits <- ncol(gmat)
gebv_mat <- matrix(rnorm(n_genotypes * n_traits, mean = 10, sd = 2),
  nrow = n_genotypes, ncol = n_traits
)
colnames(gebv_mat) <- colnames(gmat)

# Economic weights
weights <- c(10, 5, 3, 3, 5, 8, 4)

# Calculate LGSI
result <- lgsi(gebv_mat, gmat, weights, reliability = 0.7)
print(result$summary)

## End(Not run)

Description

Implements the LMSI which combines phenotypic information with molecular marker scores from statistically significant markers (Lande & Thompson, 1990). The index is I = b_y' * y + b_s' * s, where y are phenotypes and s are marker scores.

Usage

lmsi(
  phen_mat = NULL,
  marker_scores = NULL,
  pmat,
  gmat,
  G_s = NULL,
  wmat,
  wcol = 1,
  selection_intensity = 2.063,
  GAY = NULL
)

Arguments

Details

Mathematical Formulation:

The LMSI maximizes the correlation between the index $I_{LMSI} = \mathbf{b}_y^{\prime}\mathbf{y} + \mathbf{b}_s^{\prime}\mathbf{s}$ and the breeding objective $H = \mathbf{w}^{\prime}\mathbf{g}$.

Combined covariance matrices:

$\mathbf{P}_L = \begin{bmatrix} \mathbf{P} & \text{Cov}(\mathbf{y}, \mathbf{s}) \\ \text{Cov}(\mathbf{y}, \mathbf{s})^{\prime} & \text{Var}(\mathbf{s}) \end{bmatrix}$

$\mathbf{G}_L = \begin{bmatrix} \mathbf{G} \\ \mathbf{G}_s \end{bmatrix}$

where $\mathbf{P}$ is the phenotypic variance, $\text{Cov}(\mathbf{y}, \mathbf{s})$ is the covariance between phenotypes and marker scores (computed from data), $\text{Var}(\mathbf{s})$ is the variance of marker scores, $\mathbf{G}$ is the genotypic variance, and $\mathbf{G}_s$ represents the genetic covariance explained by markers.

Index coefficients:

$\mathbf{b}_{LMSI} = \mathbf{P}_L^{-1} \mathbf{G}_L \mathbf{w}$

Accuracy:

$\rho_{HI} = \sqrt{\frac{\mathbf{b}_{LMSI}^{\prime} \mathbf{G}_L \mathbf{w}}{\mathbf{w}^{\prime} \mathbf{G} \mathbf{w}}}$

Selection response:

$R_{LMSI} = k \sigma_{I_{LMSI}} = k \sqrt{\mathbf{b}_{LMSI}^{\prime} \mathbf{P}_L \mathbf{b}_{LMSI}}$

Expected genetic gain per trait:

$\mathbf{E}_{LMSI} = k \frac{\mathbf{G}_L^{\prime} \mathbf{b}_{LMSI}}{\sigma_{I_{LMSI}}}$

Value

List of class "lmsi" with components:

b_y

Coefficients for phenotypes (n_traits vector).

b_s

Coefficients for marker scores (n_traits vector).

b_combined

Combined coefficient vector [b_y; b_s] (2*n_traits vector).

P_L

Combined phenotypic-marker covariance matrix (2*n_traits x 2*n_traits).

G_L

Combined genetic-marker covariance matrix (2*n_traits x n_traits).

G_s

Genomic covariance matrix explained by markers (n_traits x n_traits).

rHI

Index accuracy (correlation between index and breeding objective).

sigma_I

Standard deviation of the index.

R

Selection response (k * sigma_I).

Delta_H

Expected genetic gain per trait (vector of length n_traits).

GA

Overall genetic advance in breeding objective.

PRE

Percent relative efficiency (if GAY provided).

hI2

Index heritability.

summary

Data frame with coefficients and metrics (combined view).

phenotype_coeffs

Data frame with phenotype coefficients only.

marker_coeffs

Data frame with marker score coefficients only.

coeff_analysis

Data frame with coefficient distribution analysis.

References

Lande, R., & Thompson, R. (1990). Efficiency of marker-assisted selection in the improvement of quantitative traits. Genetics, 124(3), 743-756.

Cerón-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Chapter 4.

Examples

## Not run: 
# Load data
data(seldata)
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate marker scores (in practice, computed from QTL mapping)
set.seed(123)
n_genotypes <- 100
n_traits <- ncol(gmat)
marker_scores <- matrix(rnorm(n_genotypes * n_traits, mean = 5, sd = 1.5),
  nrow = n_genotypes, ncol = n_traits
)
colnames(marker_scores) <- colnames(gmat)

# Simulate phenotypes
phen_mat <- matrix(rnorm(n_genotypes * n_traits, mean = 15, sd = 3),
  nrow = n_genotypes, ncol = n_traits
)
colnames(phen_mat) <- colnames(gmat)

# Economic weights
weights <- c(10, 5, 3, 3, 5, 8, 4)

# Calculate LMSI
result <- lmsi(phen_mat, marker_scores, pmat, gmat,
  G_s = NULL, wmat = weights
)
print(result$summary)

## End(Not run)

Description

Build all possible Smith-Hazel selection indices from trait combinations, with optional exclusion of specific traits.

This function systematically evaluates indices for all combinations of ncomb traits, which is useful for identifying the most efficient subset of traits for selection.

Usage

lpsi(ncomb, pmat, gmat, wmat, wcol = 1, GAY, excluding_trait = NULL)

Arguments

Value

Data frame of all possible selection indices with metrics (GA, PRE, Delta_G, rHI, hI2)

Examples

## Not run: 
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
wmat <- weight_mat(weight)

# Build all 3-trait indices
result <- lpsi(ncomb = 3, pmat = pmat, gmat = gmat, wmat = wmat, wcol = 1)

# Exclude specific traits
result <- lpsi(
  ncomb = 3, pmat = pmat, gmat = gmat, wmat = wmat,
  excluding_trait = c(1, 3)
)

## End(Not run)

Description

A synthetic dataset containing Single Nucleotide Polymorphism (SNP) marker data for the 100 maize genotypes found in maize_pheno. Designed for testing genomic selection indices (GESIM/RGESIM) and relationship matrices.

Usage

data(maize_geno)

Format

A numeric matrix with 100 rows (genotypes) and 500 columns (SNP markers):

Rows

100 genotypes, matching the Genotype column in maize_pheno.

Columns

500 simulated SNP markers, coded as 0, 1, or 2.

Source

Simulated for the selection.index package to provide reproducible examples.

Examples

data(maize_geno)
dim(maize_geno)

Description

A synthetic dataset containing multi-environment phenotypic records for 100 maize genotypes. Designed to demonstrate phenotypic selection index calculations and variance-covariance matrix estimations.

Usage

data(maize_pheno)

Format

A data frame with 600 rows and 6 variables:

Genotype

Factor representing 100 unique maize genotypes.

Environment

Factor representing 2 distinct growing environments.

Block

Factor representing 3 replicate blocks within each environment.

Yield

Numeric vector of grain yield in kg/ha.

PlantHeight

Numeric vector of plant height in centimeters.

DaysToMaturity

Numeric vector of days to physiological maturity.

Source

Simulated for the selection.index package to provide reproducible examples.

Examples

data(maize_pheno)
head(maize_pheno)

Description

Mean performance of phenotypic data

Usage

mean_performance(
  data,
  genotypes,
  replications,
  columns = NULL,
  main_plots = NULL,
  design_type = c("RCBD", "LSD", "SPD"),
  method = c("REML", "Yates", "Healy", "Regression", "Mean", "Bartlett")
)

Arguments

Value

Dataframe of mean performance analysis

Examples

mean_performance(data = seldata[, 3:9], genotypes = seldata[, 2], replications = seldata[, 1])

Description

Implements the MESIM by maximising the squared accuracy through the generalised eigenproblem combining phenotypic data with molecular marker scores. Unlike Smith-Hazel LPSI, **no economic weights are required**.

Usage

mesim(pmat, gmat, S_M, S_Mg = NULL, S_var = NULL, selection_intensity = 2.063)

Arguments

Details

Eigenproblem (Section 8.1):

$(\mathbf{T}_M^{-1}\mathbf{\Psi}_M - \lambda_M^2 \mathbf{I}_{2t})\boldsymbol{\beta}_M = 0$

where:

$\mathbf{T}_M = \begin{bmatrix} \mathbf{P} & \mathrm{Cov}(\mathbf{y},\mathbf{s}) \\ \mathrm{Cov}(\mathbf{s},\mathbf{y}) & \mathrm{Var}(\mathbf{s}) \end{bmatrix}$

$\mathbf{\Psi}_M = \begin{bmatrix} \mathbf{C} & \mathrm{Cov}(\mathbf{g},\mathbf{s}) \\ \mathrm{Cov}(\mathbf{s},\mathbf{g}) & \mathrm{Var}(\mathbf{s}) \end{bmatrix}$

Theoretical distinction:

• T_M uses phenotypic covariances: Cov(y, s) provided via S_M

• Psi_M uses genetic covariances: Cov(g, s) provided via S_Mg

• Since y = g + e, if Cov(e, s) ~= 0, then Cov(y, s) ~= Cov(g, s)

• Chapter 8.1 assumes marker scores are pure genetic predictors, so S_M can be used for both (default behavior when S_Mg = NULL)

The solution $\lambda_M^2$ (largest eigenvalue) equals the maximum achievable index heritability.

Key metrics:

$R_M = k_I \sqrt{\boldsymbol{\beta}_{M}^{\prime}\mathbf{T}_M\boldsymbol{\beta}_{M}}$

$\mathbf{E}_M = k_I \frac{\mathbf{\Psi}_M\boldsymbol{\beta}_{M}}{\sqrt{\boldsymbol{\beta}_{M}^{\prime}\mathbf{T}_M\boldsymbol{\beta}_{M}}}$

Value

Object of class "mesim", a list with:

summary

Data frame with coefficients and metrics.

b_y

Coefficients for phenotypic data.

b_s

Coefficients for marker scores.

b_combined

Combined coefficient vector.

E_M

Expected genetic gains per trait.

sigma_I

Standard deviation of the index.

hI2

Index heritability (= leading eigenvalue).

rHI

Accuracy $r_{HI}$.

R_M

Selection response.

lambda2

Leading eigenvalue (maximised index heritability).

selection_intensity

Selection intensity used.

References

Ceron-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Section 8.1.

Examples

## Not run: 
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate marker score matrices (in practice, compute from data)
S_M <- gmat * 0.7 # Cov(y, s) - phenotype-marker covariance
S_Mg <- gmat * 0.65 # Cov(g, s) - genetic-marker covariance
S_var <- gmat * 0.8 # Var(s) - marker score variance

# Most rigorous: Provide all three covariance matrices
result <- mesim(pmat, gmat, S_M, S_Mg = S_Mg, S_var = S_var)
print(result)

# Standard usage: Cov(g,s) defaults to Cov(y,s) when errors uncorrelated
result_standard <- mesim(pmat, gmat, S_M, S_var = S_var)

# Backward compatible: Chapter 8.1 simplified notation
result_simple <- mesim(pmat, gmat, S_M)

## End(Not run)

Description

Exported wrapper for missing value estimation in RCBD, Latin Square, and Split Plot designs. Calls the centralized design_stats engine for ANOVA computations.

Description

Implements the two-stage Linear Genomic Selection Index where selection is based on GEBVs at both stages with covariance adjustments due to selection effects.

Usage

mlgsi(
  Gamma1,
  Gamma,
  A1,
  A,
  C,
  G1,
  P1,
  wmat,
  wcol = 1,
  selection_proportion = 0.1,
  use_young_method = FALSE,
  k1_manual = 2.063,
  k2_manual = 2.063,
  tau = NULL,
  reliability = NULL
)

Arguments

Details

Mathematical Formulation:

Stage 1: The genomic index is $I_1 = \mathbf{\beta}_1' \mathbf{\gamma}_1$

Coefficients: $\mathbf{\beta}_1 = \mathbf{\Gamma}_1^{-1}\mathbf{A}_1\mathbf{w}_1$

Stage 2: The index uses economic weights directly: $I_2 = \mathbf{w}' \mathbf{\gamma}$

Adjusted genomic covariance matrix:

$\mathbf{\Gamma}^* = \mathbf{\Gamma} - u \frac{\mathbf{A}_1'\mathbf{\beta}_1\mathbf{\beta}_1'\mathbf{A}_1}{\mathbf{\beta}_1'\mathbf{\Gamma}_1\mathbf{\beta}_1}$

Adjusted genotypic covariance matrix:

$\mathbf{C}^* = \mathbf{C} - u \frac{\mathbf{G}_1'\mathbf{b}_1\mathbf{b}_1'\mathbf{G}_1}{\mathbf{b}_1'\mathbf{P}_1\mathbf{b}_1}$

where $u = k_1(k_1 - \tau)$

Accuracy at stage 1: $\rho_{HI_1} = \sqrt{\frac{\mathbf{\beta}_1'\mathbf{\Gamma}_1\mathbf{\beta}_1}{\mathbf{w}'\mathbf{C}\mathbf{w}}}$

Accuracy at stage 2: $\rho_{HI_2} = \sqrt{\frac{\mathbf{w}'\mathbf{\Gamma}^*\mathbf{w}}{\mathbf{w}'\mathbf{C}^*\mathbf{w}}}$

Value

List with components:

• beta1 - Stage 1 genomic index coefficients

• w - Economic weights (stage 2 coefficients)

• stage1_metrics - List with stage 1 metrics (R1, E1, rho_HI1)

• stage2_metrics - List with stage 2 metrics (R2, E2, rho_HI2)

• Gamma_star - Adjusted genomic covariance matrix at stage 2

• C_star - Adjusted genotypic covariance matrix at stage 2

• rho_I1I2 - Correlation between stage 1 and stage 2 indices

• k1 - Selection intensity at stage 1

• k2 - Selection intensity at stage 2

• summary_stage1 - Data frame with stage 1 summary

• summary_stage2 - Data frame with stage 2 summary

References

Ceron-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Chapter 9, Section 9.4.

Examples

## Not run: 
# Two-stage genomic selection example
# Stage 1: Select based on GEBVs for 3 traits
# Stage 2: Select based on GEBVs for all 7 traits

# Compute covariance matrices
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Simulate GEBV covariances (in practice, compute from genomic prediction)
set.seed(123)
reliability <- 0.7
Gamma1 <- reliability * gmat[1:3, 1:3]
Gamma <- reliability * gmat
A1 <- reliability * gmat[1:3, 1:3]
A <- gmat[, 1:3]

# Economic weights
weights <- c(10, 8, 6, 4, 3, 2, 1)

# Run MLGSI
result <- mlgsi(
  Gamma1 = Gamma1, Gamma = Gamma, A1 = A1, A = A,
  C = gmat, G1 = gmat[1:3, 1:3], P1 = pmat[1:3, 1:3],
  wmat = weights, selection_proportion = 0.1
)

print(result$summary_stage1)
print(result$summary_stage2)

## End(Not run)

Description

Implements the two-stage Linear Phenotypic Selection Index where selection occurs at two different stages with covariance adjustments due to selection effects using Cochran/Cunningham's method.

Usage

mlpsi(
  P1,
  P,
  G1,
  C,
  wmat,
  wcol = 1,
  stage1_indices = NULL,
  selection_proportion = 0.1,
  use_young_method = FALSE,
  k1_manual = 2.063,
  k2_manual = 2.063,
  tau = NULL
)

Arguments

Details

Mathematical Formulation:

Stage 1 index coefficients:

$\mathbf{b}_1 = \mathbf{P}_1^{-1}\mathbf{G}_1\mathbf{w}$

Stage 2 index coefficients:

$\mathbf{b}_2 = \mathbf{P}^{-1}\mathbf{G}\mathbf{w}$

Adjusted phenotypic covariance matrix (Cochran/Cunningham):

$\mathbf{P}^* = \mathbf{P} - u \frac{Cov(\mathbf{y},\mathbf{x}_1)\mathbf{b}_1\mathbf{b}_1'Cov(\mathbf{x}_1,\mathbf{y})}{\mathbf{b}_1'\mathbf{P}_1\mathbf{b}_1}$

Adjusted genotypic covariance matrix:

$\mathbf{C}^* = \mathbf{C} - u \frac{\mathbf{G}_1'\mathbf{b}_1\mathbf{b}_1'\mathbf{G}_1}{\mathbf{b}_1'\mathbf{P}_1\mathbf{b}_1}$

where $u = k_1(k_1 - \tau)$

Selection response: $R_1 = k_1 \sqrt{\mathbf{b}_1'\mathbf{P}_1\mathbf{b}_1}$, $R_2 = k_2 \sqrt{\mathbf{b}_2'\mathbf{P}^*\mathbf{b}_2}$

Value

List with components:

• b1 - Stage 1 index coefficients

• b2 - Stage 2 index coefficients

• stage1_metrics - List with stage 1 metrics (R1, E1, rho_H1)

• stage2_metrics - List with stage 2 metrics (R2, E2, rho_H2)

• P_star - Adjusted phenotypic covariance matrix at stage 2

• C_star - Adjusted genotypic covariance matrix at stage 2

• rho_12 - Correlation between stage 1 and stage 2 indices

• k1 - Selection intensity at stage 1

• k2 - Selection intensity at stage 2

• summary_stage1 - Data frame with stage 1 summary

• summary_stage2 - Data frame with stage 2 summary

References

Cochran, W. G. (1951). Improvement by means of selection. Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 449-470.

Cunningham, E. P. (1975). Multi-stage index selection. Theoretical and Applied Genetics, 46(2), 55-61.

Young, S. S. Y. (1964). Multi-stage selection for genetic gain. Heredity, 19(1), 131-144.

Examples

## Not run: 
# Two-stage selection example
# Stage 1: Select based on 3 traits
# Stage 2: Select based on all 7 traits

# Compute variance-covariance matrices
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

# Stage 1 uses first 3 traits
P1 <- pmat[1:3, 1:3]
G1 <- gmat[1:3, 1:3]

# Stage 2 uses all 7 traits
P <- pmat
C <- gmat

# Economic weights
weights <- c(10, 8, 6, 4, 3, 2, 1)

# Run MLPSI (default: stage1_indices = 1:3)
result <- mlpsi(
  P1 = P1, P = P, G1 = G1, C = C, wmat = weights,
  selection_proportion = 0.1
)

# Or with non-contiguous traits (e.g., traits 1, 3, 5 at stage 1):
# P1 <- pmat[c(1,3,5), c(1,3,5)]
# G1 <- gmat[c(1,3,5), c(1,3,5)]
# result <- mlpsi(P1 = P1, P = P, G1 = G1, C = C, wmat = weights,
#                 stage1_indices = c(1, 3, 5))

print(result$summary_stage1)
print(result$summary_stage2)

## End(Not run)

Description

Implements the two-stage Predetermined Proportional Gain LGSI where breeders specify desired proportional gains between traits at each stage using GEBVs.

Usage

mppg_lgsi(
  Gamma1,
  Gamma,
  A1,
  A,
  C,
  G1,
  P1,
  wmat,
  wcol = 1,
  d1,
  d2,
  U1 = NULL,
  U2 = NULL,
  selection_proportion = 0.1,
  use_young_method = FALSE,
  k1_manual = 2.063,
  k2_manual = 2.063,
  tau = NULL
)

Arguments

Details

Mathematical Formulation:

The PPG genomic coefficients are:

$\mathbf{\beta}_{P_1} = \mathbf{\beta}_{R_1} + \theta_1 \mathbf{U}_1(\mathbf{U}_1'\mathbf{\Gamma}_1\mathbf{U}_1)^{-1}\mathbf{d}_1$

$\mathbf{\beta}_{P_2} = \mathbf{\beta}_{R_2} + \theta_2 \mathbf{U}_2(\mathbf{U}_2'\mathbf{\Gamma}\mathbf{U}_2)^{-1}\mathbf{d}_2$

where proportionality constants are:

$\theta_1 = \frac{\mathbf{d}_1'(\mathbf{U}_1'\mathbf{\Gamma}_1\mathbf{U}_1)^{-1}\mathbf{U}_1'\mathbf{A}_1\mathbf{w}}{\mathbf{d}_1'(\mathbf{U}_1'\mathbf{\Gamma}_1\mathbf{U}_1)^{-1}\mathbf{d}_1}$

Covariance Adjustment:

The genetic covariance matrix $\mathbf{C}^*$ is adjusted using phenotypic PPG coefficients $\mathbf{b}_{P1} = \mathbf{P}_1^{-1}\mathbf{G}_1\mathbf{P}_1^{-1}\mathbf{d}_1$, which reflect the same proportional gain constraints as the genomic coefficients $\mathbf{\beta}_{P1}$. This ensures the adjustment reflects the actual PPG selection occurring at stage 1.

Important: When using custom U1 matrices (subset constraints), the phenotypic proxy $\mathbf{b}_{P1}$ uses the standard Tallis formula (all traits constrained), while the genomic index $\mathbf{\beta}_{P1}$ respects the U1 subset. This may cause $\mathbf{C}^*$ to be slightly over-adjusted. For exact adjustment, use U1 = NULL (default, all traits constrained). Calculating the exact restricted phenotypic proxy would require implementing the full MPPG-LPSI projection matrix method for the phenotypic coefficients.

Note: Input covariance matrices (C, P1, Gamma1, Gamma) should be positive definite. Non-positive definite matrices may lead to invalid results or warnings.

Value

List with components similar to mlgsi, plus:

• beta_P1 - PPG genomic stage 1 coefficients

• beta_P2 - PPG genomic stage 2 coefficients

• b_P1 - PPG phenotypic stage 1 coefficients (used for C* adjustment)

• theta1 - Proportionality constant for stage 1

• theta2 - Proportionality constant for stage 2

• gain_ratios_1 - Achieved gain ratios at stage 1

• gain_ratios_2 - Achieved gain ratios at stage 2

References

Ceron-Rojas, J. J., & Crossa, J. (2018). Linear Selection Indices in Modern Plant Breeding. Springer International Publishing. Chapter 9, Section 9.6.

Examples

## Not run: 
# Two-stage proportional gain genomic selection
gmat <- gen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])
pmat <- phen_varcov(seldata[, 3:9], seldata[, 2], seldata[, 1])

reliability <- 0.7
Gamma1 <- reliability * gmat[1:3, 1:3]
Gamma <- reliability * gmat
A1 <- reliability * gmat[1:3, 1:3]
A <- gmat[, 1:3]

# Desired proportional gains
d1 <- c(2, 1, 1)
d2 <- c(3, 2, 1, 1, 1, 0.5, 0.5)

weights <- c(10, 8, 6, 4, 3, 2, 1)

result <- mppg_lgsi(
  Gamma1 = Gamma1, Gamma = Gamma, A1 = A1, A = A,
  C = gmat, G1 = gmat[1:3, 1:3], P1 = pmat[1:3, 1:3],
  wmat = weights, d1 = d1, d2 = d2
)

## End(Not run)

Description

Implements the two-stage Predetermined Proportional Gain LPSI where breeders specify desired proportional gains between traits at each stage.

Usage

mppg_lpsi(
  P1,
  P,
  G1,
  C,
  wmat,
  wcol = 1,
  d1,
  d2,
  stage1_indices = NULL,
  selection_proportion = 0.1,
  use_young_method = FALSE,
  k1_manual = 2.063,
  k2_manual = 2.063,
  tau = NULL
)

Arguments

Details

Mathematical Formulation (Chapter 9.3.1, Eq 9.17):

The PPG coefficients are computed using the projection matrix method:

$\mathbf{b}_{M_1} = \mathbf{b}_{R_1} + \theta_1 \mathbf{U}_1(\mathbf{U}_1'\mathbf{G}_1\mathbf{P}_1^{-1}\mathbf{G}_1\mathbf{U}_1)^{-1}\mathbf{d}_1$