The analysis uses the ICRISAT (International Crops Research Institute for the Semi-Arid Tropics) district-level agricultural database, a long-term harmonized resource on Indian crop production. The panel covers 52 years (1966–2017), 311 districts across 20 states, and reports yield information for 23 crops. In this work, we focus on six major crops of central importance to Indian food security and the agricultural economy rice, wheat, maize, groundnut, cotton, and sugarcane.

The data are organized in a panel format with district–year as the observational unit, such that each record corresponds to a unique (District Code, Year) combination. Our consistency checks confirm that this key is unique for all 16,146 observations, indicating the absence of duplicate entries. The main characteristics of the dataset are summarized in Table 1.

The dataset reports crop-wise yield (kg/ha) together with the corresponding cultivated area, enabling consistent quantification of productivity. This coverage facilitates rigorous examination of yield variation across contrasting agro-climatic zones, cropping systems, and management practices. Descriptive statistics for the six focal crops are summarized in Table 2.

Data quality is crucial for building accurate yield prediction models, so we apply a straightforward preprocessing pipeline. All negative yield values are reset to zero, as such values are physically implausible and typically arise from data entry errors or missing-value codes. Specifically, we identified a small number of observations with yield values of −1.00 kg/ha across the six crops (fewer than 0.2% of total data). These values most likely represent placeholder codes for “data unavailable” in the original database rather than actual negative yields, which are physically impossible. Resetting these to zero follows standard data cleaning practice; sensitivity analysis confirms that results are unchanged when these observations are excluded entirely. For each crop, we cap extremely high yields at the 99.5^th percentile, which retains 99.5% of observations while limiting the influence of atypical or erroneous extremes and stabilising model performance.

We also verify temporal consistency by examining year-to-year changes at the district level and flagging cases in which inter-annual yield differences exceed five standard deviations for manual review, as such abrupt shifts are unlikely to be agronomically realistic and may indicate data quality problems. Zero-yield observations are retained but treated carefully in error metrics; in particular, the Mean Absolute Percentage Error (MAPE) is computed only for observations with yield above 100 kg/ha to avoid division by near-zero values, which is especially important for cotton where zero-yield entries are frequent. The 100 kg/ha threshold for MAPE calculation is justified on both mathematical and practical grounds: (1) division by near-zero values produces extreme, meaningless percentages (e.g., a 5 kg/ha error on 1 kg/ha yield produces 500% MAPE); (2) sub-100 kg/ha yields typically represent crop failures, abandoned plots, or negligible production where percentage error lacks practical meaning; (3) we report alternative metrics (MAE, MedAE, R²) that require no such exclusion, ensuring transparent evaluation. Supplementary sensitivity analysis using thresholds of 50, 100, and 200 kg/ha confirms that model rankings and conclusions are robust to this choice. This cleaning process yields refined distributions with improved statistical properties: for example, rice yields are reduced from extreme values above 8,000 + kg/ha to approximately 4105 kg/ha, wheat from over 10,000+ to 4,485 kg/ha, and maize from more than 15,000+ to 5,898 kg/ha. These capped values remain within agronomically realistic ranges while effectively removing data quality artefacts.

An important feature of the ICRISAT dataset is the absence of missing yield values for the six target crops across all 16,146 district–year observations. This level of completeness is uncommon in agricultural datasets and removes the need for imputation procedures that may introduce bias. The 0% missing rate for all six crops (Table 2) ensures that the models are trained exclusively on empirically observed yields rather than on statistically inferred values.

At the same time, many districts report zero yields for specific crops, reflecting the fact that not all crops are cultivated in all regions. In our analysis, we explicitly distinguish between true zeros (indicating non-cultivation in a given district–year) and missing data. These true zeros carry meaningful information about spatial and temporal cropping patterns and are therefore retained, while being treated carefully during network construction and feature engineering to avoid misinterpretation as data quality issues.

Temporal dependencies are fundamental to agricultural yield patterns. Crop yields in a given year are influenced by previous years’ productivity through soil health carryover, pest pressure buildup, farmer learning, and autocorrelated weather patterns. All mathematical definitions used in the proposed pipeline are provided in Equations 1–29. We engineer five types of temporal features for each crop:

1. Lag features: Direct yield values from previous years:

where yt−k(c) represents the yield of crop c at time t − k. These features capture immediate historical performance and provide baseline information for prediction.

2. Rolling mean features: Moving averages smooth short-term fluctuations and capture medium-term trends:

where w = 3 is the window size. We use a 3-year rolling mean to balance responsiveness to recent changes with stability against annual volatility.

3. Year-over-year change: First differences capture growth trends and productivity improvements:

This feature helps models learn whether yields are increasing, decreasing, or stable.

4. Rolling standard deviation: Volatility measures quantify yield stability:

where w = 3. High volatility may indicate vulnerability to environmental shocks or inconsistent management practices.

5. Temporal aggregates: District-level mean yields over the reference period (2008-2017) provide baseline productivity expectations for each region.

Algorithm 1 formalizes the temporal feature generation process.

Crop diversification affects yield through multiple mechanisms including risk spreading, pest and disease management, soil health, and labor distribution. We compute two diversification metrics:

1. Simpson diversity index: A commonly used measure in ecology adapted for agricultural diversification:

where pi=Ai∑j=1nAj is the proportion of area under crop i, and n is the total number of crops in the district. This index ranges from 0 (complete specialization) to nearly 1 (maximum diversification), with higher values indicating greater diversity.

2. Number of active crops: A simple count of crops with non-zero cultivation area:

where ⊮(·) is the indicator function. This metric provides an intuitive measure of cropping system complexity.

Districts with higher diversification may exhibit different yield patterns due to resource competition, complementary effects, or risk management strategies. These features enable models to account for cropping system context when predicting individual crop yields.

Network-based features will be described in detail in Section 3.3 following the network construction methodology. These features include six centrality measures derived from the district similarity network: degree, strength, closeness, betweenness, eigenvector, and clustering coefficient.

The district similarity network represents structural relationships between geographic units based on their crop yield patterns. The intuition is that districts with similar yield profiles may share common characteristics such as climate, soil quality, agricultural practices, or institutional factors, even if they are not geographically proximate.

Network Construction Procedure:

1. Reference period selection: We use the most recent decade (2008-2017) as the reference period for network construction. This 10-year window balances data richness with contemporary relevance, ensuring that the network reflects current agricultural patterns while avoiding excessive influence from historical conditions that may no longer apply.

2. District yield profile aggregation: For each district d, we compute the mean yield across the reference period for all 23 crops in the database:

where y¯d(c)=110∑t=20082017yd,t(c) represents the average yield of crop c in district d over the reference period. Missing or zero values are retained as zeros in the profile vector.

3. Similarity computation: We use cosine similarity to measure the resemblance between district yield profiles:

Cosine similarity ranges from -1 to 1, with higher values indicating greater similarity in yield patterns. This metric is scale-invariant, meaning it focuses on the pattern of yields rather than absolute magnitudes, making it suitable for comparing districts with different overall productivity levels.

4. Enhanced thresholding: Unlike previous studies that use relatively low thresholds (0.5-0.7), we employ a stringent threshold of τ = 0.80. This higher threshold ensures that only districts with genuinely similar yield patterns are connected, creating a sparser, more interpretable network:

5. Top-k pruning: To further enhance network interpretability and focus on the strongest relationships, we implement top-k pruning where each node retains only its k = 15 strongest connections. For each district d, we:

where topk(·) selects the k neighbors with highest similarity scores. This strategy prevents highly connected hubs from dominating the network while maintaining computational efficiency.

6. Network density optimization: The resulting network has 311 nodes (districts) and 2,996 undirected edges, yielding a density of 6.2%. This density falls within the optimal range (5-20%) for meaningful pattern extraction while avoiding excessive connectivity that would dilute signal.

Algorithm 2 formalizes the enhanced district network construction process.

Network Centrality Measures:

From the constructed network, we extract six centrality measures for each district, which serve as features 350 in our prediction models:

1. Degree centrality: The number of direct connections:

High-degree districts are connected to many similar regions and may represent diverse or typical cropping patterns.

2. Strength centrality: Sum of edge weights (similarity scores):

This weighted measure accounts for the strength of similarities, not just their count.

3. Closeness centrality: Inverse of average shortest path length:

High closeness indicates a district is structurally central, able to reach others through short paths.

4. Betweenness centrality: Fraction of shortest paths passing through the node:

where σ_st is the number of shortest paths from s to t, and σ_st(d) is the number passing through d. High betweenness districts act as bridges between different network communities.

5. Eigenvector centrality: Recursive measure valuing connections to well-connected nodes:

where λ is the largest eigenvalue of the adjacency matrix. This metric identifies districts connected to other influential districts.

6. Weighted clustering coefficient: Measures local network density:

where ad′d″ indicates whether an edge exists between neighbors d′ and d″. High clustering suggests the district belongs to a tightly-knit group of similar regions.

Figure 1 visualizes the district similarity network structure. (Note: This figure would be generated from the network construction code and saved during pipeline execution).

The crop co-occurrence network characterises which crops tend to be cultivated together within the same district–year observations, thereby reflecting farmers’ crop portfolio decisions that may be shaped by complementarity, risk diversification, shared resource use, or market opportunities (Lin and Schilstra, 2019).

To construct this network, we first identify, for each district–year, the set of crops with non-zero yields. For every pair of crops (c_i,c_j) that co-occur in a given district–year, we increment the corresponding edge weight w_cic_j by one. Aggregating over all observations yields an undirected, weighted network with 23 nodes (crops) and 253 edges, where edge weights represent the frequency with which crop pairs are jointly cultivated.

The co-occurrence network serves as a complementary representation to the district-level network but is not explicitly incorporated as a feature source in the present modelling framework. Rather, it is used to elucidate broader cropping-system patterns and has the potential to guide future investigations into crop–crop interaction effects.

1. Naive (Lag-1) Predictor:

The simplest baseline uses the previous year’s yield as the prediction:

For missing lag-1 values, we use the training set mean. This baseline represents the “persistence forecast” commonly used in time-series analysis and provides a minimum performance threshold.

2. Rolling Mean (3-year average):

This baseline averages the three most recent years:

This approach smooths year-to-year volatility and may perform better than lag-1 in volatile environments. For missing values, we use available recent years or fall back to the training mean.

1. Ridge Regression.

Ridge regression extends ordinary least squares by introducing an ℓ₂ penalty on the regression coefficients, which reduces variance and alleviates instability caused by multicollinearity among predictors.

In our setting, the regularization strength is fixed at α = 10.0, a value deliberately chosen to be larger than the default in order to counteract the high degree of feature redundancy in a rich, multi-source design matrix. Prior to model fitting, all predictors are transformed using RobustScaler, which centers each feature by its median and rescales it by the interquartile range:

This preprocessing step dampens the influence of extreme values and yields a more stable optimization landscape. Overall, Ridge regression serves as an interpretable linear baseline: it captures first-order effects in a transparent manner and provides a parametric reference against which the added complexity of non-linear, ensemble-based approaches can be meaningfully evaluated.

2. Random Forest.

Random Forest represents a non-parametric, ensemble-based approach that combines the predictions of many decision trees, each trained on a bootstrap sample of the data and a randomly selected subset of features. The final prediction is obtained by averaging over the individual trees.

where T_b denotes the prediction of the b-th tree and B = 300 is the number of trees in the ensemble. We increase the number of estimators from the scikit-learn default of 100 to 300 to obtain a more stable aggregate model and to reduce variance. The maximum depth of the trees is left unconstrained (max_depth = None), allowing the ensemble to discover complex, high-order interactions when supported by the data, while min_samples_leaf = 1 permits fine-grained partitioning at the leaves. A fixed random seed (random_state = 42) ensures reproducibility, and n_jobs = -1 enables the use of all available CPU cores to keep training times manageable. Because decision trees operate on the original feature scales and split locally in feature space, the Random Forest is naturally robust to outliers and heterogeneous feature magnitudes, and therefore does not require explicit feature scaling. In practice, this model captures non-linear relationships and interaction effects that are inaccessible to purely linear methods, making it a strong and relatively robust benchmark.

3. Gradient Boosting.

Gradient Boosting takes a different ensemble perspective by constructing the model in a sequential, stage-wise fashion. Instead of averaging many independent trees, it iteratively adds new trees that are trained to correct the residual errors of the current ensemble. Formally, the model at iteration m can be written as.

where h_m is the m-th base learner, fitted to the negative gradient of the loss with respect to F_m−1, and ν is the learning rate that controls the contribution of each new tree. In this study, we employ scikit-learn’s GradientBoostingRegressor with n_estimators = 100, learning_rate = 0.1, and shallow trees of depth max_depth = 3. This configuration strikes a pragmatic balance: the number of boosting stages is large enough to capture non-linear patterns, while the combination of a moderate learning rate and shallow trees provides built-in regularization. As with the other models, random_state = 42 is used to ensure replicability. Gradient Boosting is often capable of achieving high predictive accuracy by focusing successive learners on the hardest-to-predict instances; however, this flexibility also makes it more sensitive to overfitting, so careful control of hyperparameters and regularization is essential, particularly when the dataset is noisy or only moderately sized.

The overall experimental workflow, including the common feature pipeline and the split into the no_network and with_network configurations, is illustrated in Figure 2. This setup allows us to isolate and quantify the incremental predictive value of the network-based features.

Standard k-fold cross-validation is inappropriate for time-series data as it violates temporal ordering and can lead to data leakage where future information influences past predictions (Roberts et al., 2017). We employ time-series split validation where training sets contain only historical observations relative to the test set.

We construct the cross-validation folds using an expanding-window, time-aware procedure tailored to the temporal span of the dataset (1966–2017). All unique years are first ordered chronologically and the full range is partitioned into four contiguous segments of approximately equal length, defined by three temporal split points ((split₁, split₂, split₃)). Using these splits, we form three non-overlapping folds. For fold f ∈ {1, 2, 3}, the training set comprises all observations from 1966 up to and including year split_f], while the corresponding test set consists of observations from the subsequent interval (split_f, split_{f+ 1}].

In practical terms, the three folds can be read as a sequence of increasingly data-rich forecasting experiments. In Fold 1, the model is trained only on the earliest block of years and then evaluated on the immediately following period. Fold 2 enlarges this training window to cover both the early and middle years, with evaluation shifted further forward in time. Finally, Fold 3 trains on all earlier segments early, middle, and later years and assesses performance on the most recent period in the series. The exact cross-validation fold definitions (training/test periods and sample sizes) are provided in Table 3.

The choice of three folds represents a balance between evaluation robustness and data sufficiency. Each fold requires sufficient test data for reliable metric estimation (at least 10 years given district year observations), while earlier folds need adequate training data to fit complex ensemble models. We acknowledge that using only three folds provides limited statistical power for temporal stability assessment; a rolling-origin design with annual test windows would provide more granular estimates but would substantially increase computational cost. Conclusions about temporal stability should therefore be interpreted with appropriate caution, recognizing they are based on three temporal segments rather than fine-grained annual assessments.

As sketched in Figure 3, this “expanding window” design respects the natural flow of time, so that future information never leaks into the past, and it mimics the way models would actually be used in practice: fitted to whatever history is available at a given point and then asked to predict subsequent outcomes. Because each fold tests on a different portion of the timeline, the procedure also probes how stable the model is across changing conditions. At the same time, later folds benefit from longer historical records, allowing the model to learn from a richer set of patterns while the evaluation remains genuinely out-of-sample.

We employ five complementary metrics to provide comprehensive performance assessment:

1. Mean Absolute Error (MAE):

MAE provides an interpretable measure of average prediction error in the same units as yield (kg/ha). It is robust to outliers and assigns equal weight to all errors.

2. Root Mean Squared Error (RMSE):

RMSE penalizes large errors more heavily than MAE due to squaring, making it sensitive to occasional large mispredictions. It is widely used in agricultural forecasting for comparability.

3. Mean Absolute Percentage Error (MAPE):

MAPE provides scale-independent percentage error, facilitating cross-crop comparison. However, it is undefined for zero yields and unstable for near-zero values. We implement a safe version excluding observations below 100 kg/ha, particularly important for cotton with many low-yield observations:

4. Median Absolute Error (MedAE):

MedAE is highly robust to outliers and provides insight into typical prediction error, complementing the mean-based metrics.

5. Coefficient of Determination (R²):

R² measures the proportion of variance explained by the model, with values close to 1 indicating excellent fit. Negative values indicate worse performance than predicting the mean. Aggregation across folds. To ensure that MAPE is not distorted by extremely low-yield observations, we performed a sensitivity analysis using yield thresholds of 50, 100, and 200 kg/ha. The resulting MAPE values remain highly stable across thresholds, and model rankings are unchanged. See Supplementary Table 1.

For each model crop feature set combination, we report:

To determine whether the observed performance gains are systematic rather than artifacts of random variation across folds, we perform paired t-tests on the mean absolute error (MAE) values.

For each crop, we examine two comparisons. The first focuses on the added value of the network-based features. Here, we compare Gradient Boosting models trained with and without network features using a paired t-test over the per-fold MAE scores. The null hypothesis states that incorporating network features does not change the error,

while the one-sided alternative asserts that the model enriched with network information achieves lower error,

In practical terms, this test asks whether the same Gradient Boosting architecture, when supplied with network-derived covariates, consistently reduces MAE relative to an otherwise identical specification that ignores network structure. We use a significance level of α = 0.05.

The second comparison evaluates the benefit of our advanced model relative to a simple operational baseline. Specifically, we contrast the Gradient Boosting model with network features against a naïve forecasting strategy (e.g, using recent historical averages). Again, we apply a paired t-test to the MAE values obtained on each fold. The null hypothesis posits no difference in expected error,

and the alternative hypothesis formalizes the expectation that the advanced model outperforms the naïve benchmark,

This second test therefore quantifies whether the additional modelling complexity and the use of network information translate into practically and statistically meaningful improvements over a simple baseline that might be used in real-world settings.

For each of these tests, we report the t-statistic and corresponding p-value, together with the percentage reduction in MAE achieved by the better-performing model. To make the results easy to interpret at a glance, we also provide a binary significance flag indicating whether the null hypothesis is rejected at the α = 0.05 level. A p-value below 0.05 is taken as evidence against the null hypothesis: in that case, the performance difference is unlikely to be explained by random fold-to-fold variation alone and can be interpreted as a genuine, statistically supported improvement.

Model performance can change over time as the underlying data distribution evolves (concept drift). To quantify temporal stability, we examine how both error and explained variance behave across the three chronological folds.

For each model–crop combination, we first estimate a simple linear trend in MAE over folds by fitting.

where β₁ captures the direction and magnitude of the temporal trend. A positive slope (β₁> 0) means that MAE becomes larger from the earlier to the later folds, indicating that the model’s predictions are gradually getting worse over time. In contrast, a negative slope (β₁< 0) means that MAE decreases across folds, which points to an improvement in predictive performance as time progresses. Alongside the slope, we report the coefficient of determination for this regression (the “trend R²”), which summarizes how consistently performance changes across folds: higher values imply that MAE follows a clear temporal pattern rather than fluctuating idiosyncratically.

An analogous analysis is carried out for the R² metric itself. By regressing fold-wise R² values on the fold index, we obtain an R² slope that describes how the explanatory power of the model evolves over time. In this case, a positive slope indicates that the model accounts for an increasing share of the variance in later periods, while a negative slope signals that the model is gradually losing explanatory strength.

To make it easier to compare different crops and models, we define a straightforward notion of stability based on the MAE trend over time. If the absolute slope satisfies (|β₁| < 10, kg/ha) per fold, we label the model as Stable, meaning that any systematic change in error is small compared with the natural variability of the task. When the slope is larger and positive, we read this as a warning sign that performance is gradually worsening, which is consistent with concept drift or an inability of the model to keep up with changing conditions. By contrast, negative slopes, especially when accompanied by rising (R²) values, suggest that the model is not only holding up over time but is also learning from the additional information present in the later folds, and thus generalizes better to the most recent years.

For the best-performing models, we go beyond summary scores and look closely at how they behave using a set of simple but informative checks. We start with predicted-versus-observed plots, where we place the model’s yield predictions on one axis and the actual yields on the other, adding a 45-degree line to represent perfect agreement. When the model performs well, the points cluster tightly around this line across the entire range of values, indicating that it reasonably accurately captures both low and high yields.

Next, we look at residual plots, in which the residuals (ei=yi−y^i)  are plotted against the predicted values. These plots help us see whether errors are roughly random or whether there are visible patterns for example, larger errors at higher predicted yields, systematic curvature, or bands of points that might indicate bias, heteroscedasticity, or unmodelled non-linear effects.

Here we look for patterns such as increasing spread with larger predictions (heteroscedasticity), systematic curvature, or bands of points that might indicate unmodelled non-linear effects or bias. To understand the residual distribution more directly, we also look at histograms of the residuals, which help to reveal skewness, heavy tails, or a small number of extreme outliers. In addition, we use Q–Q plots comparing the empirical residual quantiles to those of a theoretical normal distribution, corresponding to the assumption.

If the model is well specified, several patterns should appear simultaneously: predicted-versus-observed points should be roughly centered on the 45-degree line; residuals should be scattered around zero with no obvious structure and roughly constant spread; the residual histogram should be close to symmetric and bell-shaped; and the Q–Q plot points should lie close to the diagonal. Noticeable departures from any of these behaviours for example, strong curvature or funnels in the residual plot, heavy tails in the histogram, or systematic bends in the Q–Q plot suggest potential misspecification, influential outliers, or violated modelling assumptions, and signal that further refinement or alternative model choices may be needed.

Algorithm 3 presents the complete model training and evaluation pipeline.

This rigorous experimental framework ensures that our findings are robust, reproducible, and scientifically sound, following best practices in machine learning evaluation for agricultural applications.

Our evaluation across the six major crops shows that the more advanced machine learning methods and Random Forest in particular deliver consistently strong predictive performance. As summarized in Table 4, the top model for each crop is a configuration that incorporates network-based features with network, underscoring the added value of the network information.

Several broad patterns stand out from Table 4. To begin with, the overall fit of the models is extremely strong. For every one of the six crops, the reported R² values are above 0.94, and for rice and sugarcane they exceed 0.98. In other words, the models capture well over 94% of the variation in observed yields, which is unusually high for large-scale agricultural yield prediction and points to near state-of-the-art performance.

The magnitude of the errors is also small when set against typical yield levels. Rice, for example, has a mean absolute error (MAE) of 35.88 kg/ha, compared with an average yield of 1,483 kg/ha, corresponding to a relative error of roughly 2.4%. Wheat shows a similar pattern, with an MAE of 54.08 kg/ha versus a mean yield of 1,489 kg/ha (about 3.6% relative error). These discrepancies are modest enough that, from a practical standpoint, the predictions lie quite close to the realized yields.

The percentage errors tell a consistent story. Mean absolute percentage error (MAPE) values range from 2.69% for rice to 8.08% for sugarcane. Even at the upper end of this interval, the model still delivers a level of accuracy that is entirely usable for planning, scenario analysis, and policy support in agricultural contexts.

Another notable feature is the stability of performance across time. Variation in the metrics across the three time-series folds is modest, and the coefficient of variation for (R²) stays low, ranging from just 0.3% for sugarcane to 6.1% for groundnut. This pattern suggests that the models are not narrowly tuned to a single time period but instead deliver a comparable level of accuracy across the different temporal segments.

Finally, the same pattern emerges when we look at model choice: in all six cases, the best-performing configuration is a Random Forest with network features (with_network). This repeated outcome suggests that Random Forest, when augmented with network-derived covariates, is particularly well suited to the structure of this problem and offers a strong default choice for district-level multi-crop yield prediction.

Table 5 provides a comprehensive comparison of all models across both feature configurations.

Rice: Rice stands out as the best-predicted crop, with the highest (R²) of 0.988 and a very small MAE of 35.88kg/ha, corresponding to only 2.4% of the mean yield. The fact that Random Forest with and without network features yield almost identical errors (35.88vs. 35.66kg/ha) suggests that, for rice, temporal and climatic signals carry most of the predictive power, while network effects play a more modest role. Gradient Boosting performs almost as well (R² = 0.987), whereas Ridge regression shows noticeably higher variability, which is consistent with the challenges posed by multicollinearity among the input features.

Wheat: Wheat also exhibits very strong performance, with (R² = 0.976) and an MAE of 54.08kg/ha (around 3.6% of the mean yield). Compared with rice, however, the variation in error across folds is slightly larger (standard deviation of 13.72kg/ha). This higher fold-to-fold variability may reflect the fact that wheat yields are more sensitive to inter-annual changes in weather, input use, or management practices, and that these factors can differ more markedly across regions.

Maize: For maize, the models achieve excellent fit (R² = 0.971) despite the very wide yield range (0–5,898kg/ha). The standard deviation of MAE (22.80,kg/ha) is higher than for rice and wheat, which is consistent with maize being cultivated under a broader set of agro-climatic and management conditions. Simple baselines perform poorly in this setting (Naive (R² = 0.474), underlining how important rich feature sets and non-linear models are for capturing maize yield variability.

Groundnut: Groundnut shows the largest error variability among the four non-perennial crops, with an MAE standard deviation of 24.12kg/ha. This pattern is consistent with groundnut’s well-known sensitivity to rainfall timing, soil moisture, and local soil characteristics. Even so, the Random Forest model achieves a strong fit with (R² = 0.946), meaning it still captures a large fraction of the yield variability. The gain over the simple baseline is especially notable: the MAE drops from 278.70kg/ha with the Naive model to 42.11kg/ha with Random Forest, making it clear that modern machine learning methods can substantially improve groundnut yield prediction.

Cotton: Cotton attains the lowest absolute MAE (13.77kg/ha), which partly reflects its lower average yield level (mean 119.82kg/ha). When expressed as a percentage, however, the errors (MAPE = 4.47%) are broadly in line with those of the other crops. Cotton is characterized by a large number of low or zero yield observations (the 75th percentile is only 202kg/ha), which makes the distribution more challenging to model. The strong performance of Random Forest in this case suggests that its flexibility is well suited to handling this mixture of zero and positive yields. We additionally evaluated performance in the low-yield regime (actual yield < 200kg/ha) to understand behavior under difficult-to-predict cases. The model tends to overpredict in this regime, consistent with training-data dominance by normal-yield observations and limited shock variables. Detailed results are in Supplementary Table 2.

Sugarcane: Sugarcane naturally exhibits the highest absolute errors, with an MAE of 129.68kg/ha, simply because its baseline yields are much larger (mean 4,484kg/ha). Even so, the fit remains excellent, with (R² = 0.986). The coefficient of variation for MAE is relatively small ((19.33/129.68 ≈ 14.9%)), indicating that the model’s performance for sugarcane is fairly consistent across the different folds. Sugarcane is also the crop that gains the most from advanced modelling: the MAE drops from 984.86kg/ha under the Naive baseline to 129.68kg/ha with the Random Forest model, an improvement of roughly 87%. This substantial reduction in error provides strong evidence that sophisticated machine learning models are well worth deploying for operational sugarcane yield forecasting.

A key question in this study is whether network-derived features add real predictive value beyond temporal and diversification cues. As shown in Table 6, comparing Gradient Boosting with and without these features leads to a somewhat surprising result the network features add little, if any, improvement in overall accuracy.

On average, the change in performance is very small. Across all six crops, the mean relative change in MAE is −0.61%, with individual effects ranging from a degradation of −4.13% for sugarcane to a modest improvement of +1.54% for cotton. These shifts are comparable in size to the natural variability observed across time-series folds and are therefore difficult to interpret as systematic gains or losses.

The crop-wise patterns are similarly muted. Rice and wheat show slight degradations in performance (−0.67% and −0.60%, respectively), maize is effectively unchanged (+0.01%), while groundnut and cotton exhibit small improvements (+0.17% and +1.54%). Sugarcane stands out as the only case with a clearly negative impact (−4.13%), suggesting that, for this crop, the added network features may introduce more noise than signal. When we repeat the comparison for Random Forest (Table 5), the picture is much the same: models with and without network features perform almost identically, reinforcing the impression that the network layer does not materially change predictive accuracy under the current feature design.

In some instances most notably sugarcane the inclusion of network features is also associated with an increase in the standard deviation of MAE across folds, which hints at overfitting or the presence of noisy, weakly informative features. Taken together, these observations run counter to the intuitive appeal of network-based approaches and motivate a closer examination of how, and to what extent, the network-derived variables are actually being used by the models. We pursue this question further through a detailed feature importance analysis.

To understand why network features contribute minimally, we examine feature importance distributions from Gradient Boosting models. Table 7 summarizes the contribution of network features to overall model importance.

The feature-importance results highlight a few simple but telling patterns.

First, network-derived features matter very little. For every crop, centrality-based network variables together account for less than 1% of total importance. The highest share appears for sugarcane at just 0.9%, while for several crops their contribution is effectively zero.

By contrast, temporal features clearly dominate. The lag-1 yield feature on its own explains between 29.8% and 48.6% of the total importance, with an average of 39.3% across crops. Other temporal variables, such as the rolling mean and lag-2 yield, also rank among the top predictors, and taken together the temporal features typically account for about 50–70% of the model’s explanatory power.

Diversification plays a secondary but visible role. The Simpson diversity index contributes between 1.6% and 4.7%, which is more than the combined contribution of all network features, but still far below that of the main temporal predictors.

There are some mild crop-specific nuances. For sugarcane, network features reach their highest relative importance (0.9%), which may reflect the influence of regional processing infrastructure captured by centrality. For maize and groundnut, network variables receive essentially zero importance, suggesting that yields are driven almost entirely by local temporal dynamics. Cotton shows a small but non-zero network contribution (0.3%), which could be linked to pest or disease pressures that spread along regional connectivity patterns.

Taken together, these findings indicate that yield prediction in this setting is overwhelmingly driven by temporal autocorrelation: what happened last year is a very strong guide to what happens this year. In comparison, the current form of the network structure adds only a weak additional signal on top of the temporal and diversification features. To verify that the learned feature importance is not fold-specific, we computed importance separately for each cross-validation fold. The top predictors show very small variability across folds, indicating stable learning rather than noise exploitation. Fold-wise stability is reported in Supplementary Table 4.

The predicted observed panels (upper left in each figure) give a compact visual summary of overall fit. For all six crops, the points lie close to the 45-degree line, which is consistent with the high R² values (0.946–0.988) reported earlier. There is no obvious tendency for the model to overpredict at low yields or underpredict at high yields; the cloud of points remains roughly centered on the reference line across the entire range, which suggests that the predictions are not systematically biased.

The plots also show that the models work well across the full spectrum of yields rather than only in some narrow band. Rice, wheat, and maize exhibit especially tight clusters, with R² values close to 0.99. Groundnut points are a little more spread out, reflecting its higher intrinsic variability. Cotton still shows a strong linear pattern despite its lower typical yields. Sugarcane, which has by far the widest yield range (up to about 12,000,kg/ha), maintains a clear, almost linear relationship between predictions and observations.

The residual panels (upper right) plot the errors ei=yi−y^iagainst the predicted values. For a well-behaved model, one expects a band of points scattered around zero with no obvious structure, and that is broadly what we see. There is no clear curvature, trend, or pattern that would point to systematic underfitting or a missing non-linear term. This supports the view that the model form is appropriate for the data at hand.

Residual spread is fairly even across most of the predicted range for rice, wheat, groundnut, and cotton, suggesting that error variance is approximately constant. Two crops merit a brief comment. For sugarcane, a small amount of extra scatter appears at the very high end of the yield range, above about 8,000,kg/ha, where a few large residuals occur. For maize, the residuals fan out slightly in the mid-range, indicating a modest increase in variance there. In all crops, a handful of large residuals (greater than about 100,kg/ha for cereals) are visible and likely correspond to years with extreme weather, local shocks, or data artefacts. These rare points do not dominate the overall pattern.

The residual histograms (lower left) give a second perspective on model fit by focusing on the distribution of errors. For each crop, the histogram is roughly bell-shaped and centered near zero, which matches the usual assumption of symmetric, approximately Gaussian noise. Mean residuals are extremely small (between about 0.0 and 0.3kg/ha), reinforcing the impression that there is no consistent tendency to overshoot or undershoot yields.

The spread of the histograms agrees with the reported RMSE values (roughly 6.6–86.0kg/ha), which indicates that our summary error metrics accurately reflect the underlying distribution rather than being driven by a few pathological cases. Deviations from a textbook normal curve do exist: in some crops the right tail is slightly longer (a few large underestimates) and the distributions show mild heavy tails. However, these departures are modest and typical for real-world agricultural data.

The Q–Q plots (lower right) compare the empirical residual quantiles to those of a theoretical normal distribution. For all crops, the points follow the diagonal line closely through the bulk of the distribution. Roughly the central 80–90% of residuals line up almost perfectly, which again supports the approximate normality assumption.

As is common, the largest deviations occur in the tails. Some points in the upper tail lie above the diagonal, indicating occasional large positive residuals (underestimation of yield), while some points in the lower tail lie below it, corresponding to large negative residuals (overestimation). These tail effects are relatively mild and confined to a small number of observations. In practice, this means that most predictions are well behaved, while a few extreme cases are harder to capture exactly. Confidence intervals based on normal errors may therefore be slightly conservative or imperfect in the extremes, but are still reasonable for day-to-day use.

Taken together, the diagnostics paint a consistent and reassuring picture. The predicted observed plots show very high predictive accuracy and little sign of systematic bias. Residuals are roughly random around zero and display only mild, crop-specific quirks, which suggests that the model structure is adequate and that major patterns in the data have been captured. The residual distributions and Q–Q plots are close enough to normal with fairly stable variance that standard statistical summaries and uncertainty estimates remain meaningful.

Small imperfections occasional outliers, slight skewness, and modest tail deviations are expected in a setting that involves weather shocks, management changes, and measurement noise. They do not materially undermine the usefulness of the models. Importantly, the diagnostics also show that the models cope well with very different yield regimes, from low-yield cotton to high-yield sugarcane. Overall, the Random Forest models with network features appear robust enough for use in operational yield forecasting, provided that users remain aware of the usual uncertainties associated with extreme events and rare outliers.