From breeding data to BreedReliabR: start here
Purpose
BreedReliabR is a downstream decision-reliability package. It does not start from a raw yield spreadsheet, SNP file, or weather table. It starts after an upstream model has produced candidate environmental-response coefficient means and their full joint uncertainty covariance.
The workflow is:
breeding data
-> upstream prediction model
-> candidate coefficient means + full joint covariance
-> TPE
-> BreedReliabR
-> fixed-k advancement analysisThis vignette uses a small genotype-specific reaction-norm model fitted with base R only to make the hand-off transparent. It is not presented as a production genomic prediction model.
1. Simulate a small multi-environment trial
set.seed(2026)
n_gen <- 8
n_env <- 6
n_rep <- 2
genotypes <- paste0("G", sprintf("%02d", seq_len(n_gen)))
environments <- paste0("E", sprintf("%02d", seq_len(n_env)))
env <- data.frame(
Environment = environments,
Heat = c(24, 27, 29, 31, 33, 35),
Water = c(72, 88, 60, 80, 55, 75)
)
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")
training_basis <- as.matrix(env[, basis_names])
rownames(training_basis) <- env$EnvironmentThe training transformation is now frozen. A real analysis must apply the same transformation to target environments rather than standardising the TPE separately.
2. Generate synthetic phenotype observations
true_coef <- cbind(
intercept = rnorm(n_gen, 5.0, 0.30),
heat_z = rnorm(n_gen, -0.25, 0.12),
water_z = rnorm(n_gen, 0.25, 0.12)
)
rownames(true_coef) <- genotypes
dat <- expand.grid(
Genotype = genotypes,
Environment = environments,
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, genotypes),
match(dat$Environment, environments),
dat$Rep),
]
rownames(dat) <- NULL
g_index <- match(dat$Genotype, genotypes)
dat$Yield <-
true_coef[g_index, "intercept"] +
true_coef[g_index, "heat_z"] * dat$heat_z +
true_coef[g_index, "water_z"] * dat$water_z +
rnorm(nrow(dat), 0, 0.25)
head(dat[, c("Genotype", "Environment", "Rep", "Yield")])
#> Genotype Environment Rep Yield
#> 1 G01 E01 1 5.519910
#> 2 G01 E01 2 5.984879
#> 3 G01 E02 1 6.184300
#> 4 G01 E02 2 5.766100
#> 5 G01 E03 1 5.060198
#> 6 G01 E03 2 5.330273At this point the data are still ordinary multi-environment trial observations. BreedReliabR has not yet been used.
3. Fit a transparent reaction-norm model outside BreedReliabR
Each genotype receives an intercept, heat-response coefficient, and water-response coefficient.
p <- length(basis_names)
X <- matrix(0, nrow = nrow(dat), ncol = n_gen * p)
for (i in seq_along(genotypes)) {
rows <- dat$Genotype == genotypes[i]
cols <- ((i - 1) * p + 1):(i * p)
X[rows, cols] <- as.matrix(dat[rows, basis_names])
}
coefficient_order <- unlist(
lapply(genotypes, function(g) paste(g, basis_names, sep = "__")),
use.names = FALSE
)
colnames(X) <- coefficient_order
fit <- lm.fit(X, dat$Yield)
stopifnot(fit$rank == ncol(X))
beta_hat <- fit$coefficients
coef_mean <- matrix(
beta_hat,
nrow = n_gen,
ncol = p,
byrow = TRUE,
dimnames = list(genotypes, basis_names)
)
sigma2 <- sum(fit$residuals^2) / fit$df.residual
coef_vcov <- sigma2 * solve(crossprod(X))
rownames(coef_vcov) <- coefficient_order
colnames(coef_vcov) <- coefficient_order
dim(coef_mean)
#> [1] 8 3
dim(coef_vcov)
#> [1] 24 24For eight candidates and three coefficients, BreedReliabR receives an 8 x 3 mean matrix and a 24 x 24 joint covariance.
In this pedagogical fixed-effect model, coef_vcov is the
joint coefficient-estimation covariance. It demonstrates the required
matrix structure; it is not a genomic prediction-error covariance and
should not be substituted for the PEV from a properly fitted genomic
model.
4. Verify ordering and covariance structure
expected_order <- unlist(
lapply(genotypes, function(g) paste(g, basis_names, sep = "__")),
use.names = FALSE
)
stopifnot(
identical(rownames(coef_vcov), expected_order),
identical(colnames(coef_vcov), expected_order),
all(is.finite(coef_mean)),
all(is.finite(coef_vcov)),
isSymmetric(coef_vcov, tol = 1e-8)
)
pev_eigen <- eigen(coef_vcov, symmetric = TRUE, only.values = TRUE)$values
pev_scale <- max(1, max(abs(pev_eigen)))
stopifnot(min(pev_eigen) >= -1e-8 * pev_scale)Candidate-major ordering means all basis coefficients for candidate 1, then all coefficients for candidate 2, and so on. Correct dimensions alone cannot prove correct ordering.
5. Specify the target population of environments
tpe_raw <- data.frame(
State = paste0("TPE", 1:4),
Heat = c(26, 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$WeightThe weights are assumptions describing the target distribution over the supplied states. They are not automatically acreage or production-area weights.
6. Preflight the hand-off
input_check <- br_check_inputs(
candidates = 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 6 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.A PASS establishes structural compatibility only. It does not validate the upstream model, covariance adequacy, biological definition of the TPE, or empirical probability calibration.
7. BreedReliabR begins here
prediction <- br_prediction(
candidates = genotypes,
coef_mean = coef_mean,
coef_vcov = coef_vcov,
basis_names = basis_names,
source = "Pedagogical genotype-specific reaction-norm model"
)
tpe <- br_tpe(
basis = tpe_basis,
weights = tpe_weights,
state_id = tpe_raw$State
)
prediction
#> BreedReliabR prediction object
#> Candidates: 8
#> Environmental basis terms: 3
#> Source: Pedagogical genotype-specific reaction-norm model
tpe
#> BreedReliabR target population of environments
#> States: 4
#> Basis terms: 3
#> Weight sum: 18. Diagnose geometric support
support <- br_support(tpe, training_basis)
support
#> BreedReliabR environmental support diagnostics
#> weight_outside_range weight_constant_mismatch weight_beyond_training_max_nn
#> 0 0 0
#> training_median_nn training_max_nn
#> 1.074186 1.495982
#> Geometric support only; these diagnostics do not establish biological correctness of the specified TPE or prediction calibration.
plot(support)
br_support() is a geometric extrapolation diagnostic. It
does not prove biological validity, prediction accuracy, or
calibration.
9. Compare environmental and prediction uncertainty
omega <- br_omega(prediction, tpe, H = c(1, 3, 5, 10))
omega$summary
#> H median_omega median_pred_var median_env_var
#> 1 1 26.987845 0.005748804 0.15514782
#> 2 3 8.995948 0.005748804 0.05171594
#> 3 5 5.397569 0.005748804 0.03102956
#> 4 10 2.698785 0.005748804 0.01551478
head(omega$candidate)
#> candidate H expected_merit pred_var env_var omega
#> 1 G01 1 5.117879 0.005748804 0.22261566 38.723823
#> 2 G01 3 5.117879 0.005748804 0.07420522 12.907941
#> 3 G01 5 5.117879 0.005748804 0.04452313 7.744765
#> 4 G01 10 5.117879 0.005748804 0.02226157 3.872382
#> 5 G02 1 4.495687 0.005748804 0.16483499 28.672920
#> 6 G02 3 4.495687 0.005748804 0.05494500 9.557640Omega is the ratio of finite-TPE environmental variance to prediction variance. It is not a probability, selection gain, genetic gain, or universal threshold.
10. Compare M0-M3 and estimate advancement probabilities
comparison <- br_compare(
prediction, tpe, H = 5, k = 3,
nsim = 1000, seed = 2026
)
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.8160000 0.8160000
#> 3 M2 5 3 1000 0.8883333 0.8883333
#> 4 M3 5 3 1000 0.8066667 0.8066667
#> probability_sum invariant_error
#> 1 3 0
#> 2 3 0
#> 3 3 0
#> 4 3 0
advancement <- br_advance(
prediction, tpe, H = 5, k = 3,
nsim = 2000, method = "M3", seed = 2026
)
advancement$candidates
#> candidate expected_merit p_advance mcse in_expected_merit_set
#> 1 G01 5.117879 1.0000 0.000000000 TRUE
#> 2 G03 4.779040 0.8655 0.007629212 TRUE
#> 3 G04 4.691907 0.5425 0.011139878 TRUE
#> 4 G05 4.667847 0.3525 0.010682784 FALSE
#> 5 G07 4.599909 0.1220 0.007318333 FALSE
#> 6 G08 4.544703 0.1085 0.006954414 FALSE
#> 7 G02 4.495687 0.0090 0.002111753 FALSE
#> 8 G06 4.132494 0.0000 0.000000000 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] 3M0 is deterministic point ranking. M1 propagates prediction uncertainty, M2 finite-TPE uncertainty, and M3 both. These are uncertainty formulations of the same prediction object, not different genomic models.
For M1-M3, p_advance is a conditional top-k inclusion
probability. Selected-set reliability is an expected retained fraction,
not the probability of recovering the complete set.
11. Examine sensitivity to H
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.8160000 0.8160000
#> 3 M2 1 3 1000 0.8340000 0.8340000
#> 4 M3 1 3 1000 0.7690000 0.7690000
#> 5 M0 3 3 1 NA NA
#> 6 M1 3 3 1000 0.8160000 0.8160000
#> 7 M2 3 3 1000 0.8703333 0.8703333
#> 8 M3 3 3 1000 0.7930000 0.7930000
#> 9 M0 5 3 1 NA NA
#> 10 M1 5 3 1000 0.8160000 0.8160000
#> 11 M2 5 3 1000 0.8883333 0.8883333
#> 12 M3 5 3 1000 0.8066667 0.8066667
#> 13 M0 10 3 1 NA NA
#> 14 M1 10 3 1000 0.8160000 0.8160000
#> 15 M2 10 3 1000 0.8930000 0.8930000
#> 16 M3 10 3 1000 0.8030000 0.8030000
#> 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 26.987845 0.005748804 0.15514782
#> 2 3 8.995948 0.005748804 0.05171594
#> 3 5 5.397569 0.005748804 0.03102956
#> 4 10 2.698785 0.005748804 0.01551478
plot(curve_results, type = "reliability")
H is the number of independent future draws from the
specified TPE whose average represents the future environmental target
under the current formulation. It is not the number of training
environments, and BreedReliabR does not automatically identify an
optimal H.
12. Export the operational hand-off
template_dir <- file.path(tempdir(), "BreedReliabR_beginner_example")
dir.create(template_dir, showWarnings = FALSE, recursive = TRUE)
candidate_coefficients <- data.frame(
candidate = rownames(coef_mean),
coef_mean,
check.names = FALSE
)
write.csv(
candidate_coefficients,
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
)These four files illustrate the operational boundary between an upstream model-fitting workflow and BreedReliabR.
13. Common mistakes
Do not supply raw SNP or phenotype tables directly to
br_prediction(). Do not replace the full joint covariance
with marginal standard errors. Do not fit a separate standardisation or
PCA to target environments. Do not interpret H as training size, Omega
as selection gain, or selected-set reliability as exact-set recovery
probability.
Increasing nsim reduces Monte Carlo error only; it
cannot correct an inadequate upstream predictive distribution or poorly
specified TPE.
- Purpose
- 1. Simulate a small multi-environment trial
- 2. Generate synthetic phenotype observations
- 3. Fit a transparent reaction-norm model outside BreedReliabR
- 4. Verify ordering and covariance structure
- 5. Specify the target population of environments
- 6. Preflight the hand-off
- 7. BreedReliabR begins here
- 8. Diagnose geometric support
- 9. Compare environmental and prediction uncertainty
- 10. Compare M0-M3 and estimate advancement probabilities
- 11. Examine sensitivity to H
- 12. Export the operational hand-off
- 13. Common mistakes