## 5 Conclusions
\(\mathbf{r_{(0b)}}\) is the vector of random effects for rows within blocks within years, with \(\mathbf{r_{(0b)}}\)\(\sim\)\(N(\mathbf{0},\)\(\sigma_{r}^{2}\mathbf{D_{r}})\), \(\mathbf{c_{(0b)}}\) is the vector of random effects for columns within blocks within years, with \(\mathbf{c_{(0b)}}\)\(\sim\)\(N(\mathbf{0},\)\(\sigma_{c}^{2}\mathbf{D_{c}})\), and \(\mathbf{e}\) is the vector of residual errors, with \(\mathbf{e}\)\(\sim\)\(N(\mathbf{0},\)\(\sigma_{c}^{2}\mathbf{D_{c}})\) and with heterogeneity defined as \(\mathbf{\Sigma}\)\(=\)\(\sigma_{e}^{2}\)(\(\mathbf{AR_{\mathbf{1}}}\)\(\otimes\)\(\mathbf{AR_{\mathbf{1}}}\)). \(\mathbf{D_{r}}\), \(\mathbf{D_{c}}\), and \(\mathbf{D_{c}}\) are block diagonal matrices where each year level has its own \(\sigma_{r}^{2}\), \(\sigma_{c}^{2}\), and \(\sigma_{e}^{2}\), respectively. All other terms were previously described. A value of \(\rho_{\text{type-B}}\) close to one (obtained from the matrix \(\mathbf{U_{0}}\)) suggests minimal genotype-by-year (GxY) interaction. The incremental Wald statistic was used to assess the significance of the year as a covariate in the model. To estimate the type-A additive genetic correlations among traits [69], a bivariate, single-environment model was fitted for each pair of traits using raw data:  
\[\mathbf{y}=\mathbf{1}\mu_{m}+\mathbf{X_{1}}\mathbf{t}+\mathbf{X_{2}}\mathbf{\beta_{(0)}}+\mathbf{Z_{1}}\bm {u_{(0)}}+\mathbf{Z_{2}}\mathbf{r_{(0b)}}+\mathbf{Z_{3}}\mathbf{c_{(0b)}}+\mathbf{e}\] (A3)  
where \(\mathbf{y}\) represents the stacked binary vector of phenotypes for two traits \(\mathbf{y}=[\mathbf{y_{1}},\mathbf{y_{2}}]\), and \(\mu_{m}\) is the bivariate version of intercept (\(\mu\)). \(\mathbf{X_{2}}\mathbf{\beta_{(0)}}\), \(\mathbf{Z_{2}}\mathbf{r_{(0b)}}\), and \(\mathbf{Z_{3}}\mathbf{c_{(0b)}}\) are stacked design effects defined as previously in the model (A2), but with \"\(\mathbf{t}\)\" now denoting the trait, \(\mathbf{u_{(0)}}\) is a stacked vector of random genetic effects within trait, with \(\mathbf{u_{(0)}}\)\(\sim\)\(N(\mathbf{0},\)\(\sigma_{g}^{2}\mathbf{G}\otimes\mathbf{U_{0}})\), \(\mathbf{G}\) is previously described, \(\mathbf{U_{0}}\) is a \(2\times 2\) (traits) vcov matrix that is unstructured (corgh), allowing for trait-specific \(\sigma_{g}^{2}\) and unique additive genetic correlations across traits (\(\rho_{\text{type-A}}\)), \(\mathbf{e}\) is a stacked vector of random residual effects, with \(\mathbf{e}\)\(\sim\)\(N(\mathbf{0},\)\(\sigma_{e}^{2}\mathbf{I}\otimes\mathbf{R_{0}})\), and \(\mathbf{R_{0}}\) is the vcov between errors that is heterogeneous (corgh), introducing residual correlations between traits and unequal \(\sigma_{e}^{2}\) across traits. To aid convergence in REML for the bivariate model, the starting values for \(\sigma_{g}^{2}\) and \(\sigma_{e}^{2}\) were derived from the univariate model (A1) for each trait. Standard errors (approximated) for \(\rho_{\text{type-B}}\) and \(\rho_{\text{type-A}}\) correlations were estimated using the Delta method [71]. For phenolic acids, flavonoids, and lipid-derived products, models (A1), (A2) and (A3) were simplified by excluding design terms \(\mathbf{X}\mathbf{\beta}\), \(\mathbf{Z_{2}}\mathbf{r_{(0)}}\), and \(\mathbf{Z_{3}}\mathbf{c_{(0b)}}\) (block, row within block and column within block). In fact, NIRS prediction failure for these traits led to missing values and partially replicated data, preventing the robust estimation of block, row and column effects.