From imported predictions to advancement
reliability: simulated maize
Decision and scope
A hypothetical maize programme must advance three of eight
unphenotyped candidates for grain yield (t/ha). Every value here is
simulated. The 12 training environments and replicated
training phenotypes illustrate a separate training population; candidate
predictions and their covariance are stipulated inputs, not a model fit
to those phenotypes. They cannot demonstrate predictive accuracy or
calibration. The data generation script is supplied in
inst/scripts/maize_demo.R in the source and is installed as
system.file("scripts", "maize_demo.R", package = "BreedReliabR").
Water deficit is a dimensionless fraction during flowering, and heat
is the number of hot days during flowering in this hypothetical
scenario. These descriptors and linear yield responses are illustrative
simplifications, not validated crop models.
1. Fit-import boundary
In a real application, fit and validate a prediction model
externally, then export coefficient means and the full joint
prediction-error covariance, including between-candidate
dependence, using the same environmental basis. Do not substitute
residual variance, omit cross-candidate covariance without
justification, or mix coefficient orders. No genomic model is fitted by
this package.
data(maize_demo)
d <- maize_demo
head(d$training_environments)
#> environment flowering_water_deficit hot_days_flowering
#> 1 TR01 0.10 2
#> 2 TR02 0.15 5
#> 3 TR03 0.20 1
#> 4 TR04 0.25 7
#> 5 TR05 0.30 3
#> 6 TR06 0.35 9
head(d$training_phenotypes)
#> genotype environment replicate yield_t_ha
#> 1 TG01 TR01 1 10.618130
#> 2 TG02 TR01 1 10.668364
#> 3 TG03 TR01 1 9.507768
#> 4 TG04 TR01 1 11.365346
#> 5 TG05 TR01 1 10.439424
#> 6 TG06 TR01 1 9.860005
d$candidates
#> candidate phenotyped
#> 1 MZ01 FALSE
#> 2 MZ02 FALSE
#> 3 MZ03 FALSE
#> 4 MZ04 FALSE
#> 5 MZ05 FALSE
#> 6 MZ06 FALSE
#> 7 MZ07 FALSE
#> 8 MZ08 FALSE
pred <- do.call(br_prediction, d$prediction_inputs)
pred
#> BreedReliabR prediction object
#> Candidates: 8
#> Environmental basis terms: 3
#> Source: Simulated fit-import inputs; not estimated from training phenotypes
Rows of coef_mean follow candidate IDs. Covariance rows
and columns follow all three coefficients of candidate 1, then all
coefficients of candidate 2, and so on. The basis is intercept,
(water_deficit - 0.35)/0.20, and
(hot_days - 7)/5. Apply these exact transformations to
training and target environments; never standardize them separately. The
uncertainty distribution is Gaussian and is assumed independent of
future environment sampling.
2. Specify the TPE
d$tpe_states
#> state_id flowering_water_deficit hot_days_flowering weight
#> 1 cool_wet 0.15 2 0.15
#> 2 typical 0.35 7 0.40
#> 3 dry 0.55 9 0.25
#> 4 hot_dry 0.60 14 0.15
#> 5 extreme 0.80 20 0.05
tpe <- br_tpe(d$tpe_basis, weights = d$tpe_states$weight,
state_id = d$tpe_states$state_id)
tpe
#> BreedReliabR target population of environments
#> States: 5
#> Basis terms: 3
#> Weight sum: 1
These fixed weights describe a hypothetical deployment mixture, with
5% mass on an extreme hot-dry state. They are assumptions, not estimates
of biological truth. Future environments are sampled independently with
replacement. H is the number of such future draws, not the number of
training trials or an effective sample size for a correlated network.
Dependence among locations/years requires a different sampling model and
is outside this demonstration.
3. Diagnose support before interpreting reliability
support <- br_support(tpe, d$training_basis)
support
#> BreedReliabR environmental support diagnostics
#> weight_outside_range weight_constant_mismatch weight_beyond_training_max_nn
#> 0.05 0 0.05
#> training_median_nn training_max_nn
#> 0.7204005 0.7204005
#> Geometric support only; these diagnostics do not establish biological correctness of the specified TPE or prediction calibration.
support$states
#> state_id weight outside_range n_outside constant_mismatch
#> cool_wet cool_wet 0.15 FALSE 0 FALSE
#> typical typical 0.40 FALSE 0 FALSE
#> dry dry 0.25 FALSE 0 FALSE
#> hot_dry hot_dry 0.15 FALSE 0 FALSE
#> extreme extreme 0.05 TRUE 2 FALSE
#> nearest_distance beyond_training_max_nn
#> cool_wet 0.2773501 FALSE
#> typical 0.4596570 FALSE
#> dry 0.2773501 FALSE
#> hot_dry 0.3602002 FALSE
#> extreme 1.4187445 TRUE
support$ranges
#> term min max sd constant
#> 1 intercept 1.00 1.0 0.0000000 TRUE
#> 2 water -1.25 1.5 0.9013878 FALSE
#> 3 heat -1.20 1.6 0.8702142 FALSE
plot(support)
The extreme state exceeds the training ranges. Distances use training
sample standard deviations in nonconstant basis columns; the intercept
contributes nothing. A change in a constant training column instead has
infinite distance. The dashed line is the largest leave-one-out
nearest-neighbour training distance, an empirical reference, not a
calibrated cutoff. Correlated/redundant descriptors can distort
Euclidean distances. Being inside each marginal range does not imply
being inside the joint training support, and proximity does not
establish model validity. Zero-weight states are still diagnosed but
contribute no weighted mass.
4. Partition uncertainty with Omega
om <- br_omega(pred, tpe, H = c(1, 3, 5, 10, 20))
om
#> BreedReliabR uncertainty diagnostic
#> H median_omega median_pred_var median_env_var
#> 1 10.3342773 0.08679504 0.8969640
#> 3 3.4447591 0.08679504 0.2989880
#> 5 2.0668555 0.08679504 0.1793928
#> 10 1.0334277 0.08679504 0.0896964
#> 20 0.5167139 0.08679504 0.0448482
head(om$candidate)
#> candidate H expected_merit pred_var env_var omega
#> 1 MZ01 1 8.8115 0.08679504 2.6023033 29.982166
#> 2 MZ01 3 8.8115 0.08679504 0.8674344 9.994055
#> 3 MZ01 5 8.8115 0.08679504 0.5204607 5.996433
#> 4 MZ01 10 8.8115 0.08679504 0.2602303 2.998217
#> 5 MZ01 20 8.8115 0.08679504 0.1301152 1.499108
#> 6 MZ02 1 8.8450 0.08679504 1.4280468 16.453092
Omega is the candidate-specific ratio of finite-environment variance
to prediction variance. Its numerator includes the
prediction-by-environment interaction in the core variance
decomposition. With fixed inputs it decreases as 1/H. It is not a
probability or a universal decision threshold. Absolute candidate
variances do not alone determine ranking reliability, which also depends
on merit gaps and joint covariance. Zero prediction variance can produce
infinite Omega; both variance components zero give Omega zero by
convention.
5. Compare all four uncertainty scenarios
comparison <- br_compare(pred, tpe, H = 5, k = 3, nsim = 2000, seed = 410)
summary(comparison)
#> 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
plot(comparison)
M0 ranks expected merit; M1 varies predictions at the mean TPE basis;
M2 varies future environments with fixed coefficient means; M3 varies
both. They use the same inputs and selection size. Seeded calls are
reproducible, but reusing a seed does not provide a paired estimator for
differences among methods. M0 reports NA set reliability because it does
not estimate uncertainty.
6. Inspect advancement probabilities and selected sets
joint <- comparison$analyses$M3
joint$candidates
#> candidate expected_merit p_advance mcse in_expected_merit_set
#> 1 MZ02 8.8450 0.6805 0.010426403 TRUE
#> 2 MZ01 8.8115 0.6165 0.010872620 TRUE
#> 3 MZ03 8.8005 0.6450 0.010699883 TRUE
#> 4 MZ05 8.6040 0.3515 0.010675855 FALSE
#> 5 MZ04 8.5085 0.2460 0.009630265 FALSE
#> 6 MZ07 8.4520 0.2340 0.009466890 FALSE
#> 7 MZ06 8.4130 0.1570 0.008134833 FALSE
#> 8 MZ08 8.2275 0.0695 0.005686376 FALSE
#> in_probability_set
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 FALSE
#> 5 FALSE
#> 6 FALSE
#> 7 FALSE
#> 8 FALSE
joint$expected_merit_set
#> [1] "MZ02" "MZ01" "MZ03"
joint$probability_set
#> [1] "MZ02" "MZ03" "MZ01"
stopifnot(abs(sum(joint$candidates$p_advance) - 3) < 1e-10)
Each draw advances exactly k candidates, so probabilities sum to k,
not one. Candidate MCSE is a binomial Monte Carlo standard error, not a
confidence interval for real-world performance. Zero observed MCSE at an
estimated probability of zero or one is not proof of certainty. Increase
nsim and repeat seeds when borderline decisions matter. Ties are
resolved deterministically by candidate ID.
Set reliability is the expected retained fraction:
the average advancement probability over the specified k-member set. It
is not the probability that every member is retained. The
probability-selected set is chosen using the same draws that estimate
its reliability, so its reported reliability can be optimistically
selected at finite nsim. This wrapper does not claim independent
evaluation or provide a set-level MCSE from marginal candidate
MCSEs.
7. Interpret reliability versus H
curve <- br_curve(pred, tpe, H = c(1, 3, 5, 10, 20), k = 3,
nsim = 2000, seed = 410)
summary(curve)
#> $reliability
#> 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
#> 17 M0 20 3 1 NA NA
#> 18 M1 20 3 2000 0.6743333 0.6743333
#> 19 M2 20 3 2000 0.9490000 0.9490000
#> 20 M3 20 3 2000 0.6681667 0.6681667
#> 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
#> 17 3 0
#> 18 3 0
#> 19 3 0
#> 20 3 0
#>
#> $omega
#> H median_omega median_pred_var median_env_var
#> 1 1 10.3342773 0.08679504 0.8969640
#> 2 3 3.4447591 0.08679504 0.2989880
#> 3 5 2.0668555 0.08679504 0.1793928
#> 4 10 1.0334277 0.08679504 0.0896964
#> 5 20 0.5167139 0.08679504 0.0448482
plot(curve)
plot(curve, type = "omega")
M1 is unchanged by H for a fixed seed. As H increases, environmental
sampling variance decreases; M3 approaches the prediction-only scenario
in the limit, but neither estimated nor exact ranking reliability must
increase monotonically in every configuration. Curves are computed at
the requested H values without smoothing or choosing a supposedly
adequate H. M0 and M1 still appear in the returned tables. Results at
different H values reuse the seed and are not independent; no
significance test for differences is supplied.
8. TPE sensitivity and research use
# A transparent alternative scenario, not an automatic correction of the TPE.
w <- d$tpe_states$weight
w[d$tpe_states$state_id == "extreme"] <- 0
tpe_alternative <- br_tpe(d$tpe_basis, w, d$tpe_states$state_id)
br_support(tpe_alternative, d$training_basis)
#> BreedReliabR environmental support diagnostics
#> weight_outside_range weight_constant_mismatch weight_beyond_training_max_nn
#> 0 0 0
#> training_median_nn training_max_nn
#> 0.7204005 0.7204005
#> Geometric support only; these diagnostics do not establish biological correctness of the specified TPE or prediction calibration.
summary(br_compare(pred, tpe_alternative, H = 5, k = 3,
nsim = 2000, seed = 410))
#> method H k nsim expected_merit_set_reliability probability_set_reliability
#> 1 M0 5 3 1 NA NA
#> 2 M1 5 3 2000 0.7113333 0.7113333
#> 3 M2 5 3 2000 0.9220000 0.9220000
#> 4 M3 5 3 2000 0.6840000 0.6840000
#> probability_sum invariant_error
#> 1 3 0
#> 2 3 0
#> 3 3 0
#> 4 3 0
Removing an unsupported state changes the scientific question; it is
justified only by deployment evidence, not by a desire to improve
apparent reliability. For a publication with real data, document data
provenance and environmental descriptors, validate the upstream
predictions out of sample, audit covariance ordering and calibration,
justify TPE weights, examine alternative plausible mixtures, assess
Monte Carlo stability, and compare predicted decisions against
independent outcomes. Passing software checks and this simulated
workflow do not establish empirical breeding utility or publication
readiness of those claims.