Prediction of Faba bean (Vicia faba L.) Yield under Artificial Shade Constraints in Moroccan Agro-fields

Faba bean (Vicia faba L.) is an important legume crop that constitutes one of the best sources of protein for human and animal nutrition. In consideration of its valuable nutritional qualities and climate adaptability, this important component of sustainable systems provides resilience in agricultural strategies and reformulation of diet (Nurmansyah et al., 2024). This crop occupies a significant place in intercropping systems (ICS) in Morocco. The most common combination for this farming practice is legumes/cereals (Tamta  et al., 2019) or legumes/wheat (Chamkhi et al., 2022). Intercropping leguminous crops is a promising way of enriching soils with biologically-fixed nitrogen, due to the legume-rhizobia symbiosis (Raji and Dörsch, 2020). However, in ICS, shading represents a potential constraint. Its effect, due to reduced light intensity, leads to changes in morphology, physiology, biomass and grain yield (Arenas Corraliza  et al., 2020).
       
In a physiological context, a study conducted by Hu  et al. (2023a) indicated that shading up to 60% caused a prominent decrease in SPAD entities in maize, leading to a considerable reduction in photosynthetic rate. Wang  et al. (2021b) also reported that 50% shading notably increased chlorophyll contents (a, b and total) in peanut leaves. Nevertheless, chlorophyll levels under shaded conditions were still lower than those in unshaded conditions. Le soja [Glycine max (L.) Merr.] cultivé à l’ombre présentait des niveaux de proline plus élevés que ceux cultivés en plein soleil, ce qui indique que les plantes à l’ombre protégeaient mieux leurs cellules des dommages (Muhammad et al., 2021).
       
According to Hu  et al. (2023b), insufficient light reduced kernel number, filling rate and weight of summer maize kernels. In addition, several studies have shown that light interception by artificial shade structures at different stages during the breeding period reduces soybean yield, mainly by influencing pod number. These results suggest that yield is sensitive to shading during flowering and pod formation (Rivest et al., 2009). For this purpose, yield prediction under this stress is imperative, as it estimates the faba bean’s economic productivity. Pre-harvest yield predictions are important for mitigating the effects of natural catastrophes and ensuring food security (Htun  et al. 2023). Yield prediction using vegetation indices has been evaluated in some studies on wheat, maize and durum (Mohammadi et al., 2024).
       
Many techniques have been developed for the selection of predictor variables such as Least absolute shrinkage and selection operator regression (Lasso) (Kumar  et al., 2021). Lasso improves prediction quality by reducing the regression coefficient compared with prediction models fitted by non-penalized maximum likelihood methods (href="#gupta_2024">Gupta et al., 2024). Tibshirani (1996) proposed Lasso, which can be used in the crop yield prediction technique. Lasso minimizes the sum of residual squares provided that the sum of the absolute value of the coefficient is less than a constant (Aravind  et al. 2022; Kumar et al. 2021) reported that elastic Net and Lasso proved to be the best models for wheat yield prediction in different localities of northwest India. Statistical methods such as stepwise regression and correlation are also frequently used in crop yield prediction (Das et al., 2018; Demirhan, 2024; Iniyan et al., 2023). Correlation studies are useful for measuring the strength and direction of relationships between different traits and grain yield (Zarei et al., 2013). Correlation is useful for revealing the magnitude and direction of the relationship between several yield component traits and grain yield (Islam et al., 2020). Stepwise regression is a statistical method for analyzing the relationship between a unique response variable or so-called dependent variable (i.e. crop yield) and independent variables or so-called controlled variables (Pourshoaib  et al., 2022). According to (Saraka and Kouassi, 2019), this method was used to assess the relative contributions of each yield component to yield by eliminating the effect of non-effective characteristics in the yield regression model.
         
Thus, yield prediction models can rapidly identify the most resistant genotypes (Skobalski  et al. 2024). Faba bean breeding programs aim to develop widely adaptable, disease-resistant, high-yielding genetic material (Mohammadi et al., 2024). The Moroccan faba bean breeding program has been in operation for more than 37 years, playing a very important role in the improvement of seed yield and quality. Up to now, 25 major faba bean varieties and 10 minor faba bean varieties have been registered in the official Moroccan catalog, all of which are exclusively inbred lines (Chetto et al., 2023). Variable selection in a multiple regression model plays an important role in prediction precision because (i) it makes the model easier to interpret, by removing redundant variables and adding no information (ii) it reduces the size of the problem to enable the algorithms to operate quicker, making it possible to handle high-dimensional data and (iii) it reduces over-fitting in the model.
       
Based on this previous work, the present study had a dual objective: (1) To determine if the yield of faba bean under shade stress can be predicted from these components. (2) To identify the yield components that correlate most strongly with yield. By identifying the most robust and easily measurable yield components. (3) Develop yield prediction models for V. faba using stepwise regression and Lasso regression. (4) Evaluate the accuracy of the field yield prediction models and select the best model for predicting bean yield under different shading treatments.
       
Moroccan plant breeding programs can employ indirect selection criteria to develop shade-tolerant varieties. These developed predictive models will contribute practically to the improvement of farming strategies in Morocco by showing which factors are important contributors to yield. These insights can guide farmers and agricultural planners in optimizing resource allocation, selecting suitable crop varieties and tailoring agronomic practices to local environmental conditions. This approach will equip breeders with valuable insights, enabling them to accurately identify and measure these varieties using precise analytical methods.

The trials were carried out at the National Institute of Agricultural Research (INRA) experimental station in Douyet, Morocco, over two consecutive growing seasons: 2020-2021 (Y1) and 2021-2022 (Y2). The site, located at 34°04’N latitude and 5^o07’W longitude, lies in the favorable agricultural zone of the Saïs plain, within the Province of Moulay Yaacoub, Fez-Meknes region, at an altitude of 416 meters. The soil is silty clay, consisting of 48.50% silt, 39.90% clay and 11.60% fine sand, predominantly dark Vertisols. Its chemical properties include a pH of 7.80, organic matter content of 3.63 mg/kg, 11.89 mg/kg of available P₂O₅ and 478.05 mg/kg of available K₂O. The region has a Mediterranean climate, with a dry season from May to October. Manual seeding was carried out on Y1 and Y2, respectively, where the area had been previously prepared by deep 3-disc plowing. Before this experiment, the land had been rotated with cereals. Then, fertilizers were applied and spread to improve soil quality, such as the base fertilizer NPK "10:30:10" (200 kg/ha). Plant-related biotic agents were intensively treated using commercialized pesticides (for pests and diseases) and cultural methods (e.g., mechanical hoeing for weeds) to prevent yield loss. Irrigation was administered only when plants were at risk of perishing. Supplemental irrigation of 30 mm was applied once after plant emergence.
      
The experiment was arranged in a split-plot design with two replications. The main plots were assigned to three shade treatments (S1, S2 and S3), while six Faba bean varieties were assigned to the subplots. Each elementary plot measured 21.6 m x 4 m, comprising six rows per variety (36 rows per plot) with 0.60 m inter-row spacing. A 3-meter-wide alley separated the replicates and treatments. Each plot covered 86.4 m², with a total experimental area of 811.8 m². The experiment was conducted on different plots within the station in each growing season. Subplots receiving the shading treatment were covered with Aluminet shading nets, made from high-density polyethylene with a metalized layer for added durability. The nets were stretched over wooden supports with metal wires and placed 1 meter above the ground. These structures shaded the plants from the top and sides, leaving the south side open for air circulation. The shading was designed to reduce the intensity of photosynthetically active radiation (PAR) by 50% and 90%, using a ZDS-10 Luxmeter (China), resulting in three treatments: no shading (S1), 50% shading (S2) and 90% shading (S3). These shading treatments were chosen based on another pot experiment conducted under olive trees’ shading conditions. Sunlight reduction was measured in about 650 Lux and 115 Lux corresponding to the S1 and S3 shade net respectively. Shading was applied from flowering to maturity, starting 81 days after seeding in Y1 and 127 days after seeding in Y2. The delay in Y2 was due to harsher weather conditions compared to Y1 (e.g., temperature increase of 1.7°C and precipitation decrease of 6.87 mm), which slowed Faba bean growth and delayed flowering. Yield per plant (YP) and its components included plant height, NTS (number of total stems per plant), NP (number of pods per plant), NS (number of seeds per plant) and WSS (100-seed weight per plant). To confirm the consistency of shade treatments, additional environmental variables (temperature, relative humidity and precipitations) were monitored throughout the experiment using specialized software (FieldClimate, METOS® by Pessl Instruments, Austria).
       
Stepwise regression is a combination of forward selection and backward elimination techniques. The model is built iteratively by adding significant variables and removing less significant ones based on predefined criteria such as p-value or Akaike Information Criterion (AIC). This regression was employed as the first modeling approach to determine the most significant variables contributing to yield prediction. The method used in this study was backward elimination, starting with a full model that included all predictors (dummy variables for shading and variety and the year variable). In each iteration, non-significant predictors (based on a significance level of 0.05) were removed, one by one, until a final model with only significant predictors was reached. The mathematical representation of the Stepwise Regression model is given as.
 

 

 
 Where,
Yi = Predicted yield (kg/ha).
b0 = Intercept.
bn= Regression coefficients for each predictor (NP, NS and WSS).
Xni = Represents the dummy variables for shading (S1, S2 and S3) and variety (V1 to V6), as well as the numeric variable for the year (Y1 and Y2).
e = Error term.
       
Lasso regression is a regularization method that performs both variable selection and regularization. It minimizes the residual sum of squares subject to the sum of the absolute values of the coefficients being less than a constant. Unlike stepwise regression, Lasso automatically shrinks less important coefficients to zero, thus retaining only the most impactful variables. The lasso regression model is expressed as.
 

Where,
l = Regularization parameter controlling the strength of the penalty.
       
The rest of the terms are the same as in the stepwise regression model.
       
To evaluate the performance of both the Stepwise and Lasso regression models, three key metrics were calculated: R-squared (R²), root mean square error (RMSE) and mean absolute percentage error (MAPE). These metrics provided a comprehensive assessment of model accuracy and goodness of fit. The R-squared (R²) statistic measures the proportion of variance in the observed yield data explained by the model was also calculated.
       
The root mean square error (RMSE) measures the square root of the average squared differences between the observed and predicted values and was calculated as.
 

  

 
Where,
n = Number of observations.
RMSE provides an absolute measure of prediction error, with smaller values indicating better model performance.
       
The mean absolute percentage error (MAPE) expresses prediction accuracy as a percentage and was calculated as.
 

   

       
All statistical analyses were performed using Python (version 3.11) programming language. The stats models package was used for implementing the Stepwise Regression model, utilizing the backward elimination approach to automatically select the most significant variables. The scikit-learn library was employed for the Lasso Regression model, with the LassoCV function used to perform cross-validation and determine the optimal regularization parameter l. Data preprocessing, including the conversion of categorical variables such as shading and variety into dummy variables, was handled using the Pandas library. For data visualization, Matplotlib and Seaborn libraries were employed to generate scatter plots that compare the predicted yields against observed yields for both regression techniques. These plots were further broken down by shading treatment to highlight differences in model performance across the different environmental conditions. Correlation analyses (Pearson type) were established using the Pandas library to reveal the relationships between yield parameters recorded in best shading conditions (S2 and S3) applied to different Faba bean varieties. In addition, a redundancy analysis (RDA) was conducted using skbio.stats.ordination.rda module to determine which predicting variables (yield parameters) were contributing the most to the precise yield prediction. All Statistical tests were established based on replicated treatments (n = 4) at a significance level of 5% (P<0.05).

Faba bean yield forecast under shading conditions
 
Stepwise regression was accordingly implemented to assess the yield prediction potential of the Faba bean under the three shading conditions (Table 1). The predicted yield was essentially based on NP, NS and WSS parameters, which had significant trends in all shading levels (Table 1). The distribution of data points (crosses) across the plot follows a diagonal trend, indicating a close alignment between the predicted and actual yields. The model’s performance was strong and accurate across the studied entities (R² = 0.84-0.88; R²_overall = 0.86; P < 0.001) (Table 1). Furthermore, an average deviation of predicted yields from the actual yields was observed as most of the predicted entities were below the observed ones (RMSE_overall = 214.8). The model had an average prediction error being less than 7% (MAPE_overall = 6.88%) reflecting its precision. Notably, the differences in the slopes of the S1, S2 and S3 levels suggest that the model captures the varying effects of shading conditions on yield, with S3 generally yielding higher predictions.


       
The lasso regression was applied to optimize the yield prediction aspect of the studied shading levels (Table 2). Similar to the previous model, yield predictions were attributed to NP, NS and WSS parameters, with statistical significance observed for all treatments (P <0.0001) (Table 2). The data points exhibit a strong diagonal alignment, signifying an excellent agreement between predicted and actual yields. The model’s performance demonstrated better accuracy and reliability across the shading levels, with R² values ranging from 0.93 to 0.96, culminating in an R²_overall of 0.95 (P < 0.001) (Table 2). The predicted yields show a highly minimal deviation from actual values, as reflected in the much lower RMSE overall value (102.6) compared to the stepwise model (Table 3). The lasso model’s prediction error was less than 3.5% (MAPE_overall = 3.28%), highlighting the model’s enhanced precision. The distinct slopes observed across shading levels (S1, S2, S3) suggest that the model accurately captures the nuanced impact of shading on yield, with S3 again yielding higher predictions overall.




       
The findings of this study offer significant insights into the impact of artificial shading on bean yield and highlight key factors influencing it. Stepwise and Lasso regression analyses demonstrated that NP, NS and WSS are crucial predictors of V. faba yield across all shading treatments. The exclusion of variables such as plant height, NTS and NP from the final models in specific treatments suggests that their influence may be limited under certain shading conditions. This could be attributed to the diminished relevance of these traits in altered light environments, as highlighted by recent research showing a significant reduction of plant height and stem characteristics of tropical tree species under high shading (Cifuentes and Moreno, 2022). In addition, our results align with those reported by (Chan Fu Wei and Molin, 2020), who found that models based on seed number effectively predict soybean yield with high accuracy (R² = 0.70). Furthermore, our findings were in line with other studies emphasizing the importance of seed number and seed weight in yield prediction (Wu et al.,  2017; Zhang et al.,  2023). Liu  et al. (2021) similarly observed that in high shading conditions, both seed number and seed weight are pivotal for accurate yield predictions, reinforcing the reliability of these predictors across varying shading scenarios.
       
The high R² values achieved in stepwise analysis demonstrate the model’s robustness in explaining the variability in Faba bean yield. In this case, the decrease in R² value from S2 to S3 reflects the increasing complexity of predicting yield under higher shading intensities. It is often explained by the increased variability and stress of plants at a higher level of shading, as will be assumed. This assumption is also strengthened by some recent studies that prove the trend of lowered predictive accuracy over specific variables with shading intensification (Pratzer et al., 2023). Unaccounted yield variability ranged from 11.2% to 16.2% with stepwise regression and from 4.2% to 6.2% with Lasso. This would prove that multiple factors, other than variables influence the yield. Further research work has to be directed towards localizing, identifying and quantifying those factors responsible for this variability due to environmental, genetic diversity and management practices.
       
Faba bean yield prediction models were validated to evaluate the performance of both stepwise and Lasso regression aspects during Y1 and Y2. Results showed that the prediction error with the Lasso approach was 2.38 and 2.09, respectively which was less than the ones obtained with the stepwise approach (Table 4), which implied the high prediction accuracy of the Lasso model. In this study, the Faba bean yield prediction model at the field scale developed by the Lasso method showed better performance and accuracy than that developed by the stepwise method. In a study conducted by Sharif et al. (2017) for modeling winter oilseed rape yield, Lasso was among the most accurate regression techniques among several others. This may be attributed to its ability to exclude variables that do not significantly contribute to the final yield. For such scenarios, where the true values of corresponding coefficients are likely near zero, regression techniques with variable selection capabilities often outperform traditional methods (Hastie et al., 2009). However, stepwise regression was less performant in this study, possibly due to the inaccurate contribution of explanatory variables. Stepwise regression is primarily designed for low-dimensional problems. When the number of observations is not substantially larger than the number of explanatory variables, high variance and overfitting can compromise the predictive power of classical approaches. These results were in line with previous studies on regional yield prediction (Cai et al., 2019; Cao et al., 2021; Liu et al., 2022). Kumar  et al. (2021) compared the performance of stepwise and Lasso regression models in wheat (Triticum aestivum L.) yield prediction using meteorological data. Their results demonstrated that Lasso regression was superior to stepwise regression in reducing data dimensionality.


       
The superior performance of the Lasso regression method highlights its potential as a valuable tool in precision agriculture, where high-dimensional datasets are common (Miller et al., 2022; Pukrongta et al., 2024). Accurate yield prediction can guide resource allocation, inform decision-making and support agricultural policies aimed at ensuring food security (Raihan, 2024). For instance, the ability to identify critical predictors from environmental, agronomic and management factors enhances the development of targeted interventions to optimize yield under variable conditions (Kumar et al., 2015). Moreover, the implications extend to sustainable agriculture practices. Farmers and stakeholders can plan more efficiently, minimize input wastage and mitigate risks associated with adverse environmental conditions by accurately predicting their crop yield.
 
Relationships of Faba bean yield prediction parameters under shading conditions
 
To delve deeper into the predictive yield dynamics, Pearson correlation analyses were accordingly conducted in the best possible shading conditions (S2 and S3) (Fig 1). Results revealed that under these conditions in Y1, Faba bean yield (YP) was positive and significantly correlated with NP, NS and WSS parameters (r = 0.35-0.72; P < 0.001), indicating that increasing these traits would increase Faba bean yield. In contrast, significant negative correlations were spotted between plant height and the YP (r = -0.39 to -0.37); P < 0.001) (Fig 1A). Furthermore, the transition from S2 to S3 levels further amplified these correlative effects.        



Regarding the second experimental year (Y2), the variables NP, NS and WSS still showed a significant positive correlation with YP in the S2 level, while a transitional shift was noticed into positive correlations between plant height and YP unlike Y1 (Fig 1B). The significant positive correlations observed between YP and NP, NS and WSS parameters in both Y1 and Y2 highlight the importance of these traits in determining yield. The increase in grain yield was attributed to an increase in the number of pods and seeds due to greater branching and more pods on the branches (Alharbi and Adhikari, 2020). The negative correlation between plant height and YP in Y1, particularly under S2 and S3 conditions suggests that taller plants may not necessarily be more productive in Faba bean. While taller plants often have a larger leaf area, excessive height can lead to increased lodging and reduced light interception at the lower canopy levels. Our results are similar to those of (Wadan and Abd el Shafi, 2014), who showed a negative correlation between plant height of wheat and grain yield. However, for barley (Hordeum vulgare L.), grain yield was positively correlated with plant height under all shading treatments (25% and 30%) (Arenas Corraliza  et al., 2020). The shift in the correlation between plant height and YP from negative in Y1 to positive in Y2 under S2 conditions suggests that the relationship between plant height and yield can be influenced by environmental factors and specific growing conditions. In Y2, taller plants may have been able to utilize the available resources more efficiently, leading to increased yield. However, further research is needed to elucidate the underlying mechanisms driving this change. The number of branches is highly variable in legumes, but is an essential determinant of grain yield (Alharbi and Adhikari, 2020), however, no correlation was reported between YG and NTS under the two shading treatments. Furthermore, the impact of these traits on yield becomes more pronounced under higher levels of shading throughout the seasons as indicated by our study. This could be attributed to the increased competition for resources, such as light and nutrients, under more shaded conditions.
       
The contribution of the most important Faba bean parameters in each shading condition was further explored using an RDA. The results revealed that the variance in shading composition could be largely explained by both axes, with 43.89% of the variance explained by RD1 and 41.56% by RD2 (Fig 2). The third shading level (S3) had 4 major contributors including NP, NS, WSS and YP with the two latter components having more significance overall. Three contributors were detected for S2 (with WSS significance), while only the NP parameter was shown to be significantly relevant to the S1 level among the two detected contributors (Fig 2). The results of the RDA offer insightful contributions to our understanding of how the Faba bean responds to varying shading conditions. Key parameters (NP, NS, WSS and YP) emerged as major contributors, with WSS and YP having greater significance. These results are consistent with the evidence showing that yield-related traits, such as seed weight and total yield, are highly sensitive to reductions in light, with plants adjusting their resource allocation to maximize reproductive success under suboptimal conditions. A study by Wang  et al. (2021) has demonstrated that high shading often induces changes in photosynthate partitioning, favoring reproductive organs over vegetative growth. This could explain the prominence of yield-related traits in the S3 shading condition. The differentiation of parameter importance across shading levels points to an adaptive response by Fab beans to varying light conditions, a concept supported by recent studies on crop plasticity (Mínguez and Rubiales, 2021). The findings reinforce the importance of NP, NS and WSS in determining yield, particularly under higher levels of shading.

In this study, the stepwise and lasso regression models were applied to predict the overall yield of faba bean varieties under 3 shading conditions over two growing seasons. Compared with the stepwise regression method, the Lasso models had higher accuracy and precision in yield prediction based on three main components (NP, NS and WSS). These entities proved to have the strongest correlations with yield, demonstrating their robustness and practical measurement in overall model prediction ability. These results can be directly applied to Moroccan plant breeding programs, providing breeders with precise information on indirect selection criteria. Ultimately, this approach will improve the efficiency and accuracy of breeding efforts aimed at improving faba bean productivity under challenging light conditions, thus contributing to the development of more resistant crop varieties. However, more research has to be established in adopting, testing and validating new yield prediction models (including the ones integrating machine learning approaches) in the field conditions. These advancements will play a pivotal role in developing resilient crop varieties tailored to the environmental constraints of Moroccan and global agroecosystems.

The authors declare that they have no conflict of interest.