--- title: "From imported predictions to advancement reliability: simulated maize" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{From imported predictions to advancement reliability: simulated maize} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4) library(BreedReliabR) ``` ## 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. ```{r import} data(maize_demo) d <- maize_demo head(d$training_environments) head(d$training_phenotypes) d$candidates pred <- do.call(br_prediction, d$prediction_inputs) pred ``` 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 ```{r tpe} d$tpe_states tpe <- br_tpe(d$tpe_basis, weights = d$tpe_states$weight, state_id = d$tpe_states$state_id) tpe ``` 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 ```{r support} support <- br_support(tpe, d$training_basis) support support$states support$ranges 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 ```{r omega} om <- br_omega(pred, tpe, H = c(1, 3, 5, 10, 20)) om head(om$candidate) ``` 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 ```{r comparison} comparison <- br_compare(pred, tpe, H = 5, k = 3, nsim = 2000, seed = 410) summary(comparison) 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 ```{r advancement} joint <- comparison$analyses$M3 joint$candidates joint$expected_merit_set joint$probability_set 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 ```{r curve} curve <- br_curve(pred, tpe, H = c(1, 3, 5, 10, 20), k = 3, nsim = 2000, seed = 410) summary(curve) 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 ```{r sensitivity} # 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) summary(br_compare(pred, tpe_alternative, H = 5, k = 3, nsim = 2000, seed = 410)) ``` 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.