翟靜靜,劉則彰,蔣思文
(華中農(nóng)業(yè)大學,湖北 武漢 430070)
母豬在養(yǎng)殖生產(chǎn)中扮演著重要的角色,在規(guī)?;B(yǎng)豬生產(chǎn)中,母豬的繁殖效率是重中之重,是養(yǎng)殖場的核心,也是決定豬場生產(chǎn)水平和經(jīng)濟效益的關鍵因素[1]。母豬繁殖性能主要指標包括總產(chǎn)仔數(shù)、產(chǎn)活仔數(shù)、斷奶仔豬數(shù)、斷奶發(fā)情間隔(Weaning to estrus interval,WEI)、PSY等[2],還包括繁殖力(如分娩率和死亡率)。
胎次也是影響母豬繁殖性能的重要因素,一般來說,繁殖性能隨著胎次數(shù)的增加而提高,在3~5胎次時達到最高水平[3,4]。在2~4胎時母豬的受胎率會隨胎次升高而增加,4胎后逐漸下降[5]。胎次主要通過影響分娩率和產(chǎn)仔數(shù)來影響母豬繁殖性能[6]。朱世平等發(fā)現(xiàn),依據(jù)胎次可將母豬繁殖性能劃分為1-2胎、3-6胎和7-8胎這三個階段,并且得出母豬繁殖的最佳胎次為3-6胎,繁殖高峰期是4胎次(P<0.05)[7]。
WEI是指母豬斷奶到再次發(fā)情配種的間隔時間。據(jù)國外的各個豬場數(shù)據(jù)分析,WEI與產(chǎn)仔數(shù)呈負相關,WEI時間越長,產(chǎn)仔數(shù)越低,而在3~5 d發(fā)情的母豬,分娩率高。因此要想提高斷奶母豬7 d內(nèi)的配種率,應縮短WEI,保持母豬的理想體況,保障產(chǎn)仔舍和配懷舍的飼養(yǎng)管理工作,減少體重損失,做好誘情工作,提高分娩成績[7]。
生產(chǎn)中把繁殖母豬發(fā)情期進行配種后沒有懷孕并重新發(fā)情的現(xiàn)象稱為返情。常規(guī)返情是受胎失敗,非常規(guī)返情是妊娠失敗,都會導致低分娩率的出現(xiàn)。母豬一般配種受胎為80%~90%,而配種未分娩的10%~20%就是返情。返情的增加,會導致配種分娩率降低,從而影響?zhàn)B殖場的經(jīng)濟效益[8,9]。返情后再配種母豬分娩率降低10%左右[10]。實際生產(chǎn)中母豬返情率一般在10%左右[11]。
時間序列法是一種統(tǒng)計分析方法,將研究對象的數(shù)量隨著時間的發(fā)生而呈現(xiàn)變動的這種特征按照時間順序進行排列,所形成的這個序列便稱為研究對象的時間序列[12]。常用的時間序列法有移動平均法、加權移動平均法和指數(shù)平均法,時間序列的平穩(wěn)性是建模必不可少的環(huán)節(jié)[13]。時間序列分為兩類:第一類是某時間序列各時期的數(shù)值是不變動的,為確定型時間序列。第二類是此時間序列各時期數(shù)值變動,充滿不確定性,為隨機型時間序列。一般情況下都是隨機型的時間序列。朱迪國等對全球豬流行性腹瀉疫情進行時間序列分析發(fā)現(xiàn),豬血液制品貿(mào)易可能是該病跨國界傳播的潛在途徑[14]。因此,在母豬繁殖生產(chǎn)數(shù)據(jù)中,可以運用時間序列對性狀的時間發(fā)展進行預測。
對于因變量為只具有兩種分類結果的屬性數(shù)據(jù),例如母豬生產(chǎn)中的母豬分娩或未分娩,就需要去建立邏輯回歸模型,這是一種經(jīng)典的分類算法[15]。Logistic模型中的參數(shù)是通過極大似然法進行估計的。邏輯回歸分析是通過一組預測變量(自變量)預測一個或多個響應變量(因變量)一種統(tǒng)計分析方法[16]。在使用Logistic模型前,需判斷是否滿足以下幾個條件:(1)因變量是二分類變量;(2)分類變量的分類是全面的且相互獨立;(3) 連續(xù)的自變量與因變量的對數(shù)轉(zhuǎn)換值之間存在線性關系;(4)自變量之間不存在共線性。利用建立的回歸模型可以對母豬性能及后期發(fā)展進行預測。Logistic模型也可以用來評估預測變量對響應變量的預期效果。Hosmer and Lemeshow運用Logistic模型預測母豬體重發(fā)育相關變量與第2胎中的非妊娠和產(chǎn)仔數(shù)的關聯(lián),發(fā)現(xiàn)是高度相關[17]。Koketsu等人利用Logistic回歸模型分析與分娩率相關的因素時發(fā)現(xiàn),胎次、分娩季節(jié)以及泌乳天數(shù)與采食量之間的交互作用會對母豬分娩率產(chǎn)生影響[18]。
研究根據(jù)豬場的生產(chǎn)記錄,時間序列和Logistic回歸模型這兩種方法,把所得的豬場數(shù)據(jù)進行匯總處理,運用時間序列去進行短期預測母豬繁殖性能高低的時間段和母豬狀態(tài),來判斷豬場中繁殖母豬在某一時間段可能會出現(xiàn)的情況;通過Logistic回歸模型進一步剖析影響因素,對影響該豬場母豬繁殖性能的因素進行調(diào)查研究,為該豬場采取相應措施提高母豬繁殖性能提供科學依據(jù),有效提高豬場經(jīng)濟效益。
此研究選擇某公司的2016-2019年7個分場的母豬繁殖性能記錄數(shù)據(jù),包括產(chǎn)仔性能、配種分娩明細及淘汰死亡等。以母豬繁殖性能為主線,對規(guī)模豬場不同生產(chǎn)指標的各種數(shù)據(jù)進行處理分析,根據(jù)對生產(chǎn)數(shù)據(jù)的分析,找到影響繁殖性能的因素,制定行之有效的解決方案,以此展開相關研究。
對數(shù)據(jù)采用Excel軟件進行統(tǒng)計,數(shù)據(jù)以平均值±標準差表示;運用GraphPad Prism 7.0作圖,使用SPSS 19.0建立邏輯回歸模型和對數(shù)據(jù)進行顯著性檢驗,P<0.05表示差異顯著,P<0.01表示差異極顯著;運用SAS 9.4軟件構建時間序列模型ARIMA進行預測。
此研究采用二元Logistic回歸分析方法對母豬胎次和年份進行了分析。母豬繁殖性能相關因素分析研究的變量與賦值見表1。根據(jù)收集到的數(shù)據(jù),將數(shù)據(jù)分為兩類情況,如受胎:1和未受胎:0;首先采用多重共線性模型篩選變量(|r|>0.7)。然后,對每個解釋變量進行單因素分析,將上述變量P<0.1的解釋變量納入logistic回歸分析模型中,運用基于最大似然估計的向前逐步回歸法(Forward: LR法)建 立logistic回歸模型,確定與母豬繁殖性能相關風險因素的依據(jù)是似然比統(tǒng)計量的P值<0.05?;貧w系數(shù)用95%置信區(qū)間(Conf idence intervals,CI)的似然比(Odds ratios,OR)表示。OR值即是相對危險度的精確估計值,如果OR大于1,變量施加的風險增加。Logistic回歸模型的公式如下:
其 中,β0為 常 量,A、B分別表示因素斷奶或配種年份、母豬斷奶或配種胎次。把胎次、年份等影響7 d斷配率、返情率、受胎率、分娩率這些指標的因素進行定義變量和賦值(見表1)。
表1 母豬繁殖性能相關因素分析研究的變量與賦值
2.1.1 時序圖檢驗
此研究收集了2016-2019.11、2016-2019.08、2016-2019.07、2016-2019.07的豬場健仔數(shù)、斷奶頭數(shù)、7 d斷配率、分娩率的數(shù)據(jù)。選取了其中80%的數(shù)據(jù)用于時間序列分析,留下20%的數(shù)據(jù)作為實際值與預測值對比。首先對原始序列進行平穩(wěn)性檢驗,分別描繪平均每頭母豬健仔數(shù)(圖1a)、斷奶頭數(shù)(圖1b)、7 d斷配率(圖1c)、分娩率(圖1d)的時序圖,從時序圖中可以看出該時間序列具有明顯的遞增趨勢性和下降趨勢性,為典型的非平穩(wěn)性序列。
2.1.2 一階差分時序圖
時間序列預處理,可知這四個指標均為不平穩(wěn)序列,ARIMA模型基于平穩(wěn)時間序列的基礎上,因此根據(jù)數(shù)據(jù)特征對數(shù)據(jù)進行差分處理,一階差分后的時序圖(圖2)顯示,各觀測值圍繞著一個常數(shù)附近上下波動,數(shù)據(jù)逐漸趨于平穩(wěn),序列表現(xiàn)平穩(wěn)的性質(zhì)。
2.1.3 模型識別
為了驗證一階差分后此序列的平穩(wěn)性,依次做出該序列的自相關圖和偏相關圖(圖3到圖6)。觀察ACF(自相關)圖和拖尾效應確定ARIMA模型的P值,觀察PACF(偏相關)和截尾效應確定ARIMA模型的q值。其中差分的次數(shù)就是模型ARIMA(p, d,q)中的階數(shù)d,因此通過最佳準則函數(shù)定價方法選擇ARIMA(0,1,1)、ARIMA(5,1,5)、ARIMA(0,1,0)、ARIMA(0,1,0)分別為健仔數(shù)、斷奶頭數(shù)、7 d斷配率、分娩率的相對最優(yōu)模型。
2.1.4 模型的預測
根據(jù)已經(jīng)建立的ARIMA模型對母豬健仔數(shù)、斷奶頭數(shù)、7 d斷配率及分娩率的數(shù)據(jù)進行預測。預測時間為9個月,健仔數(shù)的預測時間從2019年3月到11月,斷奶頭數(shù)預測時間為2018年12月到2019年8月,7 d斷配率和分娩率預測時間都是從2018年11月到2019年7月。為評估模型的預測效果,引入均方誤差(Mean square error,MSE)和平均絕對百分比誤差(Mean absolute percentage error,MAPE)兩個指標。預測的偏差較小,與實際觀察值基本吻合(表2到表5)。預測值的變化趨勢與實際值大致保持一致,實際值基本落在預測曲線95%的可信區(qū)間內(nèi),該模型擬合預測成功,具有一定的實際價值。由表6誤差評估表可見,MSE值約小,說明預測模型具有更好的精確度,MAPE<20%,可以從中判斷該模型擬合預測很成功。
2.2.1 受胎率的logistic回歸分析
原始數(shù)據(jù)為192 874條配種記錄,剔除掉胎次為0的數(shù)據(jù),剩下143 464條記錄用于后續(xù)的因素分析。由表7可知,配種年份和配種胎次會影響母豬受胎率。高胎次(≥7)(OR值:0.889,95%CI:0.810~0.975)的母豬跟低胎次(1~2)的母豬相比受胎率小于0.889倍;在2019年(OR值:1.247,95%CI:1.157~1.344)配種的母豬比2016年配種的母豬成功受胎的概率高1.247倍。
2.2.2 返情率的logistic回歸分析
收集了7個分區(qū)豬場192 874條配種記錄,剔除掉胎次為0的數(shù)據(jù),剩下143 464條記錄用于后續(xù)的因素分析。由表8可知,配種年份和配種胎次會影響母豬返情率。3~6胎 次(OR值:1.272,95%CI:1.208~1.339)的母豬比1~2胎次的母豬出現(xiàn)返情現(xiàn)象的概率高1.272倍;在2019年(OR值:0.802,95%CI:0.744~0.864)配種的母豬比2016年配種的母豬發(fā)生返情的概率小。
2.2.3 分娩率的logistic回歸分析
表2 ARIMA模型預測母豬健仔數(shù)效果表 頭
表3 ARIMA模型預測母豬斷奶頭數(shù)效果表 頭
表4 ARIMA模型預測母豬7 d斷配率效果表 %
原始數(shù)據(jù)為192 874條配種記錄,剔除掉胎次為0的數(shù)據(jù),剩下143 464條記錄用于后續(xù)的因素分析。由表9可知,配種年份和配種胎次會影響母豬分娩率。2017年(OR值:0.782,95%CI:0.747~0.817)的母豬比2016年的母豬分娩率低;高胎次(≥7)(OR值:0.851,95%CI:0.808~0.897)的母豬成功分娩的概率比低胎次(1~2)的母豬小85%。
2.2.4 7 d斷配率的logistic回歸分析
收集了7個分區(qū)豬場2016-2019年間142 478條斷奶配種記錄,剔除掉WEI>100的數(shù)據(jù),剩下138 114條記錄用于后續(xù)的因素分析。由表 10可知,斷奶年份和斷奶胎次會影響母豬7 d斷配率。2019年(OR值:1.610,95%CI:1.526~1.698)的母豬斷奶發(fā)情間隔大于7 d的發(fā)生概率是2016年的1.610倍;3~6胎次(OR值:0.939,95%CI:0.904~0.975)的母豬斷奶發(fā)情間隔大于7 d的發(fā)生概率比低胎次(1~2)的母豬低。
數(shù)據(jù)的時間序列大多是不平穩(wěn)的,數(shù)據(jù)本身具有季節(jié)趨勢,若采用簡單的模型進行分析,往往不能得到準確的預測結果。隨著對時間研究的發(fā)展,時間序列研究的方法更加多樣,其中ARIMA模型是常用的非平穩(wěn)的模型[13]。時間序列預測法對于中短期預測的效果要優(yōu)于長期預測的效果,因為時間過長發(fā)生外界因素變化的可能性會加大。世界上沒有一個模型可以做到100%的精確,只能隨著時間發(fā)展進行改善,不斷添加新的數(shù)值和因素。此研究中以健仔數(shù)、斷奶頭數(shù)、7 d斷配率和分娩率的實際數(shù)據(jù)為基礎,擬合ARIMA模型,同時對未來9個月的數(shù)據(jù)進行了預測和比較。研究結果表明健仔數(shù)在未來幾個月呈上升趨勢,與實際值相比最大偏差為0.8頭;斷奶頭數(shù)預測值具有一定的波動性,與實際值最大偏差為1.67;7 d斷配率預測值持續(xù)上升,實際值時而下降時而上升,最大偏差為11.55個百分點;分娩率預測結果呈下降趨勢,實際值具有波動,最大偏差為4.26個百分點。這幾個指標雖有一定偏差,但實際值基本落在預測曲線95%的置信區(qū)間內(nèi),說明模型對原始序列信息的提取比較充分。通過預測效果表和MAPE值(值越小,代表預測的準確度越強)來看,利用ARIMA模型對這些數(shù)據(jù)進行擬合和預測是合理可行的,但也會有一定的局限性,生產(chǎn)過程中可能會因環(huán)境變化或疾病、飼養(yǎng)管理方案改變引起差異。在短期預測當中,豬場未來一段時間的母豬生產(chǎn)情況將會高于前段時間的均值,或是有走低的趨勢,這值得生產(chǎn)者注意調(diào)整適當?shù)酿B(yǎng)殖方法。時間序列分析方法可以用在養(yǎng)豬生產(chǎn)中,可以預測未來一段時間的生產(chǎn)變化發(fā)展。
表5 ARIMA模型預測母豬分娩率效果表 %
表6 ARIMA模型誤差評估
表7 與母豬受胎率相關因素的logistic回歸結果
表8 與母豬返情率相關因素的logistic回歸結果
對受胎率、返情率、分娩率和7 d斷配率的二元logistic回歸分析,表明受胎率、返情率、分娩率和7 d斷配率均與年份和胎次有回歸關系。此研究發(fā)現(xiàn)在2017年受胎率高于2016年,返情率、分娩率和7 d斷配率低于2016年,說明2017年該場母豬配種要比2016年更好,而分娩率低可能是受胎后其他因素導致的。2018和2019年受胎率和7 d斷配率高于2016年,返情率和分娩率低于2016年,分娩率進一步降低可能和非洲豬瘟的暴發(fā)有關。而隨著胎次的增加受胎率、分娩率和7 d斷配率下降,返情率上升。此結果與朱世平等得出母豬繁殖的最佳胎次為3~6胎相反[7],說明該公司繁殖母豬的繁殖性能仍有待提高。該公司豬場年份之間具有差異,特別是2017年性能比較好,后來隨著胎次增加繁殖性能不增反減,可能是有新建場區(qū),母豬養(yǎng)殖水平和設施等還未跟上,也可能是豬場管理制度不嚴,需要從飼養(yǎng)管理方案到技術設備、母豬淘汰制度、發(fā)情鑒定和足夠的后備母豬儲備、專業(yè)的人員操作進行完善[19],細化采精、配種和分娩等記錄,提高豬場效率。
表9 與母豬分娩率相關因素的logistic回歸結果
表10 與母豬7 d斷配率相關因素的logistic回歸結果
時間序列分析方法對于母豬繁殖性能方面有一定的短期預測效果,可進一步探索用于其他生產(chǎn)指標。相較于之前的研究,此次研究根據(jù)數(shù)據(jù)更全面分析了母豬的繁殖性能。并應用logistics模型和時間序列ARIMA模型對母豬繁殖性能進行因素分析及短期預測,試驗數(shù)據(jù)量大,為生產(chǎn)提供了更具參考價值的依據(jù)。