三元肥效模型異方差診斷及其可行廣義最小二乘回歸建模研究

孔慶波，章明清，李娟，許文江，章贊德，姚建族

（1福建省農業科學院土壤肥料研究所，福州350013；2福建省亞熱帶植物研究所，福建廈門361006；3大田縣土壤肥料技術推廣站，福建大田366100；4永春縣農業技術推廣站，福建永春362200）

0 引言

【研究意義】肥料效應函數法是中國測土配方施肥技術的一個重要分支體系。但根據田間肥效試驗結果，應用普通最小二乘法（OLS）建立肥效模型的成功率明顯偏低，嚴重制約了該法的計量精確性和實用價值。因此，探討如何提高田間肥效試驗建模成功率具有意義。【前人研究進展】在常用肥效模型中，二次多項式肥效方程能較好地反映作物產量和施肥量之間的數量關系，模型計算和參數估計簡便易行，是研究和應用最多的模型種類[1-3]。但在施肥實踐中，OLS法回歸分析建立的一元肥效模型典型式僅占 60%左右，二元肥效模型典型式則只有40.2%[4-5]，而三元肥效模型的典型式更低至 23.6%[6]。針對肥效模型在應用中存在的建模成功率偏低的問題，國內外學者就肥效模型選擇與其適用性[7-9]、試驗設計與參數估計方法[10-12]、類特征肥效模型構建[13-15]及非典型式推薦施肥優化方法[12]等方面進行了深入研究，提出了許多改進措施，但至今未能得到較好地解決。當前廣泛應用的“3414”設計方案的肥效試驗結果既使采用蒙特卡洛（Monter Carlo）建模法，二元和三元模型典型式比率也僅有 56.7%和 37.3%[12]。這種成敗幾乎各半的建模結果，困擾著計量施肥研究和應用。雖然任何田間試驗都可能存在失誤從而導致失敗，但是，在不同肥料之間、肥料與其他產量限制因子之間存在著復雜的聯應關系，我們沒有理由認為其余的過半數試驗均是失誤。【本研究切入點】隨著回歸分析理論的發展和完善，人們發現，當回歸方程隨機誤差項的方差存在顯著水平差異時，普通最小二乘法回歸建模就不再適用[16]。此時其參數估計值就不再具有最小方差的優良性，模型參數顯著性的t檢驗、回歸方程顯著性的F檢驗和擬合優度R2檢驗等都將失去可靠性，同時擬合參數具有較大的方差，回歸方程預測能力也將失效，結果導致肥效模型失去實用價值。【擬解決的關鍵問題】本文利用近年來“3414”試驗設計方案完成的早稻氮磷鉀田間肥效試驗結果，探討三元二次多項式肥效模型的異方差診斷方法和可行廣義最小二乘回歸建模法（FGLS）在消除或緩解異方差危害方面的應用效果，為提高三元肥效模型建模可靠性和成功率提供一種新方法。

1 材料與方法

1.1 早稻氮磷鉀肥效試驗資料的收集整理

近年來，在水稻測土配方施肥工作中，福建在福州市、寧德市、南平市、三明市、龍巖市、莆田市、漳州的龍海市以及平和縣、泉州的南安市等水稻主產區的主要項目縣先后完成了171個氮磷鉀田間肥效試驗。這些試驗均采用“3414”設計方案[17]，即（1）N0P0K0；（2）N0P2K2；（3）N1P2K2；（4）N2P0K2；（5）N2P1K2；（6）N2P2K2；（7）N2P3K2；（8）N2P2K0；（9）N2P2K1；（10）N2P2K3；（11）N3P2K2；（12）N1P1K2；（13）N1P2K1；（14）N2P1K1。其中，“2”水平為試驗前的當地推薦施肥量，“0”水平表示不施肥，“1”水平和“3”水平的施肥量分別為“2”水平的50%和150%。本研究選擇8個代表性試驗點的田間試驗結果來探討三元二次多項式肥效模型的異方差診斷方法及其應對技術。這些試驗點基礎土樣的主要理化性狀采用常規方法[18]測定，處理（6）的施肥量采用當地推薦施肥量（表 1），試驗設計方案和田間管理與文獻[17]相同。8個試驗點的產量結果見表2。

1.2 多元線性回歸模型的異方差及其診斷

經典線性回歸分析針對模型隨機誤差項u?的一個重要假設是u?在不同觀測點具有相同方差，即誤差項u?的方差是一個等于σ2的常數[16]。但是，這種等方差假設在實際問題中往往難以得到滿足。當回歸模型存在嚴重異方差時，普通最小二乘法就喪失了實用價值[19-20]。異方差的定性診斷方法是模型殘差分析[16]。當回歸模型的誤差項具有同方差時，殘差圖上n個點的散布是隨機的，無任何規律。如果存在異方差，殘差圖上點的分布就會呈現出一定的變化趨勢，如殘差值增大或減小等。

因線性回歸模型的異方差有多種不同表現形式[19]，其定量診斷方法也有多種不同方法。本研究選擇常用的懷特（White）檢驗、帕克（Park）檢驗和冠因克-巴塞特檢驗（Koenker-Basett，KB檢驗）等，具體使用方法可參考相關專著[19]。

表1 早稻8個代表性試驗點的供試土壤主要理化性狀及其處理（6）氮磷鉀施肥量Table 1 Main physical and chemical properties of soils in 8 representative experience sites and fertilization rate in 6th treatment in early rice

表2 早稻8個代表性試驗點各施肥處理的稻谷產量Table 2 Rice grain yield of each fertilizer treatment in 8 representative experiments in early rice

1.3 可行廣義最小二乘回歸建模法

目前常用的三元二次多項式肥效模型可用矩陣形式表示為：

式中，Y?表示擬合產量；X表示肥效試驗設計矩陣，β表示肥效模型參數向量，u?表示模型擬合殘差向量。用普通最小二乘法（OLS）進行回歸分析時，模型參數的計算式為：式中，Y表示各施肥處理的試驗產量向量。針對異方差問題，數理統計學家已經提出了消除或緩解線性回歸模型異方差的許多方法[19,21]，但應用較廣泛的方法是可行廣義最小二乘法（FGLS）或稱之為可行加權最小二乘法（FWLS）。這種方法是對原回歸模型（1）式加權，使之變成一個新的不存在異方差性的模型，然后采用OLS估計其參數[19]。因此，FGLS就是對加了權重的殘差平方和實施OLS，即：

式中，Wi為給定的第i個處理產量的權數。如何獲得各處理產量的 Wi？一種可行的方法是對原模型進行OLS估計，得到隨機誤差項的近似估計量u?，以此構建權數估計量。即：，進而構造權矩陣W，然后用OLS估計新模型。記參數估計量為，即FGLS的參數估計量計算式為：

這就是原模型（1）式的可行廣義最小二乘估計，它是無偏和有效的估計量[16]。因為?*β是β的最佳線性無偏估計量，因此，FGLS的總離差平方和SST*、回歸平方和 SSR*以及殘差平方和 SSE*的計算式分別為：

進而可以計算得到FGLS回歸模型的F值和擬合優度R2等統計檢驗參數。上述的具體數學原理和計算方法可參閱相關專著[16,21-22]。

1.4 肥效模型典型性判別方法

假設三元二次多項式肥效模型的數學形式如下：

式中，N、P、K 分別表示 N、P2O5和 K2O 施肥量（kg·hm-2）；ε為擬合殘差（kg·hm-2）；Y?為擬合產量（kg·hm-2）。根據無約束最優化方法[6]，若函數f (N,P, K)在點X處的一階梯度向量g(X) （其中，X=（N,P, K））等于零向量；同時，若肥效函數的Hesse矩陣G(X )為負定，則此函數有極大值。對式（8）的三元二次多項式肥效模型而言，Hesse矩陣G(X)的各階主子式假設為G1、G2和G3，則它們的行列式值分別為：G1=2b4；G2=4b4b5-b72；G3=2(4b4b5b6+b7b8b9-b4b92-b5b82-b6b72)。因此，各類三元肥效模型的典型性判別方法如下：①若g(X* )=0，且G1＜0，G2＞0，G3＜0，則矩陣G(X )為負定，肥效模型具有全局最高產量點。同時，若g(X*)=0所代表的點X*對應的N、P、K數值均落在試驗設計的施肥量范圍內，則該肥效模型屬于典型式。②若G(X*)為負定，但點X*對應的N、P、K數值有一個或一個以上超過試驗設計施肥量，則該肥效模型屬于推薦施肥量外推的非典型式。③若g(X*)=0，而且 G1＞0，G2＞0，G3＞0，則 Hesse矩陣G (X )為正定，此函數有全局最低產量點，肥效型屬于存在最低產量點的非典型式。④若g(X*)=0，G1、G2和 G3不滿足正定和負定的條件而且不等于零，則Hesse矩陣為不定，點X*為駐點，該肥效模型屬于無最高產量點的非典型式。

1.5 肥效模型異方差檢驗和參數估計的計算機實現

三元二次多項式肥效模型的異方差檢驗、OLS和FGLS回歸建模均涉及到繁雜的數學計算，本研究采用MATLAB軟件來完成相關工作，其中，OLS法回歸診斷調用 regstats功能函數，OLS回歸分析調用regress功能函數，FGLS法則調用fgls功能函數。根據上述統計學原理，在 MATLAB軟件平臺上編寫計算程序完成異方差檢驗及其繪圖。

2 結果

2.1 OLS法回歸分析及其肥效模型典型性判別

根據表1的8個試驗點施肥量和表2相應的14個處理試驗產量結果，應用OLS進行回歸分析（表3）。結果表明，1號和2號的兩個試驗點未達到統計顯著水平，試驗資料不可用；3—8號的6個試驗點，回歸模型均達到統計顯著水平。因為通過顯著性檢驗的三元二次多項式肥效模型有典型式和3種不同類別的非典型式[6]，因而，通過顯著性檢驗的 6個回歸模型還需進一步進行典型性判別（表 3），以便確定推薦施肥方法。結果顯示，盡管3號試驗點肥效模型的一次項參數符號為正，二次項參數符號為負，滿足了植物營養的一般規律，但是，該肥效模型不存在最高產量點，邊際產量導數法推薦施肥的結果不可靠。4號試驗點有最高產量點，邊際產量導數法推薦施肥量沒有外推，但是，P的一次項參數為負數，屬于參數符號不合理的非典型式。5號試驗點有最高產量點，N、P、K的一次項和二次項參數的符號滿足植物營養的一般規律，但鉀肥的最高施肥量超過試驗設計施肥量，屬于推薦施肥結果遠外推的非典型式。6—8號試驗點的肥效模型存在最高產量點，模型一次項和二次項參數的符號正常，邊際產量導數法推薦施肥量在試驗設計施肥量范圍內，屬于典型肥效模型。

因此，采用目前常用的OLS回歸建模，1號和2號試驗點屬于未達統計顯著水平的肥效模型，3—5號試驗點建立的肥效模型屬于非典型式，只能棄之不用；只有 6—8號試驗點的肥效模型滿足植物營養學的一般肥效規律，屬于三元典型肥效模型。

2.2 肥效模型的異方差及其統計檢驗

OLS回歸建模要求模型隨機誤差項具有相同方差，否則，回歸建模結果就不可靠[16]，在應用上易對肥效試驗結果產生誤判。為此，需對表2的8個試驗資料進行回歸診斷。以各處理的試驗產量為X軸，OLS擬合殘差為Y軸，繪制圖1的殘差散點圖。結果表明，無論是統計未達顯著水平的1號試驗點肥效模型（圖1-a），還是3號試驗點的非典型肥效模型（圖1-b）或者8號試驗點的典型肥效模型（圖1-c），隨著水稻產量水平提高，其擬合殘差均呈現逐漸散開的漏斗狀分布趨勢，表明都不滿足隨機分布的假設。結果預示著14個處理產量的擬合殘差方差不相等，且隨著產量水平提高而增大，但1號試驗點隨著稻谷產量水平增加，擬合殘差的增加絕對值最大，3號試驗點則其次，8號試驗點擬合殘差絕對值的變化幅度最小。

表3 三元二次多項式肥效模型的OLS回歸分析Table 3 OLS regression analysis of ternary quadratic polynomial’s fertilizer response model

圖1的擬合殘差圖給出了三元肥效模型異方差的定性分析結果，要回答異方差的嚴重程度，還需要做定量分析（表 4）。由于異方差在線性回歸模型中可能有不同的表現形式[19]，在不同的定量檢驗方法中，只要有一種檢驗結果達到統計顯著水平，就說明回歸模型存在該種形式的異方差。結果表明，擬合殘差絕對值變化幅度最小的8號試驗點的OLS回歸模型異方差未達顯著水平，其余7個試驗點建立的OLS回歸方程至少有一種異方差模型達到統計顯著水平，即存在嚴重的異方差。

2.3 可行廣義最小二乘法對肥效模型建模結果的影響

根據表4的異方差檢驗結果，1—7號試驗點肥效模型存在統計顯著水平的異方差，表3的OLS回歸建模結果就失去了可靠性。為此，采用可行廣義最小二乘法回歸建模（表5）。

圖1 三元肥效模型擬合殘差圖Fig. 1 Model fitting residual graph of ternary fertility efficiency

表4 肥效模型異方差的統計檢驗Table 4 Statistical test of model’s heteroscedasticity

表5 可行廣義最小二乘法回歸建模結果Table 5 Modeling results using feasible generalized least squares regression

結果顯示，可行廣義最小二乘法明顯改善了肥效模型的建模效果。1號和2號試驗點的F值和擬合優度R2由OLS法的未達到顯著水平提高到極顯著水平，而且肥效模型轉化為典型式。其他6個試驗點的肥效模型F值和擬合優度R2同樣得到顯著的提升；3號、4號和5號試驗點分別由OLS法的無最高產量點的、P項參數符號不合理、推薦施肥量外推的非典型式轉化為典型式。但是，OLS法屬于典型式的6號試驗點，盡管模型存在異方差，若采用FGLS法則導致回歸模型由典型式變為非典型式，結果反而變劣。對異方差未達顯著水平差異的8號試驗點，采用FGLS回歸也能提高擬合優度，模型同樣是典型式。

因此，在8個代表性試驗點中，采用FGLS法回歸建模，顯著提高了肥效模型的擬合優度和典型式的出現幾率，明顯地提高了建模成功率。

FGLS法為什么能提高建模成功率呢？比較1號、3號和8號試驗點的兩種參數估計方法，FGLS回歸建模明顯地降低了肥效模型各個參數的方差（表 6）；表5和表3的均方誤差（MSE）結果也表明，FGLS回歸建模大幅度降低了誤差均方。結果說明FGLS大幅度提高了肥效模型的擬合精度和模型預測能力，從而提高了肥效模型的可靠性。

2.4 消除異方差后的肥效模型推薦施肥量

根據FGLS回歸建模得到的典型肥效模型，結合N、P2O5、K2O和稻谷價格分別為4.3、5.0、6.0和2.0元/kg的市場均價，應用邊際產量導數法計算推薦施肥量（表7）。OLS回歸建模能得到典型肥效模型的7號和8號試驗點，兩種建模方法所得肥效模型的推薦施肥量和預測產量水平十分相近；1—5號試驗點，由OLS回歸建模未達顯著水平或是非典型肥效模型，用FGLS回歸建模轉化成了典型肥效模型，其推薦施肥量和預測產量水平也符合現有經驗認識，尤其是，1號、3號和5號試驗點的基礎土壤堿解氮含量分別為186、213、275 mg·kg-1（表 1），接近或超過早稻高產臨界指標[17]，3個試驗點的氮肥推薦用量明顯低于常規推薦水平。結果表明，FGLS回歸建模的推薦施肥量是可靠的。

表6 OLS與FGLS回歸的參數方差比較Table 6 Parameter variance compare between the regression methods of OLS and FGLS

表7 OLS肥效模型與FGLS肥效模型推薦施肥量的比較Table 7 Compare the recommended fertilization rate between the fertilizer response model of OLS and FGLS

2.5 肥效模型異方差出現幾率及其FGLS法效果評價

為了評價FGLS對三元肥效模型的建模效果，利用福建省早稻測土配方施肥項目在全省各地主要項目縣完成的171個“3414”設計方案的氮磷鉀田間肥效試驗結果，逐個進行回歸分析，統計結果見表8。FGLS大幅度提高了三元肥效模型的擬合優度，每個試驗點的R2和F值均達到統計顯著水平，而OLS則有16%試驗點的肥效模型未達統計顯著水平。

模型典型性判別分析表明，在3種非典型肥效模型類型中，無論是OLS還是FGLS回歸建模，無最高產量點的非典型肥效模型所占比例高達40%左右，其次是參數符號不合理的非典型肥效模型所占試驗點比例大致15%，而氮磷鉀推薦施肥量屬于外推的非典型肥效模型所占比例 6%左右，兩種方法的結果差異很小；OLS回歸建模的典型肥效模型占試驗點數的27.68%，而FGLS回歸建模的典型肥效模型比例達到34.99%，比OLS回歸建模提高了7.31個百分點。

表8 不同參數估計方法對肥效模型建模結果的影響Table 8 Effects of modeling results by using different methods of estimate parameters in fertilizer response model

3 討論

3.1 多項式肥效模型產生異方差的原因

在經典回歸建模理論中，針對線性回歸方程隨機誤差項?u的一個重要假設是?u在不同觀測點具有同方差性。異方差將導致嚴重后果[16]，包括：（1）盡管OLS參數估計的無偏性仍然成立，并且基于此的產量預測也是無偏的，但參數估計值的方差不再是最小的，即最小二乘法失去了有效性，對大樣本也是如此。（2）由于異方差的影響，造成難以正確估計模型參數的標準誤，且無法辨別是正的還是負的偏差，導致回歸模型的t 統計量和F統計量失去了可靠性。（3）由于參數估計量不是有效的，從而對產量的預測也將不是有效的。二次多項式肥效模型是建立在經典正態線性回歸模型理論基礎上的，要正確建立肥效模型就必須遵守回歸分析理論的基本假設[19]。長期以來，我們在建立作物肥效模型時，大都忽略了回歸分析理論有關同方差性等條件的基本假設，這是目前建模成功率低[4-6]的重要原因之一。

總結本研究結果并結合專業特點，多項式肥效模型產生異方差的主要原因有：（1）隨著試驗設計和田間試驗管理水平的提高，肥效模型的誤差方差σ2可能會減少。反之，若田間試驗操作或管理出現失誤，可能導致部分處理出現異常產量結果，使模型誤差方差σ2出現非常數；（2）由于田間土壤肥力水平存在空間異質性，可能導致部分試驗處理出現異常產量值，導致出現異方差；（3）多項式肥效模型屬于經驗回歸模型，對某些試驗點來說可能導致模型設定不正確；（4）在多數情況下，當土壤養分肥力較低時，作物產量處于低產水平，此時通過施肥可迅速提高產量，而因施肥造成減產或平產的可能空間很小，即產量方差較小；但在高產階段，通過施肥可能繼續增產，也可能不再增產，甚至可能因過量施肥造成減產，即產量變化的可能空間較大，此時方差必然較大。從這個意義上講，依據不同施肥量和產量建立的經驗肥效模型必然帶有異方差性質。福建的171個早稻氮磷鉀肥效試驗結果表明，這種異方差性達到統計顯著水平的試驗點占到試驗總數的25%左右（表8）。因此，通過控制試驗田土壤空間異質性、提高試驗管理水平，降低試驗產量異常值的出現幾率，是消減異方差和提高建模成功率的重要途徑。

3.2 多項式肥效模型異方差的檢驗方法

針對異方差的嚴重危害性，在現代回歸分析理論中已經研究提出了多種異方差定量診斷方法[19]。本研究選擇懷特（White）檢驗、冠因克-巴塞特（KB）檢驗、帕克（Park）檢驗等3種常用的異方差定量檢驗方法，探討三元肥效模型異方差性質。結果表明，不同試驗點建立的三元二次多項式肥效模型，誤差項的異方差有不同的表現形式。由于目前還難以明確在什么試驗條件下會出現哪種形式的異方差，因而，在實用上，通常同時選用幾種不同異方差檢驗方法進行檢驗。只要有一種檢驗方法達到統計顯著水平，就說明肥效模型存在該種形式的異方差。

當前，構建異方差定量檢驗方法的共同思路是，檢驗隨機誤差項的方差與回歸模型自變量或因變量觀測值之間的相關性及其數學形式[19]，各種檢驗方法正是在這個思路下發展起來的。本文選用3種常用的檢驗方法，它們各有自己的特色和優點。其中，懷特檢驗不僅能夠檢驗異方差的存在，同時在多變量的情況下，還能夠判斷出是哪一個變量引起的異方差。帕克檢驗的優點是不但能確定有無異方差性，而且還能給出異方差的具體函數形式。冠因克-巴塞特檢驗的一個優點是即使回歸模型中的誤差項不是正態分布的，它仍能適用。因此，這三種方法都是檢驗肥效模型異方差的有效方法。

3.3 多項式肥效模型異方差的修正方法

用普通最小二乘法估計肥效模型參數時，對參數的估計是以產量擬合殘差平方和最小為條件的。在最小化過程中，對每個施肥處理的產量擬合殘差平方所提供的信息的重要程度是同等看待的，它們在決定肥效模型參數的過程中所起的作用是相同的，或者說取了相同的權數，即權數都為 1。這在同方差條件下，因不同施肥處理的作物產量水平在偏離均值程度上是相同的，這樣做是合理的。但是，在異方差條件下，由于各施肥處理的產量水平在偏離均值程度上相差很大，這樣做就未必合理。

本研究針對三元二次多項式肥效模型異方差問題，采用FGLS回歸建模。該法是對原模型加權，即：對較小 OLS擬合殘差平方的施肥處理產量賦予較大的權數，對較大OLS擬合殘差平方的施肥處理產量賦予較小的權數，使回歸模型轉化成不存在異方差性的新模型，然后再次采用OLS估計其參數。結果表明，與OLS回歸建模方法相比，FGLS方法能將絕大部分試驗點的試驗結果由 OLS方法的未達統計顯著水平轉化為達到統計顯著水平，或者將部分試驗點的非典型式模型轉化為典型肥效模型（表3和表5），明顯提高了建模成功率。

在對付線性回歸模型異方差上，除了FGLS外，“OLS+穩健標準誤”也是常用方法[22]。該法是在異方差情況下，仍然采用OLS，但對OLS估計量的標準差進行修正。與OLS估計比較，肥效模型參數估計值沒有變化，但是參數估計值的方差和標準差變化明顯。經過這種修正的模型參數標準差，模型參數的顯著性檢驗和模型產量預測等方面就變為有效的。那么，究竟是使用“OLS+穩健標準誤”還是FGLS？從肥效模型角度看，我們總是希望將OLS建模出現的非典型模型盡可能多地轉化為典型肥效模型。由于“OLS+穩健標準誤”方法的模型參數沒有變化，顯然不能實現這一目標，而FGLS方法則能把部分非典型肥效模型轉化為典型式。從這個意義上講，FGLS方法具有優勢。同時，表7的建模結果使我們注意到，OLS回歸建模是典型肥效模型的試驗資料，若采用FGLS回歸建模，可能出現非典型式，使結果變劣。因此，在實際工作中，我們也不應對異方差反應過度[19]。如果OLS方法能得到典型肥效模型，就無需考慮其他修正方法。

福建 171個早稻氮磷鉀肥效試驗資料的總結結果表明，FGLS回歸建模的三元典型肥效模型比例僅達到34.99%，只比OLS回歸建模提高7.31個百分點。究其原因，多項式肥效模型不僅存在異方差，而且還存在嚴重的多重共線性[23]問題，制約了OLS法的有效性[19]，對FGLS回歸建模也不例外。只有較好地克服了這兩個建模問題，才可能較大幅度地提高建模成功率。如何同時克服多項式肥效模型的異方差和多重共線性問題，提出最優建模策略，則有待于未來的進一步研究。

4 結論

根據171個早稻“3414”設計的氮磷鉀田間肥效試驗結果，現代回歸分析理論的可行廣義最小二乘法建模技術是消除三元二次多項式肥效模型異方差的有效方法。該法明顯減小模型各個參數的方差，改善了肥效模型的擬合精度和模型預測能力，進而提高了田間肥效試驗資料的建模成功率。

[1] VALKAMA E, UUSITALO R, TURTOLA E. Yield response models to phosphorus application: a research synthesis of Finnish field trials to optimize fertilizer P use of cereals. Nutrient Cycling in Agroecosystems,2011, 91(1): 1-15.

[2] GRWGORT M C, ZORITA M D, DARDANELLI J, Bongiovanni R G.Regional model for nitrogen fertilization of site-specific rainfed corn in haplustolls of the central Pampas, Argentina. Precision Agriculture,2011, 12(6): 831-849.

[3] 毛達如, 張承東. 推薦施肥技術中施肥模型與試驗設計的研究. 土壤通報, 1991, 22(5): 216-218.MAO D R , ZHANG C D. Fertilization model and the experiment design of fertilization recommendation technology. Chinese Journa1 of Soil Science, 1991, 22(5): 216-218. (in Chinese)

[4] 王興仁. 二元二次肥料效應曲面等產線圖在科學施肥中的位置（一）. 土壤通報, 1985, 16(1): 30-34.WANG X R. The positions of contour map of yield of binary quadratic fertilizer response curve in scientific fertilization (part 1).Chinese Journa1 of Soil Science, 1985, 16(1): 30-34. (in Chinese)

[5] 王興仁. 二元二次肥料效應曲面等產線圖在科學施肥中的位置（二）. 土壤通報, 1985, 16(2): 86-88.WANG X R. The positions of contour map of yield of binary quadratic fertilizer response curve in scientific fertilization (part 2).Chinese Journa1 of Soil Science, 1985, 16(2): 86-88. (in Chinese)

[6] 章明清, 林仁塤, 林代炎, 姜永. 極值判別分析在三元肥效模型推薦施肥中的作用. 福建農業學報, 1995, 10(2): 54-59.ZHANG M Q, LIN R X, LIN D Y, JIANG Y. Function of distinguish analysis on extreme value in recommendatory fertilization for three-fertilizer efficiency model. Fujian Journal of Agricultural Sciences, 1995, 10(2): 54-59. (in Chinese)

[7] CERRATO M E, BLACKMER A M. Comparison of models for describing corn yield response to nitrogen fertilizer. Agronomy Journal, 1990, 82(1): 138-143.

[8] 楊靖一, WADSWORTH G A, GREENWOOD D J. 三種蔬菜N肥效應曲線的比較研究. 植物營養與肥料學報, 1995, 1(1): 71-78.YANG J Y, WADSWORTH G A, Greenwood D J. Comparative study on three kinds of vegetables response to N fertilizer curve. Plant Nutrition and Fertilizer Science, 1995, 1(1): 71-78. (in Chinese)

[9] 陳新平, 周金池, 王興仁, 張福鎖, 寶德俊. 小麥-玉米輪作制中氮肥效應模型的選擇——經濟和環境效益分析. 土壤學報, 2000,37(3): 346-354.CHEN X P, ZHOU J C, WANG X R, ZHANG F S, BAO D J.Economic and environmental evaluation on models for describing crop yield response to nitrogen fertilizers at winter-wheat and summer-corn rotation system. Acta Pedologica Sinica, 2000, 37(3):346-354. (in Chinese)

[10] 吳平, 陶勤南. 氮磷鉀肥效模型的研究及其應用. 浙江農業大學學報, 1989, 15(4): 383-388.WU P, TAO Q N. Study on N, P and K fertilizer response model and its application. Journal of Zhejiang Agricultural University, 1989,15(4): 383-388. (in Chinese)

[11] 章明清. 米氏肥效方程的灰色建模法. 土壤通報, 1997, 28(1):18-19.ZHANG M Q. Gray modeling method of Mitsherlich fertilizer response equations. Chinese Journa1 of Soil Science, 1997, 28(1):18-19. (in Chinese)

[12] 章明清, 徐志平, 姚寶全, 林瓊, 顏明娟. Monte Carlo 法在多元肥效模型參數估計和推薦施肥中的應用. 植物營養與肥料學報, 2009,15(2): 366-373.ZHANG M Q, XU Z P, YAO B Q, LIN Q, YAN M J. Using Monte Carlo method for parameter estimation and fertilization recommendation of multivariate fertilizer response model. Plant Nutrition and Fertilizing Science, 2009, 15(2): 366-373. (in Chinese)

[13] COLWELL J D. The erivation of fertilizer recommendations for crops in non-uniform environment. Crop Quality and Economy, 1974:936-961.

[14] 王興仁, 陳倫壽, 毛達如, 張才東. 分類回歸綜合法及其在區域施肥決策中的應用. 土壤通報, 1989, 20(1): 17-21.WANG X R, CHEN L S, MAO D R, ZHANG C D. Classification of regression synthesis and its application for regional fertilization decision-making. Chinese Journa1 of Soil Science, 1989, 20(1): 17-21.(in Chinese)

[15] 毛達如, 張承東. 多點肥料效應函數的動態聚類方法. 北京農業大學學報, 1991, 17(2): 49-54.MAO D R, ZHANG C D. Cluster analysis of quadratic response function on the fertilizer dispersed experiment. Acta Agriculturae Universitatis Pekinensis, 1991, 17(2): 49-54. (in Chinese)

[16] MONTGOMERY D C, PECK E A, VINING G G. 線性回歸分析導論. 第5版. 王辰勇, 譯. 北京: 機械工業出版社, 2016: 91-268.MONTGOMERY D C, PECK E A, VINING G G. Introduction to Linear Regression Analysis. 5th ed. WANG C Y, trans. Beijing: China Machine Press, 2016: 91-268. (in Chinese)

[17] 李娟, 章明清, 孔慶波, 姚寶全, 顏明娟. 福建早稻測土配方施肥指標體系研究. 植物營養與肥料學報, 2010, 16(4): 938-946.LI J, ZHANG M Q, KONG Q B, YAO B Q, YAN M J. Soil testing and formula fertilization index for early rice in Fujian Province. Plant Nutrition and Fertilizer Science, 2010, 16(4): 938- 946. (in Chinese)

[18] 魯如坤. 土壤農業化學分析方法. 北京: 中國農業科技出版社,2000: 146-196.LU R K. Soil Agricultural Chemical Analysis Method. Beijing: China's Agricultural Science and Technology Press, 2000: 146-196. (in Chinese)

[19] GUJARATI D N, PORTER D C. 計量經濟學基礎 (上冊). 第5版.費劍平, 譯. 北京: 中國人民大學出版社, 2011: 313-399.GUJARATI D N, PORTER D C. Econometrics Foundation. 5th ed.FEI J P, trans. Beijing: China Renmin University Press, 2011: 313-399.(in Chinese)

[20] 尹光霞. 多元線性回歸模型中的異方差性問題. 湖北大學學報(自然科學版), 2003, 25(2): 121-125.YIN G X. Problems on heteroscedasticity in multilinear regression models. Journal of Hubei University (Natural Science Edition), 2003,25(2): 121-125. (in Chinese)

[21] 何曉群，劉文卿. 應用回歸分析. 3版. 北京: 中國人民大學出版社,2013: 158-170.HE X Q, LIU W Q. Applied Regression Analysis. 3rd ed. Beijing:China Renmin University Press, 2013: 158-170. (in Chinese)

[22] 陳強. 高級計量經濟學及Stata應用. 2版. 北京: 高等教育出版社,2015: 87-99.CHEN Q. Advanced Econometric and Stata Applications.2nd ed.Beijing: Higher Education Press, 2015: 87-99. (in Chinese)

[23] 章明清, 李娟, 孔慶波, 嚴芳. 作物肥料效應函數模型研究進展與展望. 土壤學報, 2016, 53(6): 1143-1156.ZHANG M Q, LI J, KONG Q B, YAN F. Progress and prospect of the study on crop-response-to -fertilization function model. Acta Pedologica Sinica, 2016, 53(6): 1143-1156. (in Chinese)