Genomic GxE prediction to BreedReliabR
Purpose
This vignette extends the beginner workflow to a genomic multi-environment breeding programme. Phenotyped training genotypes and completely unphenotyped advancement candidates are separated. Candidate environmental-response coefficients are predicted through genomic relationships, their full joint prediction-error covariance is constructed explicitly, and only then does BreedReliabR begin.
The transparent mixed-model calculation is a reference implementation for the hand-off. It is not intended to replace specialised genomic-model software.
1. Experimental dimensions
set.seed(2026)
n_train <- 24
n_candidate <- 8
n_total <- n_train + n_candidate
n_env <- 8
n_rep <- 2
n_markers <- 150
training_genotypes <- paste0("T", sprintf("%02d", seq_len(n_train)))
candidate_genotypes <- paste0("C", sprintf("%02d", seq_len(n_candidate)))
all_genotypes <- c(training_genotypes, candidate_genotypes)The 24 training genotypes have phenotypes and SNPs. The eight candidates have SNPs but no phenotype records.
2. Training environments and frozen basis
env <- data.frame(
Environment = paste0("E", sprintf("%02d", seq_len(n_env))),
Heat = c(24, 27, 29, 31, 33, 35, 28, 32),
Water = c(72, 88, 60, 80, 55, 75, 65, 90)
)
cor(env$Heat, env$Water)
#> [1] -0.04071977
heat_mean <- mean(env$Heat)
heat_sd <- sd(env$Heat)
water_mean <- mean(env$Water)
water_sd <- sd(env$Water)
env$heat_z <- (env$Heat - heat_mean) / heat_sd
env$water_z <- (env$Water - water_mean) / water_sd
env$intercept <- 1
basis_names <- c("intercept", "heat_z", "water_z")
p <- length(basis_names)
training_basis <- as.matrix(env[, basis_names])
rownames(training_basis) <- env$EnvironmentA real analysis should choose environmental descriptors on biological and statistical grounds. The low correlation here is only a teaching design that keeps separate response coefficients interpretable.
3. SNP markers and genomic relationship matrix
marker_frequency <- runif(n_markers, min = 0.10, max = 0.50)
M <- sapply(
marker_frequency,
function(freq) rbinom(n_total, size = 2, prob = freq)
)
M <- as.matrix(M)
rownames(M) <- all_genotypes
colnames(M) <- paste0("SNP", seq_len(n_markers))
stopifnot(all(dim(M) == c(n_total, n_markers)))
p_train <- colMeans(M[training_genotypes, , drop = FALSE]) / 2
W <- sweep(M, MARGIN = 2, STATS = 2 * p_train, FUN = "-")
relationship_denominator <- 2 * sum(p_train * (1 - p_train))
stopifnot(is.finite(relationship_denominator), relationship_denominator > 0)
K_raw <- tcrossprod(W) / relationship_denominator
K_raw <- (K_raw + t(K_raw)) / 2
rownames(K_raw) <- all_genotypes
colnames(K_raw) <- all_genotypes
ridge <- 1e-6 * max(1, mean(diag(K_raw)))
K_model <- K_raw + diag(ridge, n_total)
rownames(K_model) <- all_genotypes
colnames(K_model) <- all_genotypes
relationship_rcond <- rcond(K_model)
stopifnot(is.finite(relationship_rcond), relationship_rcond > 0)The marker simulation demonstrates the interface only. It is not intended to reproduce realistic linkage disequilibrium, family structure, population stratification, or breeding demography.
The small ridge is an explicitly specified part of this upstream reference model so that the relationship matrix can be inverted. It is not a repair of a BreedReliabR PEV.
4. Genomic covariance among environmental responses
G_beta <- matrix(
c(
0.20, 0.02, -0.01,
0.02, 0.06, -0.015,
-0.01, -0.015, 0.05
),
nrow = p,
byrow = TRUE,
dimnames = list(basis_names, basis_names)
)
G_beta_eigen <- eigen(G_beta, symmetric = TRUE, only.values = TRUE)$values
stopifnot(min(G_beta_eigen) > 0)
G_u <- kronecker(K_model, G_beta)
stopifnot(all(dim(G_u) == c(n_total * p, n_total * p)))The ordering induced by K_model %x% G_beta is
genotype-major: all three coefficients for genotype 1, then all three
for genotype 2, and so on.
5. Simulate genomic effects and training phenotypes
chol_G <- chol(G_u)
u_true <- as.vector(t(chol_G) %*% rnorm(n_total * p))
beta_true <- c(
intercept = 5.00,
heat_z = -0.20,
water_z = 0.25
)
dat <- expand.grid(
Genotype = training_genotypes,
Environment = env$Environment,
Rep = seq_len(n_rep),
KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
)
dat <- merge(dat, env, by = "Environment", sort = FALSE)
dat <- dat[
order(match(dat$Genotype, training_genotypes),
match(dat$Environment, env$Environment),
dat$Rep),
]
rownames(dat) <- NULL
stopifnot(
nrow(dat) == n_train * n_env * n_rep,
!any(dat$Genotype %in% candidate_genotypes)
)
X <- as.matrix(dat[, basis_names])
Z <- matrix(0, nrow = nrow(dat), ncol = n_total * p)
for (i in seq_len(n_train)) {
rows <- dat$Genotype == training_genotypes[i]
cols <- ((i - 1) * p + 1):(i * p)
Z[rows, cols] <- as.matrix(dat[rows, basis_names])
}
sigma_e2 <- 0.25^2
dat$Yield <- as.vector(
X %*% beta_true + Z %*% u_true +
rnorm(nrow(dat), mean = 0, sd = sqrt(sigma_e2))
)Candidates have no phenotype records. Their predictions arise through genomic covariance with the phenotyped training population.
6. Fit the transparent genomic reaction-norm model
The model is y = X beta + Z u + e, with
u ~ N(0, G_u) and e ~ N(0, sigma_e2 I). The
covariance parameters are treated as known in this controlled reference
example so the prediction/PEV hand-off can be inspected directly.
G_u_inv <- solve(G_u)
R_inv_scalar <- 1 / sigma_e2
C11 <- crossprod(X) * R_inv_scalar
C12 <- crossprod(X, Z) * R_inv_scalar
C21 <- t(C12)
C22 <- crossprod(Z) * R_inv_scalar + G_u_inv
C <- rbind(cbind(C11, C12), cbind(C21, C22))
stopifnot(all(dim(C) == c(p + n_total * p, p + n_total * p)))
C_rcond <- rcond(C)
stopifnot(is.finite(C_rcond), C_rcond > 0)
rhs <- c(
as.vector(crossprod(X, dat$Yield) * R_inv_scalar),
as.vector(crossprod(Z, dat$Yield) * R_inv_scalar)
)
theta_hat <- solve(C, rhs)A real genomic analysis normally estimates variance components. BreedReliabR can only propagate uncertainty that the upstream model supplies; it cannot restore uncertainty omitted by a plug-in PEV.
7. Define complete candidate coefficient targets
BreedReliabR requires each candidate’s complete coefficient
fixed population coefficient + genomic deviation, not just
its random BLUP deviation.
n_targets <- n_candidate * p
D <- matrix(0, nrow = n_targets, ncol = length(theta_hat))
target_names <- character(n_targets)
candidate_index <- (n_train + 1):n_total
row_id <- 1L
for (g in candidate_index) {
for (j in seq_len(p)) {
D[row_id, j] <- 1
random_position <- p + (g - 1) * p + j
D[row_id, random_position] <- 1
target_names[row_id] <- paste(all_genotypes[g], basis_names[j], sep = "__")
row_id <- row_id + 1L
}
}
rownames(D) <- target_names
candidate_coef_vector <- as.vector(D %*% theta_hat)
coef_mean <- matrix(
candidate_coef_vector,
nrow = n_candidate,
ncol = p,
byrow = TRUE,
dimnames = list(candidate_genotypes, basis_names)
)
stopifnot(all(dim(coef_mean) == c(n_candidate, p)))
round(coef_mean, 4)
#> intercept heat_z water_z
#> C01 4.8931 -0.1550 0.2324
#> C02 4.8239 -0.2313 0.1370
#> C03 5.1101 -0.0898 0.2143
#> C04 4.8765 -0.1293 0.0819
#> C05 5.0599 -0.2231 0.2064
#> C06 4.5490 -0.2865 0.2875
#> C07 4.9791 -0.1175 0.2540
#> C08 4.9058 -0.1098 0.23038. Derive the full joint candidate PEV
For the target vector t = D theta, the joint PEV is
D C^{-1} D'. The full inverse of C is not
required.
C_inv_Dt <- solve(C, t(D))
coef_vcov <- D %*% C_inv_Dt
rownames(coef_vcov) <- target_names
colnames(coef_vcov) <- target_names
stopifnot(all(dim(coef_vcov) == c(n_targets, n_targets)))
pev_scale <- max(1, max(abs(coef_vcov)))
pev_asymmetry <- max(abs(coef_vcov - t(coef_vcov)))
stopifnot(pev_asymmetry <= 1e-8 * pev_scale)
pev_eigen <- eigen(coef_vcov, symmetric = TRUE, only.values = TRUE)$values
eigen_scale <- max(1, max(abs(pev_eigen)))
stopifnot(min(pev_eigen) >= -1e-8 * eigen_scale)
expected_order <- unlist(
lapply(candidate_genotypes, function(g) paste(g, basis_names, sep = "__")),
use.names = FALSE
)
stopifnot(
identical(rownames(coef_vcov), expected_order),
identical(colnames(coef_vcov), expected_order)
)Do not silently symmetrise, bend, or diagonalise an invalid upstream covariance merely to obtain a passing BreedReliabR input.
Marginal standard errors can be useful descriptively but are not a substitute for this full matrix.
9. Define the TPE using the frozen training transformation
tpe_raw <- data.frame(
State = paste0("TPE", seq_len(4)),
Heat = c(27, 30, 33, 35),
Water = c(82, 72, 63, 56),
Weight = c(0.20, 0.30, 0.30, 0.20)
)
tpe_raw$heat_z <- (tpe_raw$Heat - heat_mean) / heat_sd
tpe_raw$water_z <- (tpe_raw$Water - water_mean) / water_sd
tpe_raw$intercept <- 1
tpe_basis <- as.matrix(tpe_raw[, basis_names])
rownames(tpe_basis) <- tpe_raw$State
tpe_weights <- tpe_raw$Weight
stopifnot(identical(colnames(coef_mean), colnames(tpe_basis)))Target states may naturally have correlated heat and water conditions. The important rule is that they are projected into the same frozen basis used by the prediction model.
10. Preflight the genomic hand-off
input_check <- br_check_inputs(
candidates = candidate_genotypes,
coef_mean = coef_mean,
coef_vcov = coef_vcov,
basis_names = basis_names,
tpe_basis = tpe_basis,
tpe_weights = tpe_weights,
training_basis = training_basis,
k = 3,
H = 5
)
input_check
#> BreedReliabR input audit
#> ========================
#>
#> Candidates 8
#> Basis coefficients 3
#> Coefficient matrix 8 x 3
#> Expected joint covariance 24 x 24
#> Supplied joint covariance 24 x 24
#>
#> PASS Candidate identifiers 8 candidate identifiers supplied.
#> PASS Candidate uniqueness Candidate identifiers are unique.
#> PASS Coefficient matrix Supplied coefficient matrix is 8 x 3 .
#> PASS Coefficient rows Coefficient rows match the candidate count.
#> PASS Coefficient values All coefficient values are finite.
#> PASS Coefficient row labels Coefficient rows match candidate order.
#> PASS Basis names 3 unique basis coefficients supplied.
#> PASS Coefficient column labels Coefficient columns match `basis_names`.
#> PASS Covariance dimensions Joint covariance dimension is 24 x 24 .
#> PASS Covariance values All covariance entries are finite.
#> PASS Covariance symmetry Covariance is symmetric within tolerance.
#> PASS Covariance PSD Covariance is positive semidefinite within tolerance.
#> PASS Downstream numerical margin Covariance is also within the tighter stochastic-sampling margin.
#> PASS Candidate-major labels Canonical candidate-major labels are consistent.
#> PASS TPE dimensions TPE contains 4 target states.
#> PASS TPE values All TPE basis values are finite.
#> PASS TPE basis labels TPE basis matches the prediction basis.
#> PASS TPE weights TPE weights are valid and sum to one.
#> PASS Training basis dimensions Training basis contains 8 states.
#> PASS Training basis values All training-basis values are finite.
#> PASS Training basis labels Training basis matches the prediction basis.
#> PASS Advancement size k k = 3 is valid.
#> PASS Future representation H Valid H value(s): 5
#>
#> Overall status: PASS
#>
#> Structural compatibility does not validate the upstream prediction model, adequacy of its covariance, biological definition of the TPE, or empirical probability calibration.The preflight audit is structural only. It does not validate genomic-model adequacy, variance-component estimation, TPE biology, or empirical calibration.
11. BreedReliabR begins here
prediction <- br_prediction(
candidates = candidate_genotypes,
coef_mean = coef_mean,
coef_vcov = coef_vcov,
basis_names = basis_names,
source = "Transparent genomic reaction-norm reference model",
metadata = list(
training_genotypes = n_train,
candidate_genotypes = n_candidate,
markers = n_markers,
environments = n_env,
candidate_phenotypes = 0
)
)
tpe <- br_tpe(
basis = tpe_basis,
weights = tpe_weights,
state_id = tpe_raw$State
)The SNP matrix, genomic relationship matrix, raw phenotypes, and
mixed-model equations do not enter BreedReliabR directly. Their
influence is carried through coef_mean and
coef_vcov.
12. Support, uncertainty, and advancement
support <- br_support(tpe, training_basis)
omega <- br_omega(prediction, tpe, H = c(1, 3, 5, 10))
comparison <- br_compare(
prediction, tpe, H = 5, k = 3,
nsim = 1000, seed = 2026
)
advancement <- br_advance(
prediction, tpe, H = 5, k = 3,
nsim = 2000, method = "M3", seed = 2026
)
support$summary
#> weight_outside_range weight_constant_mismatch weight_beyond_training_max_nn
#> 1 0 0 0
#> training_median_nn training_max_nn
#> 1 1.013575 1.287813
omega$summary
#> H median_omega median_pred_var median_env_var
#> 1 1 0.64900249 0.2210285 0.14765335
#> 2 3 0.21633416 0.2210285 0.04921778
#> 3 5 0.12980050 0.2210285 0.02953067
#> 4 10 0.06490025 0.2210285 0.01476533
summary(comparison)
#> method H k nsim expected_merit_set_reliability probability_set_reliability
#> 1 M0 5 3 1 NA NA
#> 2 M1 5 3 1000 0.5023333 0.5023333
#> 3 M2 5 3 1000 0.9343333 0.9343333
#> 4 M3 5 3 1000 0.4996667 0.4996667
#> probability_sum invariant_error
#> 1 3 0
#> 2 3 0
#> 3 3 0
#> 4 3 0
advancement$candidates
#> candidate expected_merit p_advance mcse in_expected_merit_set
#> 1 C03 4.989239 0.5640 0.011088372 TRUE
#> 2 C05 4.888837 0.4815 0.011172684 TRUE
#> 3 C07 4.831466 0.4180 0.011028962 TRUE
#> 4 C04 4.792353 0.3945 0.010928626 FALSE
#> 5 C08 4.770645 0.3765 0.010833922 FALSE
#> 6 C01 4.738952 0.3515 0.010675855 FALSE
#> 7 C02 4.677075 0.3065 0.010309165 FALSE
#> 8 C06 4.320493 0.1075 0.006926173 FALSE
#> in_probability_set
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 FALSE
#> 5 FALSE
#> 6 FALSE
#> 7 FALSE
#> 8 FALSE
sum(advancement$candidates$p_advance)
#> [1] 3p_advance is a conditional fixed-k inclusion
probability. It is not a probability of biological superiority.
Selected-set reliability is an expected retained fraction, not exact-set
recovery probability. Omega is a relative environmental-to-prediction
uncertainty diagnostic, not selection gain.
13. Sensitivity to future environmental representation
curve_results <- br_curve(
prediction, tpe,
H = c(1, 3, 5, 10),
k = 3,
nsim = 1000,
seed = 2026
)
summary(curve_results)
#> $reliability
#> method H k nsim expected_merit_set_reliability probability_set_reliability
#> 1 M0 1 3 1 NA NA
#> 2 M1 1 3 1000 0.5023333 0.5023333
#> 3 M2 1 3 1000 0.8340000 0.8340000
#> 4 M3 1 3 1000 0.4840000 0.4840000
#> 5 M0 3 3 1 NA NA
#> 6 M1 3 3 1000 0.5023333 0.5023333
#> 7 M2 3 3 1000 0.8976667 0.8976667
#> 8 M3 3 3 1000 0.4980000 0.4980000
#> 9 M0 5 3 1 NA NA
#> 10 M1 5 3 1000 0.5023333 0.5023333
#> 11 M2 5 3 1000 0.9343333 0.9343333
#> 12 M3 5 3 1000 0.4996667 0.4996667
#> 13 M0 10 3 1 NA NA
#> 14 M1 10 3 1000 0.5023333 0.5023333
#> 15 M2 10 3 1000 0.9553333 0.9553333
#> 16 M3 10 3 1000 0.4983333 0.4983333
#> probability_sum invariant_error
#> 1 3 0
#> 2 3 0
#> 3 3 0
#> 4 3 0
#> 5 3 0
#> 6 3 0
#> 7 3 0
#> 8 3 0
#> 9 3 0
#> 10 3 0
#> 11 3 0
#> 12 3 0
#> 13 3 0
#> 14 3 0
#> 15 3 0
#> 16 3 0
#>
#> $omega
#> H median_omega median_pred_var median_env_var
#> 1 1 0.64900249 0.2210285 0.14765335
#> 2 3 0.21633416 0.2210285 0.04921778
#> 3 5 0.12980050 0.2210285 0.02953067
#> 4 10 0.06490025 0.2210285 0.01476533
plot(curve_results, type = "reliability")
BreedReliabR does not automatically identify an optimal H. The curve shows how the supplied uncertainty model responds as future environmental representation changes.
14. Operational export
template_dir <- file.path(tempdir(), "BreedReliabR_genomic_example")
dir.create(template_dir, showWarnings = FALSE, recursive = TRUE)
write.csv(
data.frame(candidate = candidate_genotypes, coef_mean, check.names = FALSE),
file.path(template_dir, "candidate_coefficients.csv"),
row.names = FALSE
)
saveRDS(coef_vcov, file.path(template_dir, "candidate_joint_PEV.rds"))
write.csv(
data.frame(state = tpe_raw$State, tpe_basis, weight = tpe_weights,
check.names = FALSE),
file.path(template_dir, "target_environments.csv"),
row.names = FALSE
)
write.csv(
data.frame(environment = rownames(training_basis), training_basis,
check.names = FALSE),
file.path(template_dir, "training_environments.csv"),
row.names = FALSE
)15. Production-model boundary
The reference example treats covariance parameters as known. A production genomic workflow should estimate and validate its own variance structure and provide the appropriate full joint candidate PEV.
A practical adapter for a package such as sommer should
be labelled verified only after its candidate coefficient targets and
complete PEV agree numerically with an independently constructed
reference under the same fitted model. Marginal standard errors alone
are insufficient.
- Purpose
- 1. Experimental dimensions
- 2. Training environments and frozen basis
- 3. SNP markers and genomic relationship matrix
- 4. Genomic covariance among environmental responses
- 5. Simulate genomic effects and training phenotypes
- 6. Fit the transparent genomic reaction-norm model
- 7. Define complete candidate coefficient targets
- 8. Derive the full joint candidate PEV
- 9. Define the TPE using the frozen training transformation
- 10. Preflight the genomic hand-off
- 11. BreedReliabR begins here
- 12. Support, uncertainty, and advancement
- 13. Sensitivity to future environmental representation
- 14. Operational export
- 15. Production-model boundary