Genomic-assisted breeding has transitioned from theoretical concepts to practical applications in breeding. Genomic selection (GS) predicts genomic breeding values (GEBV) using dense genetic markers. Single-stage models predict GEBVs from phenotypic observations in one step, fully accounting for the entire variance-covariance structure among genotypes, but face computational challenges. Two-stage models, preferred for their simplicity and efficiency, first calculate adjusted genotypic means accounting for spatial variation within each environment, then use these means to predict GEBVs. However, unweighted (UNW) two-stage models assume independent errors among adjusted means, neglecting correlations among estimation errors. Here, we show that fully-efficient two-stage models perform similarly to UNW models for randomized complete block designs but substantially better for augmented designs. Our simulation studies demonstrate the impact of the fully-efficient methodology on prediction accuracy across different implementations and scenarios. Incorporating non-additive effects and augmented designs significantly improved accuracy, emphasizing the synergy between design and model strategy. Consistent performance requires the estimation error covariance to be incorporated into a random effect (Full_R model) rather than into the residuals. Our results suggest that the fully-efficient methodology, particularly the Full_R model, should be more prevalent, especially as GS increases the appeal of sparse designs. We also provide a comprehensive theoretical background and open-source R code, enhancing understanding and facilitating broader adoption of fully-efficient two-stage models in GS. Here, we offer insights into the practical applications of fully-efficient models and their potential to increase genetic gain, demonstrating a 13.80 % improvement after five selection cycles when moving from UNW to Full_R models.