Identifying the Driving Factors of Food Nitrogen Footprint in China, 2000–2018: Econometric Analysis of Provincial Spatial Panel Data by the STIRPAT Model

: This paper studies the EKC hypothesis and STIRPAT model. Based on the panel data of carbon emission intensity and other inﬂuencing factors of 30 provinces in China from 2000 to 2018, the spatial effect of per capita food nitrogen footprint (FNF) and the effect of different socio-economic factors in China were studied by using exploratory spatial data analysis and ﬁxed effect spatial Durbin model for the ﬁrst time. The results show that: (1) there is a spatial agglomeration effect and a positive spatial dependence relationship in China’s provincial per capita FNF (FNFP), which veriﬁes that the relationship between China’s FNF and economy is in the early stage of EKC hypothesis curve. (2) The driving forces of China’s FNF were explored, including Engel’s coefﬁcient of urban households (ECU), population density (PDEN), urbanization, nitrogen use efﬁciency (NUE) and technology. (3) The results show that there is a signiﬁcant spatial spillover effect of FNFP. The ECU and NUE can reduce the regional FNFP, and can slow down the FNFP of surrounding provinces. (4) Policy makers need to formulate food nitrogen emission reduction policies from the food demand side, food consumption side and regional level.


Introduction
Over the past decades, China's economy has been developing at an amazing speed and has been widely concerned. By 2020, China's GDP reached CNY 101.6 trillion, and its total economic volume has reached a new level of CNY 1 trillion [1]. The gap between urban and rural development continues to narrow, the income growth of rural residents is faster than that of urban residents, and the urbanization exceeds 60% [1]. The rapid growth of the economy and the acceleration of urbanization are accompanied by the trend characteristics of high nutrition, high protein, and diversification in food consumption of residents [2]. The dietary structure of residents has changed from plant-based food consumption to animal and plant-based food, a staple food consumption, to non-staple food [2][3][4]. The process of food consumption and production is both sides of the relationship between social supply and demand, independent of each other, but also an organic whole. The choice of food consumption directly determines the input and output of nitrogen in the process of food production, and food consumption is the main driving force of food production [5]. As the basis element of maintaining agricultural production, nitrogen plays an irreplaceable role in green agricultural development [6]. However, excessive accumulation of reactive nitrogen in the environment will cause a series of negative effects on human health and ecosystem [7]. The European nitrogen assessment identified five major threats to nitrogen pollution: water quality, air quality, greenhouse gas balance, ecosystem and biodiversity, and soil quality [8]. These impacts impede the progress towards sustainable development goals because they affect human health, resource management, livelihoods and the economy [9][10][11][12][13]. In recent years, nitrogen pollution management methods have changed [14], including new ideas of consumption and production, in order to seriously solve the problem of nitrogen [15][16][17][18][19][20][21][22][23]. Existing studies have shown that the loss of reactive nitrogen from food production and consumption is the most important factor causing China's nitrogen footprint, which has become one of the major environmental problems facing China [24]. In order to keep within limits man-made nitrogen emissions, realize the sustainable development of ecological environment and human health, it is very important to reveal the driving factors of nitrogen emissions.
In the field of environmental economics, there is abundant literature on the relationship between economic development and environmental quality, which can be better understood by means of EKC and STIRPAT which is a stochastic regression model for measuring the impact of population, affluence and technological changes on the environment. EKC Hypothesis is an important theory of empirical research on the relationship between economic growth and environment [25,26]. Since Grossman and Krueger (1991) [27] initiated this work, EKC Hypothesis revealed that the relationship between environmental deterioration and per capita income is an inverted U-shape, and the EKC Hypothesis has aroused great attention. Panayotou (1993) [28] confirmed the validity of the EKC Hypothesis. At present, the EKC Hypothesis, as a classical environmental economic theory, is widely used in the study of carbon emissions [29][30][31].
On the other hand, the STIRPAT model [32][33][34] is one of the most famous models in the literature of environmental economics. It is a mathematical generalization of the traditional IPAT (Impact = Population + Affluence + Technology) model proposed by Ehrlich and Holden (1971) [35]. It is a random version derived from the IPAT model. It overcomes the shortcomings of the assumption that the fixed proportion of environmental factors changes in the IPAT theoretical model [36]. The main goal of STIRPAT model is to estimate the impact of driving forces on ecology [37], to predict the non-proportional and nonmonotonic functional relationships among the factors affecting the ecosystem [34,38], and to explain the driving forces of SO 2 [39], carbon dioxide (CO 2 ) [40] and PM 2.5 [41], and other pollutant emissions. The advantage of the STIRPAT model is that it can evaluate the effect of driving factors [34], estimate the coefficient of each variable, and modify the influencing factors [42]. This model studies the effects of population, affluence and technology on carbon dioxide emission level [43,44]. Different scholars extended the STIRPAT model to incorporate additional factors into the traditional STIRPAT environmental impact model, enabling the extended STIRPAT model to study the impact of urbanization [45], industrial structure [46], trade openness [47] on greenhouse gas emissions, and other environmental consequences, which provides policy basis for energy conservation and emission reduction. When discussing the feasibility of EKC model, the literature is relatively rich. Due to the space problem, this paper only summarizes the different studies of nitrogen oxide (NOx) EKC in nitrogen emissions, in addition, discusses the research progress of random influence EKC based on STIRPAT theory. Table 1 shows the review of different studies on NOx EKC, which is summarized from literature, country, time frame, research methods, variables and conclusions. Among the environmental indicators used in EKC research internationally, CO 2 emissions are still the preferred indicators, because it is the most common in empirical literature [48]. Nevertheless, the use of NOx as a pollution indicator [28,[49][50][51][52][53][54][55][56] has been recorded in earlier and recent literatures, even though it is slightly reduced compared with CO 2 . The EKC hypothesis of NOx has been studied in OECD [54,66,71], Asia [48,68,70], Europe [55,57,58,64], America [63], developing countries, and developed countries [28,56,69], and most of them have been verified.
Similarly, many influencing variables of NOx emission have been identified. GDPP income is still the most widely used economic indicator to date, followed by per capita income [48,52,55]. Export [48,60,64], urbanization [48,56,61,62,70], foreign direct investment [50,61,62,65], and industrial structure [48,51,61,69] have been widely studied as influential variables in the research model, and are considered to have a significant impact on the relationship between NOx emissions and economic growth. The research conclusions reflect that most of the NOx emissions and the economy/income are in an inverted U-shape or an inverted N shape, which verifies the existence of EKC.
Based on the STIRPAT theory framework and EKC hypothesis, the current research mainly focuses on the following fields: The impact of income on carbon emissions in the African continent [72], the determinants of CO 2 emissions in Tianjin [73], the decoupling of the income-environment dynamic relationship between CO 2 and air pollutants at the Italian sector level [74], the driving factors of CO 2 emissions in high-, middle-and lowincome countries [75], the impact of Australia's main driving forces on the environment [76], the determinants of CO 2 emissions in Chinese cities [77], the driving forces of energy-related CO 2 emissions in China [78], the impact of human driving forces on the ecological footprint of Gannan Pastoral Area [79], and the impact of international trade on CO 2 emissions [80], etc. There are also few literatures based on the combination of STIRPAT and EKC to study the impact of urbanization [70], exports and foreign direct investment [61] on NOx emissions. The results of the study confirmed the existence of EKC within the framework of STIRPAT.
Nowadays, it is rare to conduct research on NOx EKC hypothesis from a spatial econometrics perspective [65,70] of available literature. Spatial effect is the key factor to evaluate the impact of economic development on environmental conditions [81,82], and the spatial correlation of data is the inherent feature of many environmental disciplines [70]. EKC finally presents the uncertainty results such as positive linear relationship, inverted Ushaped relationship, U-shaped relationship, N-shaped relationship, and no relationship [62]. When there is a spatial relationship, ignoring the spatial control in the econometric model will lead to biased estimation [83,84]. The impact of NOx is more local or regional than other pollutants in the world. It can be expected that it will be one of the pollutants more likely to be realized by EKC hypothesis [58]. One direction of future research may be to include spatial aspect (i.e., geographical proximity), so as to reveal possible regional differences and the sources of different EKC models [66].
Many literatures have adopted different environmental deterioration variables. The environmental variables listed in Table 1 above are mainly single substance pollutants, such as CO 2 , SO 2 , CH 4 and N 2 O, discuss the relationship between environmental deterioration and socio-economic development. The use of a single substance pollutant measurement indicator will cause errors in the measurement results due to insufficient comprehensiveness [85]. In fact, the overall concept of pollution is composed of different components. The performance of these components may be different, and it is not easy to merge into a single measurement standard [52]. The standard can be used for reference to the construction method of the Comprehensive Environmental Index (CEI) [86]. This study uses FNF as an environmental quality indicator. Based on the EKC model and STIRPAT model, this study uses the FNF as an environmental quality indicator, and selects the panel data of the FNF and related regional development factors of 30 provinces in China from 2000 to 2018 to make empirical estimation through spatial econometric models, then analyzes whether there is a spatial autocorrelation relationship between FNFs in various regions of China, and estimates the spatial spillover effects of various factors by establishing a spatial panel model of FNF and each influential factor.

Calculation of Food Nitrogen Footprint
According to the definition of NF proposed by Leach et al. (2012) [87], FNF in this study can be divided into food consumption nitrogen footprint (FCNF) and food production nitrogen footprint (FPNF). The calculation method of FPNF and FCNF is based on the adjusted version of N-calculator. FNF reflects the sum of N losses associated with food production and food consumption. FPNF refers to the loss of all N from food production to food consumption. FCNF refers to the amount of food nitrogen consumed by human beings, assuming that the nitrogen excreted by adults is the same as that in food described by Leach et al. (2012) [87].
The FNF PC is calculated according to the following formula: In Equation (1), FNF PC is theFNFP, FCNF i is the per capita FCNF of item i food, FPNF i is the per capita FPNF of item i food, and i refers to all kinds of food (n = 8).

FCNF Calculation
The per capita FCNF is calculated by the following formula: In Equation (2), FCNF PC is the per capita FCNF, C i is the per capita food consumption of item i food, and N i is the nitrogen content of item i food.
In order to estimate FCNF, we assume that nitrogen does not accumulate in the human body, which means that all consumed nitrogen is lost in the sewage flow [87]. Due to the lack of available data on the effectiveness or any temporal variation of nitrogen removal from wastewater treatment plants, this study did not consider wastewater treatment using nitrogen removal technology [88], and all nitrogen consumed by human beings was released into the environment in the form of human waste. The level of nitrogen content of different foods is found by multiplying the food protein content by the conversion coefficient of 0.16 [87]. The nitrogen content of 8 items of different foods is shown in Table 2 [4,[89][90][91].

FPNF Calculation
Virtual Nitrogen Factor (VNF) is introduced into the NF accounting of food production. VF refers to any nitrogen produced in the process of food production that will not be directly consumed by human beings, including all nitrogen lost by nitrogen fertilizer application in farmland, livestock breeding and food processing [87]. In order to define the boundary and avoid repeated calculation, in the calculation of VF, the nitrogen loss caused by energy consumption in food production is generally excluded, and this part of nitrogen loss is calculated into the energy NF. VNFs are calculated by dividing total nitrogen loss by total available nitrogen. The per capita FPNF is calculated as follows: In Equation (3), FPNF PC is the per capita FPNF, and VNF i is the VNF of food item i. We determined the VNF in this study based on regional characteristics, the time span, main food consumption, and data availability ( Table 2).

Environmental Kuznets Curve Hypothesis
EKC was originally an empirical hypothesis, which described an inverted U-shaped curve between economic development and environmental quality. Various indicators of environmental quality degraded with economic growth at the initial stage, and after reaching the threshold stage, environmental degradation began to decrease [25]. Since the EKC hypothesis is based on the relative importance of scale effect, composition effect, and technique effect, structural modeling should be performed to identify them. However, some simplified models have been used to test the validity of EKC hypotheses in empirical literature. The following general simplified model is used to test the EKC hypothesis in the empirical literature [92].  (4), Y is the dependent variable representing environmental pollutants, X is the independent variable used as an economic income indicator, Z is the vector of other control variables that can affect Y, α is the intercept, β i is the estimated regression coefficient, and e is the error term. In this paper, the environmental pressure dependent variable Y is measured by FNFP, and the economic output variable X is measured by GDPP. Many studies use the logarithmic form of the simplified model described above [60,[93][94][95][96]. Select the most appropriate function form for the data and have a higher explanatory power within the data range [97]. The significance of model parameter β i is estimated and tested. Dinda (2004) [92] pointed out that the following five results may occur: i. If β 1 = β 2 = 0, then there is no relationship between Y and X.
ii. If β 1 > 0 and β 2 = 0, then there is a monotonic increasing or linear relationship between Y and X.
iii. If β 1 < 0 and β 2 = 0, then there is a monotonic decreasing relationship between Y and X.
iv. If β 1 > 0, β 2 < 0, then there is an inverted U-shaped relationship between Y and X, thus, the EKC hypothesis is valid.

The Extended STIRPAT Model
The STRIPAT model is used as the theoretical basis of this study. Firstly, it is used to verify the existence of EKC between GDPP and FNFP, and secondly, to estimate the driving factors of per capita FNF in China. As a model widely used in ecological environment research, the main idea of STRIPAT is that environmental impact is a function of Population (P), Affluence (A), and Technology (T). The stochastic form of STIRPAT modified by Dietz and Rose [34], which is expressed as follows: Given this, we establish the following panel data model:  (6), I is the environmental impact, and the FNFP is selected to measure the environmental impact, the PDEN is selected to represent the population P, the GDPP is selected to represent the affluence A, the ratio of the number of patent applications to the total population is selected to represent the technology T. b, c and d are the elastic coefficients of lnP, lnA and lnT, respectively. α is a constant term, and ε is an error term.
Through logarithmic transformation of the variables in Equation (4), a logarithmic regression form for estimating and testing hypotheses is obtained. Combining with Equation (6), we deduce the existence verification model of EKC under the framework of STIRPAT model, which is expressed as follows: In order to achieve the research purpose, it is necessary to consider many driving factors that affect the FNFP as independent variables. Based on the explanation in Table 1, this paper constructs an empirical model after the extended STIRPAT model, as shown in Equation (8). Table 3 summarizes the selection, interpretation, and description of variables.

. Global Spatial Autocorrelation Analysis
The global Moran's I index can be used to measure the spatial distribution pattern of FNF and to test the spatial autocorrelation of FNF in China. The calculation formula of global Moran's I index is as follows: In Equation (9), I is the global Moran index; x i and x j are the FNF PC of the i and j provincial spatial units respectively; n is the number of selected provincial units; x is the mean value of FNF PC ; W ij is the spatial weight matrix. This study adopts 3 forms: adjacency weight matrix, geological distance weight matrix, and economic distance weight matrix.

Local Spatial Autocorrelation
Local Moran's I index indicates whether there is a high-value or low-value concentration of FNF and its neighboring provinces. The local Moran's I index is calculated as follows: In Equation (9), S is the standard deviation of FNF PC of each province; x i , x j , n, x, W ij have the same meaning as x i , x j , n, x, W ij in global spatial autocorrelation analysis Equation (9).

Spatial Panel Econometric Model
If exploratory spatial data analysis does find the spatial autocorrelation of the data, spatial econometric analysis should be used. Spatial econometric analysis was first proposed by Anselin (1988) [83]. The model in this study is based on the spatial panel data model proposed by Elhorst (2015) [98], which combines the advantages of spatial model and panel data model, and considers the spatial effect of driving factors ofFNFP.
Generally speaking, there are three spatial panel data models [99]: spatial lag model (SLM), spatial error model (SEM), and spatial Durbin model (SDM). The general form of these spatial panel data models [98] is as follows: Based on the empirical model of extended STIRPAT (Equation (8)) and the general form of spatial panel data model (Equation (11)), the spatial Durbin panel model of this study is constructed. The specific expression (Equation (12)) is as follows: In Equation (12), W is a standardized non-negative spatial weighting matrix, which represents the adjacency relationship between regions. ρ is the regression coefficient of the dependent variable on the spatial autocorrelation, which indicates the impact of provincial FNF on other provinces. θ is the coefficient of the spatial lag term of the independent variable, which represents the spatial spillover effect of the independent variable on neighboring provinces. The more significant θ is, the stronger the spatial interaction of the independent variables is. β is used to measure the effect of independent variables on the change of FNF in the province. α i and γ t are spatial and time fixed effects, respectively, and ε is an unobservable random error.
In Equation (12), when θ is equal to 0, The model is Spatial Lag Model (SLM). SLM only contains the spatial lagged term of the dependent variable, which is expressed as: In Equation (12), when θ is equal to −ρβ, the model is the Spatial Error Model (SEM). SEM only contains the spatial dependence of random errors, which is expressed as:

Test of Spatial Panel Model
In estimating, spatial econometric models need to be screened using appropriate methods and procedures, mainly using LM tests and (Robust) LM tests to test whether spatial factors are included, and to establish a spatial model. Then which form of spatial panel model is adopted by Wald and LR tests [100]. In the process of a test, the following criteria can be selected: According to the selection conditions of SLM and SEM established by Anselin et al. (2008) [101], the corresponding OLS model is established for panel data, and the residual items are analyzed and judged by LM test and (Robust) LM test, and the form of the model is determined according to the significant level of the Lagrangian Multiplier contained in the LM and Robust LM tests. If the Lagrangian multiplier of the lag model is more significant than that of the error model, and the Robust-LMLag of the lagged model is significant, but the Robust-LMError of the error model is not significant, the SLM model is selected. If the Lagrangian multiplier of the error model is more significant than that of the lag model, and Robust-LMErrror of the error model is significant, but the Robust-LMLag of the lag model is not significant, the SEM model is selected.
Lesage and pace (2009) [84] believe that the loss caused by ignoring the hysteresis of error term itself will be less than the loss caused by ignoring the common spatial effect of independent variable and dependent variable. On this basis, Elhorst (2014) [102] proposed that SDM model should be given priority in model estimation, and at the same time, LR and Wald tests were used to consider whether SDM should be converted into SLM model or SEM model. The effects achieved by the LR and Wald tests are approximately equivalent. According to the above conversion conditions of spatial model, the null assumption of SDM conversion to SLM and SEM is expressed as: θ = 0 and θ = −ρβ. Under the LM test condition, if the model passes LR test or Wald test, H 0 is refused, then SDM model can be considered.

Data Sources
The data of FNF covering a period of 19 years, from 2000 to 2018, involving the per capita consumption data of 8 items food, including grain, vegetables, melons and fruits, animal meat, poultry meat, aquatic products, eggs and milk. Geographically, it covers 30 provinces in China. Due to data availability, this study does not include Tibet, Taiwan, Hong Kong, and Macao. Specific data sources: (1) urban per capita food consumption data.  Table 4 shows the descriptive statistical results of the main indicators used in this paper.

Spatial Distribution of Food Nitrogen Footprint in China
In order to explore the spatial effects of FNF in China, we calculated the spatial statistics of Moran's I across the region, represented by the Moran scatter diagram, with examples in 2000, 2010 and 2018. As shown in Figure 1, the Moran scatter plot shows the global Moran's I-index, Z-statistic, and p-value. The horizontal axis represents theFNFP, and the vertical axis represents spatial lag ofFNFP. Each value represents the corresponding province of China. The points in the first quadrant H-H (High High) and the third quadrant L-L (Low Low) indicate that the FNFP of the province and neighboring provinces are higher and lower, respectively. The points in the second quadrant L-H (Low High) and the fourth quadrant H-L (High Low) indicate that the FNFP of the province is significantly lower and higher than that of adjacent provinces, respectively. The slope of the point regression line in the scatter plot is actually the value of Moran's I. It can be seen from Figure 1 that China's FNFP has a significant spatial dependence, with Moran's I index passing the 1% significance test in 2000 and 2010 and reaching the 5% significance level in 2018. The positive value of Moran's I reveals the positive spatial autocorrelation of FNFP. The higher Moran's I value, the more obvious the spatial correlation of FNFP. In addition, according to our calculations for all 19 years during the study period, the Moran's I values ranged from 0.138 to 0.351, and all values passed the 1% or 5% significance test, which indicates that the positive spatial autocorrelation has long-term stability [104].

Selection of Spatial Panel Model
Through the above exploratory spatial data analysis, the existence of spatial autocorrelation of FNFP was confirmed. In accordance with the steps proposed by (Elhorst, 2014a) [105], firstly, without considering the spatial effect, OLS was used to estimate the non-spatial model to test the relationship between economic growth and FNFP. Secondly, regarding panel data analysis, we examined four types of fixed effects: pooled OLS, spatial fixed effects, time fixed effects, and two-way fixed effects. According to the estimation results, the test of application space model is given in the results (Table 5). Table 5. Estimation results of food nitrogen footprint EKC using panel data without spatial effect considered.  First of all, we tested the relationship between spatial error and spatial lag of panel data. In the absence of spatial factors, we used a Lagrange multiplier test for four models (pooled OLS, spatial fixed effects, time fixed effects, and two-way fixed effects) of ordinary panel data. The no-spatial error effect and no-spatial lag effect of LM and Robust LM were tested by the spatial toolkit provided by Elhorst (2014a) [105]. Results (Table 5) show that in the mixed OLS, spatial fixed, temporal fixed, and double fixed effects, the probabilities of the LM test no spatial lag are significant at 1% level except for Robust LM test no spatial error of time fixed effects with at 5%, and therefore four types of fixed refused the null hypothesis of no spatial error relationship of the model. The probability of the LM test with no spatial error is significant at the 1% level. Pooled OLS fails to pass the Robust LM test no spatial error, and the spatial fixed effects, time fixed effects and two-way fixed effects pass the significant levels of 1%, 5% and 5% of Robust LM test no spatial lag, respectively. The goodness of fit of the pooled OLS, spatial fixed effects, time fixed effects and two-way fixed effects was compared from R 2 , the goodness of fit of spatial fixed effect model is greatly improved to 0.8037, which is far higher than the other three fixed effects with R 2 lower than 0.5.

Pooled OLS Spatial Fixed Effects Time Fixed Effects Two-Way Fixed Effects
Next, in the comparison of fixed effect model and random effect model, Hausman test results show that in the case of 9 degrees of freedom, the test value is 20.39, which passes the 5% significance level test and refuses the random effect model. Therefore, the spatial fixed effect model will be selected in this study.
Once again, in order to further test the fitting effect of SEM and SLM, the SDM is introduced. The Spatial Durbin model contains both spatial error term and spatial lag term. Wald and LR tests are used to estimate and compare the three spatial econometric models.
Finally, SEM, SLM and SDM were respectively constructed for the panel data. Due to the spatial correlation, the basic assumption model is no longer applicable, so we use the quasi-maximum likelihood method (QMLE) for statistical test, and the test results are shown in Table 6. 158.06 *** Notes: *, **, *** indicate that statistics are significant at the 10%, 5%, and 1% level of significance, respectively. Standard errors z in parenthesis.
As can be seen from Table 6, SLM passes the 1% significance test of Wald spatial lag test and LR spatial lag test, it rejects the null hypothesis H 0 : θ = 0, thus rejects the null hypothesis that "the model will degenerate into SLM". SEM passes the 1% significance level test of Wald spatial error test and LR spatial error test, the null hypothesis H 0 : θ + ρβ = 0 is rejected, thereby rejecting the null hypothesis that "the model will degenerate into a spatial error model". In conclusion, SDM is more suitable to describe the panel data of FNFP at provincial level in China.

Results of EKC Validation and Influencing Factors Analysis
From the estimated results of the fixed effects model (Table 5) and the SDM model with fixed effects (Table 6), the linear term of GDPP is significantly negative, and the quadratic term of GDPP is significantly positive, which means that the relationship between FNFP and economic development is not the inverted U-shaped assumed by the classical EKC, but a U-shaped relationship.
From the estimation results of the SDM model with fixed effects (Table 6), it can be seen that PDEN, technology, urbanization and ECU have a positive impact on the growth of FNFP at a significant level of 1%. In the case of other factors unchanged, each 1% increase of these four factors can promote the FNFP increase by 0.2859%, 0.0264%, 0.1862% and 0.3389%, respectively. As urban income increases, the ECU will decrease, leading to a decline in the FNF. NUE has a significant positive effect on the increase of FNFP at a 5% level, with an elasticity coefficient of 0.0439. Strengthening the NUE and reducing the amount of nitrogen fertilizer application per unit of grain yield will slow down the FNF. The elastic coefficients of foreign trade and industrial structure are small, and has not passed the significance test.
According to the spatial lag coefficient, the elasticity coefficient of ρ at 5% significant level is 0.1216, which indicates that the panel data of FNFP in China have a positive spatial dependence, and the provincial FNFP is spatially clustered. The coefficients of W*lnOPEN, W*lnNUE, and W*lnECU in the spatial lag terms of explanatory variables passed the 1% significance level test, and the values were −0.0588, −0.1331, and −0.7414, respectively. This suggested that the explanatory variables in the SDM model with fixed effects also had spatial correlation, which means that foreign trade, NUE and ECU in neighboring provinces every increase an average of 1%, then the FNFP in this province changed considerably by −0.0588%, −0.1331% and −0.7414% in the circumstances of other variables keeps invariant. Moreover, W*lnPDEN reached the 5% significance level, its elasticity coefficient was −0.3979. W*lnURB passed the 10% significance level test with a coefficient of −0.1748.

Results of Spatial Effect Analysis
In fact, when estimating the spatial panel model, there are not many explanations for the coefficients of explanatory variables, that is, the coefficients of explanatory variables are not of great significance. What really needs to be explained are the direct and indirect effects, namely spatial spillover effects. Spatial spillover effect refers to the change of explanatory variables in adjacent regions caused by the change of explanatory variables in one region. The direct effect refers to the influence of the change of a certain factor in a certain area on the local area, and the indirect effect refers to the influence of the change of a certain factor in a certain area on its adjacent area. This indirect effect is transmitted through spatial interaction, and gradually decreases with the increase of the distance between spatial units [106].
In order to further analyze the effects of spatial interaction, this paper decomposed the influence of driving factors on FNF into direct effect and indirect effect according to the method provided by Lesage and Pace (2009) [84]. The effect of the explanatory variable on the FNFP of a local area is a direct effect. The impact of the explanatory variable on the FNFP in other regions is an indirect effect (i.e., the spatial spillover effect of influencing factors). The sum of the two is the total effect. Table 7 shows the effect decomposition results of SDM estimation based on fixed effects. The first and second columns of Table 7 show the direct and indirect effects of each explanatory variable on FNFP. The third column is the total effect of each explanatory variable on FNFP. The results show that: Table 7. Decomposition of spatial effects. (1) The total and direct effects of GDPP on FNFP were negative. At the same time, the indirect effect of GDPP had a negative impact. The significance test of direct effect reached the level of 1%, the significance test of total effect passed the level of 10%, and the indirect effect did not reach the level of significance test. It shows that under the condition of other variables unchanged, GDPP can directly reduce the FNFP of local provinces.

Variable Direct Effect Indirect Effect
(2) The total and indirect effects of PDEN on FNFP were negative, while the direct effects were positive. The direct effect and indirect effect reached the significance level of 1% and 5%, respectively, and the total effect did not pass the significance test. This means that the increase of PDEN may worsen the local nitrogen emissions, but has a mitigating effect on the FNFP of the surrounding provinces.
(3) The direct effect of technology development on FNFP is positive, with a significance level of 1%, while the indirect effect and the total effect do not pass the significance level test. This suggests that technology development has not reduced nitrogen emissions in the region.
(4) Foreign trade has a negative effect on the FNFP, the significant level is 1%, which illustrates that foreign trade can reduce the FNFP. The direct effect is positive without the significance level, the indirect effect is negative without the significance level of 1%. This implies that foreign trade is conducive to improving the FNFP in neighboring provinces. (5) The direct effect of urbanization on FNFP is positive, and the significance level was 1%. while the indirect effect and total effect did not pass the significance level test. This makes clear that the development of urbanization stimulates nitrogen emissions in the region.
(6) The NUE has a negative effect on the FNFP, with a significant level of 5%, indicating that the FNFP could be reduced with the increase of NUE. The direct effect was positive, and the indirect effect was negative. The direct effect and indirect effect reached the significant level of 5% and 1%, respectively. This indicates that the increase of NUE may worsen the local nitrogen emission, but the increase in the FNFP of the surrounding provinces is not obvious. (7) The total effect of the ECU on the FNFP is negative, with a significance level of 1%, which indicates that ECU can reduce FNFP. The direct effect was positive, and the indirect effect was negative. The direct effect and indirect effect reached the significant level of 5% and 1%, respectively. This explains that the reduction of ECU may improve the local nitrogen emission, but can increase the FNFP of surrounding provinces. (8) The total effect, direct effect, and indirect effect of industrial structure on the FNFP did not pass the significance level test.

Discussion
The empirical results of this study provide strong evidence for the spatial correlation of food-related nitrogen emissions and the spatial spillover effect of FNF in China. The highly significant global Moran's I index and local Moran's I index proposed HH and LL spatial aggregation patterns of FNF. On the basis of the extended STIRPAT theoretical framework and spatial econometric analysis approach, this study empirically verified that economic development, PDEN, technology, urbanization process, ECU and NUE are the driving forces of FNF in China.
Considering that the GDP-FNF model with quadratic term of GDPP has the best fitting, the valuable findings of this study are that there is a U-shaped EKC in the relationship between GDP and FNFP, which is consistent with the research conclusion of Och (2017) [48], but our conclusion is inconsistent with the previous similar research results on NOx emissions. Brajer et al. (2011) [52] and Luo et al. (2014) [51] believe that there is an inverted U-shaped EKC in the relationship between China's NOx emissions and GDPP, while GE et al. (2018) [70] verify that there is an inverted N-shaped EKC in this relationship. We found this inconsistency caused by two main reasons: (1) we applied the extended STIRPAT theory to the explained variable as the FNFP, which is a composite index rather than a single substance pollutant environmental index, and the explanatory variable is relevant factors that has a significant contribution to the FNFP, while Their research does not control the potential spatial autocorrelation between regional emissions, even if their samples are from cities that are usually spatially related. The existence of spatial relationship provides a potential explanation for the instability of EKC parameters, and the behavior of neighboring provinces will affect the behavior of the region itself.
In the study of influencing factors for NOx emission, the main research methods currently are OLS regression model, unit root, cointegration test, and Granger causality test. With the in-depth research of econometric methods, vector error correction model (VECM), generalized Moment estimation (GMM), fully modified least squares method (FMOLS), generalized least squares estimator (GLS) and cross-section dependence test (Cross-section dependence test), co-integration autoregressive distribution lag (ARDL) boundary test, and other methods are generally widely used in the study of NOx emission. The environmental indicator of this study is FNF, which is a composite environmental indicator to quantify the reactive nitrogen released during the food life cycle and its impact on the environment. The above-mentioned methods are applied to study the influencing factors of FNFP in China's provinces, only discuss the development and change of FNF in the provinces themselves, while ignoring the spatial correlation and spatial characteristics of FNF between each province. Ge et al. (2018) [70] found that there was a significant inverted N-shaped EKC between NOx emission and urbanization by using spatial econometric model. Using exploratory spatial data analysis and fixed effect SDM, this paper further concludes that China's inter-provincial FNFP presents spatial agglomeration effect, and has a significant positive spatial dependence on the FNFP of neighboring provinces. Therefore, this paper not only analyzes the influencing factors of provincial FNFP development, but also fully explores the spatial spillover effect of provincial FNFP development elements on surrounding provinces through spatial econometric method.
According to the analysis of the influencing factors of FNFP, the ECU has the greatest impact on FNFP. In the context of the acceleration of urbanization in China and the improvement of people's living standards, the decrease of ECU will improve the status of FNF, cause changes in nutritional diet, increase the proportion of plant-based foods, and reduce the ratio of animal-based foods. Another major factor is population growth leading to an increase in the FNFP. Population puts great pressure on food security, resulting in nitrogen emissions at both ends of food production and consumption. Technology development has not dropped the FNF, which indicates that technology has not decreased nitrogen emissions in food production, transportation and storage, in particular, the fertilization process of food production has not yet reflected green production.
The basic findings of this study refute the existence of EKC hypothesis of China's FNFP, thus emphasizing that the relationship between China's FNF and economic development is still in the early stage of EKC hypothesis. Among many influencing factors, ECU and use efficiency have been found to reduce the FNFP. Government decision-makers need to consider three policy-level measures: (1) At the food production side, the Chinese government has put forward the Action Plan for Zero Growth in the Use of Fertilizers by 2020 [107], and implemented such technical paths as optimizing the structure of chemical fertilizer use, reducing the use of nitrogen fertilizer, applying organic fertilizer, improving fertilization methods, promoting precision fertilization, and steadily improving NUE. (2) At the food consumption side, policy makers advocate the social fashion of Practicing Strict Economy, and Opposing Waste, and repeatedly emphasize the need to stop food and beverage waste. At the same time, the Opinions on Promoting Food and Beverage Economy and Opposing Food and Beverage Waste by Standardization [108] issued by the government proposes to vigorously carry out the promotion activities of food and beverage saving standards.
(3) At the regional level, the change of regional FNF and the spatial spillover effect of socio-economic factors on the FNF of neighboring provinces should be fully considered. When formulating relevant policies for the prevention and control of nitrogen pollution, Government functional sector should take into account the whole cooperation between provinces and regions. With the implementation of the policy and the development of economy, the turning point of FNFP will come in advance.
Due to the lack of relevant statistical data and parameters, it is more difficult to estimate the FNF and the error of the results, which contributes to the uncertainties of our study. Strengthening the research of food VNF based on the local situation in China is conducive to improving the accuracy of FNF estimation and reducing the nitrogen waste in each process from the perspective of food life cycle.
The influence of FNF on the environment of non-point pollution may be one of the directions of future nitrogen footprint research. In this paper, FNF and application of nitrogen fertilizer are analyzed. However, the nitrogen loss and waste in both ends of food production and consumption need further study. The main sources of non-point pollution caused by these nitrogen losses are analyzed. Then, the ecological effect of nitrogen footprint on the watershed environment is better evaluated.

Conclusions
In this study, the modified N-Calculator method proposed by Leach et al. (2012) [87] was used to calculate the FNFP of 30 provinces in China from 2000 to 2018. The extended STIRPAT theory was used as the research framework, the EKC idea was incorporated, and the exploratory spatial data analysis method was applied to study the spatial autocorrelation of China's FNFP. On this basis, the spatial Durbin panel model with fixed effect was adopted to study the influence and spatial effect of socio-economic factors on provincial FNF in China, and the following conclusions were drawn: First, this study first verified the quantitative relationship between China's economy and FNFP. Compared with traditional econometric methods, spatial econometric techniques have never been used to explore the relationship between FNF and economy. Due to the spatial correlation of FNF to adjacent areas, the parameters estimated by the spatial panel data model are more reliable than those of the traditional panel model. The results show that the relationship between economy and FNFP is different from the current inverted U-shaped EKC hypothesis between economy and NOx, but appears as a U-shaped EKC hypothesis curve, which clarifies that the relationship between China's FNF and economy is in the early stage of EKC hypothesis curve.
Second, in the study on the influencing factors of the FNFP, it is verified that the main driving forces of China's FNFP are the ECU, PDEN, urbanization, NUE, and technology.
The improvement of ECU and nitrogen utilization efficiency will slow down the FNFP in China. Population and urbanization have increased the FNFP in China, and technology does not slow down the FNF in China. It is proved that there is a positive spatial dependence of the FNFP, indicating that there is a significant spatial spillover effect in the FNFP. By means of spatial partial differential decomposition of the estimated results, ECU and nitrogen utilization efficiency can reduce the FNFP in the region, and can slow down the FNFP in the surrounding provinces.
Third, the paper provides a framework for the selection of spatial econometric model of FNF, LM test results of panel data show that the no fixed and spatial fixed panel model are the choice of this study. Further spatial LR and Wald spatial lag (error) test results show that SDM cannot be simplified as SAR or SEM form, SDM is more suitable for this research, the research on the spatial correlation problem of SDM can more fully explain the space effect, make up space lag panel overlook the defects of explanatory variables space effect, therefore, research on spatial correlation problem, can give priority to the SDM model. In the study of spatial correlation, SDM can more comprehensively explain the spatial effect and make up for the defect that the spatial lag panel ignores the spatial effect of explanatory variables. Therefore, SDM model can be given priority in the study of spatial correlation.

Policy Recommendations
First, the ECU, the PDEN, the Urbanization, the NUE, and the Technology are the main factors for the growth of FNFP in the province. By guiding healthy and scientific eating habits, reducing the proportion of meat and fish food in the diet structure, stopping food waste, and promoting the pattern of polite, saving, green, and low-carbon consumption, it is an effective way to guide the rapid growth of FNFP in each province. Furthermore, reducing nitrogen fertilizer application level per unit grain yield, and improving NUE, which are also effective suggestions to reduce the growth of FNFP in China's provinces. It is suggested to implement a strict nitrogen reduction policy, accelerate the development and promotion of new fertilizers such as slow-release fertilizer, water-soluble fertilizer, and biological fertilizer, and promote the replacement of inorganic nitrogen fertilizer by organic fertilizer in proportion. Second, the significant spatial spillover effect of FNFP suggests that policy makers, especially local governments, should not only pay attention to the local nitrogen emission, but also consider the impact of nitrogen emissions on neighboring provinces. The Chinese government should formulate policies at the national level to break the administrative boundaries of local governments and encourage local governments to actively participate in promoting food nitrogen emission reduction and economic growth. In addition, policy makers should formulate food nitrogen emission reduction plans at both the food production and consumption sides, advocate green production and green consumption, and strengthen the application and promotion of technology in nitrogen reduction.

Conflicts of Interest:
The authors declare no conflict of interest.