1. Introduction

1.1. Background

With the global population growth and the acceleration of urbanization, the demand for food has been steadily increasing. The security of food supply is essential for national economic stability and sustainability, which led to the importance of an efficient and stable Grain Supply Chain (GSC) to provide sufficient and safe food, stabilize market prices, and stimulate agricultural production. The GSC encompasses various stages, beginning with farmers and growers, moving on to producers and food companies, and then wholesalers and retailers [1]. Agriculture and grain production serve as the foundational components of the GSC. Previous research has pointed out that the GSC of China is facing severe challenges, including environmental pollution, the irrational use of fertilizers and pesticides, and a lack of breeding technologies [2]. In addition, grain production is also constrained by land and water resources, affecting the stability and adequacy of the GSC.

1.2. Related Research

Both domestically and internationally, scholars have analyzed GSC management and agricultural productivity through various approaches and models. Some scholars assert that the source of geographical diversity, cash flow management, and suppliers are the most important enablers for the ripple effects and the systemic interrelationships in GSC using Pareto analysis, total interpretive structural modeling, and matrice d'impacts croisés multiplication appliquée à un classement (MICMAC) methods, showing the complex relationships between variables, sustainability, and food security [3]. An impaired decision model under the government-subsidized producer framework was used to examine how government subsidies affect GSC loss reduction and member benefits [4]. Some scholars have used factor analysis, entropy weight method-TOPSIS, and grey relational analysis to show that Guangdong's GSC has low external risks of low yields and food reserves. The main elements determining higher internal risks include food production, food enterprise risk resistance, and the consumer price index for food [5]. This study emphasizes food production's importance to GSC integrity. However, in the field of food production research, factor analysis is one of the most common methods for categorizing the influencing factors of grain production into comprehensive input factors and disaster factors, ranking the contribution of each, and proposing corresponding suggestions [6]. An extensive range of nonlinear models, such as Gaussian Process Regression and Gradient Boosting Decision Trees, have been developed to construct a food production forecasting model for Chongqing, demonstrating the area of cultivation, agricultural machinery, labor, and disaster-affected area substantially impact food production [7]. Using spatial autocorrelation and a Geodetector model, some researchers developed a model for the quality of rural production spaces in Songnen Plain and examined the influencing factors and spatial agglomeration characteristics. Although various factors impacted the rural production environment, the main influences on the quality were the purchasing power of production resources, infrastructure, and agricultural mechanization [8]. Geographically weighted regression has been used to compare grain yield factors in Henan and Heilongjiang provinces. The two provinces were significantly different in terms of the influencing factors, which included input factors, socioeconomic factors, and the geographical environment in this comparative study [9].

Currently, research concerning the GSC mainly emphasizes technological advancements, governmental assistance, and risk evaluation. Nonetheless, there is limited research addressing grain production in certain regions, particularly in Heilongjiang Province, which has flat landscapes, rich soil, and significant organic matter content, making it the largest crop production base in China. However, with the uncertainties of the economic structure and climate change, crop production in this region encounters various problems, such as land degradation and inadequate agricultural infrastructure, which pose serious challenges for the GSC [10]. The existing literature predominantly emphasizes macro-level discussions, from the perspectives of national policies and industrial dynamics, addressing market supply and demand relationships, agricultural production models, and so forth. The research lacks empirical research on the micro-level factors, specifically how the individual stages of grain cultivation are influenced by external factors, making it difficult to assess Heilongjiang Province's grain production and growth potential.

1.3. Objection

This study aims to close the gap by building a Multiple Linear Regression (MLR) model within the framework of Stepwise Regression analysis focusing on the factors impacting Grain Yield in Heilongjiang Province from 1992 to 2022 within the context of the GSC. It provides empirical analysis and strategic recommendations for the grain production inputs, the stability of the GSC, and the agricultural sustainable development in Heilongjiang Province.

2. Methodology

2.1. Data

This study focuses on eight independent variables covering the years from 1992 to 2022 and aims to explore the elements affecting Grain Yield in Heilongjiang. The data was derived from the China Statistical Yearbook, the China Agricultural Yearbook, and the Heilongjiang Provincial Statistical Yearbook. The name and definition of each variable are illustrated in Table 1.

Table 1: Variable Definition

2.2. Models

In this study, Stepwise Regression Analysis and the conventional MLR model were conducted to perform the Grain Yield model.

Stepwise Regression Analysis is a method for building Linear Regression Models that assesses the statistical significance of each independent variable in the model. Using a stepwise approach, this process entails the method of an ideal predictive model, with variables being added and removed according to the statistical importance and contribution to the overall performance. The Stepwise Regression method utilizes the Coefficient of Multiple Determination (R²) and the significance of Regression Coefficients as its criteria. Specifically, the dependent variable was analyzed in each independent variable separately. The primary regression equation was used by the regression equation incorporated with the most significant independent variable, followed by other less significant independent variables sequentially. If the outcome improves the overall model performance by an increase in R² while other parameters remain significant, then the variable was kept, otherwise eliminated. This procedure was repeated until there were no considerable changes to R². The flowchart of Stepwise Regression Analysis Variable Selection is shown in Figure 1 below.

Figure 1: Stepwise Regression Analysis Variable Selection Flowchart [11]

Suppose that the MLR model consists of one dependent variable (Y) and k independent variables, which can be represented as follows:

Y = β₀ + β₁X₁ + β₂X₂ +… + β_KX_K + μ (1)

In this MLR model, k represents the number of independent variables from one to eight. Y is the dependent variable, X₁, X₂...X₈ are the independent variables, while β₀, β₁, β₂…β₈ are the regression coefficients, and μ represents the residual. While the residual μ shows the percentage of Y variation that cannot be determined by the current independent factors, the regression coefficients show the degree to which the independent variables affect the dependent variable.

2.3. Research Framework

This study has following steps:

Step 1: Descriptive Statistical Analysis. The objective of the first step is to identify the overall characteristics, trends, outliers, and potential correlations of all variables in the raw data. This process serves as a foundation for further studies concentrating on the relationships between Grain Yield and various factors.

Step 2: Benchmark Regression Model (BRM) Analysis. This step shows the basic connections between variables through a simple linear regression equation, which helps to build later models and to simplify the overall model. If the BRM shows inadequate fit, the study moves on to Step 3.

Step 3: Variable Selection Based on Stepwise Regression. The purpose of this step is to resolve Multicollinearity issues that are present in the BRM and to improve the model fit by eliminating independent variables that do not contribute to the dependent variable.

Step 4: MLR Analysis and model validation. This final step aims to assess the importance of each parameter in the regression model by utilizing the tests of economic and statistical significance, Residual Normality Analysis, Homoscedasticity Analysis, Autocorrelation Analysis, and Multicollinearity Analysis to ensure that all variables are significant and uncorrelated, thereby confirming the rationality and validity of the regression model.

3. Results

3.1. Descriptive Statistics of the Data

Descriptive Statistical Analysis of all variables was conducted in R, as shown in Table 2.

Table 2: Descriptive Statistics of the Data

3.2. Benchmark Regression

Table 3 shows the BRM statistic results of the dependent variable against the eight independent variables.

Table 3: Benchmark Regression Results

From the economic point of view, the regression coefficients suggest a positive correlation between Grain Yield in Heilongjiang and several factors, including Grain Sown Area, Effective Irrigated Area, Pesticide Usage, Total Power of Agricultural Machinery, and Value Added in Primary Industry. Alternatively, Grain Yield exhibits a negative correlation with Disaster-Affected Areas, Employment in Primary Industry, and Per Capita Regional GDP. Normally, the labor input can facilitate the management of crop cultivation and maintain large-scale agricultural production to increase the total Grain Yield. Thus, Employment in the Primary Industry should be positively correlated with the Grain Yield. Usually, a greater Per Capita Regional GDP implies advanced and proficient financial resources for agriculture research and development (R&D) of the area, leading to a surge in the Grain Yield. Hence, there should be a positive relationship between Per Capita Regional GDP and Grain Yield. However, the findings from the BRM challenge these assumptions, revealing the possible econometric issues within the model that need additional studies and validations.

From the statistical point of view, while the overall significance of the model is sufficient, the significance of the regression coefficients for individual independent variables is relatively weak. Only the regression coefficients for the independent variables X1 and X3 are significant, as shown by their corresponding p-values from the T-statistics, both of which are less than 0.05. Therefore, the next step is using a Stepwise Regression method to determine the important variables.

3.3. Multiple Linear Stepwise Regression

Following the approach outlined in Section 2.3, R² is used as the c^riteria to prioritize the significance of independent variables. The sequence of X1, X6, X2, X8, X5, X7, X4, and X3 is arranged from the most influential to the least, with X1 being the most influential because its regression equation has the highest R² value of 0.97778. Therefore, the regression equation of X1 was selected as the primary regression equation. Eventually, additional independent variables were added individually, with both R² and the significance of the regression coefficients, to determine which inadequate significance variables should be removed. The last selected independent variables are X1, X2, X3, X6, and X8. A MLR analysis is given on these five variables, with the results shown in Table 4.

Table 4: Multiple Linear Regression Results

The MLR model is expressed as follows: Y = 1738.02164 + 0.44847*X1 + 0.19329*X2 – 0.0408*X3 – 0.05332*X6 + 0.13628*X8

Grain Yield was positively correlated with Grain Sown Area, Effective Irrigated Area, and Value Added in Primary Industry, but negatively correlated with Disaster-Affected Area and Employment in Primary Industry. The regression coefficients show the estimated influence of each independent variable on the Grain Yield when all other variables remain the same. For example, for one unit increase in the Grain Sown Area, the Grain Yield will be increased by 44.847 million tons if all other variables remain constant. The intercept of the regression equation was 1738.02164, indicating the value was not affected by the independent variables.

The F-test is used to appraise the overall significance of the regression model. If the p-value of the F-statistic is below the predetermined significance level of 0.05, the Null Hypothesis may be rejected, meaning there are significant differences among the samples. In this model, the overall linear relationship between the dependent variable and independent variables was significant due to the p-value being significantly less than 0.05. The T-test is used to assess the significance of individual independent variables. Most of the regression coefficients of independent variables were significant, as calculated by the p-value of the T-statistic. The p-value for independent variables X1, X2, and X3 were significantly less than 0.05, implying that Grain Sown Area, Effective Irrigated Area, Disaster-Affected Area, and Employment in Primary Industry have a substantial impact on Grain Yield. R² is used as the primary indicator to examine the accuracy and dependability of the model. It quantifies the amount of variability in the dependent variable by the regression model. A higher R² means a stronger correlation between the model and the data. With the R² value of 0.9987, the model indicates that 99.78% of the variance of Grain Yield can be explained by independent variables.

3.4. Model Validation

In the following section, a comprehensive validation is conducted for the model in Section 3.3. The intention is to seek the satisfaction of the model’s assumptions, thereby improving the model’s accuracy and dependability.

3.4.1. Residual Normality Analysis

Figure 2: Normal quantile-quantile (Q-Q) plot (Photo credit: Original)

Many statistical inference methods, such as the T-test and F-test, rely on the assumption of normality. A significant deviation of the errors from a normal distribution may lead to issues such as biases and overfitting within the model. Figure 2 presents the Q-Q plot for the MLR model. The error points were very near to the diagonal line, which means that the random errors were normally distributed, particularly in the central region. The error points at the tails of the two ends, corresponding to 1994, 1999, and 2007 respectively, were deviated from the overall normal trend.

3.4.2. Homoscedasticity Analysis

Table 5: Homoscedasticity Analysis

It is crucial to determine whether the variance of the random errors stays constant. Suppose there was no Homoscedasticity in this model that makes the estimated regression coefficients, not the Best Linear Unbiased Estimator (BLUE), meaning that the validity of the T-test and F-test were seriously compromised. The White Test was employed using Eviews to identify the presence of Heteroscedasticity in the model. The Null Hypothesis of the White Test was Homoscedasticity. From Table 5, it can be observed that the p-value is 0.5264 greater than 0.05, indicating the Null Hypothesis cannot be rejected and there was no evidence of Heteroscedasticity in random errors.

3.4.3. Autocorrelation Analysis

Table 6: Autocorrelation Analysis

If the random errors were Autocorrelation, it may lead to an overestimation of the T-value, causing the Ordinary Least Squares (OLS) estimator to be ineffective. The Durbin-Watson Test was employed using R to assess autocorrelation, with the results presented in Table 6. The Null Hypothesis of the Durbin-Watson Test was no Autocorrelation, meaning there was no correlation among random errors. The p-value of the Durbin-Watson Test is 0.074 greater than 0.05, demonstrating the Null Hypothesis cannot be rejected and there was no evidence of Autocorrelation in random errors.

3.4.4. Multicollinearity Analysis

Table 7: Multicollinearity Analysis

Multicollinearity happens when two or more independent variables have a strong linear relationship. It is achieved by calculating the Variance Inflation Factor (VIF) for each independent variable. When VIF value exceeds 10, it may be considered indicative of Multicollinearity leading to a regression coefficient instability. Thus, Multicollinearity must be addressed to improve model performance. The VIF value shown in Table 7 indicates a significant enhancement in Multicollinearity compared to the BRM.

4. Discussion

4.1. Reason Analysis

These results indicate that Heilongjiang's Grain Sown Area, Effective Irrigated Area, Disaster-Affected Area, Employment in Primary Industry, and Value Added in Primary Industry have a significant impact on Heilongjiang's Grain Yield from 1991 to 2022.

As an irreplaceable resource, land resources are the most important factor in food production which is consistent with the above hypothesis. When all other factors remain constant, the yield will increase as the Grain Sown Area and Effective Irrigated Area increase. It shows that Grain Yield depends on environmental conditions.

Similarly, Disaster-Affected Areas such as climate anomalies, floods, and droughts have greater significance when compared to other seemingly important predictors. It causes the risk of grain yield reduction resulting in an inverse correlation between the Disaster-Affected Area and Grain Yield. This illustrates that protecting the environment and preventing the land resources from being Natural disasters are also important.

However, the negative coefficient associated with Employment in Primary Industry contradicts the common perception that "the more the labor, the higher the Grain Yield." This is consistent with the “Law of Diminishing Marginal Returns” theory as labor continues to move away from the peak, it leads to a decrease in Grain Yield when other conditions remain constant [12]. Moreover, past research confirmed the educational level of Employment in Primary Industry is also an important factor for the quality and yield of agricultural production due to their inadequate knowledge of cultivation and operational skills for modern equipment. Therefore, Employment in Primary Industry may have a negative correlation with the yield [13]. In this instance, the Songhua River Plain and the Sanjiang Plain are the primary grain-producing regions of Heilongjiang province. Excessive labor input of these huge plains could cause inefficient production. It becomes crucial to enhance agricultural mechanization and automation.

Although the significance of the coefficient for the Value Added in Primary Industry is low after variable selection, it’s still been retained which suggests a positive correlation between Value Added in Primary Industry and Grain Yield. The Value Added in Primary Industry could serve as an indicator of agricultural economic growth. An increase in the Grain Yield of the Primary Industry signifies a strong demand for agricultural products that positively influences the Grain Yield.

4.2. Recommendations

From the results of the empirical analysis and the challenges faced by agricultural production in the Heilongjiang region, the following recommendations are proposed to boost the Grain Yield.

4.2.1. Stabilize Grain Sown Area and Increase Effective Irrigated Area

The government should motivate farmers to cultivate large-scale rice fields, increase investments in water conservation, and ensure adequate irrigation resources to facilitate soil fertility. Furthermore, efforts should be made to improve the crop rotation and fallow land system. In recent years, the Chinese government has implemented policies and regulations aimed at promoting the use of land fallow in regions that experienced land degradation or ecological fragility. A land fallow trial program was initiated in 2016 [14]. Heilongjiang has formed a rotation model with corn and soybeans as the main crop, and potatoes as well as other grains as the secondary crops, which has gained acceptance by most farmers. By implementing the crop rotation and fallow land system, integrating of comprehensive management with planting and soil treatment could improve the ecological environment and ensure that resources are fully utilized, thus securing long-term sustainable grain production [15].

4.2.2. Strengthen Disaster Monitoring and Reduce Affected Areas

The government can strengthen the monitoring and early warning system to establish disaster resilience prevention. Additionally, it provides farmers with post-disaster relief measures in assisting agricultural reconstruction.

4.2.3. Reduce Ineffective Labor Inputs and Promote Agricultural Innovation

The government can minimize inefficient labor inputs while guaranteeing farmers steady work, supporting agricultural automation developments to further improve Grain Yield and production efficiency.

4.2.4. Invest in the Primary Sector and Promote Collaboration across Industries

The government should continue its funding in agriculture and associated fields, such as management techniques, infrastructures, and technology R&D. Supporting the growth of the digital sector and industry-to-industry cooperation will increase market demand and create a desirable market environment. This approach is beneficial to the stability of the GSC and the enhancement of the Grain Yield.

4.3. Limitations

Stepwise Regression has led to the retention of five independent variables and there is a high residual, suggesting that unobserved random variables may be the determining factor in Grain Yield. It is not feasible to conduct a comprehensive study for all elements due to the underlying limitations of both the model and the data. Therefore, further data collection is required to determine how supplementary factors affect the Grain Yield. Aside from the above-mentioned policy recommendations, the enrichment of Grain Yield and agricultural sustainability in Heilongjiang can be achieved by reinforcing agricultural product quality and safety, the growth of eco-friendly farming, as well as the market framework. By adopting these strategies, the Grain Yield of Heilongjiang is anticipated to grow ahead towards a better efficiency and sustainable future in national food security and economic advancement.

5. Conclusion

This study discovered the key factors of Grain Sown Area, Effective Irrigated Area, Disaster-Affected Area, Employment in Primary Industry, and Value Added in Primary Industry may influence Grain Yield. An increase in the Grain Sown Area, Effective Irrigated Area, and Value Added in Primary Industry may have a positive effect on Grain Yield, whereas an increase in the Disaster-Affected Area and Employment in Primary Industry may lead to a decrease in Grain Yield. These findings provide empirical analysis and recommendations for the crop production inputs as well as the GSC production stability in the Heilongjiang region. Although Stepwise Regression did not solve the Multicollinearity problem, the regression coefficient for each independent variable was more significant, and the Multicollinearity was substantially refined compared to the BRM. Therefore, the Stepwise Regression for variable selection was proven to be effective. The comprehensiveness of this study is constrained by the limitations of the model and the data. The influential factors were primarily based on production conditions and market demand. To further explore the possibilities of facilitating grain production stability and the interactions with the GSC, among other segments, future research should aim to establish models that are more closely aligned with the data while investigating additional factors affecting Grain Yield. This would contribute to the development of a more efficient, secure, and flexible GSC system to satisfy the growing demand.