Genomic Prediction

This document introduces genomic prediction in Mandala through a genomic relationship matrix (GRM). The central idea is that mandala_gp() is a small convenience wrapper around mandala() for genomic prediction models.

Under the hood, the model still uses the ordinary Mandala syntax:

  1. prepare or provide a GRM;
  2. include GM(geno, GRM) in the random model;
  3. predict tested and untested genotypes from the fitted genomic model.

The examples use the same small MET phenotype dataset as the introduction document and add the shared example GRM sim_GRM_1000.rds.

library(mandala)

Data

The phenotype data are plot-level MET observations. The GRM contains the larger simulated genotype pool; mandala_grm_prep() automatically aligns it to the genotypes observed in the data.

pheno_path <- system.file("extdata", "fullrep_MET_n1000.csv", package = "mandala", mustWork = TRUE)
grm_path <- system.file("extdata", "sim_GRM_1000.rds", package = "mandala", mustWork = TRUE)

met <- read.csv(pheno_path, stringsAsFactors = FALSE, check.names = FALSE)
for (v in c("geno", "env", "rep", "block", "row", "col")) met[[v]] <- factor(met[[v]])
met$yld <- as.numeric(met$yld)

GRM <- readRDS(grm_path)

c(
  n_records = nrow(met),
  n_genotypes = length(unique(met$geno)),
  n_environments = length(unique(met$env)),
  grm_dimension = nrow(GRM)
)
#>      n_records    n_genotypes n_environments  grm_dimension 
#>           1000             50             10           1000

Genomic Prediction with mandala()

This model uses genomic relationships for the genotype main effect, plus a homogeneous genotype-by-environment term and block-within-replication term. The only extra preparation is aligning the GRM to the phenotype data with mandala_grm_prep().

G_list <- mandala_grm_prep(
  GRM = GRM,
  data = met,
  geno_col = "geno",
  verbose = FALSE
)

fit_gm <- mandala(
  fixed = yld ~ env,
  random = ~ GM(geno, GRM) + geno:env + env:rep:block,
  data = met,
  matrix_list = G_list,
  verbose = FALSE
)

summary(fit_gm)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~GM(geno, GRM) + geno:env + 
#>     env:rep:block, data = met, matrix_list = G_list, verbose = FALSE)
#> 
#> Variance Components:
#>      component   estimate  std.error   z.ratio bound %ch
#>   GM(geno,GRM) 0.54540814 0.12620151  4.321724     P  NA
#>       geno:env 0.07574029 0.02727827  2.776580     P  NA
#>  env:rep:block 0.08046126 0.02323121  3.463499     P  NA
#>       R.sigma2 0.48906366 0.03213805 15.217590     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0539983 0.4258882 16.563030
#>     envY1_L2 -0.2672079 0.2010568 -1.329017
#>     envY1_L3 -0.7769428 0.2010437 -3.864547
#>     envY1_L4 -1.1048010 0.2010409 -5.495405
#>     envY1_L5 -0.4548158 0.2010549 -2.262148
#> 
#> Converged: TRUE  |  Iterations: 10
#> 
#> Model Notes:
#> - Model size class: complex (n=1000, q=610, MME dim=620).
#> - Selected prediction SEs are available through mandala_predict().
#> - Selected fixed-effect tests are available with mandala_fixed_tests(type = 'selected').
#> 
#> Random Effects (BLUPs) [first 5]:
#>        random level   estimate std.error    z.ratio
#>  GM(geno,GRM)    G1 -0.8081274 0.4365360 -1.8512275
#>  GM(geno,GRM)   G10  0.1119708 0.4362316  0.2566774
#>  GM(geno,GRM)   G11  1.3988170 0.4336945  3.2253512
#>  GM(geno,GRM)   G12  0.3443569 0.4358125  0.7901494
#>  GM(geno,GRM)   G13  1.1030691 0.4365890  2.5265616
#> 
#> logLik: -1225.159   AIC: 2458.318   BIC: 2477.909   logLik_Trunc: -315.410
pred_gm <- mandala_predict(fit_gm, "geno", verbose = FALSE)
head(pred_gm)
#>   geno predicted_value std_error
#> 1   G1        5.786939 0.1562591
#> 2  G10        6.740107 0.1556430
#> 3  G11        8.031981 0.1553483
#> 4  G12        6.950213 0.1552926
#> 5  G13        7.736765 0.1553067
#> 6  G14        6.554666 0.1561896

The genomic relationship can also be used for genotype-by-environment effects. This keeps the same GRM, but allows the genomic genotype signal to vary by environment.

fit_gm_env <- mandala(
  fixed = yld ~ env,
  random = ~ GM(geno, GRM) + GM(geno, GRM):env + env:rep:block,
  data = met,
  matrix_list = G_list,
  verbose = FALSE
)

summary(fit_gm_env)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~GM(geno, GRM) + GM(geno, 
#>     GRM):env + env:rep:block, data = met, matrix_list = G_list, 
#>     verbose = FALSE)
#> 
#> Variance Components:
#>         component   estimate  std.error    z.ratio bound %ch
#>      GM(geno,GRM) 0.55363275 0.12613741  4.3891241     P  NA
#>  GM(geno,GRM):env 0.01544381 0.02496859  0.6185297     P  NA
#>     env:rep:block 0.07782892 0.02297782  3.3871322     P  NA
#>          R.sigma2 0.55102082 0.03064919 17.9783183     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0538731 0.4325091 16.309192
#>     envY1_L2 -0.2674258 0.2169121 -1.232877
#>     envY1_L3 -0.7784810 0.2169276 -3.588667
#>     envY1_L4 -1.1101087 0.2168907 -5.118287
#>     envY1_L5 -0.4530669 0.2169387 -2.088455
#> 
#> Converged: TRUE  |  Iterations: 15
#> 
#> Model Notes:
#> - Model size class: complex (n=1000, q=610, MME dim=620).
#> - Selected prediction SEs are available through mandala_predict().
#> - Selected fixed-effect tests are available with mandala_fixed_tests(type = 'selected').
#> 
#> Random Effects (BLUPs) [first 5]:
#>        random level   estimate std.error    z.ratio
#>  GM(geno,GRM)    G1 -0.8204456 0.4363476 -1.8802572
#>  GM(geno,GRM)   G10  0.1187481 0.4360939  0.2722995
#>  GM(geno,GRM)   G11  1.4070947 0.4338527  3.2432544
#>  GM(geno,GRM)   G12  0.3427839 0.4357695  0.7866174
#>  GM(geno,GRM)   G13  1.1098259 0.4364918  2.5426042
#> 
#> logLik: -1229.034   AIC: 2466.069   BIC: 2485.660   logLik_Trunc: -319.285
pred_gm_env <- mandala_predict(fit_gm_env, "geno", verbose = FALSE)
head(pred_gm_env)
#>   geno predicted_value std_error
#> 1   G1        5.789891 0.1639351
#> 2  G10        6.731705 0.1631657
#> 3  G11        8.023645 0.1628169
#> 4  G12        6.956366 0.1627007
#> 5  G13        7.725547 0.1627762
#> 6  G14        6.557935 0.1638938

Genomic Prediction with mandala_gp()

mandala_gp() is a convenience wrapper for the same model family. It can use the same prepared GRM list and stores the genomic information needed by helper functions for tested, untested, and cross-validation predictions.

gp_fit <- mandala_gp(
  data = met,
  fixed = yld ~ env,
  geno_col = "geno",
  random = ~ GM(geno, GRM) + geno:env + env:rep:block,
  matrix_list = G_list,
  verbose = FALSE
)

summary(gp_fit$fit)
#> Model statement:
#> mandala(fixed = yld ~ env, random = ~GM(geno, GRM) + geno:env + 
#>     env:rep:block, data = met, matrix_list = `<prepared matrix_list>`, 
#>     method = "sparse", engine = "auto")
#> 
#> Variance Components:
#>      component   estimate  std.error   z.ratio bound %ch
#>   GM(geno,GRM) 0.54540814 0.12620151  4.321724     P  NA
#>       geno:env 0.07574029 0.02727827  2.776580     P  NA
#>  env:rep:block 0.08046126 0.02323121  3.463499     P  NA
#>       R.sigma2 0.48906366 0.03213805 15.217590     P  NA
#> 
#> Fixed Effects (BLUEs) [first 5]:
#>       effect   estimate std.error   z.ratio
#>  (Intercept)  7.0539983 0.4258882 16.563030
#>     envY1_L2 -0.2672079 0.2010568 -1.329017
#>     envY1_L3 -0.7769428 0.2010437 -3.864547
#>     envY1_L4 -1.1048010 0.2010409 -5.495405
#>     envY1_L5 -0.4548158 0.2010549 -2.262148
#> 
#> Converged: TRUE  |  Iterations: 10
#> 
#> Model Notes:
#> - Model size class: complex (n=1000, q=610, MME dim=620).
#> - Selected prediction SEs are available through mandala_predict().
#> - Selected fixed-effect tests are available with mandala_fixed_tests(type = 'selected').
#> 
#> Random Effects (BLUPs) [first 5]:
#>        random level   estimate std.error    z.ratio
#>  GM(geno,GRM)    G1 -0.8081274 0.4365360 -1.8512275
#>  GM(geno,GRM)   G10  0.1119708 0.4362316  0.2566774
#>  GM(geno,GRM)   G11  1.3988170 0.4336945  3.2253512
#>  GM(geno,GRM)   G12  0.3443569 0.4358125  0.7901494
#>  GM(geno,GRM)   G13  1.1030691 0.4365890  2.5265616
#> 
#> logLik: -1225.159   AIC: 2458.318   BIC: 2477.909   logLik_Trunc: -315.410

Predict Tested Genotype Means

pred_g <- mandala_gp_predict(gp_fit, "geno", verbose = FALSE)
head(pred_g)
#>   geno predicted_value std_error
#> 1   G1        5.786939 0.1562591
#> 2  G10        6.740107 0.1556430
#> 3  G11        8.031981 0.1553483
#> 4  G12        6.950213 0.1552926
#> 5  G13        7.736765 0.1553067
#> 6  G14        6.554666 0.1561896

pred_ge <- mandala_gp_predict(gp_fit, "geno:env", verbose = FALSE)
head(pred_ge)
#>   geno   env predicted_value std_error
#> 1   G1 Y1_L1        6.050564 0.3047670
#> 2   G1 Y1_L2        5.734578 0.3047672
#> 3   G1 Y1_L3        5.479284 0.3039156
#> 4   G1 Y1_L4        5.277320 0.3047260
#> 5   G1 Y1_L5        5.768992 0.3038358
#> 6   G1 Y2_L1        6.130175 0.3038194

Predict Tested and Untested Genotypes

The example GRM contains many genotypes with no phenotype records in this small vignette dataset. With include_untested = TRUE, Mandala returns the tested predictions from the fitted model and projected predictions for untested genotypes present in the full GRM.

pred_all <- mandala_gp_predict(
  gp_fit,
  classify_term = "geno",
  include_untested = TRUE,
  GRM = GRM,
  verbose = FALSE
)

summary(pred_all)
#> 
#> Genomic prediction carried out successfully
#> -------------------------------------------
#>                                                item value
#>                           Total genotypes predicted  1000
#>       Tested genotypes with phenotypic observations    50
#>  Untested genotypes without phenotypic observations   950
head(pred_all[, c("geno", "predicted_value", "std_error", "gebv", "source")])
#>   geno predicted_value std_error        gebv source
#> 1   G1        5.786939 0.1562591 -0.80812741 tested
#> 2  G10        6.740107 0.1556430  0.11197081 tested
#> 3  G11        8.031981 0.1553483  1.39881704 tested
#> 4  G12        6.950213 0.1552926  0.34435694 tested
#> 5  G13        7.736765 0.1553067  1.10306909 tested
#> 6  G14        6.554666 0.1561896 -0.04512281 tested
tail(pred_all[, c("geno", "predicted_value", "std_error", "gebv", "source")])
#>       geno predicted_value std_error        gebv   source
#> 995   G995        6.555601        NA -0.09479182 untested
#> 996   G996        6.559796        NA -0.09059622 untested
#> 997   G997        6.618727        NA -0.03166578 untested
#> 998   G998        6.615381        NA -0.03501141 untested
#> 999   G999        6.594937        NA -0.05545487 untested
#> 1000 G1000        6.671627        NA  0.02123497 untested

For advanced workflows, mandala_gp_predict_untested() remains available when only the explicit GRM projection table is required.

Genotype-Holdout Cross-Validation

This quick example masks whole genotypes, fits the model on the remaining genotypes, and predicts the held-out genotypes from the GRM.

set.seed(2026)

gp_cv <- mandala_gp_cv(
  data = met,
  fixed = yld ~ env,
  geno_col = "geno",
  random = ~ GM(geno, GRM) + geno:env + env:rep:block,
  GRM = GRM,
  k = 5,
  repeats = 1,
  verbose = FALSE
)

summary(gp_cv)
#> 
#> Mandala GP cross-validation summary
#> -----------------------------------
#>  repeats k  fold_by mean_pred_ability sd_pred_ability mean_rmse mean_mae n_folds_ok
#>        1 5 genotype         0.7259711       0.1239532  6.591855 6.576399          5
gp_cv$by_fold
#>   repeat_id fold n_train n_test n_eval pred_ability     rmse      mae status
#> 1         1    1      40     10     10    0.7105344 6.469911 6.459938     ok
#> 2         1    2      40     10     10    0.5576787 6.440720 6.428987     ok
#> 3         1    3      40     10     10    0.6757847 6.677241 6.661173     ok
#> 4         1    4      40     10     10    0.8817879 6.575837 6.557245     ok
#> 5         1    5      40     10     10    0.8040697 6.795567 6.774654     ok
plot(
  gp_cv$by_fold$fold,
  gp_cv$by_fold$pred_ability,
  pch = 16,
  xlab = "Fold",
  ylab = "Prediction ability",
  main = "Genotype-holdout GP cross-validation"
)
abline(h = gp_cv$summary$mean_pred_ability, col = "gray40", lty = 2)
legend(
  "bottomright",
  legend = paste0("mean = ", round(gp_cv$summary$mean_pred_ability, 3)),
  bty = "n"
)

Relationship-Matrix Diagnostics

mandala_relationship_diagnostics(GRM[levels(met$geno), levels(met$geno)])
#> 
#> Mandala relationship-matrix diagnostics
#>   Matrix             : Relationship matrix 
#>   Size               : 50 x 50 
#>   Diagonal           : min 0.97275 | mean 1.00912 | max 1.05171 
#>   Off-diagonal       : min 0.09397 | mean 0.28709 | max 0.7491 
#>   Diagonal flags     : 0 outside [0.8, 1.2]
#>   Duplicate flags    : 0 pairs with standardized relationship > 0.95

Selected Literature

Genomic prediction uses dense marker information to predict genetic merit through genomic relationships or marker effects. The examples above emphasize GBLUP-style modeling through a GRM, which is closely related to ridge-regression and genomic best linear unbiased prediction. Genotype-holdout cross-validation is commonly used to estimate empirical prediction ability.