Maize Yield Prediction at an Early Developmental Stage Using Multispectral Images and Genotype Data for Preliminary Hybrid Selection

Abstract

Assessing crop production in the field often requires breeders to wait until the end of the season to collect yield-related measurements, limiting the pace of the breeding cycle. Early prediction of crop performance can reduce this constraint by allowing breeders more time to focus on the highest-performing varieties. Here, we present a multimodal deep learning model for predicting the performance of maize (Zea mays) at an early developmental stage, offering the potential to accelerate crop breeding. We employed multispectral images and eight vegetation indices, collected by an uncrewed aerial vehicle approximately 60 days after sowing, over three consecutive growing cycles (2017, 2018 and 2019). The multimodal deep learning approach was used to integrate field management and genotype information with the multispectral data, providing context to the conditions that the plants experienced during the trial. Model performance was assessed using holdout data, in which the model accurately predicted the yield (RMSE 1.07 t/ha, a relative RMSE of 7.60% of 16 t/ha, and R² score 0.73) and identified the majority of high-yielding varieties, outperforming previously published models for early yield prediction. The inclusion of vegetation indices was important for model performance, with a normalized difference vegetation index and green with normalized difference vegetation index contributing the most to model performance. The model provides a decision support tool, identifying promising lines early in the field trial.

1. Introduction

Breeding crop varieties for improved performance in the field is costly and can encompass several years, as the individuals must be selected from large populations grown across multiple environments and years [1]. The prediction of end-of-season traits in the field can bypass the major time limitation imposed by the plant’s long life cycle, by allowing preliminary selection of the most promising individuals based on the plant phenotype at early developmental stages [2,3]. This would enable researchers to select a subset of the varieties to focus on, allowing the resources to be employed on collecting additional phenotyping measurements, sequencing the genome of selected individuals, or defining crossing plans for seed production [4]. High-throughput phenotyping (HTP) allows the monitoring of crop development and measuring plant traits under field or greenhouse conditions, using a wide range of sensors deployed on ground platforms, uncrewed aerial vehicles, or via satellite. The richness of HTP datasets enables breeders to explore features at a scale that was not previously feasible. For example, canopy temperature and a normalized difference vegetation index (NDVI) were shown to increase the accuracy of wheat (Triticum aestivum) genomic selection by 70%, using best linear unbiased prediction models (BLUP) [5]. In addition, plant spectral features collected using thermal, RGB, and multispectral sensors can be used to identify water-stressed individuals [6,7], for disease identification and severity estimation [8,9,10,11], and for lodging quantification [12,13]. For early yield prediction, studies using soybean (Glycine max) and wheat suggest that grain yield, seed size and plant maturity date have a high correlation to spectral features from crops at an early developmental stage [2,3]. In sorghum (Sorghum bicolor), it was observed that dynamic changes in the spectral features throughout the plant’s early development were better for differentiating plant performance than images collected at later stages. This might be due to biomass accumulation and canopy density at later stages, which decrease the sensitivity of vegetation indices [4]. However, defining hand-crafted feature descriptors is a laborious process, and potentially leads to errors as the expert is required to make many decisions as to which aspects of the data better describe the target feature. In addition, the feature descriptors may not be robust enough to be applied to other datasets, restricting their future use.

Deep learning models build multilevel hierarchical representations of the relationship between features, allowing the model to automatically extract features and learn from the datasets [14,15]. Convolutional neural networks (CNN) are the most widely used deep learning architecture for crop yield prediction, followed by deep neural networks (DNN) [16,17,18,19,20]. Recently, a study used MODIS satellite data for the large-scale prediction of maize (Zea mays) and soybean yield. The proposed deep CNN outperformed six other machine and deep learning models and predicted the yield across 2208 counties in the USA with an error of 8.74% [21]. In another study, a deep CNN trained with growers’ field images (RGB and NDVI), was used for predicting wheat and barley (Hordeum vulgare) yield one month after sowing, with an 8.8% of error [22]. Despite the successful application of deep learning for yield prediction, it has not been effectively explored for early yield prediction in a field-trial context where there are few replicate plots from the same hybrid, providing reduced representations of the variety’s phenotypic plasticity in comparison to a grower’s field, where the same variety is sown throughout.

Because plant performance is a result of the individual’s genetic and environmental conditions, the addition of complementary data types, such as weather, soil, genotype, and field management, can enrich the features that will represent the plant varieties for the model, potentially generating more accurate yield predictions within the field-trial context. Several studies have used multiple data types as input for machine learning models for crop yield prediction in the field [23,24]; however, using hand-crafted features may prevent the model from capturing inter- and cross-modality. In most multimodal deep learning models, each deep learning module specializes in a single data type depicting one aspect of the phenomenon, which are then concatenated and used to inform the prediction [25,26]. Multimodal deep learning models have been used for crop yield prediction, incorporating historical yield to weather and genotype data [27], multi-layer soil data and satellite images for yield variation prediction within a wheat field, with around 80% accuracy [28], and fusing canopy structure, temperature and texture for the yield prediction of soybean lines, with an error of 15.9% and R² of 0.72 [29]. Multimodal deep learning can assist in the difficult task of early yield prediction as it accounts for the interaction of many factors that can dramatically influence plant performance.

The objective of this study is to assess the effectiveness of using multimodal deep learning models for early yield prediction of maize hybrids grown under field trial conditions for three consecutive years, using genotype, field management and spectral data collected by UAV. Specifically, the main goals of this study are: (1) to assess if a deep learning model based on tabular features can outperform established machine learning methods for yield prediction; (2) evaluate the impact of using multimodal deep learning for early yield prediction in comparison to unimodal models, based on single-data type (i.e., genotype and field management data); (3) measure the potential contribution of vegetation indices to model performance. As a result, we present a decision support tool for breeders that can assist the selection of maize hybrids at an early developmental stage.

2. Materials and Methods

2.1. Experimental Conditions and Dataset Preparation

Data was collected and made available by the Genomes to Fields (G2F) initiative (www.genomes2fields.org, accessed 25 March 2021) [30], and comprises three consecutive field trials conducted in 2017, 2018 and 2019 in College Station, TX, USA (doi: 10.25739/4ext-5e97, 10.25739/96zv-7164 and 10.25739/d9gt-hv94, accessed on 25 March 2021). The trials were designed using randomized complete blocks with at least two replicates per condition; the plants were sown in two-row plots, with 0.76 m row spacing and 7 m plot length. Commercial check varieties were sown around the trial to form the border. Grain yield was collected from the two-row plots using a plot combine harvester; the harvested yield is used as a ground-truth for model training. Approximately 1500 plots were sown per year, using 1113 maize hybrid lines obtained from the combination of 1040 parental lines. In 2017, maize hybrids were grown with different sowing dates, with individuals sown at the optimal dryland period on 3 March, and others sown at delayed planting on 6 April, referred to as P1 and P2, respectively. The plants grown in 2017 were also subjected to two fertilizer treatments: F1 being the optimal fertilization (using 70–75 gallons per acre of 32:0:0) and F2, reduced fertilization (43 gallons per acre of 32:0:0), whereas, in 2018–19, the plants were managed at uniform optimal conditions. The treatments were named P1F1 (optimal planting and fertilization), P1F2 (optimal planting and reduced fertilization), and P2F1 (delayed planting with optimal fertilization).

Multispectral images were collected using a Micasense RedEdge sensor (Blue: 460–480 nm, Green: 550–560 nm, Red: 660–670 nm, Red Edge: 710–720 nm, and NIR: 840–860 nm) deployed on a UAV platform. The images were collected at 110 m height, around midday; a minimum of twelve ground control points were used to increase GPS accuracy.

2.2. Multispectral Image Processing and Preparation as Model Input

Orthomosaics were assembled using Metashape (v1.6, Agisoft); the reflectance was calibrated for the datasets 2017–18 using the reflectance panel provided by the sensor manufacturer and DLS, the orthomosaics presented an average ground resolution of 7.4 cm/pix. The shapefiles were drawn using the plotshpcreate R library [31] and, using a custom cut_plots_using_shapefile.sh script, ran within a singularity container [32] with Gdal (v3.2). The GeoTiff plots were converted to NumPy arrays, straightened to a vertical position in the frame, and normalized to 0–1 range by dividing the pixel value by 32768, which corresponds to 100% of the reflectance in 16-bit integer files, using a custom python script (msi_processing.ipynb). To prepare the images as input for the deep learning model, the plots were cut in half and stacked horizontally to fit the expected input in the deep learning model architecture, generating a square image of 40 × 40 pixels. The halved plots were randomly rotated and flipped, before being stacked using msi_utils.py. Eight vegetation indices (VIs) were calculated into separate image bands and added to the multispectral image (described in

Table 1. Vegetation indices, included in the plot images; the band number refers to the order in which the VI appears after the original five spectral bands (red, green, blue, NIR and red edge).

Band Number	Identification	Equation
6	Normalized Difference Vegetation Index (NDVI)	$\frac{R_{840}-R_{668}}{R_{840}+R_{668}}$
7	Normalized Difference Vegetation Index Red Edge (NDVI-RE)	$\frac{R_{840}-R_{717}}{R_{840}+R_{717}}$
8	Normalised Difference Red Edge Red (NDRE-R)	$\frac{R_{717}-R_{668}}{R_{717}+R_{668}}$
9	Enhanced Normalized Difference Vegetation Index (ENDVI)	$\frac{\left(\left(R_{840}+R_{560}\right)-\left(2\times R_{475}\right)\right)}{\left(\left(R_{840}-R_{560}\right)+\left(2\times R_{475}\right)\right)}$
10	Canopy Chlorophyll Content Index (CCCI)	$\frac{\left(\frac{R_{840}-R_{717}}{R_{840}+R_{717}}\right)}{\left(\frac{R_{840}-R_{668}}{R_{840}+R_{668}}\right)}$
11	Green Normalized Difference Vegetation Index (GNDVI)	$\frac{R_{840}-R_{560}}{R_{840}+R_{560}}$
12	Green Leaf Index (GLI)	$\frac{\left(2\times R_{560}\right)-R_{475}-R_{668}}{\left(2\times R_{560}\right)+R_{475}+R_{668}}$
13	Optimized Soil-Adjusted Vegetation Index (OSAVI)	$\frac{R_{840}-R_{717}}{R_{840}+R_{717}+0.16}\times \left(1+0.16\right)$

) using msi_utils.py (Figure 1); these VIs were selected based on their reported capacity to estimate plant traits. The code used in this study is publicly available through Github, along with the singularity containers to run the scripts; links can be found in the data availability statement.

The experimental field data included parental lines, seed source (stock), fertilizer application, planting date, and days after sowing, and were processed to be used as training data for the tabular machine learning models. The categorical features were encoded using tab_processing.ipynb that employs the scikit-learn (0.23.2) label encoder for turning categorical features into numerical vectors and one-hot encoders, and the fastai library (v2.1.5) was used for constructing the embedding layer for each feature [33]. The categorical columns, mapped to numerical values, were normalized to fit the range 0–1. There were 5.41% maize plots with missing values from a total of 4487 plots; the majority were filled using the value from their replicate, and 1.18% plots were deleted as the replicate values were also missing. The datasets from 2017, 2018 and 2019 were combined and split into training and holdout sets (90:10). The training dataset was further split into training and validation sets (70:30). The holdout dataset contained data from the three years that were not used during model development but were subject to the same technical and environmental conditions. Results in the holdout dataset in this case indicate the model’s performance for predicting early yield under different environmental conditions.

2.3. Machine Learning Models

2.3.1. Machine Learning Module for Tabular Data

The machine learning modules for yield prediction, based on tabular data, were trained using parental lines, seed source (stock), fertilizer application (optimal or reduced), planting (optimal or late), and the number of days after sowing. Due to the impact of hyperparameters on the performance of machine learning models, the hyperparameters for the random forest, XGBoost, and tabular deep neural network (tab-DNN) were defined after conducting one thousand trials with different hyperparameter combinations, using the Optuna (v2.3) framework. Optuna searches the user-defined hyperparameter space and efficiently discards values with low performance, returning the best hyperparameter set encountered across the trials [34]. For random forest, we searched the minimum sample leaf, the number of estimators and minimal sample split values; XGBoost searched for the optimal learning rate, the number of estimators, maximum depth and minimum child weight; and the DNN searched for the learning rate, optimization function, and the number of layers. Two encoding methods were tested for random forest and XGBoost: ordinal encoding and one-hot-encoding using scikit-learn (0.23.2). The DNN was implemented using fastai v2.1.5 [33] and PyTorch v1.7 [35]. Categorical features were encoded using entity embeddings [36]. Embedding layers learn from the feature distribution in the whole dataset, clustering categorical variables with higher similarity and exposing the intrinsic relationship between them; they are an effective way to map features in an embedding space [36]. The output from the embedding layers was concatenated to the continuous features and was used as input for the DNN (Figure 2). The implementation of the machine learning models for tabular data used tabular_models.ipynb.

2.3.2. Deep Learning Module for Multispectral Imagery

The spectral deep neural network (sp-DNN) for yield prediction was trained using the processed multispectral images, plus the calculated vegetation indices bands, resulting in input images with 13 spectral bands. The sp-DNN uses a modified version of the residual neural network architecture with 18 layers (ResNet18) to process the multispectral images. The ResNet architecture employed comprises optimized convolutional layers, and a self-attention layer [37,38]. The self-attention mechanism helps model the relationship between distant regions of the input image, computing an attention weight matrix that indicates the most relevant regions for the prediction [39]. The hyperparameter search library Optuna (v2.3) [34] was used to determine the sp-DNN architecture: the number of layers, activation function, learning rate, and if attention layers were applied. The general model architecture is shown in Figure 2, and the model script is available at msi_model.ipynb.

2.3.3. Multimodal Deep Learning Framework

The multimodal deep learning model was composed of the tab-DNN and sp-DNN modules described in Section 2.3.1. and Section 2.3.2., followed by a fusion module with two linear layers and ReLU as an activation function. The weights from the last layer from tab-DNN and sp-DNN are concatenated and fed as input for the fusion module. Gradient blending was employed to address the calibration of multiple loss functions for the modules [40], with scaling weights depending on how long each module takes to converge (tabular module 0.1, spectral module 0.55, and fusion module 0.35). Scaling the modules helps to prevent overfitting when the modules are trained simultaneously, as different input data and module architectures are likely to converge at different times [40]. Each module outputs its yield prediction, allowing a comparison of performance against the fusion module prediction (Figure 3). Multimodal model performance was compared, using pre-trained modules against training the architecture from scratch with co-learning. The pre-trained multimodal architecture employed previously trained tab-DNN and sp-DNN as modules, with a further 60 epochs of co-training to fine-tune the fusion module, whereas the multimodal architecture was trained from scratch using initialized tab-DNN and sp-DNN modules, and they were trained along with the fusion module. Scripts used for the implementation of this model are multimodal_model.py, multimodal_utils.py, and multimodal_model.ipynb.

2.4. Evaluation Metrics

Model performance was evaluated using fivefold stratified cross-validation (scikit-learn, 0.23.2). The stratification was performed by year, to ensure that the training and validations datasets were balanced. In stratified cross-validation, the dataset is split between training and validation sets for N consecutive folds, each fold having an equal proportion of samples from the three years. One fold is used once as validation, while the remaining four folds are used for training; this step is repeated until each fold was used as validation once. The metrics employed to measure performance were R² (coefficient of determination), root mean squared error (RMSE) and relative RMSE (RMSE %), using the following equations:

$$\mathrm{RMSE}=\sqrt{{{\displaystyle \sum }}_{i=1}^n\frac{\begin{array}{c}{\left(yi-\widehat{y}\right)}^2\end{array}}{n-1}}$$

(1)

$$\mathrm{RMSE}\%=\frac{RMSE}{ymax-ymin}\times 100$$

(2)

$${\mathrm{R}}^2=1-\frac{{{\displaystyle \sum }}_i{\left(yi-\widehat{y}i\right)}^2}{{{\displaystyle \sum }}_i(yi-\mathrm{y}\text{¯})^2}$$

(3)

where y and ŷ are the observed and predicted yield, respectively, for each sample i, and the total number of samples is indicated by n, whereas ymax–ymin are used to indicate the range of observed yield, and ȳ represents the mean of the observed data. The contribution of the VIs to the multimodal prediction performance was assessed using permutation importance, in which the position of the pixels are randomized in the selected VI, breaking the structured association between the feature and prediction [41,42,43]. To measure the VI importance, we compared the difference in performance on RMSE and R², between the prediction with the randomized VI and the intact VI. A custom python script was used to perform feature assessment: permutation_nb.ipynb.

3. Results

3.1. Maize Grain Yield Distribution

Harvested yield distribution varied across years due to the different seasonal growing conditions. The impact of environmental conditions in each year can be estimated by comparing the performance of nineteen check varieties grown across the three years, under optimal fertilization and planting date conditions (Figure 4a), where the plants grown in 2017 have a higher overall performance than plants grown in the other years. The 2017 trial presents two discernible groups with mean yield values of 10 and 6 t/ha⁻¹, composed of hybrids subjected to P1F1 and P2F1 treatments (Figure 4b). In the 2019 trial, all hybrids were subjected to P1F1 and displayed a consistently low yield compared to previous years, with values ranging from 2 to 11 t/ha⁻¹ although the majority of the plots yielded between 6 to 9 t/ha⁻¹ (Supplementary File S1). In total, 1113 maize hybrid lines were sown from a combination of 1040 unique parent lines (Figure 4c). The under-representation of some parental lines may decrease the accuracy of yield prediction for these hybrids, as there will be a smaller dataset to learn from.

3.2. Machine Learning Hyperparameter Optimization

For the random forest model, the minimum sample leaf was the most important hyperparameter, with 0.91 importance, followed by the number of estimators that define the size of the estimator forest. The optimal number of estimators can vary widely depending on the dataset employed, usually achieving a threshold beyond which there is no increase in accuracy that can justify the computational cost [44]. For XGBoost, the learning rate had the most impact on model performance, with 0.75 importance, followed by the number of estimators. Random forest and XGBoost-optimized model hyperparameters are available in the script tabular_models.ipynb. The tab-DNN and sp-DNN had the same hyperparameter importance order, with the learning rate as the most influential hyperparameter, followed by the optimization function, the number of layers, and the attention layer. The optimized tab-DNN model was composed of a 100-layer network, with ReLU as the activation function, RMSE as the loss function, ranger as the optimization function, a learning rate of 1e-3, batch normalization, dropout (0.5), and an embeddings dropout (0.5). The selection of the optimization function and the number of layers was more relevant to the sp-DNN. The sp-DNN was composed of a ResNet with 18 layers, a self-attention layer, a learning rate of 1e-3, and using Adam as the optimizer [45]. A hyperparameter importance plot is available in Supplementary File S2.

3.3. Prediction of Yield Based on Crop Metadata

Random forest demonstrated the best overall performance for predicting grain yield based on the genotype and crop management data, followed by the tab-DNN (

Table 2. Performance comparison of different machine learning models, using fivefold cross-validation. Values are presented with their standard variation. The best model’s result is highlighted in bold.

Model	Feature Encoding	Validation Dataset			Holdout Dataset
		RMSE	RMSE %	R2	RMSE	RMSE %	R2
Random Forests	Ordinal encoding	1.48 ± 0.05	10.46 ± 0.37	0.53 ± 0.04	1.48	10.59	0.49
	One-hot encoding	1.46 ± 0.08	10.35 ± 0.54	0.55 ± 0.05	1.43	10.11	0.53
XGBoost	Ordinal encoding	2.17 ± 0.05	15.40 ± 0.38	0.00 ± 0.06	2.28	15.77	−0.14
	One-hot encoding	2.21 ± 0.08	15.67 ± 0.57	−0.04 ± 0.05	2.12	15.02	−0.03
DNN	Embeddings	1.67 ± 0.04	11.87 ± 0.33	0.41 ± 0.02	1.77	12.57	0.27

). The machine learning models presented between 0.11–0.27% of relative RMSE variation when predicting yield in the validation and holdout datasets, suggesting that the models were not overfitting the training dataset and performed similarly on unseen data. The three models presented a lower R² score when predicting the yield from the 2019 hybrids. The categorical feature-encoding method had little effect on the prediction accuracy of the validation and holdout sets, with one-hot encoding causing no improvement in the random forest or XGBoost performance. In addition, training the models with one-hot encoded data took longer, as was also observed in a previous study [46]. The range of the predicted grain yield is often more narrow than the actual yield distribution (Figure 5; Supplementary File S3). Although the year feature was hidden from the models, we can observe that random forest and the tab-DNN presented higher error rates for predicting the 2018 and 2019 yields. XGBoost performed worse than other models, grouping hybrids into clusters and providing a single yield value per group, as shown in Figure 5. To assess the capacity of the models to identify high-performance hybrids, each plot was placed into categories based on their performance and labeled as high- (above 10 t/ha), medium- (between 6 and 9.99 t/ha), and low-yield plots (below 6 t/ha). Sorting the hybrids into categories showed that random forest and tab-DNN were able to identify approximately 76% of the high-performance hybrids, but were incapable to identify more than 38% of the low-performance plots (Figure 5G–I).

3.4. Spectral Feature Deep Learning Model Performance

The sp-DNN was trained using thirteen spectral bands, including the eight calculated VIs. The trained model predicted grain yield with 8.99% relative RMSE in the holdout dataset, as shown in

Table 3. Yield prediction performance of sp-DNN using multispectral features. The validation dataset values were obtained through fivefold cross-validation and presented with standard deviation.

	RMSE	RMSE %	R2
Validation dataset	1.46 ± 0.14	10.35 ± 1.01	0.55 ± 0.08
Holdout dataset	1.27	8.99	0.63

. The model predictions accurately reflect the treatment groupings P1F2 and P2F1, with no changes in the image capture and reflectance calibration. The capacity of the model to group P1F2 and P2F1 predictions shows that the spectral features were indicative of the plant conditions. In addition, the bimodal tendency observed in the ground-truth yield is accentuated in the predicted yield (Figure 6). The model R² score was higher for the prediction of 2017 and 2018 hybrids, although 2017 presented a significantly higher error rate than other years (p-value > 0.05), but no significant difference in prediction error was observed between most treatments except for P1F1, compared to P1F2. The sp-DNN model successfully identified the majority of high-yielding varieties (above 10 t/ha), surpassing the tabular machine learning models and tab-DNN (Figure 6E).

3.5. Multimodal Deep Learning Model Performance

Multimodal co-learning enables sharing complementary information from different modalities, influencing model performance even if one of the modalities is absent during test and validation [25]. In our case, we used co-learning to train the modules and used the weights to feed the fusion module. In this framework, each module outputs an individual prediction, based solely on its data type, being completely independent of the sister modality used during training. Co-learning improved the performance of both modules, as shown in

Table 4. Performance comparison of multimodal models with late feature fusion. Validation dataset results were obtained with fivefold cross-validation and presented with standard deviation. The weighted prediction was obtained by scaling each module’s prediction and summing the predicted yield. The best-performing models are highlighted in bold.

Multimodal Framework		Validation Dataset			Holdout Dataset
		RMSE	RMSE %	R2	RMSE	RMSE %	R2
Feature Fusion	Tabular module	1.29 ± 0.20	9.15 ± 1.40	0.64 ± 0.11	1.53	10.84	0.47
	Spectral module	1.27 ± 0.17	9.02 ± 1.22	0.65 ± 0.10	1.14	8.06	0.71
	Fusion module	1.16 ±0.05	8.24 ± 0.38	0.71 ± 0.02	1.17	8.27	0.69
	Weighted prediction	1.13 ± 0.04	8.00 ± 0.31	0.73 ± 0.016	1.07	7.60	0.73
Feature Fusion with pre-trained modules	Tabular module	1.29 ± 0.21	9.10 ± 1.50	0.64 ± 0.12	2.16	15.31	−0.07
	Spectral module	1.27 ± 0.18	8.96 ± 1.31	0.65 ± 0.11	1.19	8.44	0.67
	Fusion module	1.14 ± 0.04	8.09 ± 0.25	0.72 ± 0.02	1.22	8.63	0.66
	Weighted prediction	1.12 ± 0.02	7.90 ± 0.16	0.74 ± 0.01	1.21	8.55	0.66

. The tab-DNN module benefitted substantially and increased the correlation of the prediction in the holdout dataset from 0.27 R² when trained independently to 0.47 R² under the co-learning scheme. The greatest prediction performance was obtained by scaling the predictions of the tabular, spectral and fusion modules, using the weights defined during training (0.5, 0.35, 0.6), achieving 7.60% of relative RMSE (equivalent to 1.22 t/ha⁻¹ RMSE) in the holdout dataset, and around 8% of relative RMSE in the validation dataset. No significant difference in performance was observed using pre-trained modules or training them from scratch, there was also no significant difference in the error rate of predictions for plots grown under any treatment or year (p > 0.05).

The maize plots were placed into categories based on their performance, labeled as high- (above 10 t/ha), medium- (between 6 and 9.99 t/ha), and low-yield plots (below 6 t/ha). Sorting the hybrids into categories showed that the multimodal model was able to identify around 75–86% of the high- and medium-performance hybrid plots and 36% of the low-performance plots, which is expected due to the low RMSE (Figure 7). Further investigation suggests that the majority of underestimated plots were from 2017, whereas 2019 had the most plots expected to perform better than they did. Comparing the prediction error of hybrids sown on multiple years or their replicates shows a small prediction error (Supplementary File S3).

3.6. Effect of Spectral Features on the Prediction of End-of-Season Traits

The impact of different VIs in the multimodal model predictions was measured using feature permutation. The performance of the multimodal model was negatively affected by removing the VI bands as shown in Figure 8. Removing all proposed VIs increased the relative RMSE from 7.60%, obtained with the complete dataset, to 10.84%. The correlation between spectral features and grain yield also decreased substantially from 0.72 to 0.45 R². GNDVI and NDVI were the most influential VIs, with the randomization of the values causing the most impact on the model performance and feature correlation to grain yield prediction. The other VIs had a minor effect on the performance, and the absence of either of them decreased the accuracy of prediction. There was little difference between the prediction performance using RGB- or MSI-only images; however, the majority of VIs used the red edge or near-infrared bands, which may influence the decision to employ a multispectral sensor instead of RGB.

4. Discussion

Plant selection at an early stage is challenging due to the impact of weather variability and the time needed for the plant to show discernible traits; however, studies in maize and soybean crops have shown that early phenotypic traits, collected with HTP platforms, can be used as an indicator of yield under optimal conditions [47,48]. In maize, a previous study using the 2017 G2F dataset identified loci associated with controlling plant height at early development, which showed a significant impact on yield prediction [48]. In our study, the use of early phenotypic data and genotype showed varied results, depending on the machine learning model employed, with the multimodal model presenting a higher prediction accuracy. Assessing the tabular machine learning models independently, Random forest outperformed XGBoost and tab-DNN in this dataset, indicating its effectiveness for crop yield prediction using a single type of data. The random forest method has been extensively applied for crop yield prediction [4,49,50,51,52,53], with previous studies showing it to outperform other machine learning algorithms in selected datasets [50,53]. XGBoost had the worst performance among the machine learning models tested, placing the samples into yield groups. This behavior suggests that the tabular dataset had an insufficient number of features for XGBoost to separate low and high-yielding hybrids, as previous studies showed XGBoost as competent when dealing with larger datasets [54,55]. Tab-DNN was selected as a module for the multimodal architecture, due to its potential for improvement within the co-learning framework. When independently trained using the tabular dataset, tab-DNN showed a higher error rate for the P1F1 treatment, which included samples grown in the three years without the data year being indicted in any feature, which may have caused the seasonal variability to impact prediction accuracy. In P2F1, the yield was consistently lower than other treatments because the hybrids were sown later in the season and experienced heat stress. A study using the same dataset hypothesized that delayed planting in that region of Texas, US, may increase photosynthetic activity, leading to increased plant height and lower yield [48].

Regarding the sp-DNN module, the bimodal tendency observed in the ground-truth yield data is accentuated in the predicted yield (Figure 6), with a larger number of plots being predicted to have a higher yield. The sp-DNN module trained with spectral data from the crop’s early stages presented a higher prediction accuracy with R² 0.63 than the machine learning algorithms based on genotype and crop management (R² 0.53). Several factors may have contributed to this result, as the spectral features are capable of indicating the plant’s health status and biomass coverage [56,57]. Both tab-DNN and sp-DNN modules showed a tendency for more conservative predictions toward the mean, predicting a narrower yield range than ground-truth data.

The multimodal model was able to predict yield across three consecutive years, with a relative RMSE of 7.6% and R² of 0.73. The modules consistently presented a low R2 score for the early yield prediction of 2019 hybrids; this may be due to the parental lines presenting a different performance in previous years, which may introduce a bias to the model. It outperformed previous models based on spectral features for early yield prediction, for both soybean (an RMSE of 9.82 and R² of 0.68) [3] and maize subjected to phosphorus treatments (R² between 0.15~0.82) [58]. The higher accuracy may be due to the capacity of multimodal deep learning to incorporate different data into the prediction. Multimodal deep learning models have been used for crop yield prediction using spectral features from later development stages, where plant growth has a higher correlation with the final yield. Nonetheless, our model had comparable performance, with two multimodal models for soybean yield prediction using images from plant reproductive stages (a relative RMSE of 15.9% and R² of 0.72) [29], and for cotton-yield prediction using images from flowering and shortly before harvest (a relative RMSE of 6.8% and R² of 0.94) [59]. Multimodality deep learning models can support early yield prediction, due to the complementary information shared within the model [25]. Training the spectral and tabular features together has improved both module predictions, increasing tab-DNN performance on the holdout from R² 0.27 to 0.47, and the sp-DNN R² from 0.63 to 0.71 in the holdout dataset. The improvement in accuracy caused by multimodal data is corroborated by previous studies; it was reported that adding spectral features improved the prediction of grain yield in wheat by 70% [5], and in maize, it increased by 0.1 in R² [60]. Here, we demonstrate that a model using two types of input data shows higher accuracy compared to single-data models, one that is valuable when trained using data from young plants with few measurements to inform the model. The multimodal model presented here showed high accuracy (relative RMSE 7.6%) for early yield prediction, being suitable as a decision support tool that can inform breeders of the most productive hybrids. However, it is important to notice that using the model for early yield prediction of another species or under different treatment may require fine-tuning, using the current model as a base for transfer learning [61]. Further training of the model can be achieved by using private or other publicly available datasets [62]. In addition, the multimodal framework can be expanded with the inclusion of new specialized modules for weather, soil and other sensor data, as previous studies show their impact on model accuracy [17,63,64].

Identifying the components that lead to deep learning model prediction is a useful step in verifying model coherence, reducing noise, and improving the overall model [65]. The contribution of the spectral components for the model prediction can be used to extract insights, as the spectral features can capture different aspects of the plant phenotype. In this study, VIs had a high impact on the prediction, causing a substantial decrease in model performance when using RGB or multispectral images alone. RGB images performed better than using multispectral images alone, which was also observed in another study assessing the early growth of maize hybrids under phosphorus treatments [58]. However, the majority of VIs employed the red edge or near-infrared bands, which may influence the decision toward collecting the raw data, using a multispectral sensor instead of RGB. Our observation of the importance of NDVI and GNDVI to predict crop traits is similar to previous studies using VIs. For example, a study using 33 different VIs reported that NDVI, NDRE and GNDVI had the strongest correlation with maize yield [49]. Other studies also demonstrate NDVI correlation with grain yield and moisture, suggesting its use as a secondary trait for plant selection during trials [58,66,67]. For wheat, NDVI provided spatial separation between the high- and low-yielding regions within a plot during the anthesis, tillering, and seedling stages [68]. NDVI and GNDVI are usually related to biomass and tend to saturate at later growth stages [57]. However, since the aim was to predict yield at an early stage, including NDVI and GNDVI was advantageous. Contrary to previous studies, OSAVI only had a minor contribution to yield prediction at an early stage, even though this index is commonly employed for prediction in sparse vegetation as it corrects soil interference in the reflectance [69,70]. The minor effect of OSAVI and other VIs may be explained due to redundancy in the information provided by the other indices, which may partially compensate for the absence of one or another VI. Nonetheless, even OSAVI and other low-contributing VIs increased the model performance to some extent, and it is advantageous to include them in the training data.

5. Conclusions

Early selection of the highest performing crop plants in a field trial can accelerate breeding, reducing the cost and duration of the breeding cycle. Here, we demonstrate a clear advantage of using multiple data types for early yield prediction, and how training specialized modules together can help improve prediction accuracy. Although previous studies have shown that adding extracted spectral features can improve genotype selection accuracy, using a multimodal deep learning algorithm was not previously explored. The multimodal model developed in this study achieved high performance for early yield prediction and could provide a decision support tool for breeders. However, it is important to notice that it was trained using single-location data, and it would require further training with new datasets to be applicable to different locations. In addition, the model accuracy can be improved by regularly updating the model after each season, as the new data can show variation in the hybrids’ performance or integrate new hybrids into the training set. This approach is commonly used in machine learning to ensure that the training dataset does not deviate much from the target dataset we are trying to predict. Further research could explore whether the inclusion of other data types such as weather, soil and genetic information can contribute to improving model accuracy, as they play an important role in plant performance.