## 3 Results and Discussion  
Figure 3 presents the distribution of ground truth maize and soybean yields for the used samples across the years 2019 to 2022. Maize yield data indicated relatively consistent distributions across the years, with median yields ranging between 13 and 14 t/ha, with a slight increase in median and minimum values in 2021 and 2022. Soybean yields had higher interannual variability in both median and spread than maize. Notably, yield data for 2021 had the highest median soybean yield, while 2019 and 2022 demonstrated broader distributions and a heavier lower tail, reflecting greater yield variability and a higher incidence of low-yielding observations. These patterns suggest heterogeneous interannual fluctuations between the two crops, potentially driven by crop-specific responses to climatic and management factors [73]. The results of the spatial autocorrelation test of input crop  
yield values per dataset using Moran's I are presented in Table 1, indicating a moderately low spatial autocorrelation for all the used datasets. While the evaluated machine learning algorithms efficiently model nonlinear relationships and can tolerate correlated predictors [74], they are not inherently resistant to the effects of spatial autocorrelation. While [PERSON] and [PERSON] [74] argued that RF, SVM, and artificial neural networks do not always produce higher prediction accuracy results with the presence of spatial autocorrelation in training and test datasets, prediction accuracy results in this study were likely affected by the moderately low spatial autocorrelation in all the datasets. Across all years and both crops, RF consistently outperformed other models in both cross-validation and final model fit, indicating its superior ability to capture nonlinear relationships in derived phenological indicators for yield prediction (Table 2). The optimal hyperparameters of all the used machine learning models are presented in Table A1. Most notably, the limited performance of linear and kernel-based models like MARS and SVM suggests the importance of nonlinearity and interactions in modeling crop yield potential. While RF's cross-validation R2 was moderately high (up to 0.409 in 2019), its final fit R2 of 0.898 suggests that the full training set provided a substantial benefit to model learning, possibly due to reduced variance and increased data, which was consistent across all four study years. As the SVM, MARS, and BRNNs did not produce high differences in the accuracy assessment metrics produced by cross-validation and final model fit, the data leakage caused by spatial autocorrelation in ground truth crop yield samples was unlikely to be the reason for this occurrence. Instead, the flexibility of RF and its tendency to model complex nonlinear relationships can increase the risk of overfitting, as a previous study noted that RF generally produced a much higher degree of overfitting in comparison to similar models [75]. Although repeated cross-validation provides a high level of resistance \begin{table}
\begin{tabular}{c c c c c} \hline \hline \multirow{2}{*}{**Year**} & \multicolumn{2}{c}{**Maze**} & \multicolumn{2}{c}{**Soybean**} \\  & **Moran’s I** & **_p_-**Value** & **Moran’s I** & **_p_-**Value** \\ \hline
2019 & 0.260 & \textless{}0.001 & 0.256 & \textless{}0.001 \\
2020 & 0.157 & \textless{}0.001 & 0.163 & \textless{}0.001 \\
2021 & 0.132 & \textless{}0.001 & 0.160 & \textless{}0.001 \\
2022 & 0.346 & \textless{}0.001 & 0.285 & \textless{}0.001 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: The results of the spatial autocorrelation test of the input crop yield values, per dataset, using Moran’s I.  
Figure 3: Violin plots representing the value distribution of the maize and soybean ground truth yield data used in the study. to overfitting in most cases, the large difference in R\({}^{2}\) values suggests the necessity for external validation using independent datasets in future studies.