Parameterization and Calibration of Wild Blueberry Machine Learning Models to Predict Fruit-Set in the Northeast China Bog Blueberry Agroecosystem

: Data deficiency prevents the development of reliable machine learning models for many agroecosystems, especially those characterized by a dearth of knowledge derived from field data. However, other similar agroecosystems with extensive data resources can be of use. We propose a new predictive modeling approach based upon the concept of transfer learning to solve the problem of data deficiency in predicting productivity of agroecosystems, where productivity is a nonlinear function of various interacting biotic and abiotic factors. We describe the process of building metamodels (machine learning models built and trained on simulation data) from simulations built for one agroecosystem (US wild blueberry) as the source domain, where the data resource is abundant. Metamodels are evaluated and the best metamodel representing the system dynamics is selected. The best metamodel is re-parameterized and calibrated to another agroecosystem (Northeast China bog blueberry) as the target domain where field collected data are lacking. Experimental results showed that our metamodel developed for wild blueberry achieved 78% accuracy in fruit-set prediction for bog blueberry. To demonstrate its usefulness, we applied this calibrated metamodel to investigate the response of bog blueberry to various weather conditions. We found that an 8% reduction in fruit-set of bog blueberry is likely to happen if weather becomes warmer and wetter as predicted by climate models. In addition, southern and eastern production regions will suffer more severe fruit-set decline than the other growing regions. Predictions also suggest that increasing commercially available honeybee densities to 18 bees/m 2 /min, or bumble bee densities to 0.6 bees/m 2 /min, is a viable way to compensate for the predicted 8% climate induced fruit-set decline in the future.


Introduction
Bog blueberry (Vaccinium uliginosum L.), is a Eurasian and North American flowering plant species in the genus Vaccinium within the heath family [1,2]. The commercial value of bog blueberries comes from their antioxidant content, which is the highest of 40 common fruits and vegetables [3]. Bog blueberries are widely distributed and managed in the highlands of Northern China [4] as a major agroecosystem supporting the local rural economy. As of 2019, farmers in three provinces in Northeast China manage 3.5 million hec-well as sharing suitable climatic conditions [27]. The pioneering research on wild blueberry has encouraged us to investigate the possibility of calibrating the predictive models built and validated on wild blueberry for the prediction of bog blueberry fruit-set.
Our primary goal in this research, therefore, was to test the hypothesis that predictive models successfully used for wild blueberry can be adapted to bog blueberry fruit-set prediction with adjusted parameter estimation and validation on a relatively small dataset. Although the two species are very similar, we still are not clear about whether a non-reciprocal pollen compatibility relationship across different genotypes exists for bog blueberry and what the proportion of plants that are self-compatible is in this species [28]. These relationships are critical for estimating the pollination efficiency and consequently predicting fruit-set [15]. We hypothesized that the pollination simulation model built for wild blueberry captures the main structure of pollination dynamics in bog blueberry. This hypothesis is based upon the similar population reproductive structure that is shared by the two blueberry species [15,28,29]. Both species are primarily heterogamous, meaning that outcrossing is essential, but they also share autogamy (a proportion of self-compatible individuals within populations). Based upon this similar reproductive and genetic structure, the metamodels based upon the wild blueberry simulation model should be able to be used to identify bog blueberry specific patterns in pollination from the limited bog blueberry data that we had available for modeling. Once we successfully transferred the metamodel into bog blueberry fruit-set prediction, we applied it as a tool to evaluate the response of Northeast China's bog blueberry agroecosystem to weather factors, which have been regarded as the most important abiotic variables affecting blueberry fruit-set (a proxy to yield) [11,12]. In addition, we also investigated the viable management strategies that could likely offset the negative effects of climate change on bog blueberry productivity.

Overview
Our study combines field observation and computer modeling ( Figure 1) to produce reliable predictions of bog blueberry's fruit-set and utilize them to predict the productivity of Northeast China's bog blueberry agroecosystem in response to various weather conditions and predicted climate change. A simulated dataset previously generated by the wild blueberry pollination simulation model [15] was employed to develop several machine learning algorithms. We use the terms "machine learning algorithm" and "metamodel" interchangeably in this paper throughout the several steps of our analytical approach: feature selection, model formulation and selection. Once the best metamodel with the highest prediction accuracy was selected, its parameters were then estimated by being calibrated to the bog blueberry fruit-set data collected from 6 fields in Northeast China in 2015 and 2016. Then the calibrated metamodel was validated against a second independent bog blueberry fruit-set data collected at the same 6 fields from 2017 to 2019, which had not been used for metamodel construction. After validation, the metamodel was used for three experiments to predict bog blueberry fruit-set in Northeast China under several weather conditions. The geographical pattern of fruit-set change, as well as, potential management strategies to mitigate negative effects of climate change were also investigated as model applications.

Datasets and Simulation Data
Data used for metamodel development was downloaded from a publicly accessible wild blueberry fruit-set dataset [30], see Mendeley Data: https://data.mendeley.com/datasets/p5hvjzsvn8/1, accessed on 25 August, 2021, which was generated by the Wild Blueberry Pollination Simulation Model [15]. The dataset contains 777 records of simulated wild blueberry fruit-set, each of which is an average of 100 replications of simulation experiments. These simulations were conducted by varying the following parameters: the average clone size (CS) in square meters within a wild blueberry field; pollinator density (common taxa) in bees per square meter per minute for Honeybee (HB), Bumble bee (BB), Andrena (AD) and Osmia (OS); the highest point (MaxUTR), the lowest point (MinUTR) and the average (AvUTR) of the upper range of the daily air temperature along the blueberry production season; the highest point (MaxLTR), the lowest point (MinLTR), the average (AvLTR) of the lower range of the daily air temperature during the blueberry production season; the total number of rainy days (RD, or precipitation are used interchangeably in this paper) during the bloom season when the daily precipitation is higher than 2.5 cm; and the average number of days that rain (AvRD) during the bloom season. Table  1 summarizes the range of parameter values that the simulations covered. These ranges represent typical wild blueberry spatial traits, bee pollinator density and four possible climatic trends around current weather conditions, which were Warm and Dry, Warm and Wet, Cool and Dry, Cool and Wet, respectively. In other words, the simulation data contain patterns that reveal the relationships between fruit-set and the features of plant, pollinator, and weather conditions.  Figure 3) in six locations (see Section 2.1 Study area) [25]. Fruit-set was calculated as the total number of viable fruits in one clone (genet or genetically unique plant) at harvest divided by the number of flower buds before bloom in the same clone. Due to uncertainties of the long duration of these field observations, sample sizes in each year and at each field site varied from 58 to 300. Unfortunately, no plant spatial traits or bee visitation rates had been recorded together with the fruit-set observations. In addition to fruit-set, historical weather data including precipitation and air temperature from 2015 to 2019 were collected from the China National Meteorological Center (http://data.cma.cn, accessed on 25 August, 2021) at the geographical location of each field site. The six field observations on bog blueberry fruit-set are included in Appendix A.

Feature Selection
In the process of developing a metamodel for fruit-set prediction, there may be features that are irrelevant or of minor importance. The contribution of some of the types of features can be determined to be minimal as compared to the most important features obtained as the result of feature selection. The total 13 features which were obtained from the publicly available dataset [30] of the Wild Blueberry Simulation Model have a different unit and value range and therefore, we normalized them into a common scale (e.g., between 0 and 1) without distorting differences in the ranges of values. We then used five machine learning methods that have their own built-in feature selection function: Decision Tree (DT), Random Forest (RF), Gradient Boost Decision Tree (GBDT), AdaBoost, and Extreme Gradient Boosting (XGBoost). All five methods were employed to select the subset of features having the highest coefficient of determination (R 2 ). Statistically, the coefficient of determination (R 2 ) explains the proportion of variance in the predicted result that is explained by the features. We iteratively tested the coefficient of determination of the five machine learning methods on all subsets of features, in which the number of features in a subset increased from 1 to 13 [31]. It was observed that the XGBoost method achieved the highest R 2 value (0.929) among the five machine learning algorithms when the number of features reached eight (8) (Figure 4). The eight (8) feature subset selected by XGBoost were CS, HB, BB, AD, OS, MaxUTR, MinUTR, and RD. These eight (8) features were used for developing the metamodels.

Model Development
The Python 3.6 software development environment was employed to develop nine metamodels from the simulated wild blueberry dataset [30], in which wild blueberry fruitset was the dependent (response) variable, and the eight (8) selected features (see Section 2.4) were used as independent (explanatory) variables.

Multiple Linear Regression (MLR)
MLR is a statistical modeling technique that involves predicting a numeric value given multiple independent variables. MLR models have been used extensively in the agricultural research field [32] to develop predictive models that assume a linear relationship between more than one explanatory variable and a response variable by fitting an additive linear equation to observed data. In this study, the goal of using multiple linear regression is to model the linear relationship between the eight (8) selected features and wild blueberry fruit-set, which is given by the general form as in Equation (1): where is wild blueberry fruit-set in the ith sample, i.e., the dependent variable; corresponds to the independent variables (CS, HB, BB, AD, OS, MaxUTR, MinUTR, and RD); is the equation intercept; is the corresponding linear coefficient of the k th independent variable; and is random error. The training process is responsible for determining the most likely coefficients for the MLR model from learning patterns in the dataset (see Section 2.3).

Support Vector Regression (SVR)
The Support Vector Machine (SVM) algorithm was first introduced by Vladimir Vapnik and his colleagues [33] and is a class of supervised machine learning algorithms for classification problems. Later on, Drucker et al. [34] proposed Support vector regression based on the concept of Vapnik's support vectors. The main purpose of SVM regression is to represent complex relationships through nonlinear mapping. Residual error is minimized by adding a hyperplane and maximizing the distance between predicted and observed values. Initially, the variables are modeled from the original space to a high-dimensional feature space by a kernel function, where the kernel functions can be (linear, polynomial, Gaussian, etc.) and depend on the relationship between the explanatory and response variables. Then, a linear model is built in the derived feature space to minimize error, where it becomes linearly separable [35]. The main hyperparameter fine-tuning of the SVM model includes a kernel function, cost parameter(C), gamma (γ) and the impact of regularization. In this study, the performance of different kernel functions was compared, and finally a Gaussian kernel function was selected. This significantly reduces the risk of overfitting and improves generalization.

K-Nearest Neighbor (KNN)
Initially, k-NN was used only for classification. However, over the past few decades, this model has also been used for both classification and regression modeling. The k-NN algorithm is a nonlinear instance-based machine learning method [36] which is based on the distance of the predictor variables to the nearest training group known to the model [37]. When k-NN is applied to solve regression problems, the value of the response is calculated as a weighted sum of the responses of all the k neighbors, in which the weight is inversely proportional to the distance from the input record [38]. The neighbors' distance can be calculated by Euclidean, Manhattan, and Minkowski distance formulas. However, in this study, the Minkowski distance formula was chosen by comparing the performance of different distance measurement metrics, given by Equation (2): where is the number of nearest neighbors, and are the distance between two points, and is a real value between 1 and 2.

Decision Trees (DT)
Decision Trees (DTs) are non-parametric algorithms that fall under the category of supervised machine learning algorithms and can handle large and complex datasets effectively without a complex parameter structure [39]. Decision tree regression assesses the features of an object and trains a model in the tree structure to predict data to produce meaningful results in the future. In this study, we used the C4.5 algorithm [40] to train the DT model for predicting wild blueberry fruit-set.

Random Forest (RF)
Random Forest (RF) first introduced by Breiman [41], is one of the most widely employed ensemble of machine learning methods which create multiple regression trees that are generated by a large set of decision trees for computing regression models [42]. When compared to using a single decision tree which often creates an unstable model, RF makes predictions by combining predictions from several decision trees using Bootstrap aggregation or a Bagging technique [43]. Bootstrapping involves random sampling of data with replacement and has been proved to effectively reduce and control variance (overfitting) of a predictive model. Random forest training includes training for each decision tree on a randomly selected subset of features and data. Then, the final prediction result is obtained by majority vote or averaging the outputs from each of the sub-models. This can significantly improve the predictability of the RF model in terms of accuracy and generalizability. With respect to handling noisy data, Random forests are more robust. In addition, Random Forests perform well at capturing tabular data with numeric features and maintaining nonlinear interactions between the response variable and the predictors [44].

Adaptive Boosting (AdaBoost)
AdaBoost regression is a machine learning meta-algorithm proposed by [45] that begins with fitting a regressor on the training dataset and then adds additional replications of the original regressor on the same dataset, but with performance improvement based on error information collected from its predecessors [46]. In this way, a final model based upon a strong learning basis can be expected. To improve performance, the algorithm can be used in conjunction with many other learning algorithms.

Gradient Boosting Decision Tree (GBDT)
Jerome Friedman [47] in 1999 first introduced GBDT. This method is a popular machine learning algorithm that optimizes the predicted value of a model in a series of steps during the learning process. With the aim of minimizing the loss function, every single iteration of the decision tree involves adjusting the values of the coefficients, weights, or biases employed on each of the input variables used to predict the target value. There are four main benefits to applying GBDT for predictive model development: (1) feature selection is inherently performed during the learning process; (2) not prone to use of collinear or identical features; (3) models are relatively easy to interpret, and (4) it is easy to specify different loss functions. Hence, in order to solve predictive problems in both classification and regression domains, GBDT is among the most widely used methods used in machine learning.

Artificial Neural Networks (ANN)
Artificial Neural Networks represents a network model for emulating a biological neural system which consists of input, hidden, and output layers [48]. The input layer directly receives the data, whereas the output layer produces the required output. The layers between the input and output layers are hidden layers where the intermediate computation takes place. The process of creating a neural network begins with the perceptron. The perceptron receives the inputs which correspond to the independent variables, and the neurons in the hidden layer multiplies the inputs by weights based upon multiple passes of the data. The results are then transferred to the output layer via an activation function (such as relu, tanh, sigmoid) [49,50]. In this study, we adopted a three-layer perceptron ANN having the activation function "ReLu" (because of its computational efficiency) [31] in the hidden layers, which uses Adam as the gradient descent method. The activation function used in our study is the non-linear rectified linear unit (ReLu) function, which has output 0 for input less than 0 and raw output otherwise. It is mathematically expressed in Equation (3).
In recent years, many agricultural researchers have been using XGBoost to predict crop yields [30]. The eXtreme Gradient Boosting (XGBoost) method is an ensemble machine learning algorithm based on the principles of the gradient boosting decision tree method. With this method boosted trees are efficiently constructed and can automatically operate parallel computation 10 times faster than gradient boosting [51]. It supports supervised machine learning algorithms, including classification, regression and ranking. For better performance, XGBoost provides three additional features compared to GBDT. First, the weights of each new tree can be scaled down by given constant leaves which reduces the influence of a single tree on the final score. Second, introduction of a regularized loss function avoids the problem of overfitting. Third, the Taylor expansion method is used to approximate the loss function thereby speeding up the process of optimization [52]. The predicted values of the tree ensemble model are computed using the formulae in Equations (4) and (5).
Equation (4) represents the regression tree space where is the input, is the output, the tree structure, and the number of leaves in the tree are identified as . Each matches the independent tree structure q and the weight , where represents the score of the leaf. This predicted value can be evaluated by: where, Ω( ) = + ‖ ‖ is the loss function which is the difference between the predicted value and the response value y . The entity Ω represents the complexity of the model and is a regularization term that has the function of smoothing the weight to avoid overfitting.

Model Evaluation and Selection
We employed a widely used statistical approach, cross-validation, to address the problem of overfitting, in which we built models using 5-fold cross validation methodology. The five folds were divided into training (four-folds) and testing (one-fold) sets. The cross-validation was repeated in n rounds or trials. In each round the wild blueberry simulation dataset was partitioned into a complementary subset. The model estimation was conducted on one subset (the training set), while the validation was conducted on the other subset (the validation or test set). The wild blueberry simulation dataset was partitioned differently in different rounds, and the validation results were averaged over the n rounds as an estimation of the model's prediction accuracy.
In this study, to examine the prediction accuracy of the nine metamodels for wild blueberry fruit-set predictions, four different statistical evaluation metrics were employed. First, the coefficient of determination (R 2 ) was used, which is defined as the proportion of the variance in the response variable that is explained by the independent variables. Second, we used the mean absolute error (MAE), which is defined as the absolute mean difference between the actual fruit-set and the predicted fruit-set . Third, we used the mean average percentage error (MAPE), defined as the prediction accuracy as an average of ratios, each of which is the prediction error − divided by the actual value . Fourth, the root mean squared error (RMSE) was used and is defined as a measure of the quadratic mean difference between predicted and actual fruitset. The mean absolute error (MAE), mean average percentage error (MAPE) and root mean squared error (RMSE) were computed using the Equations (6)-(8); respectively: where is the number of samples, is the observed fruit-set and is the model-predicted fruit-set. The coefficient of determination (R 2 ) was computed according to [53] and is given by Equation (9): where , and n are as defined above, and is the mean observed wild blueberry fruit-set. The metamodel with the lowest prediction error and the highest R 2 was selected as the best metamodel.

Model Calibration and Validation
When the best metamodel was selected as the candidate for wild blueberry fruit-set prediction, it needed to be further calibrated to fit the reproductive dynamics of bog blueberry before specific predictions could be made for bog blueberry. The best metamodel trained on the wild blueberry simulation dataset was regarded as the best representation of the relationships between fruit-set and the eight variables by learning the basic pattern of pollination dynamics produced by the Wild Blueberry Pollination Simulation Model. We hypothesized that the selected metamodel has the potential to describe a similar pattern of pollination dynamics in the bog blueberry agroecosystem [30]. This hypothesis was tested in two steps.
Step one was based upon exploring and estimating the best parameter set of the best metamodel for bog blueberry.
Step two was based upon comparing the predictions with bog blueberry fruit-set not used previously in order to examine the goodness of fit.
The best metamodel of wild blueberry prediction had eight parameters (or features, including clone size, bee taxon-specific densities and weather) that determine fruit-set. However, the prediction only needs weather data from bog blueberry fields as input. Other parameters such as clone size, densities of different bee taxa are unknown and need to be predetermined before fruit-set prediction can be made for the bog blueberry agroecosystems. Therefore, to calibrate the best metamodel to bog blueberry, it was necessary to find optimal values for the subset of parameters that (1) minimized the difference between the predicted bog blueberry fruit-set and the observed of bog blueberry fruit-set on six bog blueberry fields in years 2015 and 2016; and simultaneously (2) maximized the total coefficient of determination (R 2 , see Equation (9)) of these parameters that explain the variance of the model. The first objective can be mathematically formulated as: where , ( ) is a metric function that evaluates the prediction error, which was specified as the RMSE in equation 8; M and N are the number of fields (i.e., (6)) and the number of years (i.e., (2)) being predicted, respectively; θ denotes the following subset of parameters selected for optimization from the metamodel: Clone size (CS), density of Honeybees (HB), Bumble bees (BB), Andrena bees (AD) and Osmia bees (OS). The second objective can be formulated as: where , ( ) is the coefficient of determination (R 2 ) of θ that explains the variance from the metamodel, which is defined in Equation (9). Finally, we looked for the set of values for the 5 parameters θ that optimize both objectives, which can be formulated as a multiobjective problem: A multi-objective optimization method [54] was employed to search the best parameter set for the metamodel that is expected to achieve the highest prediction accuracy and interpretability for bog blueberry fruit-set. To minimize the random effect, we replicated the optimization algorithm 100 times, each of which was considered to have reached convergence when both of the objectives did not make further improvement in values smaller than 10 −6 upon successive iterations.
Once the best values of the five parameters (except for the three weather relevant parameters) were found for the metamodel, its predictions were validated on the mean fruit-set observed in the same six bog blueberry fields in the years 2017-2019. Our conceptual approach was that if the validation results were acceptable, it indicated that the metamodel had successfully learned the difference in pollination patterns between wild blueberry and bog blueberry. The metamodel could then be regarded as successfully calibrated to bog blueberries and would be ready for making realistic field predictions.

Prediction Applications
We conducted three prediction experiments to evaluate how different weather conditions expected in the near future (5~10 years) might affect bog blueberry fruit-set in northeastern China. Specifically, we used the calibrated metamodel, i.e., XGBoost, to make predictions for bog blueberry fruit-set in the 151 growth regions with several combina-tions of temperature and precipitation variation while keeping other parameters unchanged. Then the variation in predicted fruit-set of these regions were analyzed to determine if a geographical pattern existed. Finally, if specific weather conditions caused a decline in fruit-set, several management efforts which were available to be adjusted in the metamodel, such as increasing commercial honeybee or bumble bee densities, were tested for potential mitigation of the fruit-set decline under the stress of weather factors.

Bog Blueberry Fruit-Set in Various Weather Conditions during Bloom
According to a recent climate research report, the most likely climate change direction in Northeast China in the near future (5~10 years) would be higher air temperatures with more rainy days during the bog blueberry production season. These conditions may pose a threat to bog blueberry's fruit-set. The primary goal of this experiment was to estimate bog blueberry's fruit-set under higher temperature and precipitation conditions. However, the opposite weather conditions were also considered in this experiment to elucidate the full effects of a highly variable fluctuation in weather on bog blueberry crop production. Specifically, we conducted a full factorial prediction experiment, in which the two factors (i.e., MaxUTR and RD) were systematically varied on seven levels (decreased or increased by 5%, 15% and 25% with respect to the baseline scenario indicating the recent weather conditions between 2017 and 2019). Therefore, there were 7*7 = 49 combinations of predictions made. The predicted fruit-set of the 151 bog blueberry growth regions in Northeast China were plotted and projected onto a GIS map with ArcGIS [55], from which the variation of predicted bog blueberry fruit-set can be observed along a gradient of weather change.

Geographical Pattern of Bog Blueberry Fruit-Set Change
We were particularly interested in whether the change in weather during bloom evenly or irregularly affects bog blueberry fruit-set across the 151 growth regions in Northeast China. If the effects were irregular, we expected to observe geographical patterns in response. For example, are specific parts of the growing regions more likely affected by weather change than others? It is important to evaluate the potential for bog blueberry production along a geographical gradient under a changing climate and provide insight into management and conservation. Specifically, we used the predictions generated in the previous experiment to analyze the changes in correlation between fruit-set and the longitudes and latitudes of the predicted regions when the air temperature positively or negatively 5%, 15% and 25% deviated from the current climate baseline (2017-2019). We also conducted a similar analysis for scenarios in which the number of rainy days during bloom were changed in the same proportional pattern.

Bee Density Enhancement for Compensating Bog Blueberry Fruit-Set Decline
Some weather conditions may negatively affect pollinator foraging activities, which could cause a decline in bog blueberry fruit-set. We were especially interested in whether some factors that are relevant to bog blueberry production can be manipulated to increase fruit-set. Hopefully these efforts can compensate or mitigate a decline in fruit-set caused by inappropriate weather conditions. Of the three groups of eight features (parameters of the XGBoost model), weather (RD, MaxUTR and MinUTR), plant spatial traits (CS) and bee density (HB, BB, OS and AD); only bee density is subject to human intervention. Commercial honeybee and bumble bees are available in many regions in Northeast China. Therefore, it is feasible to increase the density of honeybees or bumble bees without reducing the density of wild bees for the purpose of enhancing the pollination service in a field. By exploring the parameter space of the XGBoost model, we aimed at finding the optimal density for honeybees or bumble bees that can independently compensate for fruit-set decline. Specifically, we made two sets of predictions with the metamodel to fit the relationships between bee density and bog blueberry fruit-set while keeping other factors unchanged in accordance with the weather conditions where fruit-set was negatively affected. Then we quantitatively estimated the percentage of extra honeybee or bumble bee investment necessary to mitigate the negative effect of climate change.

Metamodel Formulation and Selection
Nine metamodels were constructed and their performances in terms of prediction ability were evaluated on the simulated wild blueberry pollination dataset with a fivefold cross-validation method. The three best metamodels that had the lowest MAE, MAPE, and RMSE and the highest R 2 were: RF (Random Forest), GBDT (Gradient Boosting Decision Tree), and XGBoost (eXtreme Gradient Boosting) ( Table 2). The remaining six metamodels were not considered further in future metamodeling steps. Even though RF and GBDT had almost the same level of prediction error as XGBoost, we selected XGBoost as the best metamodel and used it for bog blueberry fruit-set prediction because it had a higher R 2 (0.932) than RF (0.922) and GBDT (0.927). The XGBoost metamodel had the lowest unexplained variance in wild blueberry fruit-set among the three-best performed metamodels. The highest R 2 achieved by XGBoost has given the best model generalization capability in dealing with unexpected patterns in validation data not used to construct the metamodels. The scatter plots ( Figure 5) illustrate how accurate the fruit-set in the validation data can be predicted by different metamodels trained in the training set of the wild blueberry pollination dataset.  In addition to model prediction performance, feature importance in the three best metamodels were also evaluated to understand the influence of potential predictive factors on blueberry fruit-set. The importance of features was estimated by the coefficient statistics between each feature and fruit-set and then standardized to 100% as a relative importance measure. Results showed that CS (clone size), RD (Raining days), and Max-UTR (the highest point of the upper range of the daily air temperature along the blueberry product season) were the three most important factors influencing blueberry fruit-set across the three best metamodels ( Figure 6). It was not surprising that CS was the most important factor affecting blueberry fruit-set due to the outcrossing nature of blueberry (i.e., differential outcrossing success depending upon clone genetics of sire and recipient [28]) and bee foraging behavior, which is confirmed by previous empirical [6] and simulation research results [15]. However, both wild blueberry and bog blueberry are wild species that are not planted, so farmers have no control over this feature. The two weather parameters, precipitation and air temperature showed the next highest influence on blueberry fruit-set.

Metamodel Calibration and Validation
The XGBoost metamodel outperformed the remaining eight metamodels in previous training and evaluation steps on the wild blueberry pollination dataset. It was therefore selected as the candidate model to predict bog blueberry fruit-set since it had learned the basic dynamics of pollination and fruit-set for blueberry species. The XGBoost model was then calibrated to the empirical bog blueberry fruit-set dataset collected at the six field sites (CBS, GHS, JMS, MH, TH and YC) in years 2015 and 2016. The calibration process is a multi-objective optimization [55] that simultaneously minimizes prediction error and maximizes interpretability by finely tuning the parameters. Figure 7A,C,E,G,I shows the predicted versus the observed bog blueberry fruit-set for the six fields in years 2015 and 2016 without calibration, i.e., the XGBoost metamodel trained on wild blueberry data were directly used in bog blueberry fruit-set prediction. After calibration, the average accuracy of bog blueberry fruit-set prediction increased from 0.658 to 0.975, as shown in Figure  7B,D,F,H,J, in which the prediction accuracy is defined as the coefficient of determination (R 2 ) of the linear regression between predicted and observed fruit-set. The visualizations in Figure 7 also shows that specific ranges of parameters, such as clone size and bee density, are also the optimal settings for bog blueberry. Table 3 presents the means and 95% confidence intervals of the optimal parameter set for the XGBoost model that had been calibrated to bog blueberry observations. The optimal parameter set obtained on the bog blueberry dataset showed a slight difference in contrast with the original version trained and tested on the wild blueberry dataset. Clone size shrunk from 18.78 to 15.32; honeybee density increased from 0.417 to 1.026; bumble bee density decreased from 0.282 to 0.232, etc. These variations of model parameters likely reflect the difference of plant physiology, phenology, and genetics in vegetative and reproductive growth between wild blueberry and bog blueberry. Table 3. The optimal parameter set of the best metamodel found by the multi-objective optimization algorithm [54] when XGBoost was calibrated to the observed fruit-set of the six bog blueberry field sites (CBS, GHS, JMS, MH, TH, and YC) in the years 2015 and 2016.  The calibrated XGBoost model with the optimal parameter set was validated on bog blueberry fruit-set data collected at the six field sites (CBS, GHS, JMS, MH, TH and YC), but in different years between 2017 and 2019. Table 4 shows that 14 out of 18 bog blueberry fruit-set predictions made by XGBoost fell into the 95% confidence interval of the field observations, which achieved an overall 78% prediction accuracy. If we consider the 4 missed predictions (Wilcoxon Signed Rank test, p < 0.05) further, the worst case (CBS in year 2018) had only a 4% prediction error between the lower bound of the 95% confidence interval (0.504) and the XGBoost prediction (0.482); while the best case (MH in year 2018) had less than 1% prediction error between the lower bound of the 95% confidence interval (0.509) and the XGBoost prediction (0.504). This evidence suggests a successful model validation if one takes into account the general acceptance rate of model prediction [56]. Therefore, the XGBoost metamodel can be adapted and trusted in bog blueberry fruit-set predictions.

Bog Blueberry Fruit-Set Predictions
The first experiment revealed that under the current weather conditions (baseline scenario (T + R)) in which the variation of air temperature and precipitation was zero, the 151 bog blueberry growth regions had a mean fruit-set of 0.403 and a median of 0.374. The highest fruit-set prediction recorded was 0.521 in Oroqin county (50. 6° N, 123.73° E). The lowest fruit-set prediction was 0.373 and recorded in multiple regions including Horqin county (46.36° N, 121.13° E). As expected, higher temperature significantly drove the fruitset decline in most regions (ANOVA, F = 6.360, df = 603, p < 0.001). The rise of air temperature alone by 5, 15 and 25% caused 3.5% (mean = 0.389), 4.2% (0.386) and 4.2% (0.386) decline in bog blueberry fruit-set in the 151 regions ( Figure 8A). The lowest fruit-set recorded was 0.317 in Sunwu county (49.25° N, 127.2° E) in both scenarios where air temperature was raised by 15 and 25%. Surprisingly, fruit-set in currently warmer weather did not drop by an equivalent percentage as the increase in air temperature. We also did not find further fruit-set decline when the air temperature in these regions was already 15% higher than current weather conditions (Tukey-Kramer HSD test, p = 0.866). This indicates a nonlinear relationship between bog blueberry fruit-set and possible warmer weather conditions. The opposite direction in weather change, i.e., when the air temperature in these regions declined by 5, 15 and 25%, promoted bog blueberry fruit-set by 7.7% (0.434), 21.8% (0.491) and 18.9% (0.479) in the 151 regions (ANOVA, F = 200.723, df = 603, p < 0.001, Figure 8A). If the weather became dryer and the air temperature remained unchanged, the average fruit-sets of the 151 bog blueberry growth regions were increased by 4.2, 13.9, and 20.8% when the number of rainy days decreased by 5, 15, and 25% ( Figure 8B), which were significantly higher than those fruit-sets associated with the baseline weather conditions (Turkey-Kramer HSD, p < 0.001). The highest fruit-set recorded was 0.528 in Manchuria county (49. 35° N, 117.19° E). On the contrary, if the weather was more wet, the average fruit-sets in the 151 bog blueberry growing regions showed 0.2, 4.4, and 8.2% decline in contrast with the baseline scenario when the number of rainy days increased by 5, 15, and 25% ( Figure 8B). It appears that only a 5% precipitation increase in the future will not (Turkey-Kramer HSD, p = 0.999) cause a decline in bog blueberry fruit-set. The lowest fruit-set in wetter weather conditions was 0.317 in 28 counties.
Overall, increased air temperature and/or precipitation posed a larger magnitude of threat to bog blueberry fruit-set than that the benefit of increased fruit-set brought about by an opposite weather direction of decreased air temperature or precipitation (see Appendix A, Figure A1). The nonsymmetrical effects of temperature or precipitation variation indicated that the most likely trend of future weather may cause severe damage to bog blueberries as we initially expected. Most bog blueberry production regions in Northeast China will suffer fruit-set declines if the weather continues to warm or become more wet, a model prediction of future climate change.
Predictions also provided evidence of interactions between air temperature and precipitation (ANOVA, df = 36, F = 15.821, p < 0.001). The highest fruit-set was 0.562 in Manchuria county (49.35° N, 117.19° E) under the cooler and dryer weather conditions (air temperature was 15% lower, and precipitation was 25% lower than the baseline situation). This record was also the highest fruit-set value in our predictions for all possible combinations of changes in temperature and precipitation. It suggests that cooler and dryer weather might be ideal for bog blueberry reproduction and harvestable yields. The lowest fruit-set was recorded as 0.317 in more than 60 counties under the warmer and wetter weather scenario (air temperature and precipitation were both increased 25% over baseline), which is much worse than the negative effect caused by increased air temperature or precipitation alone. This is a warning signal to future bog blueberry production in Northeast China, again as warmer more rainy springs are to be expected. Bog blueberry fruit-set predictions in the 151 regions under different weather variation scenarios are geographically shown in Figures 9 and 10, in which fruit-set of counties were categorized into three levels (less than 0.4, between 0.4 and 0.5, higher than 0.5).
The second experiment found that overall, bog blueberry fruit-set is positively correlated with latitude of the production regions. This is true for both field observations and metamodel predictions. Northern regions had higher fruit-set than southern ones. But no apparent trend was detected in this relationship when air temperatures varied while precipitation was static (see Appendix A, Table A1, p = 0.595). This finding suggests that warmer or cooler weather in the near future is not likely to change the geographical pattern of the locations of the highest bog blueberry productivity in Northeast China. Surprisingly, a strong trend in the correlations between latitude and the change in fruit-set was detected when precipitation was varied (Table A2, R2 = 0.833, p < 0.001). The higher the precipitation level, the stronger the correlation between latitude and fruit-set. The trend indicated that if the level of precipitation goes up in the future, the northern bog blueberry production regions (such as Mohe, Tahe, and Oroqin counties) will have an advantage in achieving high bog blueberry fruit-set compared to the southern regions (such as Changbaishan, Shulan and Tonghua counties). However, if weather becomes dryer, the southern regions will benefit more from optimal weather growing conditions and the difference in fruit-set levels between northern and southern regions will shrink.
The decline in bog blueberry fruit-set due to climate induced increases in air temperature and rainy days can be estimated at 8%, from 0.403 down to 0.37 ( Figure 8B). According to the third prediction experiment, the 8% decline in fruit-set can be compensated by intentionally increasing honeybee density from the current 1.0 bees/m 2 /min to 18 bees/m 2 /min ( Figure 11A). If all other factors are fixed and only honeybee density is increased, the 8% of bog blueberry fruit-set that has to be compensated for requires much higher honeybee density than a scenario where all factors are fully interactive and able to be manipulated. However, due to the recognized high efficiency of bumble bees as pollinators for blueberry, increasing commercial bumble bee density alone from the current 0.232 bees/m 2 /minute to 0.568 bees/m 2 /minute ( Figure 11B) can enhance bog blueberry fruit-set by 8%. One should note that even though the more efficient bumble bee can offset climate induced fruit-set decline, a bottleneck exists since our predictions show that it is almost impossible to have higher bog blueberry fruit-set than 0.60 by only introducing bumble bees (bumble bee density higher than 0.65 bees/m 2 /minute) because predicted proportion fruit set increases at a decreasing rate with increasing bumble bee forager density.

Modelling Approach and Interpretation
This research contributed to our scientific knowledge with a novel modeling approach to predicting the fruit-set of Northeast China's bog blueberry based upon a very limited amount of data for model development. As a major agroecosystem in Northeast China; bog blueberry's yield, fruit quality and economic stability are important to the local rural economy. Machine learning models that capture the complex nonlinear relationships between fruit-set and various interacting biotic and abiotic factors are essential for developing adaptive bog blueberry management strategies. However, due to the lack of systematic field experiments and observations on bog blueberries in Northeast China, it was unlikely that in the next few years agricultural scientists could directly develop and train reliable machine learning models for prediction purposes. We therefore utilized previously developed and validated computer simulation models based on US wild blueberry production as a reliable data generator to feed nine machine learning models. The best machine model with the highest prediction performance from the previous step was parameterized and calibrated on a small bog blueberry data set (less than 1500 samples). After calibration, the metamodel was used for bog blueberry fruit-set prediction. This approach allows machine learning models to be adapted to new prediction tasks even though the available data for retraining is limited. The theoretical foundation for this approach is well-known and referred to as transfer learning [57]. The approach provides a general machine learning framework where a model trained (aka knowledge gained) in one domain of interest can be applied (aka transferred) to a new but similar domain where data is usually limited. In recent years, transfer learning has been widely used in deep learning models for crop disease detection [58][59][60], but it is not common for non-deep learning models [61], which may need more fine-tuning with multi-objective optimization techniques [54]. However, transfer learning would work with any type of machine learning algorithms where the base features are the same [62] with different magnitudes of interpretability [63]. In our research, the best machine learning model trained on the wild blueberry simulation dataset was regarded as the best representation of the relationships between fruit-set and the eight predictor variables by learning the basic pattern of pollination dynamics of multiple blueberry species. Although bog blueberry may have diverse genotypes (e.g., the proportion of self-compatible genotypes, [64]) that could cause different pollination dynamics, these differences can still be explained by a metamodel with reparameterization on a small amount of bog blueberry data. The difference of decomposed feature-level variance before and after metamodel calibration (the visualizations in Figure  7) might shed light on this point. Unlike the deep neural network-based transfer learning approach [65], in which only the very few last layers are replaced and retrained while front layers are kept as storage of common knowledge between the source and the target domain [57], our approach uses a similar concept of knowledge storage, but with a different way of retraining [62]. In our approach, the basic patterns of pollination dynamics are captured and stored in the machine learning models (metamodels) trained on vast amounts of simulation data. However, we did not split these metamodels, instead, the difference in pollination dynamics between different blueberry species was offset by adjusting a partial parameter set. Specifically, clone size (CS) and the density of different bee species (OS, BB, AD, and HB) were chosen for adjusting with the help of a multi-objective optimization method. In doing so, the best value of these features in the metamodels maximize the likelihood of predicting bog blueberry fruit-set. The merits of this technique are threefold. First, directly adjusting parameters (not for model manipulation) instead of the structure of a metamodel for retraining is helpful in some scenarios where the features learned are not multiscale. Splitting the model for retraining (e.g., like deep neural networks do) may not work. Second, directly adjusting parameters is straightforward and can be done very quickly, because searching an optimal set by means of a strong multi-objective optimization procedure has much lower computational complexity than training a deep neural network until it converges [54]. Third, since the number of parameters in our approach (four) is much less than the number of hyperparameters in a deep neural network (usually thousands), the retraining process requires much less data [63], which implies that our approach is more robust and universal than many transfer learning situations where data insufficiency is common. However, we acknowledge that our method, like all metamodeling approaches based on simulation data, is restricted by how effectively the simulation model can be exploited. If the simulation experimental design is not representative of the system or the number of experiments conducted is not adequate, information loss is expected and missed patterns might bias prediction results.

Practical Implications of Model Prediction
Our metamodel predictions showed that warmer, wetter weather is detrimental to optimal fruit set. As a measure of potential yield, a decline in fruit set has direct implications on profitability and sustainability of bog blueberry production in Northeast China. The current geographic distribution of bog blueberry is the cool, temperate (circumpolar) regions of the northern hemisphere [66]. This current climate envelope of cool, dry weather suggests that bog blueberry will not be tolerant to higher temperature and increased precipitation [67]. Our metamodel prediction experiments focused on reproduction and fruit set. Based upon wild blueberry, some of the mechanisms of decreased fruit set due to climate warming and increased precipitation are: reduced activity of bees during rainy days, reduction of the longevity of stigma viability during warm springs, faster rate of progression through the bloom period by the flowering clones (genetically unique individual plants), asynchrony between the emergence of wild bees and the early blooming clones, and increased blossom loss due to infection by Monilinia vaccinii-corymbosi (Reade) (the causal fungal organism of mummy berry disease) [15,26,68]. Parkinson and Mulder [2] also support our hypothesis that bog blueberry fruit set will be affected in a similar way that wild blueberry fruit set is by warming. In Alaska, USA, they found that bog blueberry exhibited a decline in fruit set in warmer lowland sites relative to cooler upland sites.
However, the loss in bog blueberry productivity due to climate warming and increased precipitation may be even greater than the losses that we modeled. Increased precipitation can result in other plant diseases due to fungi that reduce photosynthetic area and overall vigor of the plant and may make plants susceptible to other pathogens [26]. Graae et al. [69] support our hypothesis that climate change may increase fungal disease in bog blueberry as it does in wild blueberry. They found increased seed infection by fungi under warmer temperatures. Also in wild blueberry, Tasnim et al. [13] found that already occurring climate change resulting in warming and increased potential evapotranspiration is having negative effects on wild blueberry photosynthetic rates and overall plant growth rates. We expect that these same responses will be observed in the sibling species, bog blueberry. One of the few experiments with bog blueberry under changing weather conditions showed that nutrient uptake under warming declines in bog blueberry [70]. Therefore, the predicted decline in bog blueberry fruit set due to climate change may be only one component of a more complex suite of detrimental responses as a result of future climate change.

Conclusions
In our research, a new predictive modeling approach was introduced based upon the concept of transfer learning to solve the problem of data deficiency in predicting productivity of agroecosystems. Our paper delineates the process of building machine learning models (metamodels) from simulations built for one agroecosystem, where data resources are abundant; evaluating and selecting the best metamodel representing the dynamics of the system; and re-parameterizing and calibrating the metamodel to another agroecosystem where crop species-specific data is scare. To test our proposal, we used the pollination dynamics of the wild (lowbush) blueberry agroecosystem in the US as the source domain and the bog blueberry agroecosystem in Northeast China as the target domain. We have shown that knowledge gained in one crop system by modeling and simulation can be used for other crop systems sharing similar traits, which could save considerable investment in field observation and data collection [71]. Another contribution of our work is that it provides a feasible path for packaging knowledge and making it available across the scientific and grower community, which is often a neglected aspect of sustainable solution development [72]. In addition to the methodological exploration, we also applied our approach to predicting the response of productivity (i.e., fruit-set, as a proxy to yield) of Northeast China's bog blueberry agroecosystem to different weather conditions, which are regarded to be the most important abiotic factor affecting bog blueberry yield. Our application showed that in the most likely climate change scenario (i.e., higher temperatures with more precipitation), a decrease in fruit-set of Northeast China bog blueberries, the worst-case scenario being an 8% reduction was expected. The geographical pattern also revealed that the southern and eastern production regions will suffer more severe fruit-set decline than the rest of the bog blueberry growing regions. We suggest that while controlling blueberry clone size and wild bee populations is unlikely to be successful, increasing commercially available honeybee density from the current 1 bee/m 2 /minute to 18 bees/m 2 /minute, or commercially available bumble bee density from the current 0.2 bee/m 2 /minute to 0.6 bees/m 2 /minute, might be a viable way to compensate for the predicted 8% climate induced fruit-set decline in the near future.  Acknowledgments: We thank the three anonymous reviewers for their careful reading of our manuscript and their insightful comments, which help us to improve the manuscript. We are also grateful to several field workers and graduate students for data collection and preprocessing.