基于混合回歸分析的荷斯坦奶牛氮排泄量預測

董瑞蘭 孫國強 于光輝

(青島農業大學動物科技學院，青島 266109)

畜牧養殖對環境造成的污染問題已引起全球的廣泛關注[1]。奶牛氮排泄量與飼糧的粗蛋白質(CP)含量和消化率密切相關[2]。在奶業生產中，為了提高產奶量和乳品質而提高飼糧CP的含量，導致70%～80%的攝入氮通過糞便和尿液排泄到水體、陸地和大氣中[3-4]，進一步造成環境的惡化。我國農業農村部開展的全國農業污染源普查中，畜禽養殖污染物的排放是基于動物的產排污系數核算的[5]。這種方法操作起來方便，但是糞污養分含量隨飼料的來源和組成不同而發生變化，用固定的系數估計變化的排放數值是不準確的。因此，通過影響氮排泄量的指標開發預測模型是準確評估養分排放的有效方法[6]，有助于提高糞污養分預測的準確性。而且，鑒于尿液中含有的氮比糞中的氮更易揮發[7]，通過排泄路徑(尿液和糞便)預測氮排放量更有意義。

關于奶牛的氮排泄量已有較多預測方程報道，有根據乳中尿素氮(milk urea nitrogen，MUN)[8-9]、氮攝入量(nitrogen intake，NI)[10-12]、瘤胃降解蛋白(rumen degradable protein，RDP)[13]預測尿氮(urinary nitrogen，UN)的一元線性方程；也有根據NI、干物質采食量(dry matter intake，DMI)和產奶量(milk yield，MY)組合[14]，或MUN、CP和DMI組合[15]預測UN的多元線性方程。有根據NI[10-11]或DMI[12]預測糞氮(fecal nitrogen，FN)的一元線性方程；也有根據NI和DMI[14]或有機物采食量(organic matter intake，OMI)組合[16-17]預測FN的二元方程。總氮(total nitrogen，TN)可根據MY[13]或NI[11-12,18]單獨預測或根據NI和DMI組合[14]，體重(body weight，BW)、CP與MY組合[18]以及BW、CP與DMI組合[19]預測。Dong等[20]用荷斯坦奶牛數據庫評估了上述預測模型，發現預測的準確性普遍偏低。因此，本研究旨在用荷斯坦奶牛的氮排泄數據庫，構建準確的UN和FN預測模型，并驗證準確性，為奶牛生產提供依據。

1 材料與方法

1.1 建立數據庫

荷斯坦奶牛氮排泄數據庫(包括32個發表研究的99個觀察值)來自Dong等[20]。采用文獻發表的荷斯坦奶牛氮排泄數據構建數據庫。在線檢索的數據庫包括中國知網(CNKI，http://www.cnki.net/)、Science Direct(http://www.sciencedirect.com/)和Journal of Dairy Science(http://www.journalofdairyscience.org/)。確立相關文獻研究的標準是：1)研究對象是荷斯坦奶牛；2)從已出版文獻的數據庫檢索出在中國進行的研究，代表現代奶牛養殖的實際條件，特點是使用全混合日糧(TMR)飼養和采用自動擠奶設備；3)需提供荷斯坦奶牛氮排泄的輸出變量信息，包括UN、FN、MUN、TN和尿氮占總氮的比例(UN/TN)；4)需提供關于飼糧的CP含量、奶牛的NI、DMI、BW和MY及乳成分等特征輸入變量信息；5)若研究中某些變量未報道，但可通過其他已知變量計算的也包括在數據庫中。例如，飼糧的CP含量可通過DMI和NI計算獲得，TN為UN與FN的總和，UN/TN為UN與TN的比值。剔除觀察值變量缺失的5個研究[21-25]，剩余27個研究的89個觀察值對應的所有變量均已知，即如果文獻的觀察值提供了關于NI、DMI、CP、BW、粗飼料比例(forage proportion，FP)、MY、乳氮(milk nitrogen，MN)、MUN、氮的全消化道表觀消化率(apparent total tract nitrogen digestibility，TTND)、乳氮占氮攝入量的比例(MN/NI)、氮攝入量占干物質采食量的比例(NI/DMI)、飼料利用效率(feed efficiency，FE，MY與DMI的比值)、氮沉積量(nitrogen retention，NR)、UN、FN、TN、UN/TN等所有變量信息，則被包含在數據庫中。如果有觀察值的某個變量未知則不被包含在數據庫中。

數據庫經JMP統計軟件馬氏離群分析(Mahalanobis outlier analysis，MOA)認定為極值的數據，其學生化殘差在-2.5～2.5的范圍之外，被鑒定為離群異常值。通過分析，共計8個觀察值被鑒定為離群異常值，剔除離群異常值后，最后剩余27個研究的81個觀察值(隨機選擇2/3數據用于建模，1/3數據用于分析模型預測值與實際觀察值之間的回歸關系)進行最終的分析。

1.2 構建預測模型

1.2.1 篩選最佳預測變量

模型的最佳預測變量采用SAS 8.1軟件Corr過程程序的皮爾森相關系數(r)確定[1]。依次分析輸入變量(NI、DMI、CP、MUN、BW、TTND、NI/DMI、MY、MN、MN/NI、FP、NR、FE)與輸出變量(UN、FN、TN、UN/TN)之間的相關關系。

1.2.2 自變量和因變量之間的回歸關系

經1.2.1中篩選出的最佳預測變量作為自變量，將奶牛的氮排泄量(UN、FN、TN、UN/TN)作為因變量，分析自變量和因變量之間的一元和多元回歸關系。

進行多元回歸關系統計分析時，采用逐步回歸分析法(step wise selection，SWS)剔除不顯著(P>0.05)的變量，篩選出最終保留在模型中的自變量。多個變量之間的共線性用方差膨脹因子(variance inflation factors，VIF)檢驗，它是衡量每個自變量與模型中其他變量之間緊密相關的指標。預測變量若具有較大的VIF(>2.5)，可導致模型發生較大的膨脹誤差。當確定存在共線性時，預測氮排放量時皮爾森相關系數最小的變量即被剔除。

1.2.2.1 一元線性模型

自變量和因變量之間的一元回歸關系用以下方程分析：

y=bx+a。

式中：x指因變量；y指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN；a是固定常數；b指對應變量的系數。

1.2.2.2 多元線性模型

自變量和因變量之間的多元回歸關系用以下方程分析：

y=b1(x1)+b2(x2)+b3(x3)+…+bn(xn)+a。

式中：x指因變量；y指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN；a是固定常數；b1，b2，b3，…，bn指各對應變量的系數。

1.2.3 模型的擬合效果

用SAS 8.1軟件分析均方根誤差(root mean square error，RMSE)和阿凱克信息論準則(Akaike’s information criterion，AIC)值。RMSE和AIC值是模型擬合好壞的指標，其值越小代表模型擬合效果越好[1]。

1.3 模型預測值與實際觀察值之間的回歸關系分析

將模型輸出的預測值作為自變量，實際觀察值作為因變量，繪圖分析模型預測值與實際觀察值之間線性回歸關系的顯著性，獲得決定系數(R2)，同時確定回歸線和等線(y=x)的偏離程度。

1.4 數據的統計與分析

用SAS 8.1軟件的混合回歸模型(mixed regression model，MRM)過程程序分析自變量和因變量之間的一元和多元回歸關系，P<0.05為回歸關系顯著，P<0.01為回歸關系極顯著。通常認為，模型預測值與實際觀察值之間R2越大，則獲得的預測模型越準確。在用SAS 8.1軟件程序分析時，不同的研究作為隨機變量也包括在模型中，用于消除不同研究之間的試驗誤差。模型預測值和實際觀察值之間的線性關系用以下方程分析：

y=bx+a。

式中：x指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN的模型預測值；y指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN的實際觀察值；a是固定常數；b指對應變量的系數。

2 結果與分析

2.1 建模數據庫

由表1可知，奶牛的BW、CP、DMI、TTND、NI平均值分別為565 kg、15.4%、19.20 kg/d、69.8%和470 g/d，變動范圍分別為407～755 kg、11.0%～19.2%、11.57～28.10 kg/d、58.9%～78.0%和214～720 g/d。奶牛的UN和FN的變動范圍較大，介于58.2～329.0 g/d、81.0～193.0 g/d，平均值分別為168.0和141.0 g/d。奶牛的TN變動范圍為139.2～481.3 g/d。UN/TN平均值為0.54，變動范圍為0.39～0.69。

表1 基于動物和飼糧因子以及氮輸入和輸出變量的建模數據庫描述性統計

2.2 奶牛氮排泄量與動物和飼糧因子之間的皮爾森相關系數

由表2可知，UN分別與NI、DMI、CP、MUN、TTND、NI/DMI、MY和MN呈顯著或極顯著正相關關系(P<0.05或P<0.01)。FN分別與NI、DMI、CP、NI/DMI、MY、MN和FE呈極顯著正相關關系(P<0.01)，而與TTND和FP呈顯著或極顯著負相關關系(P<0.05或P<0.01)。TN分別與NI、DMI、CP、MUN、NI/DMI、MY和MN呈極顯著正相關關系(P<0.01)。UN/TN分別與DMI、MUN、TTND和FP呈顯著或極顯著正相關關系(P<0.05或P<0.01)，而與NR呈極顯著負相關關系(P<0.01)。

表2 奶牛氮排泄量與動物和飼糧因子之間的皮爾森相關系數

2.3 一元線性預測模型

2.3.1 UN預測模型

預測UN的模型見表3(模型1～6)。輸入變量NI、CP、NI/DMI、MY、MN/NI、MN對回歸關系都有顯著影響(P<0.05)。根據RMSE和AIC值可知，基于單一變量的UN模型擬合效果為：MN/NI

表3 奶牛尿氮、糞氮和總氮的一元線性預測模型

2.3.2 FN預測模型

預測FN的模型見表3(模型7～12)。輸入變量NI、MY、MN、DMI、NI/DMI、CP對回歸關系都有顯著影響(P<0.05)。當單獨用NI(模型7)預測FN時，其RMSE和AIC值最低。但單獨用CP(模型12)預測FN時，其RMSE和AIC值最高。用NI預測的FN與實際觀察值之間的比較見圖1-b，回歸分析顯示，基于NI的回歸模型預測值和實際觀察值之間的R2為0.84(P<0.01)。

2.3.3 TN預測模型

預測TN的模型見表3(模型13～18)。輸入變量NI、MY、MN、NI/DMI、CP、FP對回歸關系都有顯著影響(P<0.05)。當單獨用NI(模型13)預測TN時，其RMSE和AIC值最低。但單獨用FP(模型18)預測FN時，其RMSE和AIC值最高。用NI預測的TN與實際觀察值之間的比較見圖1-c，回歸分析顯示，基于NI的回歸模型預測值和實際觀察值之間的R2為0.93(P<0.01)。

圖1 尿氮(a)、糞氮(b)、總氮(c)和尿氮占總氮的比例(d)的實際觀察值和模型預測值之間的關系

2.3.4 UN/TN預測模型

UN/TN預測模型見表3(模型19～22)。輸入變量TTND、MN/NI、CP、FE對回歸關系都有顯著影響(P<0.05)。當單獨用TTND(模型19)預測UN/TN時，其RMSE和AIC值最低。單獨用TTND預測的UN/TN與實際觀察值之間的比較見圖1-d，回歸分析顯示，基于TTND的回歸模型預測值和實際觀察值之間的R2為0.90(P<0.01)。

2.4 多元線性模型

2.4.1 UN預測模型

預測UN的多元線性回歸模型見表4(模型23和24)。在多個輸入變量中，CP、DMI與NI/DMI組合以及DMI與NI/DMI組合對UN多元線性模型有極顯著影響(P<0.01)。模型23的RMSE和AIC值與模型24相比較低。用模型23預測的UN與實際觀察值之間的比較見圖2-a，回歸分析表明，模型23的預測值和實際觀察值之間的R2為0.90(P<0.01)。

表4 奶牛尿氮、糞氮和總氮的多元線性預測模型

2.4.2 FN預測模型

預測FN的多元線性回歸模型見表4(模型25和26)，在多個輸入變量中，MY與NI組合以及MY與DMI、FP組合對FN多元線性模型有極顯著影響(P<0.01)。模型25的RMSE和AIC值較模型26低。用模型25預測的FN與其實際觀察值之間的比較見圖2-b，回歸分析表明，模型25的預測值和實際觀察值之間的R2為0.87(P<0.01)。

圖2 尿氮(a)、糞氮(b)、總氮(c)和尿氮占總氮的比例(d)的實際觀察值和模型預測值之間的關系

2.4.3 TN預測模型

預測TN的多元線性回歸模型見表4(模型27～29)，在多個輸入變量中，NI與NI/DMI組合、DMI與NI/DMI組合以及DMI與MY組合對TN多元線性模型有極顯著影響(P<0.01)。模型27的AIC和RMSE值是3個TN預測模型中最低的。模型29的AIC和RMSE值是3個TN預測模型中最高的。用模型27預測的TN與其實際觀察值之間的比較見圖2-c，回歸分析顯示，模型27的預測值和實際觀察值之間的R2為0.92(P<0.01)。

2.4.4 UN/TN預測模型

預測UN/TN的多元線性回歸模型見表4(模型30和31)。NI、TTND、FP與NR組合以及NI、TTND與NR組合對UN/TN多元線性模型有極顯著影響(P<0.01)。模型31的RMSE和AIC值最低。用模型31預測的UN/TN與實際觀察值之間的比較見圖2-d，回歸分析顯示，模型31的預測值和實際觀察值之間的R2為0.89(P<0.01)。

3 討 論

3.1 建模數據庫

數據庫涵蓋了奶牛泌乳初期、中期和后期各產乳階段。但需指出的是，由于研究使用了發表文獻數據，收集的數據是各試驗處理的觀察值平均值，而不是動物個體的試驗數據，所以本研究未考慮動物之間的個體差異。與美國1990—1995年泌乳奶牛數據庫[12]DMI范圍10.3～28.2 kg/d相比，荷斯坦奶牛的最低和最高DMI數值與之接近。將荷斯坦奶牛的氮排泄數據庫與Spek等[15]建立的北歐和北美數據庫相比，中國NI的平均值低于北歐(470 g/d vs. 485 g/d)和北美(470 g/d vs. 637 g/d)。同樣，UN的平均值也低于北歐(168 g/d vs. 185 g/d)和北美(168 g/d vs. 212 g/d)；FN的平均值也低于北歐(141 g/d vs. 159 g/d)和北美(141 g/d vs. 223 g/d)。

3.2 UN預測模型

用單一變量預測UN時，基于NI(模型1)預測的準確性(AIC值=510.2，R2=0.91)最高。關于NI和UN之間存在顯著的相關關系已有較多報道[10-12]。Castillo等[10]和Kebreab等[11]研究發現，NI和UN之間呈指數關系，原因可能為在氮平衡為正的情況下，UN不是直接測定，而是用物料守恒法計算差值獲得，由此導致了UN的偏高估計。荷斯坦奶牛排泄的UN與NI之間呈現正相關線性關系，而且用NI可準確預測奶牛的UN，該結果與Kebreab等[11]和Reed等[12]的結論一致。本研究中UN和NI之間的回歸關系斜率與Reed等[12](0.34 vs. 0.33)和Kebreab等[11](0.34 vs. 0.38)報道的范圍接近，而且本研究的回歸關系常數項值更低(4.8 vs. 12.0和4.8 vs. 20.0)。

當用多個變量預測UN時，模型23提供了較好的預測方法(R2=0.90)。Spek等[15]用DMI和CP對歐洲和北美奶牛的UN進行預測，發現在輸入變量DMI和CP的基礎上增加第3個變量MUN，比單獨用MUN或CP以及同時用MUN和CP的組合預測得更準確。本研究中，模型23的回歸關系常數項與Spek等[15]報道的一致，而且AIC值更低(488.3 vs. 988.0)。

3.3 FN預測模型

用單一變量NI、MY、MN、DMI、NI/DMI、CP預測FN時，模型7的擬合效果和準確性最高。Castillo等[10]基于英國的奶牛數據發現，FN與NI存在線性函數關系(FN=0.21NI+52.3)。此方程的斜率與本研究的模型7一致，但是Castillo等[10]建立的模型的常數項值約為本研究中模型7的1.2倍。Vasconcelos等[26]也報道了奶牛FN與NI之間存在線性關系。Kebreab等[11]報道了奶牛FN與NI之間也存在類似的線性關系，但FN的斜率為0.28，截距數值為10。Schuba等[27]研究發現NI和FN之間呈二次關系，即FN=117.3+3.24×NI2-0.004 3×NI，但其AIC值高達1 509，說明FN和NI之間二次關系擬合效果低。

用多個變量預測荷斯坦奶牛排泄的FN時，基于NI和MY的FN預測模型(模型25)擬合性略高，預測效果較好(R2=0.87)。將多元線性模型(模型25)與一元線性模型(模型7)相比，發現多元線性模型的預測值和實際觀察值之間的R2較高，而常數項值也低于一元線性模型。

3.4 TN預測模型

將NI、CP、MY、MN、FP與NI/DMI作為單一輸入變量時，基于NI(模型13)預測TN的準確性最高(AIC值=505.6，R2=0.93)。用多個變量預測TN時，模型27預測的準確性較高。將基于NI的TN預測模型(模型13)分別與基于NI的UN預測模型(模型1)和基于NI的FN預測模型(模型7)相比，模型13的擬合效果優于模型1和模型7。

3.5 UN/TN預測模型

考慮到UN在TN中占較大比例，建立了UN/TN預測模型。與MN/NI(AIC值=-157.2，模型20)、CP(AIC值=-139.8，模型21)或FE(AIC值=-148.3，模型22)相比，將TTND(AIC值=-158.0，模型19)作為單一輸入變量的模型預測效果最好。通過UN/TN與動物或飼糧因子的多變量回歸關系及逐步回歸分析發現，變量TTND與UN/TN的相關系數最高。因此，TTND被證明是預測奶牛UN/TN的最佳預測變量，該發現與在肉牛中[1]的研究結果一致。Huhtanen等[14]將CP作為奶牛UN/TN的預測變量，發現二者存在顯著的線性關系。但本研究結果發現用TTND預測UN/TN比CP效果更好，原因是TTND與CP密切相關。飼糧CP和氮的消化率越高，就會導致越多的氮以瘤胃氨或小腸氨基酸的形式被吸收。那么，超過奶牛營養需要的剩余氮量就會進入尿液，從而增加了UN的比例。需要注意的是，本研究僅基于發表文獻進行建模分析，最好通過試驗獲得數據進一步驗證其預測性。

4 結 論

① 基于NI的一元UN、FN和TN回歸模型以及基于TNND的一元UN/TN回歸模型的RMSE和AIC值最低，一元回歸模型擬合效果好。

② 基于DMI、CP與NI/DMI組合的多元UN回歸模型，基于NI與MY組合的多元FN回歸模型，基于NI與NI/DMI組合的多元TN回歸模型的擬合效果好。