BreedReliabR: Methodology and empirical
validation
Ashutosh Sawarkar
2026-10-01
This vignette describes an offline, model-agnostic workflow. The
included code uses the package’s simulated maize example; it does not
reproduce or refit the G2F analyses. The external model fit, coefficient
means, and complete prediction-error covariance must be supplied by the
user’s validated workflow.
Workflow
external model fitting
-> coefficient means + full PEV
-> br_prediction()
-> br_tpe()
-> br_support()
-> br_omega()
-> br_compare() / br_advance()
-> br_curve()
-> external validation against independent outcomes
The candidate-major coefficient order is fixed as candidate 1’s basis
terms, then candidate 2’s terms, and so on. The prediction covariance
must retain cross-candidate and cross-coefficient terms.
br_prediction() checks dimension, symmetry, and
positive-semidefinite behavior, but it cannot verify the upstream
model’s biological or statistical validity.
Import the model outputs and specify a TPE
The simulated example has a compact pre-specified prediction object
and five hypothetical TPE states. For a real analysis, replace these
objects with coefficients and a full PEV from the external fit, without
changing the package input order.
library(BreedReliabR)
data(maize_demo)
prediction <- do.call(br_prediction, maize_demo$prediction_inputs)
tpe <- br_tpe(maize_demo$tpe_basis,
weights = maize_demo$tpe_states$weight,
state_id = maize_demo$tpe_states$state_id)
prediction
#> BreedReliabR prediction object
#> Candidates: 8
#> Environmental basis terms: 3
#> Source: Simulated fit-import inputs; not estimated from training phenotypes
tpe
#> BreedReliabR target population of environments
#> States: 5
#> Basis terms: 3
#> Weight sum: 1
Weights define an assumed target distribution. Equal weights are
useful when they match the target question, but they are not
automatically production-area weights.
Check environmental support and uncertainty scale
support <- br_support(tpe, maize_demo$training_basis)
support$summary
#> weight_outside_range weight_constant_mismatch weight_beyond_training_max_nn
#> 1 0.05 0 0.05
#> training_median_nn training_max_nn
#> 1 0.7204005 0.7204005
omega <- br_omega(prediction, tpe, H = c(1, 3, 5, 10, 21))
omega$summary
#> H median_omega median_pred_var median_env_var
#> 1 1 10.3342773 0.08679504 0.89696401
#> 2 3 3.4447591 0.08679504 0.29898800
#> 3 5 2.0668555 0.08679504 0.17939280
#> 4 10 1.0334277 0.08679504 0.08969640
#> 5 21 0.4921084 0.08679504 0.04271257
Support is geometric: ranges and nearest-neighbour distances can flag
extrapolative states, but cannot establish biological validity or
calibration. Omega(H) compares finite-TPE environmental variance with
prediction variance. It is a relative diagnostic, not a probability or a
universal threshold.
Compare M0–M3 and calculate advancement probabilities
comparison <- br_compare(prediction, tpe, H = 5, k = 3,
nsim = 2000, seed = 410)
comparison$summary
#> method H k nsim expected_merit_set_reliability probability_set_reliability
#> 1 M0 5 3 1 NA NA
#> 2 M1 5 3 2000 0.6743333 0.6743333
#> 3 M2 5 3 2000 0.8573333 0.8573333
#> 4 M3 5 3 2000 0.6473333 0.6473333
#> probability_sum invariant_error
#> 1 3 0
#> 2 3 0
#> 3 3 0
#> 4 3 0
head(comparison$candidates)
#> method H candidate expected_merit p_advance mcse in_expected_merit_set
#> 1 M0 5 MZ02 8.8450 1 0 TRUE
#> 2 M0 5 MZ01 8.8115 1 0 TRUE
#> 3 M0 5 MZ03 8.8005 1 0 TRUE
#> 4 M0 5 MZ05 8.6040 0 0 FALSE
#> 5 M0 5 MZ04 8.5085 0 0 FALSE
#> 6 M0 5 MZ07 8.4520 0 0 FALSE
#> in_probability_set
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 FALSE
#> 5 FALSE
#> 6 FALSE
M0 is deterministic point ranking and does not estimate reliability.
M1 propagates prediction uncertainty, M2 propagates finite-TPE
uncertainty while holding coefficients fixed, and M3 propagates both.
The results are conditional on the supplied predictive distribution,
TPE, H, and fixed-k rule.
For a selected set, reliability is the expected retained fraction
k^-1 * sum(p_i^A); it is not the probability of exact-set
recovery. The candidate-level Monte Carlo standard error describes
simulation error only. The probability-selected set is selected and
evaluated using the same draws, so its reliability estimate can be
optimistic at finite simulation size.
Interpret the H curve
curve <- br_curve(prediction, tpe, H = c(1, 3, 5, 10), k = 3,
nsim = 2000, seed = 410)
curve$summary
#> method H k nsim expected_merit_set_reliability probability_set_reliability
#> 1 M0 1 3 1 NA NA
#> 2 M1 1 3 2000 0.6743333 0.6743333
#> 3 M2 1 3 2000 0.7106667 0.7106667
#> 4 M3 1 3 2000 0.5945000 0.5945000
#> 5 M0 3 3 1 NA NA
#> 6 M1 3 3 2000 0.6743333 0.6743333
#> 7 M2 3 3 2000 0.8041667 0.8041667
#> 8 M3 3 3 2000 0.6383333 0.6383333
#> 9 M0 5 3 1 NA NA
#> 10 M1 5 3 2000 0.6743333 0.6743333
#> 11 M2 5 3 2000 0.8573333 0.8573333
#> 12 M3 5 3 2000 0.6473333 0.6473333
#> 13 M0 10 3 1 NA NA
#> 14 M1 10 3 2000 0.6743333 0.6743333
#> 15 M2 10 3 2000 0.9161667 0.9161667
#> 16 M3 10 3 2000 0.6561667 0.6561667
#> 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
plot(curve, type = "reliability")
plot(curve, type = "omega")
H is the number of future independent TPE draws under the reference
formulation. br_curve() computes each requested H
separately and does not impose monotonicity on estimated
reliability.
External validation and empirical findings
For external validation, freeze the upstream fit, TPE, predictions,
probabilities, and selected sets before inspecting outcomes. Construct
an independent phenotype reference, then score candidate-level events
with a candidate-level prevalence baseline. Keep uncertainty in the
external ranking separate from Monte Carlo error in predicted
probabilities.
The project’s prediction-frozen 2024 validation and forward
year-held-out 2022 analysis found poor absolute calibration; M2 had the
worst Brier performance; M3’s Brier loss was slightly below M1’s in both
settings; and predicted selected-set reliability was substantially
overoptimistic. All probability methods performed worse than
prevalence-only forecasts, including after the corrected Stage L3a
candidate-level denominator audit. Performance alignment was stronger
descriptively for candidates with historical phenotype records, without
supporting a causal interpretation. The 2022 analysis uses a later
authoritative G2F release for a year-held-out analysis; it is not an
exact reconstruction of the separate official 2022 competition test
design.
These empirical results are not a claim of general calibration.
Correctly propagating the specified uncertainties does not repair an
inadequate upstream predictive distribution or establish improved
genetic gain. Full provenance, scripts, and the preserved
pre-correction/corrected L3a audit trail are kept in the separate
reproducibility archive; raw G2F data are not included here.