李 娟 章明清? 許文江 孔慶波 姚寶全

（1 福建省農(nóng)業(yè)科學(xué)院土壤肥料研究所，福州 350013）

（2 福建省亞熱帶植物研究所，福建廈門 361000）

（3 福建省農(nóng)田建設(shè)與土壤肥料技術(shù)推廣總站，福州 350003）

肥料效應(yīng)函數(shù)是實(shí)現(xiàn)計(jì)量施肥的主要技術(shù)手段之一，其中，二次多項(xiàng)式是研究和應(yīng)用最多的肥效模型［1-4］。但在施肥實(shí)踐中，應(yīng)用普通最小二乘法（OLS）建立的二元或三元二次多項(xiàng)式肥效模型出現(xiàn)大量的非典型式［5-7］，諸如“最高產(chǎn)量”不高、“最佳產(chǎn)量”不佳、模型一次項(xiàng)系數(shù)為負(fù)數(shù)或二次項(xiàng)系數(shù)為正數(shù)、推薦施肥量遠(yuǎn)外推等令人困惑的異常結(jié)果。針對(duì)肥效模型在應(yīng)用中存在的問題，國(guó)內(nèi)外學(xué)者就肥效模型選擇與其適用性［8-10］、試驗(yàn)設(shè)計(jì)與參數(shù)估計(jì)方法［11-13］、類特征肥效模型構(gòu)建［14-16］及非典型式推薦施肥優(yōu)化方法［13］等方面進(jìn)行了深入研究，提出了諸多改進(jìn)措施，但對(duì)肥效模型設(shè)定本身是否存在問題及其應(yīng)對(duì)方法則鮮見研究報(bào)道［4］。當(dāng)前廣泛應(yīng)用的“3414”設(shè)計(jì)方案的田間試驗(yàn)結(jié)果，既使采用蒙特卡洛（Monter Carlo）建模法［13］，二元和三元二次多項(xiàng)式典型肥效模型的比例也僅有56.7%和37.3%，這種狀況嚴(yán)重制約了肥效模型的計(jì)量精確性和應(yīng)用價(jià)值。

初步研究［4］表明，二次多項(xiàng)式肥效模型存在嚴(yán)重的多重共線性，制約了OLS法回歸建模的有效性。為此，本研究利用福建省近年來 “3414”試驗(yàn)設(shè)計(jì)方案完成的早稻氮磷鉀田間肥效試驗(yàn)結(jié)果，探討三元二次多項(xiàng)式肥效模型多重共線性診斷方法，以及主成分回歸（PCR）技術(shù)與其在消除多重共線性危害上的應(yīng)用效果，旨在為提高氮磷鉀田間肥效試驗(yàn)的建模成功率提供一種新方法。

1　材料與方法

1.1　早稻氮磷鉀“3414”田間肥效試驗(yàn)資料的收集整理

近年來，福建省在水稻測(cè)土配方施肥工作中，在福州市、寧德市、南平市、三明市、龍巖市、莆田市、漳州的龍海市以及平和縣、泉州的南安市等水稻主產(chǎn)區(qū)，土肥技術(shù)力量較強(qiáng)的主要項(xiàng)目縣的早稻上先后完成了171個(gè)氮磷鉀田間肥效試驗(yàn)。這些試驗(yàn)均采用“3414”設(shè)計(jì)方案［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。其中，“0”水平表示不施肥，“2”水平為試驗(yàn)前的當(dāng)?shù)氐租浲扑]施肥量，“1”水平和“3”水平的施肥量分別為“2”水平的50%和150%。根據(jù)肥效模型建模過程中可能出現(xiàn)的未達(dá)統(tǒng)計(jì)顯著水平、典型肥效模型和三種不同形式的非典型肥效模型［7］等幾種類別，本研究選擇7個(gè)代表性早稻試驗(yàn)結(jié)果（選擇依據(jù)主要是：根據(jù)普通最小二乘法建立三元肥效模型可能出現(xiàn)的幾種情形，在兼顧主要稻田土壤類型基礎(chǔ)上，分別選擇了OLS法未通過顯著性檢驗(yàn)的1個(gè)代表性試驗(yàn)點(diǎn)、OLS法產(chǎn)生的3種不同類型非典型式的3個(gè)代表性試驗(yàn)點(diǎn)、OLS法為典型式的3個(gè)代表性試驗(yàn)點(diǎn)）來探討三元二次多項(xiàng)式肥效模型的多重共線性診斷方法及其主成分回歸技術(shù)，并利用171個(gè)早稻氮磷鉀田間肥效試驗(yàn)結(jié)果評(píng)價(jià)主成分回歸的應(yīng)用效果。這7個(gè)試驗(yàn)點(diǎn)涵蓋了灰泥田、黃泥田和灰沙田等福建稻田的主要土壤類型，采用常規(guī)方法［18］測(cè)定的土壤理化性狀（表1）也大致反映了各肥力指標(biāo)的主要變化范圍。7個(gè)代表性試驗(yàn)點(diǎn)N2P2K2處理的施肥量見表1，相關(guān)試驗(yàn)產(chǎn)量見表2。供試水稻品種采用當(dāng)?shù)卮竺娣e種植的良種，試驗(yàn)設(shè)計(jì)方案和田間管理措施與文獻(xiàn)［17］相同。

表1　早稻7個(gè)代表性試驗(yàn)點(diǎn)土壤主要理化性狀及N2P2K2處理施肥量Table 1　Main physical and chemical properties of the soils in the seven representative experiment sites and fertilizer application rate of their Treatment N2P2K2

表2　早稻7個(gè)代表性試驗(yàn)點(diǎn)各施肥處理的稻谷產(chǎn)量Table 2 Early rice yields in the seven representative trial sitesrelative to treatment(kg hm-2)

1.2　多項(xiàng)式肥效模型的多重共線性診斷方法

假設(shè)三元二次多項(xiàng)式肥效模型的數(shù)學(xué)形式如下：

式中，N、P、K分別表示N、P2O5和K2O施肥量，kg hm-2；ε為擬合殘差，kg hm-2；Y為擬合產(chǎn)量，kg hm-2。在應(yīng)用OLS回歸分析時(shí)，根據(jù)“3414”設(shè)計(jì)方案14個(gè)施肥處理的氮、磷、鉀施肥量與式（1）模型對(duì)應(yīng)的9個(gè)回歸變量N、P、K、N2、P2、K2、NP、NK和PK可組成一個(gè)試驗(yàn)設(shè)計(jì)矩陣X：

回歸變量： N P K N2P2K2NP NK PK

矩陣X共有14行和9列，每一行表示一個(gè)試驗(yàn)處理，每一列表示一個(gè)回歸變量。在回歸分析的經(jīng)典理論中，假設(shè)式（1）模型的9個(gè)回歸變量在兩兩之間不存在線性相關(guān)［19］。但在實(shí)際問題中，這些回歸變量之間總是會(huì)存在一定程度的相關(guān)性。在數(shù)學(xué)上，如果存在某些常數(shù)c1、c2…c9與相應(yīng)各回歸變量的線性組合等于常數(shù)c0這一關(guān)系近似成立，則表示變量間存在多重共線性。因此，一個(gè)常用但不完全合適的共線性度量方法是樣本相關(guān)系數(shù)平方即r2［20］。精確共線性對(duì)應(yīng)于r2=1，非共線性對(duì)應(yīng)于r2=0。當(dāng)r2越接近于1，近似共線性就越強(qiáng)。

度量多重共線性嚴(yán)重程度的一個(gè)重要指標(biāo)是方陣XTX的條件數(shù)［19-20］，即：

式中，k為條件數(shù)，λmax、λmin分別為方陣XTX的最大特征值和最小特征值。一般應(yīng)用經(jīng)驗(yàn)認(rèn)為，若k＜100，則多重共線性的程度很??；若100≤k≤1000，則存在中等強(qiáng)度的多重共線性；若k＞1000，則存在嚴(yán)重的多重共線性。

多重共線性的另一個(gè)重要診斷指標(biāo)是方差膨脹因子（VIF）［20］，即：

式中，(XTX)-1表示方陣XTX的逆矩陣，其中，X矩陣中的數(shù)據(jù)需事先進(jìn)行標(biāo)準(zhǔn)化處理；diag表示取矩陣的對(duì)角線元素，j表示第j個(gè)回歸系數(shù)（j=1，2…9）。模型中各個(gè)回歸系數(shù)的VIF度量了回歸變量之間的相關(guān)性對(duì)該項(xiàng)回歸系數(shù)方差的影響，一個(gè)或更多個(gè)較大的VIF就表明存在多重共線性。一個(gè)經(jīng)驗(yàn)法則［20］認(rèn)為，當(dāng)VIF超過5或10時(shí)，就表明回歸模型存在嚴(yán)重的多重共線性。

1.3　多項(xiàng)式肥效模型的主成分回歸分析

當(dāng)多項(xiàng)式肥效模型存在嚴(yán)重多重共線性時(shí)，一個(gè)有效解決辦法是重構(gòu)回歸變量，使之成為無線性相關(guān)的新變量。主成分分析為重構(gòu)回歸變量提供了可靠的數(shù)學(xué)方法，它將三元二次多項(xiàng)式肥效模型的9個(gè)回歸變量轉(zhuǎn)換成互不相關(guān)的9個(gè)主成分，然后選擇合適的幾個(gè)主成分進(jìn)行OLS回歸分析。

因此，從試驗(yàn)處理的原始施肥量數(shù)據(jù)出發(fā)，計(jì)算試驗(yàn)設(shè)計(jì)矩陣X的協(xié)方差矩陣S，然后求取S矩陣的特征值和特征向量。對(duì)式（1）而言，特征值共有9個(gè)，從大到小排列為λ1、λ2…λ9，與之對(duì)應(yīng)的特征向量為q1、q2…q9，每個(gè)q均為長(zhǎng)度為9的列向量。設(shè)向量x= (N, P, K, N2, P2, K2, NP, NK, PK)T，則9個(gè)主成分的矩陣Z表達(dá)式為：

描述某個(gè)回歸變量Xj與某個(gè)主成分Zi相關(guān)程度的樣本相關(guān)系數(shù)計(jì)算式［21-22］為：

r的平方即為決定系數(shù)。若r2(Xj,Zi)數(shù)值大，就表明第i個(gè)主成分主要反映了第j個(gè)回歸變量的信息，其結(jié)果可用于解釋各個(gè)主成分的專業(yè)含義。

因此，應(yīng)用中可根據(jù)各個(gè)主成分對(duì)試驗(yàn)設(shè)計(jì)矩陣的方差貢獻(xiàn)大小，選用適當(dāng)?shù)闹鞒煞謧€(gè)數(shù)結(jié)合14個(gè)試驗(yàn)處理產(chǎn)量結(jié)果進(jìn)行OLS回歸分析，并進(jìn)行相應(yīng)的統(tǒng)計(jì)檢驗(yàn)，這種方法稱為主成分回歸。

1.4　多重共線性診斷和主成分回歸的計(jì)算機(jī)實(shí)現(xiàn)

本研究采用MATLAB軟件進(jìn)行多重共線性診斷和主成分回歸分析。在多重共線性診斷中，回歸變量間的相關(guān)系數(shù)計(jì)算調(diào)用corrcoef功能函數(shù)；方陣XTX的條件數(shù)、逆矩陣和VIF計(jì)算分別調(diào)用cond、inv等功能函數(shù)；X矩陣的協(xié)方差矩陣S則調(diào)用cov功能函數(shù)。在主成分回歸中，主要調(diào)用了pca和regress等功能函數(shù)。上述計(jì)算的數(shù)學(xué)原理和有關(guān)功能函數(shù)使用方法參閱相關(guān)專著［21-22］。

2　結(jié) 果

2.1　多項(xiàng)式肥效模型的多重共線性診斷

根據(jù)表2的試驗(yàn)產(chǎn)量和表1的相關(guān)處理施肥量，利用OLS對(duì)式（1）進(jìn)行回歸建模，然后對(duì)各個(gè)回歸模型進(jìn)行統(tǒng)計(jì)檢驗(yàn)和模型典型性判別［7］，結(jié)果見表3。在7個(gè)代表性肥效模型中，1號(hào)試驗(yàn)點(diǎn)未達(dá)到統(tǒng)計(jì)顯著水平，其余6個(gè)試驗(yàn)點(diǎn)的肥效模型均滿足統(tǒng)計(jì)檢驗(yàn)要求。但是，2號(hào)試驗(yàn)點(diǎn)的P項(xiàng)系數(shù)為負(fù)數(shù)，3號(hào)試驗(yàn)點(diǎn)的推薦施肥量外推，4號(hào)試驗(yàn)點(diǎn)無最高產(chǎn)量點(diǎn)，均屬于非典型式，僅有第5、6、7號(hào)三個(gè)試驗(yàn)點(diǎn)的肥效模型滿足了植物營(yíng)養(yǎng)的一般肥效規(guī)律，屬于典型肥效模型。

1號(hào)至4號(hào)試驗(yàn)點(diǎn)的OLS回歸建模，出現(xiàn)了未通過顯著性檢驗(yàn)或者非典型式，是試驗(yàn)失誤造成的還是回歸方法不合理？為了明確這類問題，對(duì)表2的7個(gè)試驗(yàn)資料應(yīng)用式（2）和式（3）進(jìn)行多重共線性診斷。以2號(hào)試驗(yàn)點(diǎn)為例，與文獻(xiàn)［4］的表2結(jié)果相同，9個(gè)回歸變量的方差膨脹因子均遠(yuǎn)大于10，平均值達(dá)31.90；試驗(yàn)設(shè)計(jì)矩陣X的施肥量數(shù)據(jù)未經(jīng)標(biāo)準(zhǔn)化處理時(shí)，矩陣條件數(shù)達(dá)7628，即使對(duì)各處理施肥量進(jìn)行標(biāo)準(zhǔn)化轉(zhuǎn)換，其條件數(shù)也達(dá)548.8，均顯示肥效模型存在強(qiáng)烈的多重共線性。

進(jìn)一步對(duì)9個(gè)回歸變量進(jìn)行簡(jiǎn)單相關(guān)性分析，表明N與N2、P與P2、K與K2之間的相關(guān)系數(shù)高達(dá)0.942 9，N與NP或NK之間、P與NP或PK之間、K與NK或PK之間的相關(guān)系數(shù)達(dá)0.738 8，N2與NP或NK之間、P2與NP或PK之間、K2與NK或PK之間的相關(guān)系數(shù)達(dá)0.732 2，存在高度線性相關(guān)。因試驗(yàn)設(shè)計(jì)方案和肥效模型設(shè)定相同，其他6個(gè)試驗(yàn)點(diǎn)的試驗(yàn)設(shè)計(jì)矩陣條件數(shù)、方差膨脹因子、自變量相關(guān)系數(shù)等的相關(guān)結(jié)果與此一致，不再贅述。因此，三元二次多項(xiàng)式肥效模型存在嚴(yán)重的多重共線性，且這種共線性主要是由模型設(shè)定本身造成的。

表3　三元二次多項(xiàng)式肥效模型的OLS回歸分析和典型性判別Table 3　OLS regression analysis and typicality discrimination of TQPMs (ternary quadratic polynomial models for fertilizing effect)

2.2　回歸變量的主成分分析及其重構(gòu)

上述診斷表明，三元二次多項(xiàng)式肥效模型存在嚴(yán)重的多重共線性。因此，從重構(gòu)回歸變量入手來消除這種共線性危害就成為提高建模成功率的關(guān)鍵。根據(jù)主成分分析原理和式（4），從試驗(yàn)設(shè)計(jì)矩陣的原始數(shù)據(jù)出發(fā)，對(duì)肥效模型回歸變量進(jìn)行主成分分析，將式（1）肥效模型的9個(gè)回歸變量按照不同的線性組合重新劃分為無線性相關(guān)的9個(gè)主成分，并與X矩陣特征值對(duì)應(yīng)從大至小排列。同樣以2號(hào)試驗(yàn)點(diǎn)為例，表4的分析表明，具有最大特征值的第1主成分對(duì)矩陣X的方差貢獻(xiàn)率達(dá)73.82%，前7個(gè)主成分累加方差貢獻(xiàn)率達(dá)99.9%以上。

表4　第2號(hào)試驗(yàn)點(diǎn)設(shè)計(jì)矩陣X的9個(gè)主成分及其特征值Table 4 Nine principal components and their eigenvalues of the matrix X designed for Field Trial Site No.2

為了明確9個(gè)主成分在反映氮磷鉀施肥效應(yīng)上的差異，根據(jù)式（5）分別計(jì)算表4的9個(gè)主成分與肥效試驗(yàn)設(shè)計(jì)矩陣X的N、P、K、N2、P2、K2、NP、NK、PK等9個(gè)回歸變量的決定系數(shù)r2（表5）。根據(jù)r2值大小結(jié)合式（1）肥效模型，第1主成分主要反映了回歸變量N效應(yīng)，部分反映了NP交互效應(yīng)；第2主成分主要反映了回歸變量K效應(yīng)，部分反映了回歸變量N效應(yīng)；第3主成分主要反映了回歸變量P效應(yīng)，部分反映了回歸變量N效應(yīng)；等等。其他試驗(yàn)點(diǎn)的主成分特點(diǎn)亦有與此類似的結(jié)果，不再贅述。

表5 第2號(hào)試驗(yàn)點(diǎn)各回歸變量與9個(gè)主成分的決定系數(shù)（r2）Table 5　Determination coefficient (r2) between various regression variables and 9 principal componentfor field trial site No.2

2.3　主成分回歸的肥效模型及其推薦施肥量

根據(jù)表4各個(gè)主成分的N、P、K、N2、P2、K2、NP、NK、PK等回歸變量的線性組合系數(shù)和試驗(yàn)設(shè)計(jì)矩陣X的每一行對(duì)應(yīng)數(shù)值，分別計(jì)算14個(gè)處理的各自得分值，組成新的設(shè)計(jì)矩陣。由于前8個(gè)或前7個(gè)主成分分別反映了X矩陣方差信息的99.99%和99.90%以上，因此，新設(shè)計(jì)矩陣由前8個(gè)主成分或前7個(gè)主成分組成。

利用新設(shè)計(jì)矩陣與對(duì)應(yīng)的14個(gè)處理產(chǎn)量，進(jìn)行OLS回歸并轉(zhuǎn)化成式（1）函數(shù)形式（表6）。由于“3414”試驗(yàn)的處理數(shù)為14個(gè)，總自由度為13；采用前8個(gè)主成分或前7個(gè)主成分進(jìn)行主成分回歸，就相當(dāng)于回歸自由度分別為8和7；因而，建立的三元肥效模型的誤差自由度分別為5和6。結(jié)果表明，1號(hào)、2號(hào)和3號(hào)試驗(yàn)點(diǎn)所得三元肥效模型被轉(zhuǎn)化為滿足植物營(yíng)養(yǎng)學(xué)一般肥效規(guī)律的典型式；4號(hào)試驗(yàn)點(diǎn)采用8個(gè)主成分進(jìn)行回歸分析時(shí)，其結(jié)果與OLS回歸相同均屬于無最高產(chǎn)量點(diǎn)的非典型式，但采用7個(gè)主成分回歸即可得到三元典型肥效模型。5號(hào)和6號(hào)試驗(yàn)點(diǎn)的OLS回歸分析得到典型式，采用前8個(gè)或前7個(gè)主成分進(jìn)行回歸分析仍然得到典型式。但是，7號(hào)試驗(yàn)點(diǎn)采用前8個(gè)或前7個(gè)主成分回歸，肥效模型均由OLS的典型式變成了非典型式，結(jié)果反而變劣。

表6 三元二次多項(xiàng)式肥效模型的主成分回歸和典型性判別Table 6　Principal component regression and typicality discrimination of TQPMs

根據(jù)主成分回歸建立的典型肥效模型以及每千克 N4.30元、P2O55元、K2O6.0元和稻谷2.0元的市場(chǎng)均價(jià)，用邊際產(chǎn)量導(dǎo)數(shù)法計(jì)算推薦施肥量（表7）。由于選用的前7個(gè)或前8個(gè)主成分綜合反映了試驗(yàn)設(shè)計(jì)矩陣X的99.9%以上的變異方差，1號(hào)至4號(hào)試驗(yàn)點(diǎn)主成分回歸模型的最高施肥量和經(jīng)濟(jì)施肥量及其預(yù)測(cè)產(chǎn)量，其數(shù)值大小均在試驗(yàn)設(shè)計(jì)施肥量和試驗(yàn)所得產(chǎn)量范圍內(nèi)，且與已有的研究結(jié)果［17］基本一致。尤其是5號(hào)和6號(hào)試驗(yàn)點(diǎn)，主成分回歸的模型推薦施肥量與OLS回歸的模型推薦施肥量差異極小甚至相同，預(yù)測(cè)早稻產(chǎn)量也基本一致，表明主成分回歸的推薦施肥是可靠的。

表7　三元二次多項(xiàng)式肥效模型不同建模方法的推薦施肥量比較Table 7　Comparison in recommended application rate between TQPMs different in modeling method

2.4　主成分回歸對(duì)三元肥效模型的建模效果

表6的代表性例子表明，主成分回歸能使某些肥效試驗(yàn)結(jié)果由OLS法的未通過顯著性檢驗(yàn)或非典型式轉(zhuǎn)化為典型肥效模型，但也可能使OLS法的典型式變成非典型式。為了評(píng)價(jià)主成分回歸對(duì)三元肥效模型的建模效果，利用近年來福建省早稻測(cè)土配方施肥工作在主要項(xiàng)目縣完成的 “3414”設(shè)計(jì)的氮磷鉀田間肥效試驗(yàn)資料，逐個(gè)進(jìn)行回歸分析和比較，統(tǒng)計(jì)結(jié)果見表8。

在171個(gè)肥效試驗(yàn)中，OLS回歸建模的三元典型肥效模型的比例占試驗(yàn)點(diǎn)總數(shù)的27.5%，而采用前8個(gè)主成分或前7個(gè)主成分回歸建模，典型肥效模型的比例則分別提高至41.5%和43.3%，為OLS法的1.5倍和1.6倍。其中，主成分回歸大幅度降低了模型系數(shù)符號(hào)不合理的非典型式，從OLS法的14.6%下降至幾乎為零；模型無最高產(chǎn)量點(diǎn)的非典型式也從OLS法的36.8%分別下降至32.2%和21.1%。但是，主成分回歸并非利用設(shè)計(jì)矩陣X的100%方差信息，其建模結(jié)果使推薦施肥量外推的非典型式個(gè)數(shù)反而有所增加。

表8　最小二乘法回歸（OLS）和主成分回歸（PCR）對(duì)三元二次多項(xiàng)式肥效模型建模結(jié)果的影響Table 8　Effects of modeling using OLS and PCR on TQPMs

建模過程還表明，如同7號(hào)試驗(yàn)點(diǎn)，主成分回歸有時(shí)也會(huì)失效，使OLS法的典型式變?yōu)榉堑湫褪?，甚至造成模型未能通過顯著性檢驗(yàn)。因此，從實(shí)用角度出發(fā)，如果主成分回歸僅針對(duì)OLS法的非典型式和統(tǒng)計(jì)不顯著的肥效模型，即：OLS結(jié)合前8個(gè)主成分的主成分回歸或者前7個(gè)主成分的主成分回歸的建模策略，結(jié)果使典型肥效模型的比例分別可提高至48.5%和55.6%，明顯優(yōu)于單獨(dú)使用OLS或者主成分回歸。

3　討 論

3.1　非典型肥效模型的產(chǎn)生原因及其多重共線性危害

經(jīng)典線性回歸分析理論假設(shè)回歸模型的各個(gè)自變量之間無多重共線性［19-20］。二次多項(xiàng)式肥效模型在數(shù)學(xué)形式上是一種非線性模型，人們通常不重視這種模型是否違背了無多重共線性的假設(shè)。故這種模型中是否存在共線性、共線性將產(chǎn)生什么危害、如何處理共線性等方面問題，均鮮見相關(guān)研究報(bào)道［4］。本研究通過方陣XTX條件數(shù)和方差膨脹因子［20］等常用指標(biāo)進(jìn)行多重共線性診斷，表明這種模型存在嚴(yán)重的多重共線性問題。

王興仁［5-6］和章明清等［7］曾報(bào)道，二元二次多項(xiàng)式肥效模型的典型式比例僅占40.2%，三元肥效模型的典型式則低至23.6%。人們對(duì)多元肥效模型出現(xiàn)大量非典型式的現(xiàn)象至今仍然困惑不解，許多人由此直接歸因于田間試驗(yàn)失敗或者對(duì)肥效模型的應(yīng)用價(jià)值失去信心。研究表明，多重共線性給多項(xiàng)式肥效模型的OLS回歸建模和推薦施肥準(zhǔn)確性帶來嚴(yán)重危害，具體表現(xiàn)在：①在完全共線性下方陣XTX不可逆，正則方程無解，無法建立肥效模型。②在近似共線性下，雖然可以得到OLS參數(shù)估計(jì)量，但此時(shí)lXTXl ≈ 0，引起方陣XTX主對(duì)角線上的元素顯著增大，從而使參數(shù)估計(jì)值的方差膨脹，難以得到精確的參數(shù)估計(jì)值，例如3號(hào)和4號(hào)試驗(yàn)點(diǎn)出現(xiàn)了非典型肥效模型。③如果模型中兩個(gè)回歸變量具有高度線性相關(guān)性，那么它們中的一個(gè)變量可以由另一個(gè)變量近似表征，這時(shí)回歸變量前的參數(shù)并不反映各自與作物產(chǎn)量之間的結(jié)構(gòu)關(guān)系，僅能反映它們對(duì)產(chǎn)量的共同影響，結(jié)果使參數(shù)失去了應(yīng)有的肥效含義。同時(shí)，回歸系數(shù)可能變得對(duì)產(chǎn)量數(shù)據(jù)的微小變化異常敏感，甚至越過零值，表現(xiàn)出系數(shù)值的正負(fù)號(hào)出現(xiàn)反常現(xiàn)象，導(dǎo)致參數(shù)估計(jì)值的專業(yè)含義不合理，例如2號(hào)試驗(yàn)點(diǎn)的P項(xiàng)本應(yīng)為正數(shù)，結(jié)果卻為負(fù)數(shù)。④由于方差膨脹，模型顯著性檢驗(yàn)失效，例如1號(hào)試驗(yàn)點(diǎn)出現(xiàn)未通過顯著性檢驗(yàn)的假象。⑤膨脹的方差使產(chǎn)量預(yù)測(cè)區(qū)間過度變大，模型的預(yù)測(cè)功能失效。

因此，現(xiàn)在基本可以明確，非典型肥效模型大多是由多重共線性引起的，肥效模型設(shè)定本身是產(chǎn)生多重共線性的主要原因。

3.2　三元二次多項(xiàng)式肥效模型的主成分回歸

多項(xiàng)式統(tǒng)計(jì)模型在農(nóng)學(xué)學(xué)科中出現(xiàn)大量非典型式的現(xiàn)象，在其他學(xué)科諸如經(jīng)濟(jì)學(xué)、社會(huì)學(xué)、心理學(xué)等學(xué)科也同樣大量存在。主成分回歸經(jīng)過近幾十年的發(fā)展，已成為這些學(xué)科對(duì)付多重共線性的有效技術(shù)［24］。主成分回歸的基本思路是利用主成分分析提取X矩陣的主成分，然后選取方差貢獻(xiàn)率達(dá)到符合要求的幾個(gè)主成分，計(jì)算各處理對(duì)應(yīng)回歸變量的得分值，組成新的設(shè)計(jì)矩陣，最后利用新設(shè)計(jì)矩陣與各處理的試驗(yàn)產(chǎn)量進(jìn)行OLS回歸建模。本研究表明，這種建模方法能使1號(hào)至4號(hào)試驗(yàn)點(diǎn)由OLS法未達(dá)統(tǒng)計(jì)顯著水平或非典型式轉(zhuǎn)化為典型式。

在近年來的測(cè)土配方施肥工作中，福建省在早稻上完成了眾多“3414”設(shè)計(jì)的田間肥效試驗(yàn)。本研究收集整理了來自土肥技術(shù)力量較強(qiáng)的主要項(xiàng)目縣完成的171個(gè)試驗(yàn)資料。針對(duì)這些試驗(yàn)結(jié)果，選用前8個(gè)主成分或前7個(gè)主成分進(jìn)行回歸建模，三元典型式出現(xiàn)比例分別提高至41.5%和43.3%。尤其是前7個(gè)主成分回歸建模使模型系數(shù)符號(hào)不合理的非典型式出現(xiàn)比例下降至零，同時(shí)大幅度降低了無最高產(chǎn)量點(diǎn)的非典型式比例（表8），表明了主成分回歸技術(shù)在提高三元肥效模型建模成功率方面具有重要的應(yīng)用價(jià)值。

對(duì)比研究表明，因各施肥處理的氮磷鉀施肥量所采用的量綱相同，數(shù)據(jù)差異不顯著，試驗(yàn)設(shè)計(jì)矩陣的原始數(shù)據(jù)是否進(jìn)行了標(biāo)準(zhǔn)化處理，對(duì)主成分回歸結(jié)果基本無影響；若將全部9個(gè)主成分均納入主成分回歸，其結(jié)果與OLS法完全一致，達(dá)不到改善建模效果的作用；雖然前5個(gè)主成分就能反映X矩陣方差變異的99.87%，但采用6個(gè)或者少于6個(gè)主成分進(jìn)行回歸分析，所建立的171個(gè)早稻三元二次多項(xiàng)式肥效模型均未能通過顯著性檢驗(yàn)。因此，對(duì)三元肥效模型而言，主成分回歸只能選取前8個(gè)主成分或者前7個(gè)主成分。

在消除或緩解多重共線性危害方面，數(shù)理統(tǒng)計(jì)學(xué)家已經(jīng)提出了諸多解決方法［19-20，24］，諸如剔除不重要自變量、擴(kuò)大樣品容量，或者采用嶺回歸、主成分分析和偏最小二乘法等有偏估計(jì)方法來提高參數(shù)估計(jì)的穩(wěn)健性。顯然，剔除不重要自變量和擴(kuò)大樣品容量的方法不適合肥效模型的構(gòu)建。在有偏估計(jì)方法中，嶺回歸估計(jì)的質(zhì)量取決于偏倚系數(shù)的選取，但目前尚無標(biāo)準(zhǔn)的決策準(zhǔn)則［24］，僅能憑經(jīng)驗(yàn)判斷。偏最小二乘回歸是通過迭代算法在X矩陣中逐步提取主成分，使其盡可能多地反映X矩陣的變異方差，同時(shí)又兼顧對(duì)因變量（產(chǎn)量）的解釋能力。然而，其迭代算法復(fù)雜，不易使用。相比而言，主成分回歸具有扎實(shí)的統(tǒng)計(jì)學(xué)基礎(chǔ)，建模思路簡(jiǎn)潔明了。該法在目前已有MATLAB、SAS和Stata等諸多成熟軟件可供使用，是較簡(jiǎn)便和實(shí)用的肥效模型參數(shù)估計(jì)方法，易于被廣大應(yīng)用人員理解和掌握。

3.3　三元二次多項(xiàng)式肥效模型的建模策略

盡管主成分回歸具有許多優(yōu)點(diǎn)，但它畢竟是一種有偏估計(jì)。受回歸偏倚的影響，有時(shí)會(huì)將諸如7號(hào)試驗(yàn)點(diǎn)的OLS回歸典型式變成了非典型式。該法在提取主成分時(shí)僅考慮了試驗(yàn)設(shè)計(jì)矩陣，而不考慮各個(gè)主成分與各處理產(chǎn)量的相關(guān)性，而且，對(duì)什么樣的試驗(yàn)數(shù)據(jù)會(huì)造成主成分回歸失效？目前在統(tǒng)計(jì)學(xué)上并無明確的答案［24-25］。鑒于此，從實(shí)用角度出發(fā)，主成分回歸僅針對(duì)OLS法未通過顯著性檢驗(yàn)肥效模型和非典型式肥效模型，而對(duì)OLS回歸得到典型式的肥效模型就不再做主成分回歸。這個(gè)策略使三元典型式的比例進(jìn)一步提高至55.6%（表8），為單獨(dú)使用OLS法的2.0倍。

4　結(jié) 論

對(duì)三元肥效模型而言，主成分回歸是解決多項(xiàng)式肥效模型多重共線性危害的一種有效方法，OLS建模法結(jié)合前7個(gè)主成分的主成分回歸是三元二次多項(xiàng)式肥效模型的最佳建模策略。

［ 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 ］ Gregoret M C，Zorita M D，Dardanelli J，et al.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 ］ 王興仁，陳新平，張福鎖，等. 施肥模型在我國(guó)推薦施肥中的應(yīng)用. 植物營(yíng)養(yǎng)與肥料學(xué)報(bào)，1998，4（1）：67—74 Wang X R，Chen X P，Zhang F S，et al. Application of fertilization model for fertilizer recommendation in China（In Chinese）. Plant Nutrition and Fertilizer Science，1998，4（1）：67—74

［ 4 ］ 章明清，李娟，孔慶波，等. 作物肥料效應(yīng)函數(shù)模型研究進(jìn)展與展望. 土壤學(xué)報(bào)，2016，53（6）：1143—1156 Zhang M Q，Li J，Kong Q B， et al. Progress and prospect of the study on crop—response—to —fertilization function model（In Chinese）. Acta Pedologica Sinica，2016，53（6）：1143—1156

［ 5 ］ 王興仁. 二元二次肥料效應(yīng)曲面等產(chǎn)線圖在科學(xué)施肥中的位置（一）. 土壤通報(bào)，1985，16（1）：30—34 Wan X R. The positions of contour map of yield of binary quadratic fertilizer response curve in scientific fertilization（part 1）（In Chinese）. Chinese Journa1 of Soil Science，1985，16（1）：30—34

［ 6 ］ 王興仁. 二元二次肥料效應(yīng)曲面等產(chǎn)線圖在科學(xué)施肥中的位置（二）. 土壤通報(bào)，1985，16（2）：86—88 Wan X R. The positions of contour map of yield of binary quadratic fertilizer response curve in scientific fertilization（part 2）（In Chinese）. Chinese Journa1 of Soil Science，1985，16（2）：86—88

［ 7 ］ 章明清，林仁塤，林代炎，等. 極值判別分析在三元肥效模型推薦施肥中的作用. 福建農(nóng)業(yè)學(xué)報(bào)，1995，10（2）：54—59 Zhang M Q，Lin R X，Lin D Y，et al. Function of distinguish analysis on extreme value in recommendatory fertilization for three—fertilizer efficiency model（In Chinese）. Fujian Journal of Agricultural Sciences，1995，10（2）：54—59

［ 8 ］ Cerrato M E，B1ackmer A M. Comparision of models for describing corn yield response to nitrogen fertilizer.Agronomy Journal，1990，82（1）：138—143

［ 9 ］ 楊靖一，Wadsworth G A，Greenwood D J. 三種蔬菜N肥效應(yīng)曲線的比較研究. 植物營(yíng)養(yǎng)與肥料學(xué)報(bào)，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（In Chinese）. Plant Nutrition and Fertilizer Science，1995，1（1）：71—78

［10］ 陳新平，周金池，王興仁，等. 小麥—玉米輪作制中氮肥效應(yīng)模型的選擇——經(jīng)濟(jì)和環(huán)境效益分析. 土壤學(xué)報(bào)，2000，37（3）：346—354 Chen X P，Zhou J C，Wang X R，et al. Economic and environmental evaluation on models for describing crop yield response to nitrogen fertilizers at winter—wheat and summer—corn rotation system（In Chinese）. Acta Pedologica Sinica，2000，37（3）：346—354

［11］ 吳平，陶勤南. 氮磷鉀肥效模型的研究及其應(yīng)用. 浙江農(nóng)業(yè)大學(xué)學(xué)報(bào)，1989，15（4）：383—388 Wu P，Tao Q N. Study on N，P and K fertilizer response model and its application（In Chinese）.Journal of Zhejiang Agricultural University，1989，15（4）：383—388

［12］ 章明清. 米氏肥效方程的灰色建模法. 土壤通報(bào)，1997，28（1）：18—19 Z h a n g M Q. G r a y m o d e l i n g m e t h o d o f Mitsherlich fertilizer response equations（In Chinese）. Chinese Journa1 of Soil Science，1997，28（1）：18—19

［13］ 章明清，徐志平，姚寶全，等. Monte Carlo 法在多元肥效模型參數(shù)估計(jì)和推薦施肥中的應(yīng)用. 植物營(yíng)養(yǎng)與肥料學(xué)報(bào)，2009，15（2）：366 —373 Zhang M Q，Xu Z P，Yao B Q，et al. Using Monte Carlo method for parameter estimation and fertilization recommendation of multivariate fertilizer response model（In Chinese）. Plant Nutrition and Fertilizing Science，2009，15（2）：366—373

［14］ C o l w e l l J D. T h e d e r i v a t i o n o f f e r t i l i z e r recommendations for crops in non—uniform environment：Fertilizer. Crop Quality and Economy，1974：936—961

［15］ 王興仁，陳倫壽，毛達(dá)如，等. 分類回歸綜合法及其在區(qū)域施肥決策中的應(yīng)用. 土壤通報(bào)，1989，20（1）：17—21 Wang X R，Chen L S，Mao D R，et al. Classification of regression synthesis and its application for regional fertilization decision—making（In Chinese）. Chinese Journa1 of Soil Science，1989，20（1）：17—21

［16］ 毛達(dá)如，張承東. 多點(diǎn)肥料效應(yīng)函數(shù)的動(dòng)態(tài)聚類方法.北京農(nóng)業(yè)大學(xué)學(xué)報(bào)，1991，17（2）：49—54 Mao D R，Zhang C D. Cluster analysis of quadratic response function on the fertilizer dispersed experiment（In Chinese）. Acta Agriculturae Universitatis Pekinensis，1991，17（2）：49—54

［17］ 李娟，章明清，孔慶波，等. 福建早稻測(cè)土配方施肥指標(biāo)體系研究. 植物營(yíng)養(yǎng)與肥料學(xué)報(bào)，2010，16（4）：938—946 Li J，Zhang M Q，Kong Q B，et al. Soil testing and formula fertilization index for early rice in Fujian Province（In Chinese）. Plant Nutrition and Fertilizer Science，2010，16（4）：938— 946

［18］ 魯如坤．土壤農(nóng)業(yè)化學(xué)分析方法．北京：中國(guó)農(nóng)業(yè)科學(xué)技術(shù)出版社，2000：146—196 Lu R K. Analytical methods for soil and agro—chemistry（In Chinese）. Beijing：China Agricultural Science and Technology Press，2000：146—196

［19］ 何曉群，劉文卿. 應(yīng)用回歸分析. 第3版. 北京：中國(guó)人民大學(xué)出版社，2013：158—170 He X Q，Liu W Q. Applied regression analysis（In Chinese）. 3rd ed. Beijing：China Renmin University Press，2013：158—170

［20］ Montgomery D C，Peck E A，Vining G G. 王辰勇，譯. 線性回歸分析導(dǎo)論. 第5版. 北京：機(jī)械工業(yè)出版社，2016：91—268 Montgomery D C，Peck E A，Vining G G. Wang C Y. trans. Introduction to linear regression analysis（In Chinese）. 5th ed. Beijing：China Machine Press，2016：91—268

［21］ 袁志發(fā)，宋世德. 多元統(tǒng)計(jì)分析. 第2版. 北京：科學(xué)出版社，2009：219—250 Yuan Z F，Song S D. Multivariate statistical analysis（In Chinese）. 2nd ed. Beijing：Science Press，2009：219—250

［22］ 薛毅，陳立萍. 實(shí)用數(shù)據(jù)分析與MATLAB軟件. 北京：北京工業(yè)大學(xué)出版社，2015：189—274，377—394 Xue Y，Chen L P. Practical data analysis and MATLAB software（In Chinese）. Beijing：Beijing University of Technology Press，2015：189—274，377—394

［23］ Gujarati D N，Porter D C，費(fèi)劍平，譯. 計(jì)量經(jīng)濟(jì)學(xué)基礎(chǔ)（上冊(cè)）. 第5版. 北京：中國(guó)人民大學(xué)出版社，2011：313—399 Gujarati D N，Porter D C. Fei J P. trans. Econometrics foundation（the first volume）（In Chinese）. 5th ed. Beijing：China Renmin University Press，2011：313—399

［24］ 王惠文，王劫，黃海軍. 主成分回歸的建模策略研究.北京航空航天大學(xué)學(xué)報(bào)，2008，34（6）：661—664 Wang H W，Wang J，Huang H J. Modeling strategy of principle component regression（In Chinese）.Journal of Beijing University of Aeronautics and Astronautics，2008，34（6）：661—664

［25］ 葉宗裕. 主成分回歸能消除多重共線性嗎? 統(tǒng)計(jì)與信息論壇，2012，27（3）：15—20 Ye Z Y. Can the principal components regression eliminate the malticollinearity（In Chinese）?Statistics & Information Forum，2012，27（3）：15—20