Wheat breeding

library(makemyprior)

Model

We use a model from Hem et al. (2021); a genomic wheat breeding model with three genetic effects: additive, dominance and epistasis The model is: yi = μ + ai + di + xi + εi, i = 1, …, 100,

where * μ is an intercept with a 𝒩(0, 10002) prior, * a = (a1, …, a100) ∼ 𝒩100(0, σa2A) is the additive effect, * d = (d1, …, d100) ∼ 𝒩100(0, σa2D) is the dominance effect, * x = (x1, …, x100) ∼ 𝒩100(0, σa2X) is the epistasis effect, and * εi is the residual effect with variance σε2.

The covariance matrices A, D and X are computed from the single nucleotide polymorphism (SNP) matrix with thousands of genetic markers.

We use the dataset wheat_breeding in makemyprior, which consists of indexes for the effects and precision matrices for each of the genetic effects, and present three priors.

We scale the precision matrices so the corresponding covariance matrices have typical variance equal to 1 with scale_precmat.

wheat_data_scaled <- wheat_data
wheat_data_scaled$Q_a <- scale_precmat(wheat_data$Q_a)
wheat_data_scaled$Q_d <- scale_precmat(wheat_data$Q_d)
wheat_data_scaled$Q_x <- scale_precmat(wheat_data$Q_x)

formula <- y ~
  mc(a, model = "generic0", Cmatrix = Q_a, constr = TRUE) +
  mc(d, model = "generic0", Cmatrix = Q_d, constr = TRUE) +
  mc(x, model = "generic0", Cmatrix = Q_x, constr = TRUE)

We do not carry out inference, as it takes time and will slow down the compilation of the vignettes by a lot, but include code so the user can run the inference themselves.

Prior 1

We do not want to say anything about the total (also called phenotypic) variance. We have expert knowledge saying that we want shrinkage towards the residual effect to avoid overfitting, and the heritability, which is the amount of phenotypic variance attributed to the genetic effects, is around 0.25. The additive, dominance and epistasis effects get about (85, 10, 5)% of the genetic variance according to the expert. The expert is quite sure about these choices, and we use a concentration parameter c, saying how much of the density mass of the prior such that Prob(logit(1/4) < logit(ω*) − logit(m) < logit(3/4)) = c for a variance proportion ω*.

This information can be implemented as a prior as follows:

prior1 <- make_prior(formula, wheat_data_scaled, prior = list(
  tree = "s1 = (d, x); s2 = (a, s1); s3 = (s2, eps)",
  w = list(s1 = list(prior = "pcM", param = c(0.67, 0.8)),
           s2 = list(prior = "pcM", param = c(0.85, 0.8)),
           s3 = list(prior = "pc0", param = 0.25))))

prior1
#> Model: y ~ mc(a, model = "generic0", Cmatrix = Q_a, constr = TRUE) + 
#>     mc(d, model = "generic0", Cmatrix = Q_d, constr = TRUE) + 
#>     mc(x, model = "generic0", Cmatrix = Q_x, constr = TRUE)
#> Tree structure: d_x = (d,x); a_d_x = (a,d_x); eps_a_d_x = (eps,a_d_x)
#> 
#> Weight priors:
#>  w[d/d_x] ~ PCM(0.67, 0.8)
#>  w[a/a_d_x] ~ PCM(0.85, 0.8)
#>  w[eps/eps_a_d_x] ~ PC1(0.75)
#> Total variance priors:
#>  V[eps_a_d_x] ~ Jeffreys'
plot_prior(prior1) # or plot(prior1)

plot_tree_structure(prior1)

Note that we do not fit the model in this vignette, as it takes some time.

posterior1 <- inference_stan(prior1, iter = 15000, warmup = 5000,
                            chains = 1, seed = 1)

plot_posterior_stan(posterior1, param = "prior", prior = TRUE)

For inference with INLA:

posterior1_inla <- inference_inla(prior1)
plot_posterior_stdev(posterior1_inla)

Prior 2

We have a good intuition on the absolute magnitude of the residual and genetic effects, but since the genetic effects are confounded, we want to use a prior saying that about 80% of the genetic variance is additive effect, and shrinkage towards the additive, and and be ignorant about the division of dominance and epistasis. The residual and genetic effects are both assumed to not be much larger than 3.

prior2 <- make_prior(formula, wheat_data_scaled, prior = list(
  tree = "s1 = (d, x); s2 = (a, s1); (eps)",
  w = list(s1 = list(prior = "dirichlet"),
           s2 = list(prior = "pc1", param = c(0.8))),
  V = list(s2 = list(prior = "pc", param = c(3, 0.05)),
           eps = list(prior = "pc", param = c(3, 0.05)))))

prior2
#> Model: y ~ mc(a, model = "generic0", Cmatrix = Q_a, constr = TRUE) + 
#>     mc(d, model = "generic0", Cmatrix = Q_d, constr = TRUE) + 
#>     mc(x, model = "generic0", Cmatrix = Q_x, constr = TRUE)
#> Tree structure: d_x = (d,x); a_d_x = (a,d_x); (eps)
#> 
#> Weight priors:
#>  w[d/d_x] ~ Dirichlet(2)
#>  w[a/a_d_x] ~ PC1(0.8)
#> Total variance priors:
#>  sqrt(V)[a_d_x] ~ PC0(3, 0.05)
#> Independent variance priors:
#>  sigma[eps] ~ PC0(3, 0.05)
plot_prior(prior2)

plot_tree_structure(prior2)

Again we do not fit the model in this vignette, as it takes some time.

posterior2 <- inference_stan(prior2, iter = 15000, warmup = 5000,
                            chains = 1, seed = 1)

plot_posterior_stan(posterior2, param = "prior", prior = TRUE)
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] makemyprior_1.2.2 rmarkdown_2.29   
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.9        utf8_1.2.4        generics_0.1.3    lattice_0.22-6   
#>  [5] digest_0.6.37     magrittr_2.0.3    evaluate_1.0.1    grid_4.4.2       
#>  [9] fastmap_1.2.0     jsonlite_1.8.9    Matrix_1.7-1      promises_1.3.0   
#> [13] fansi_1.0.6       scales_1.3.0      jquerylib_0.1.4   cli_3.6.3        
#> [17] shiny_1.9.1       rlang_1.1.4       visNetwork_2.1.2  munsell_0.5.1    
#> [21] splines_4.4.2     withr_3.0.2       cachem_1.1.0      yaml_2.3.10      
#> [25] tools_4.4.2       dplyr_1.1.4       colorspace_2.1-1  ggplot2_3.5.1    
#> [29] httpuv_1.6.15     buildtools_1.0.0  vctrs_0.6.5       R6_2.5.1         
#> [33] mime_0.12         lifecycle_1.0.4   htmlwidgets_1.6.4 MASS_7.3-61      
#> [37] shinyjs_2.1.0     pkgconfig_2.0.3   pillar_1.9.0      bslib_0.8.0      
#> [41] later_1.3.2       gtable_0.3.6      glue_1.8.0        Rcpp_1.0.13-1    
#> [45] xfun_0.49         tibble_3.2.1      tidyselect_1.2.1  sys_3.4.3        
#> [49] knitr_1.49        farver_2.1.2      xtable_1.8-4      htmltools_0.5.8.1
#> [53] maketools_1.3.1   labeling_0.4.3    compiler_4.4.2

References

Hem, Ingeborg Gullikstad, Maria Lie Selle, Gregor Gorjanc, Geir-Arne Fuglstad, and Andrea Riebler. 2021. “Robust modeling of additive and nonadditive variation with intuitive inclusion of expert knowledge.” Genetics, January.