Machine Learning and Deep Learning Approaches for Crop Yield Prediction: A Comprehensive Analysis of Model Performance and Computational Efficiency

Unpredictable climate conditions and population growth have exacerbated the global challenge of ensuring food security [1-4]. In this study, yield (hg/ha) is considered the target variable, representing crop output per unit land area and distinguishing it from total production. We utilized temperature, rainfall, and pesticide usage data for crop yield prediction (CYP) across 101 countries over 21 years. Traditionally, CYP relied on statistical models that assume linearity and rigorous mathematical distributions. Consequently, they are often inadequate for capturing the complex, nonlinear, and dynamic interdependencies common in agronomic systems [5]. However, this scenario has dramatically shifted toward data-driven approaches, particularly machine learning (ML) and deep learning (DL). For example, Random Forest (RF), Support Vector Machines (SVM), and Gradient Boosting (GB) algorithms in ML have provided encouraging results in recent studies [6-8]. These methods can handle high-dimensional, heterogeneous data and capture nonlinear feature interactions. DL models, particularly recurrent architectures like Long Short-Term Memory (LSTM) networks, are designed to capture temporal dependencies in sequential data; however, their effectiveness in tabular CYP datasets remains uncertain. These capabilities are crucial for processing diverse environmental and temporal data [9-13]. Despite these advancements, several limitations remain. Only a limited number of studies have benchmarked a comprehensive set of algorithms, particularly with respect to computational resource requirements. Though DL models are increasingly applied, their performance is rarely benchmarked across different architectural recurrent neural network (RNN) variants (e.g., RNNs vs. LSTMs vs. attention-based models), and their computational implications are often overlooked.

To confront these gaps, we propose a multifaceted approach that includes carefully chosen preprocessing techniques, Exploratory Data Analysis (EDA) to gain more insights into crop yield, and systematic comparison of model efficiency. Unlike previous studies examining only a few models, we comprehensively evaluate 14 ML and 3 DL RNN variants for CYP using a multiyear, multiregional dataset. We selected these models to represent diverse learning paradigms, from traditional statistical methods to cutting-edge deep learning architectures. A key contribution of this study is the joint evaluation of predictive accuracy and computational efficiency under a unified experimental framework, ensuring a fair comparison across models.

Our study aims to address these challenges with the following objectives:

• To benchmark a broad set of ML and DL models for CYP, including ensemble methods and recurrent architectures.

• To evaluate models’ performance based on prediction metrics (MSE, MAE, R²) and deployment-oriented metrics (training time, inference time, and memory usage).

• To propose AgroStackNet and benchmark its performance against the other models.

The remainder of this paper is organized as follows: Section 2 presents the related work, Section 3 describes the materials and methods, Section 4 details the methodology and proposed AgroStackNet framework, Section 5 presents the results and discussion, and Section 6 concludes the study with future research directions.

3.1 Experimental setup

We conducted this study on a computing system with an Intel Core i9-12900 processor and an NVIDIA GeForce RTX 3060 GPU running on Windows 11. The details of the implementation platform are delineated in Table 1. The process flow of the proposed multifaceted CYP framework is depicted in Figure 1.

Table 1. Implementation platform details

Figure 1. Multifaceted crop yield prediction (CYP) model

3.2 Dataset description

We utilized the CYP dataset from Kaggle, which compiles the Pesticides & Yield dataset from Food and Agriculture Organization (FAO) and the Rainfall & Temperature dataset from the World Data Bank [32-34]. It includes 28,242 records from 1990 to 2013 across 101 countries, including Russia, Ecuador, American Samoa, Gambia, China, India, and other geographically diverse countries. Table 2 presents the dataset attributes and their relevance to CYP.

Table 2. Crop yield prediction (CYP) dataset attributes

3.3 Data preprocessing

We began the preprocessing pipeline with the removal of irrelevant attributes (e.g., unnamed columns), followed by verification of missing values. No missing entries were observed. Outliers were handled using the Winsorization technique, where extreme values in numerical features were limited using a 5% threshold. The Interquartile Range (IQR) and Z-score methods were employed for exploratory identification of outliers; however, Winsorization was applied as the final transformation method. Tfigurehe effect of this transformation is illustrated in Figure 2(a)-(c).

To address skewness in feature distributions, a logarithmic transformation of the form $\ln (1+x)$ was applied as shown in the Eq. (1).

log_value=\ln (1+value)   (1)

Prior to transformation, feature values were adjusted to ensure non-negativity. The distributions before and after transformation are shown in Figure 2(d)-(e).

(a)

(b)

(c)

(d)

(e)

Figure 2. (a) Winsorization Before/After – Rainfall, (b) Winsorization Before/After – Pesticides, (c) Winsorization Before/After – Temperature, (d) Log Transformation Before/After – Pesticides, (e) Log Transformation Before/After – Temperature

Categorical features were encoded using label encoding to convert them into numerical representations suitable for machine learning models; however, no ordinal interpretation was assumed during analysis.

Feature scaling was performed in two sequential stages. First, Min-Max scaling was applied to normalize features within the range [0,1]. Subsequently, Z-score standardization was applied to ensure zero mean and unit variance across features.

All preprocessing steps, including transformations and scaling, were applied exclusively to the training data and subsequently applied to the test data to prevent data leakage.

These preprocessing steps collectively enhanced data quality and improved model robustness.

3.4 Exploratory crop yield data analysis

After careful data preprocessing, we performed correlation analysis as part of EDA to explore relationships between the agricultural features and crop yield, which is illustrated in Figure 3. The crop yield has a correlation of -0.11 with average temperature and -0.06 with rainfall, which implies that these two environmental factors, among others, are weakly correlated with crop yield. The negative correlation of Area with rainfall (-0.23) and pesticide use (-0.31) indicates an association where regions with lower rainfall exhibit relatively higher pesticide usage; however, no causal relationship can be inferred.

Figure 3. Heatmap of correlation coefficients for agricultural variables

A weak positive correlation (0.10) between pesticide usage and crop yield suggests a weak linear association with crop yield. Similarly, the negligible correlation (0.03) between temperature and pesticide usage indicates minimal linear association between temperature and pesticide usage. A moderate correlation (0.31) between temperature and rainfall reflects a climatic trend where higher temperatures are associated with increased precipitation.

Following the correlation analysis, we generated a density map (Figure 4) to better understand crop yield distribution across various countries concerning average annual rainfall. Most of the reported yields hover around the 500-1,500 mm average rainfall, and there are high yield outliers spread across various levels of rainfall. It is important to note that over 2,000 mm regions of annual rainfall have few recorded data points, hinting at possible issues like waterlogging and nutrient-deficient soils. Among all countries, Belgium has the highest yield record (~501,412 hg/ha), which may be attributed to advanced agricultural practices, efficient irrigation systems, and better resource management. These outcomes suggest that rainfall alone does not fully explain variations in crop yield.

Figure 4. Crop yield distribution by rainfall and country

Table 3. High-yield records by crop item and country

Table 3 presents high-yield records across different crop types and countries, highlighting dominant contributors to global production. We noticed India’s significant role in producing substantial yields in cassava and sweet potatoes, constituting 25% and 18% of the recorded global share, respectively. This might be due to 14% of the dataset being of India.

Among the cultivated crops, potatoes exhibited the highest average yield (Figure 5), followed by cassava, sweet potatoes, and yams. The top ten potato yields and the corresponding features across different countries and years are illustrated in Table 4. As per the table analysis, Belgium topped in potato yield in 2011, 2004, and 2002. These instances were correlated with high pesticide usage (5,740-9,204 tonnes) and low rainfall (847 mm). During 2007-2013, New Zealand was able to sustain yields at a moderate level (~478,154-495,751 hg/ha) with moderate pesticide usage (~5,086 tonnes) and high rainfall (1,732 mm).

Figure 5. Average crop yield by crop type

Table 4. Top ten potato yield by country, item, year, pesticides usage, and average rainfall

To investigate variations in crop yield among different countries, we created a boxplot shown in Figure 6. It also illustrates that Belgium is at the top in crop yield, reinforcing its highly productive agricultural practices.

Figure 6. Crop yield distribution by country

We visualized the crop yield trend over the years for the top fifteen countries by yield using the line plot shown in Figure 7. We noticed that countries like Denmark, Netherlands, and Germany have consistently high yields, while others like Egypt and Jamaica show lower values. Some countries display stable trends, while others experience significant fluctuations. These fluctuations may be influenced by climate change, rainfall, or pesticide usage.

Figure 7. Temporal analysis of crop yield for the top 15 countries

We created a spline interpolation plot (Figure 8) and noticed a general upward trend in crop yield from 1990 to 2013.

Figure 8. Average crop yield over time (1990-2013)

While these exploratory analyses provide valuable insights, the weak linear relationships observed indicate that CYP requires models capable of capturing complex nonlinear interactions. This motivates the use of ML to improve the accuracy of agricultural predictions.

We scrupulously evaluated the ML models using 10-fold CV on 80% of the training data and unseen test sets. The mean performance results of the 10-fold CV and test dataset performance are summarized in Tables 7 and 8, respectively. All models were evaluated under identical experimental conditions to ensure fair comparison. This includes the same train-test split, consistent preprocessing pipeline, identical cross-validation strategy, and standardized evaluation metrics.

Table 7. Performance comparison of machine learning models (10-fold CV)

5.1 Performance evaluation using CV

To ensure robustness and generalization, 10-fold cross-validation was performed. This evaluation aimed to evaluate the stability of the models across several partitions of the data and detect early signs of overfitting before test set evaluation.

RF, with MSE of 1.17 × 10⁷ and R² of 0.9985, exhibited the best performance during CV and was followed closely by BR (MSE = 1.30 × 10⁷, R² = 0.9984). These two models maintained their high predictive power throughout the different validation folds. XGBst and CatBst also had similarly competitive R² values (0.994) but had comparatively higher MAE values. This suggests that these models make slightly less accurate predictions than ensemble-based models.

Conversely, linear models like LR, Ridge, and Lasso showed the lowest predictive performance, with an R² value of 0.166. This implies that linear methods are unable to handle complex nonlinear data. In addition, AdaBst and GB had a mediocre performance with R² of 0.67 and 0.93, respectively. However, their increased MAE values exposed poor accuracy.

After CV, we did not eliminate underperforming models outright. Instead, we retained all models for further assessment on the test set to validate their generalization ability on unseen datasets and verify whether their CV results were consistent.

5.2 Model evaluation on the test set

The models’ performances on the test dataset largely mirrored CV results, with RF providing the best accuracy with R² of 0.9986 and MSE of 1.08 × 10⁷. BR performed similarly with R² of 0.9984 and MSE of 3.98 × 10⁷, which is additional evidence for its stability. Of the boosting-based models, near the top was XGBst (R² = 0.995, MSE = 8.12 × 10⁷). LGBM and CatBst, despite their high validation scores, underperformed on the test set, indicating potential overfitting caused by dataset distribution shifts. A key observation was the failure of linear models to generalize, as reflected in their poor performance on the test set (R² = 0.16, MSE = 668.72 × 10⁷ for LR). Although AdaBoost achieved moderate performance in both CV (R² = 0.67) and test evaluation (R² = 0.72), it remained significantly inferior to other ensemble models.

Generally, the similarity between CV and test results for top-performing models (RF, BR, and XGBst) indicates reliability. On the other hand, all the linear models underperformed during both tests, confirming their incompatibility with the complex structure of the dataset.

5.3 Uncertainty quantification

In addition to generating point predictions from models, it is crucial, especially for critical applications in CYP, to quantify the confidence and reliability of these outputs. To achieve this, we performed uncertainty quantification using two methods: Gaussian Process Regression (GPR) and bootstrap resampling.

5.3.1 Predictive uncertainty with gaussian process regression

GPR quantifies prediction uncertainty by providing mean predictions and associated predictive variances used to construct confidence intervals. A Gaussian Process (GP) is a collection of random variables such that any finite number has a joint Gaussian distribution. It is fully specified by its mean function m(x) and covariance (kernel) function $k\left(x, x^{\prime}\right)$ as specified in Eq. (8):

$f(x) \sim \mathrm{GP}\left(\mathrm{m}(\mathrm{x}), \mathrm{k}\left(\mathrm{x}, \mathrm{x}^{\prime}\right)\right)$   (8)

where, $m(x)=\mathrm{E}[\mathrm{f}(\mathrm{x})]$ and $\mathrm{k}\left(\mathrm{x}, \mathrm{x}^{\prime}\right)=\operatorname{Cov}\left(\mathrm{f}(\mathrm{x}), \mathrm{f}\left(\mathrm{x}^{\prime}\right)\right)$.

5.3.2 Posterior predictive distribution for uncertainty quantification

Given training data $D=\left\{\left(x_i, y_i\right)\right\}_{i=1}^N$ (inputs $x_i$, observed outputs $y_i$), and assuming a GP prior, the posterior predictive distribution for a new test point $\mathrm{x}_*$ is also Gaussian as given in Eq. (9):

$\mathrm{p}\left(\mathrm{f}\left(x_*\right) \mid \mathrm{X}, \mathrm{y}, x_*\right)=\mathrm{N}\left(\mu\left(x_*\right), \sigma^2\left(x_*\right)\right.$ (9)

where the posterior mean $\mu\left(x_*\right)$ and variance $\sigma^2\left(x_*\right)$ are given in Eqs. (10) and (11):

$\mu\left(x_*\right)=k_{* f}^T\left(K_{f f_{+}} \sigma_n^2 I\right)^{-1} y$   (10)

$\sigma^2\left(x_*\right)=\mathrm{k}_{* *}-k_{* f}^T\left(K_{f f+\sigma_n^2} I\right)^{-1} \mathrm{k}_{* f}$   (11)

Here, $K_{f f}$ is the training data covariance matrix, $k_{* f}$ is the covariance vector between the test point and training data, $\mathrm{k}_{* *}$ is the test point variance, and $\sigma_n^2$ is the noise variance.

We observed that 95.37% of test points fell within our study’s 95% confidence intervals, as illustrated in Figure 9 This outcome suggests that the model’s uncertainty estimates are well-calibrated.

Figure 9. Gaussian process regression predictions with 95% confidence intervals

5.3.3 Performance uncertainty via bootstrap resampling

In bootstrap resampling, we repeatedly resampled the test set for 2,000 bootstrap iterations to re-evaluate the models on these resamples. Confidence intervals were calculated based on the models’ performance in all the iterations to quantify their variability and robustness. The results of this bootstrap-based uncertainty analysis on the test dataset are presented in Table 8.

Table 8. Performance comparison of machine learning models with 2000 bootstrap-based uncertainty estimates (Test dataset)

5.4 Computational efficiency and deployment considerations

One of the notable objectives of this study is to determine the efficacy of the models not only based on performance evaluation metrics but also on their computational efficiency, such as training time, inference time, and memory usage. The results highlight a trade-off between the accuracy and computational efficiency of the models.

5.4.1 Training time

RF took the longest training period, with 15.46 s in cross-validation and 18.92 s when trained on the entire training dataset before the test evaluation. Its high computational cost was justified by its superior accuracy. On the other hand, the boosting models trained significantly faster, with XGBst and LGBM completing their training within 0.28 and 0.27 s during CV, while training time before final testing took 0.34 s and 0.31 s, respectively. These results indicate that boosting models are better suited for situations where frequent retraining is needed, whereas RF is optimal in situations where accuracy comes first. Though BR models have slightly lower accuracy, these models are viable for deployment due to minimal computational requirements. The linear models like Ridge, Lasso, and LR were trained simultaneously, each taking under 0.05 s. However, these models are not viable for deployment due to their much lower accuracy.

5.4.2 Inference time

Model inference speed is imperative, especially in scenarios that require real-time solutions. XGBst and LGBM provided the fastest inference times during CV, requiring 0.010 and 0.016 s, respectively, while on the test dataset, they needed 0.011 and 0.017 s. These metrics definitively demonstrate their efficiency in time-constrained scenarios.

RF showed slower inference times, needing 0.411 s during CV and 0.437 s during test evaluation. BR did almost as well but was better in terms of a lower inference time of 0.193 s in CV and 0.206 s in test evaluation, and thus, it stands out as a good balance.

Linear models had negligible inference times, taking under 0.002 s, but had no predictive ability. These results further confirm that XGBst and LGBM are optimal in speed and accuracy. RF has been proven useful in cases where more latency is acceptable in return for better performance.

5.4.3 Memory consumption

Memory efficiency is another important factor when deploying a model, especially when working with limited resources. RF used the most memory, at 126.93 MB for CV and 132.11 MB for test evaluation, which could be problematic for deployment in low-resource environments. However, BR reduced memory usage to 76.42 MB for CV and 81.87 MB for test evaluation with good performance.

Boosting models showed greater efficiency as XGBst required 104.34 MB for CV and 107.82 MB for test evaluation. Nonetheless, the memory efficiency of LGBM made it the most computationally efficient, using only 1.71 MB and 1.96 MB, for CV and test evaluation, respectively. These findings show that LGBM is optimal for memory-limited deployments.

Lastly, while linear models consumed less than 1 MB, their poor performance limited their practicality. In sum, these results point to LGBM outperforming other models for memory efficiency, followed by XGBst, while RF and BR models had considerably higher memory requirements to achieve accuracy. Bar plots illustrating the performance comparison of models in 10-fold CV and test dataset are shown in Figures 10 and 11 respectively.

Figure 10. Performance comparison of models in 10-fold CV

Figure 11. Performance comparison of models in test dataset

5.4.4 Trade-offs and deployment suitability

Balancing accuracy and computational efficiency is crucial for model selection. RF remains the best choice for high-accuracy applications where computational cost is secondary. However, BR provides a viable alternative with reduced resource demands. XGBst and LGBM offer lower training times, minimal inference latency, and efficient memory usage for real-time applications. Among them, LGBM is the most memory-efficient. In environments with strict resource constraints, LGBM is ideal due to its balance of efficiency and accuracy. Evaluating computational efficiency across both CV and test datasets ensures a well-rounded understanding of each model’s feasibility for deployment. These results confirm that tree ensemble models consistently outperform linear models when predicting crop yield. This emphasizes the need for nonlinear modeling approaches for large-scale agricultural data.

Moreover, the gap between CV and test performance for some boosting models (like CatBst and LGBM) raises concerns about overfitting, which indicates that their robustness could be enhanced through specific dataset tuning in future efforts. This study evaluates all the candidate models to choose the most suitable ones for real-world deployment based on predictive accuracy, generalization ability, and computational efficiency rather than prematurely removing weaker candidates.

5.5 Residual analysis of top performing machine learning models

We created residual plots to interpret the nature and distribution of errors in prediction for the top 6 models (RF, BR, XGBst, CatBst, KNN, and LGBM). These plots illustrated the differences between predicted and actual values (residuals) and predicted values. For an unbiased model, residuals should exhibit a random scatter around zero with consistent variance. The model has not fully captured underlying nonlinear relationships if the residual plot has persistent curved patterns. An inconsistent spread of residuals suggests varying error variance with predicted values. A systematic shift of residuals consistently above or below the zero line for specific predictions indicates model bias.

We can observe a broadly consistent and uniform scatter around the zero line in the residual plots of the models in Figure 12. This indicates the efficacy of these models in capturing underlying relationships in crop yield data with minimized systematic biases and improved modeling of nonlinear relationships. These plots reinforce confidence in models’ reliability and generalization ability.

Figure 12. Residual plots of top 6 models by mean R²

5.6 Performance analysis of recurrent neural network variants

This study was conducted to probe how depth, hidden unit size, and attention mechanisms influence predictive performance, training cost, and resource utilization of all the RNN variants (simple RNN, LSTM, and Attention-LSTM models). The performance of the best-performing RNN variants within each variant category is presented in Table 9 for a broader context. Among the simple RNN configurations, the three-layer RNN with 64 neurons and Tanh activation (R² = 0.9916, MSE = 6.73 × 10⁷) demonstrated the best performance. This performance was followed closely by the two-layer RNN with 128 neurons and Tanh activation (R² = 0.9896, MSE = 8.30 × 10⁷). We also noticed that the predictive performance and computational efficiency of a two-layer RNN with 128 neurons and ReLU were lower than those of RNN_2L_128N_tanh. These results emphasize that the optimal performance of any configuration depends on layers, neurons, and activation functions rather than just the depth of the network. This also indicates that Tanh activation was more effective in capturing long-term dependencies, as evidenced by the superior performance of Tanh-activated RNNs compared to their ReLU counterparts. Furthermore, smaller RNNs came with competitive inference speeds and a noticeable trade-off in accuracy.

Table 9. Top 5 models per architecture category based on R² Score

Among all the RNN variants, the LSTM models significantly outperformed their basic RNN counterparts. Even smaller LSTM configurations, such as LSTM_2L_128N and LSTM_3L_64N, maintained high accuracy (R² > 0.996) and reasonable inference times around ~0.79–0.93s.

We also witnessed that integrating attention mechanisms into LSTMs can improve the performance of LSTM. The Attention-LSTM with 3 layers and 128 neurons achieved high accuracy, R²=0.9965, MSE = 2.77 × 107, and MAE = 2255.42 with a moderate inference time (0.8744 s).

5.7 Comparison with traditional machine learning models

Though the predictive performance of LSTM and Attention-LSTM is commendable, ML models, such as RF (R² = 0.9986) and BR (R² = 0.9984), outperformed them in accuracy and significantly in training time.

The LSTM, with three layers and 128 neurons, achieved the best performance across all architectures. However, its training time (1039.78 s) is significantly higher than that of RF (14.73 s) and XGBst (0.15 s).

In terms of resource consumption, deeper RNNs were computationally more expensive. Interestingly, the RNN architectures utilized memory on par with other ML models and considerably less than RF. For instance, the highest memory usage of RNN is 5.89 MB, but for RF, it is 127.55 MB despite using GPU acceleration. The longer training times and larger memory footprints remained a limitation of all the RNN variants compared with ML models.

However, inference times across DL and ML models were relatively comparable. For example, the two-layer RNN with 64 neurons and Tanh achieved an inference time of 0.6617 s, close to that of RF and KNN, having 0.3832s and 0.3586s, respectively. It suggests that once RNNs are trained, they can perform at speeds comparable to traditional ML models in time-sensitive applications.

Although RNN variants required more time and memory to train, optimized configurations achieved performance close to top ML models while maintaining real-time inference capabilities. These observations can also be drawn from the scatter plots shown in Figure 13.

Figure 13. Performance evaluation of recurrent neural network (RNN) models

The results from this study confirm that:

• LSTM and Attention-LSTM architectures outperform simple RNN variants.
• Tanh activation is superior to ReLU in deeper RNNs in this study.
• Model depth, neuron count, and activation functions substantially impact performance and resource requirements.

Attention mechanisms can enhance the robustness and performance of the models.

5.8 Analysis of AgroStackNet hybrid ensemble models

Out of all the evaluated configurations, the five best-performing AgroStackNet hybrids are presented in Table 10. The model combining XGBst, CatBst, RF, and BR with Ridge as the meta-learner achieved the lowest test MSE (1.0777 × 10⁷) and the highest R² score (0.998653), indicating excellent generalization.

Table 10. Top AgroStackNet hybrid models and their performance metrics

Other configurations that used Ridge or BRidge as meta-learners also performed comparably well. However, they had a minimal trade-off in training time and memory usage. Overall, the Ridge consistently performed effectively across the best ensembles.

These results demonstrate that carefully assembled hybrid ensembles can significantly increase prediction accuracy and maintain efficiency. This approach eliminates the need to rely on overly complex model stacks.

5.9 Comparative insights with prior research

To contextualize our findings, we benchmarked our models’ performance against prior studies conducted on the same FAO-based crop yield dataset [28-31]. Table 11 summarizes the results of selected models from these studies, provides concise methodological descriptions, and presents the results of our corresponding implementations. We can observe from the table that our approach consistently achieves superior performance compared to prior studies and were made possible by rigorous preprocessing and advanced outlier handling techniques. Also, it can be noticed that these earlier works did not investigate the computational efficiency of the models, leaving a critical gap in CYP. This study evaluates computational efficiency as a core objective while achieving better predictions to address this gap. For example, our RF model achieved an R² of 0.9986 and an MAE of 433.02. This performance not only significantly surpasses prior RF implementations (with R² values ranging from 0.9674 to 0.9837 in prior studies) but also outperforms the overall best result from any prior model (R² = 0.9852, MAE = 2113.05 [30]). The slight improvements in R² score may seem modest, but they can significantly impact decisions in large-scale agricultural forecasting, such as those related to resource allocation, food security, and crop insurance schemes.

Table 11. Comparison with the prior studies