某豬場(chǎng)母豬繁殖性能預(yù)測(cè)和logistic回歸分析

翟靜靜，劉則彰，蔣思文

（華中農(nóng)業(yè)大學(xué)，湖北 武漢 430070）

母豬在養(yǎng)殖生產(chǎn)中扮演著重要的角色，在規(guī)?；B(yǎng)豬生產(chǎn)中，母豬的繁殖效率是重中之重，是養(yǎng)殖場(chǎng)的核心，也是決定豬場(chǎng)生產(chǎn)水平和經(jīng)濟(jì)效益的關(guān)鍵因素[1]。母豬繁殖性能主要指標(biāo)包括總產(chǎn)仔數(shù)、產(chǎn)活仔數(shù)、斷奶仔豬數(shù)、斷奶發(fā)情間隔（Weaning to estrus interval,WEI）、PSY等[2]，還包括繁殖力（如分娩率和死亡率）。

胎次也是影響母豬繁殖性能的重要因素，一般來(lái)說(shuō)，繁殖性能隨著胎次數(shù)的增加而提高，在3～5胎次時(shí)達(dá)到最高水平[3,4]。在2～4胎時(shí)母豬的受胎率會(huì)隨胎次升高而增加，4胎后逐漸下降[5]。胎次主要通過(guò)影響分娩率和產(chǎn)仔數(shù)來(lái)影響母豬繁殖性能[6]。朱世平等發(fā)現(xiàn)，依據(jù)胎次可將母豬繁殖性能劃分為1-2胎、3-6胎和7-8胎這三個(gè)階段，并且得出母豬繁殖的最佳胎次為3-6胎，繁殖高峰期是4胎次（P＜0.05）[7]。

WEI是指母豬斷奶到再次發(fā)情配種的間隔時(shí)間。據(jù)國(guó)外的各個(gè)豬場(chǎng)數(shù)據(jù)分析，WEI與產(chǎn)仔數(shù)呈負(fù)相關(guān)，WEI時(shí)間越長(zhǎng)，產(chǎn)仔數(shù)越低，而在3～5 d發(fā)情的母豬，分娩率高。因此要想提高斷奶母豬7 d內(nèi)的配種率，應(yīng)縮短WEI，保持母豬的理想體況，保障產(chǎn)仔舍和配懷舍的飼養(yǎng)管理工作，減少體重?fù)p失，做好誘情工作，提高分娩成績(jī)[7]。

生產(chǎn)中把繁殖母豬發(fā)情期進(jìn)行配種后沒(méi)有懷孕并重新發(fā)情的現(xiàn)象稱為返情。常規(guī)返情是受胎失敗，非常規(guī)返情是妊娠失敗，都會(huì)導(dǎo)致低分娩率的出現(xiàn)。母豬一般配種受胎為80%～90%，而配種未分娩的10%～20%就是返情。返情的增加，會(huì)導(dǎo)致配種分娩率降低，從而影響?zhàn)B殖場(chǎng)的經(jīng)濟(jì)效益[8,9]。返情后再配種母豬分娩率降低10%左右[10]。實(shí)際生產(chǎn)中母豬返情率一般在10%左右[11]。

時(shí)間序列法是一種統(tǒng)計(jì)分析方法，將研究對(duì)象的數(shù)量隨著時(shí)間的發(fā)生而呈現(xiàn)變動(dòng)的這種特征按照時(shí)間順序進(jìn)行排列，所形成的這個(gè)序列便稱為研究對(duì)象的時(shí)間序列[12]。常用的時(shí)間序列法有移動(dòng)平均法、加權(quán)移動(dòng)平均法和指數(shù)平均法，時(shí)間序列的平穩(wěn)性是建模必不可少的環(huán)節(jié)[13]。時(shí)間序列分為兩類：第一類是某時(shí)間序列各時(shí)期的數(shù)值是不變動(dòng)的，為確定型時(shí)間序列。第二類是此時(shí)間序列各時(shí)期數(shù)值變動(dòng)，充滿不確定性，為隨機(jī)型時(shí)間序列。一般情況下都是隨機(jī)型的時(shí)間序列。朱迪國(guó)等對(duì)全球豬流行性腹瀉疫情進(jìn)行時(shí)間序列分析發(fā)現(xiàn)，豬血液制品貿(mào)易可能是該病跨國(guó)界傳播的潛在途徑[14]。因此，在母豬繁殖生產(chǎn)數(shù)據(jù)中，可以運(yùn)用時(shí)間序列對(duì)性狀的時(shí)間發(fā)展進(jìn)行預(yù)測(cè)。

對(duì)于因變量為只具有兩種分類結(jié)果的屬性數(shù)據(jù)，例如母豬生產(chǎn)中的母豬分娩或未分娩，就需要去建立邏輯回歸模型，這是一種經(jīng)典的分類算法[15]。Logistic模型中的參數(shù)是通過(guò)極大似然法進(jìn)行估計(jì)的。邏輯回歸分析是通過(guò)一組預(yù)測(cè)變量（自變量）預(yù)測(cè)一個(gè)或多個(gè)響應(yīng)變量（因變量）一種統(tǒng)計(jì)分析方法[16]。在使用Logistic模型前，需判斷是否滿足以下幾個(gè)條件：(1)因變量是二分類變量；(2)分類變量的分類是全面的且相互獨(dú)立；(3) 連續(xù)的自變量與因變量的對(duì)數(shù)轉(zhuǎn)換值之間存在線性關(guān)系；(4)自變量之間不存在共線性。利用建立的回歸模型可以對(duì)母豬性能及后期發(fā)展進(jìn)行預(yù)測(cè)。Logistic模型也可以用來(lái)評(píng)估預(yù)測(cè)變量對(duì)響應(yīng)變量的預(yù)期效果。Hosmer and Lemeshow運(yùn)用Logistic模型預(yù)測(cè)母豬體重發(fā)育相關(guān)變量與第2胎中的非妊娠和產(chǎn)仔數(shù)的關(guān)聯(lián)，發(fā)現(xiàn)是高度相關(guān)[17]。Koketsu等人利用Logistic回歸模型分析與分娩率相關(guān)的因素時(shí)發(fā)現(xiàn)，胎次、分娩季節(jié)以及泌乳天數(shù)與采食量之間的交互作用會(huì)對(duì)母豬分娩率產(chǎn)生影響[18]。

研究根據(jù)豬場(chǎng)的生產(chǎn)記錄，時(shí)間序列和Logistic回歸模型這兩種方法，把所得的豬場(chǎng)數(shù)據(jù)進(jìn)行匯總處理，運(yùn)用時(shí)間序列去進(jìn)行短期預(yù)測(cè)母豬繁殖性能高低的時(shí)間段和母豬狀態(tài)，來(lái)判斷豬場(chǎng)中繁殖母豬在某一時(shí)間段可能會(huì)出現(xiàn)的情況；通過(guò)Logistic回歸模型進(jìn)一步剖析影響因素，對(duì)影響該豬場(chǎng)母豬繁殖性能的因素進(jìn)行調(diào)查研究，為該豬場(chǎng)采取相應(yīng)措施提高母豬繁殖性能提供科學(xué)依據(jù)，有效提高豬場(chǎng)經(jīng)濟(jì)效益。

1 材料與方法

1.1 數(shù)據(jù)來(lái)源

此研究選擇某公司的2016-2019年7個(gè)分場(chǎng)的母豬繁殖性能記錄數(shù)據(jù)，包括產(chǎn)仔性能、配種分娩明細(xì)及淘汰死亡等。以母豬繁殖性能為主線，對(duì)規(guī)模豬場(chǎng)不同生產(chǎn)指標(biāo)的各種數(shù)據(jù)進(jìn)行處理分析，根據(jù)對(duì)生產(chǎn)數(shù)據(jù)的分析，找到影響繁殖性能的因素，制定行之有效的解決方案，以此展開(kāi)相關(guān)研究。

1.2 數(shù)據(jù)統(tǒng)計(jì)分析

對(duì)數(shù)據(jù)采用Excel軟件進(jìn)行統(tǒng)計(jì)，數(shù)據(jù)以平均值±標(biāo)準(zhǔn)差表示；運(yùn)用GraphPad Prism 7.0作圖，使用SPSS 19.0建立邏輯回歸模型和對(duì)數(shù)據(jù)進(jìn)行顯著性檢驗(yàn)，P＜0.05表示差異顯著，P＜0.01表示差異極顯著；運(yùn)用SAS 9.4軟件構(gòu)建時(shí)間序列模型ARIMA進(jìn)行預(yù)測(cè)。

1.3 Logistic模型的構(gòu)建

此研究采用二元Logistic回歸分析方法對(duì)母豬胎次和年份進(jìn)行了分析。母豬繁殖性能相關(guān)因素分析研究的變量與賦值見(jiàn)表1。根據(jù)收集到的數(shù)據(jù)，將數(shù)據(jù)分為兩類情況，如受胎：1和未受胎：0；首先采用多重共線性模型篩選變量（|r|＞0.7）。然后，對(duì)每個(gè)解釋變量進(jìn)行單因素分析，將上述變量P＜0.1的解釋變量納入logistic回歸分析模型中，運(yùn)用基于最大似然估計(jì)的向前逐步回歸法（Forward: LR法）建 立logistic回歸模型，確定與母豬繁殖性能相關(guān)風(fēng)險(xiǎn)因素的依據(jù)是似然比統(tǒng)計(jì)量的P值＜0.05。回歸系數(shù)用95%置信區(qū)間（Conf idence intervals，CI）的似然比（Odds ratios，OR）表示。OR值即是相對(duì)危險(xiǎn)度的精確估計(jì)值，如果OR大于1，變量施加的風(fēng)險(xiǎn)增加。Logistic回歸模型的公式如下：

其 中，β0為 常 量，A、B分別表示因素?cái)嗄袒蚺浞N年份、母豬斷奶或配種胎次。把胎次、年份等影響7 d斷配率、返情率、受胎率、分娩率這些指標(biāo)的因素進(jìn)行定義變量和賦值（見(jiàn)表1）。

表1 母豬繁殖性能相關(guān)因素分析研究的變量與賦值

2 結(jié)果與分析

2.1 時(shí)間序列分析預(yù)測(cè)

2.1.1 時(shí)序圖檢驗(yàn)

此研究收集了2016-2019.11、2016-2019.08、2016-2019.07、2016-2019.07的豬場(chǎng)健仔數(shù)、斷奶頭數(shù)、7 d斷配率、分娩率的數(shù)據(jù)。選取了其中80%的數(shù)據(jù)用于時(shí)間序列分析，留下20%的數(shù)據(jù)作為實(shí)際值與預(yù)測(cè)值對(duì)比。首先對(duì)原始序列進(jìn)行平穩(wěn)性檢驗(yàn)，分別描繪平均每頭母豬健仔數(shù)（圖1a）、斷奶頭數(shù)（圖1b）、7 d斷配率（圖1c）、分娩率（圖1d）的時(shí)序圖，從時(shí)序圖中可以看出該時(shí)間序列具有明顯的遞增趨勢(shì)性和下降趨勢(shì)性，為典型的非平穩(wěn)性序列。

2.1.2 一階差分時(shí)序圖

時(shí)間序列預(yù)處理，可知這四個(gè)指標(biāo)均為不平穩(wěn)序列，ARIMA模型基于平穩(wěn)時(shí)間序列的基礎(chǔ)上，因此根據(jù)數(shù)據(jù)特征對(duì)數(shù)據(jù)進(jìn)行差分處理，一階差分后的時(shí)序圖（圖2）顯示，各觀測(cè)值圍繞著一個(gè)常數(shù)附近上下波動(dòng)，數(shù)據(jù)逐漸趨于平穩(wěn)，序列表現(xiàn)平穩(wěn)的性質(zhì)。

2.1.3 模型識(shí)別

為了驗(yàn)證一階差分后此序列的平穩(wěn)性，依次做出該序列的自相關(guān)圖和偏相關(guān)圖（圖3到圖6）。觀察ACF（自相關(guān)）圖和拖尾效應(yīng)確定ARIMA模型的P值，觀察PACF（偏相關(guān)）和截尾效應(yīng)確定ARIMA模型的q值。其中差分的次數(shù)就是模型ARIMA（p， d，q）中的階數(shù)d，因此通過(guò)最佳準(zhǔn)則函數(shù)定價(jià)方法選擇ARIMA（0，1，1）、ARIMA（5，1，5）、ARIMA（0，1，0）、ARIMA（0，1，0）分別為健仔數(shù)、斷奶頭數(shù)、7 d斷配率、分娩率的相對(duì)最優(yōu)模型。

2.1.4 模型的預(yù)測(cè)

根據(jù)已經(jīng)建立的ARIMA模型對(duì)母豬健仔數(shù)、斷奶頭數(shù)、7 d斷配率及分娩率的數(shù)據(jù)進(jìn)行預(yù)測(cè)。預(yù)測(cè)時(shí)間為9個(gè)月，健仔數(shù)的預(yù)測(cè)時(shí)間從2019年3月到11月，斷奶頭數(shù)預(yù)測(cè)時(shí)間為2018年12月到2019年8月，7 d斷配率和分娩率預(yù)測(cè)時(shí)間都是從2018年11月到2019年7月。為評(píng)估模型的預(yù)測(cè)效果，引入均方誤差（Mean square error，MSE）和平均絕對(duì)百分比誤差（Mean absolute percentage error，MAPE）兩個(gè)指標(biāo)。預(yù)測(cè)的偏差較小，與實(shí)際觀察值基本吻合（表2到表5）。預(yù)測(cè)值的變化趨勢(shì)與實(shí)際值大致保持一致，實(shí)際值基本落在預(yù)測(cè)曲線95%的可信區(qū)間內(nèi)，該模型擬合預(yù)測(cè)成功，具有一定的實(shí)際價(jià)值。由表6誤差評(píng)估表可見(jiàn)，MSE值約小，說(shuō)明預(yù)測(cè)模型具有更好的精確度，MAPE＜20%，可以從中判斷該模型擬合預(yù)測(cè)很成功。

2.2 因素分析

2.2.1 受胎率的logistic回歸分析

原始數(shù)據(jù)為192 874條配種記錄，剔除掉胎次為0的數(shù)據(jù)，剩下143 464條記錄用于后續(xù)的因素分析。由表7可知，配種年份和配種胎次會(huì)影響母豬受胎率。高胎次（≥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個(gè)分區(qū)豬場(chǎng)192 874條配種記錄，剔除掉胎次為0的數(shù)據(jù)，剩下143 464條記錄用于后續(xù)的因素分析。由表8可知，配種年份和配種胎次會(huì)影響母豬返情率。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模型預(yù)測(cè)母豬健仔數(shù)效果表 頭

表3 ARIMA模型預(yù)測(cè)母豬斷奶頭數(shù)效果表 頭

表4 ARIMA模型預(yù)測(cè)母豬7 d斷配率效果表 %

原始數(shù)據(jù)為192 874條配種記錄，剔除掉胎次為0的數(shù)據(jù)，剩下143 464條記錄用于后續(xù)的因素分析。由表9可知，配種年份和配種胎次會(huì)影響母豬分娩率。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個(gè)分區(qū)豬場(chǎng)2016-2019年間142 478條斷奶配種記錄，剔除掉WEI＞100的數(shù)據(jù)，剩下138 114條記錄用于后續(xù)的因素分析。由表 10可知，斷奶年份和斷奶胎次會(huì)影響母豬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)的母豬低。

3 小結(jié)與討論

數(shù)據(jù)的時(shí)間序列大多是不平穩(wěn)的，數(shù)據(jù)本身具有季節(jié)趨勢(shì)，若采用簡(jiǎn)單的模型進(jìn)行分析，往往不能得到準(zhǔn)確的預(yù)測(cè)結(jié)果。隨著對(duì)時(shí)間研究的發(fā)展，時(shí)間序列研究的方法更加多樣，其中ARIMA模型是常用的非平穩(wěn)的模型[13]。時(shí)間序列預(yù)測(cè)法對(duì)于中短期預(yù)測(cè)的效果要優(yōu)于長(zhǎng)期預(yù)測(cè)的效果，因?yàn)闀r(shí)間過(guò)長(zhǎng)發(fā)生外界因素變化的可能性會(huì)加大。世界上沒(méi)有一個(gè)模型可以做到100%的精確，只能隨著時(shí)間發(fā)展進(jìn)行改善，不斷添加新的數(shù)值和因素。此研究中以健仔數(shù)、斷奶頭數(shù)、7 d斷配率和分娩率的實(shí)際數(shù)據(jù)為基礎(chǔ)，擬合ARIMA模型，同時(shí)對(duì)未來(lái)9個(gè)月的數(shù)據(jù)進(jìn)行了預(yù)測(cè)和比較。研究結(jié)果表明健仔數(shù)在未來(lái)幾個(gè)月呈上升趨勢(shì)，與實(shí)際值相比最大偏差為0.8頭；斷奶頭數(shù)預(yù)測(cè)值具有一定的波動(dòng)性，與實(shí)際值最大偏差為1.67；7 d斷配率預(yù)測(cè)值持續(xù)上升，實(shí)際值時(shí)而下降時(shí)而上升，最大偏差為11.55個(gè)百分點(diǎn)；分娩率預(yù)測(cè)結(jié)果呈下降趨勢(shì)，實(shí)際值具有波動(dòng)，最大偏差為4.26個(gè)百分點(diǎn)。這幾個(gè)指標(biāo)雖有一定偏差，但實(shí)際值基本落在預(yù)測(cè)曲線95%的置信區(qū)間內(nèi)，說(shuō)明模型對(duì)原始序列信息的提取比較充分。通過(guò)預(yù)測(cè)效果表和MAPE值（值越小，代表預(yù)測(cè)的準(zhǔn)確度越強(qiáng)）來(lái)看，利用ARIMA模型對(duì)這些數(shù)據(jù)進(jìn)行擬合和預(yù)測(cè)是合理可行的，但也會(huì)有一定的局限性，生產(chǎn)過(guò)程中可能會(huì)因環(huán)境變化或疾病、飼養(yǎng)管理方案改變引起差異。在短期預(yù)測(cè)當(dāng)中，豬場(chǎng)未來(lái)一段時(shí)間的母豬生產(chǎn)情況將會(huì)高于前段時(shí)間的均值，或是有走低的趨勢(shì)，這值得生產(chǎn)者注意調(diào)整適當(dāng)?shù)酿B(yǎng)殖方法。時(shí)間序列分析方法可以用在養(yǎng)豬生產(chǎn)中，可以預(yù)測(cè)未來(lái)一段時(shí)間的生產(chǎn)變化發(fā)展。

表5 ARIMA模型預(yù)測(cè)母豬分娩率效果表 %

表6 ARIMA模型誤差評(píng)估

表7 與母豬受胎率相關(guān)因素的logistic回歸結(jié)果

表8 與母豬返情率相關(guān)因素的logistic回歸結(jié)果

對(duì)受胎率、返情率、分娩率和7 d斷配率的二元logistic回歸分析，表明受胎率、返情率、分娩率和7 d斷配率均與年份和胎次有回歸關(guān)系。此研究發(fā)現(xiàn)在2017年受胎率高于2016年，返情率、分娩率和7 d斷配率低于2016年，說(shuō)明2017年該場(chǎng)母豬配種要比2016年更好，而分娩率低可能是受胎后其他因素導(dǎo)致的。2018和2019年受胎率和7 d斷配率高于2016年，返情率和分娩率低于2016年，分娩率進(jìn)一步降低可能和非洲豬瘟的暴發(fā)有關(guān)。而隨著胎次的增加受胎率、分娩率和7 d斷配率下降，返情率上升。此結(jié)果與朱世平等得出母豬繁殖的最佳胎次為3～6胎相反[7]，說(shuō)明該公司繁殖母豬的繁殖性能仍有待提高。該公司豬場(chǎng)年份之間具有差異，特別是2017年性能比較好，后來(lái)隨著胎次增加繁殖性能不增反減，可能是有新建場(chǎng)區(qū)，母豬養(yǎng)殖水平和設(shè)施等還未跟上，也可能是豬場(chǎng)管理制度不嚴(yán)，需要從飼養(yǎng)管理方案到技術(shù)設(shè)備、母豬淘汰制度、發(fā)情鑒定和足夠的后備母豬儲(chǔ)備、專業(yè)的人員操作進(jìn)行完善[19]，細(xì)化采精、配種和分娩等記錄，提高豬場(chǎng)效率。

表9 與母豬分娩率相關(guān)因素的logistic回歸結(jié)果

表10 與母豬7 d斷配率相關(guān)因素的logistic回歸結(jié)果

時(shí)間序列分析方法對(duì)于母豬繁殖性能方面有一定的短期預(yù)測(cè)效果，可進(jìn)一步探索用于其他生產(chǎn)指標(biāo)。相較于之前的研究，此次研究根據(jù)數(shù)據(jù)更全面分析了母豬的繁殖性能。并應(yīng)用logistics模型和時(shí)間序列ARIMA模型對(duì)母豬繁殖性能進(jìn)行因素分析及短期預(yù)測(cè)，試驗(yàn)數(shù)據(jù)量大，為生產(chǎn)提供了更具參考價(jià)值的依據(jù)。