中國荷斯坦牛乳房炎風險評估模型建立及預測分析

史良玉，李文龍，唐永杰，米思遠，肖 煒，劉 林，張 毅，王雅春，俞 英*

（1.中國農業大學動物科學技術學院，北京 100193；2.北京市畜牧總站，北京 100029；3.北京奶牛中心，北京 100192）

奶牛乳房炎不僅使奶牛產奶量顯著下降，造成嚴重經濟損失，還可導致乳汁成分發生改變，使牛奶的營養價值和食用價值顯著下降[1-2]。奶牛發生乳房炎時，由于炎性反應程度不同，乳汁成分的改變程度也不盡相同。根據乳房或乳汁有無肉眼可見的變化，奶牛乳房炎分為臨床乳房炎和隱性乳房炎2 類。

目前，測定乳房炎反應的指標是奶牛群體改良（Dairy Herd Improvement，DHI）測定記錄中的乳汁體細胞數（Somatic Cell Count，SCC）。DHI 是一套針對奶牛泌乳性能及乳成分的完整奶牛生產性能記錄體系。目前國際上奶牛隱性乳房炎乳SCC 的判斷標準多為20 萬/mL～50 萬/mL。馬裴裴等[3]研究表明，北京地區中國荷斯坦牛群宜以10 萬/mL～50 萬/mL SCC 作為隱性乳房炎標準。為克服SCC 在統計分析中的不足，通常將其轉換成遵循正態分布的體細胞評分（Somatic Cell Score，SCS）的形式[4]。

及早發現奶牛乳房炎有利于患病奶牛的治療，進而可以減少與乳房炎相關的經濟損失[5]。奶牛乳房炎風險預測已有多種模型，其中Logistic 回歸、深度學習、隨機森林等模型的準確率及預測價值并無顯著差異[6]。奶牛乳房炎風險評估方面的研究主要集中在分析風險因素上，而在奶牛乳房炎預測方面鮮有研究。Logistic 回歸通過分析一組預測變量（自變量），預測一個或多個響應變量。這種統計分析方法可用于評估預測變量對響應變量的預期效果。本研究根據已獲得的DHI 記錄建立Logistic 回歸模型，利用Logistic 回歸對中國荷斯坦牛DHI 數據中各測定指標進行分析，篩選出對不同程度奶牛乳房炎的預測具有統計學意義的風險因素；同時構建Logistic 風險評估模型，以實現提前發現具有乳房炎高發病風險的奶牛，進而在生產中提前防范，防止奶牛乳房炎的發生。

1 材料與方法

1.1 實驗材料 來自北京市奶牛中心1998—2016 年的北京地區DHI 數據，共包括121 個牛場78 386 頭中國荷斯坦牛的1 967 310 條測定記錄。數據內容包括個體號、胎次、產犢日期、測定日期、產奶量、乳脂量、乳脂率、乳蛋白量、乳蛋白率和SCC 等。

1.2 實驗方法

1.2.1 數據整理 利用R 3.4.1 軟件對DHI 數據進行整理。針對獲得的原始DHI 記錄，刪除SCC 缺失的牛只記錄，由于2009—2013 年DHI 記錄中未包括原始SCC值，故將其刪除。Santman 等[7-9]認為，產犢后前4 d的SCC 值會偏高，不能作為標準，故在本研究中予以剔除。篩選同一頭奶牛連續2 個月都具有DHI 的記錄，并根據其中第2 個月的SCC 將中國荷斯坦牛分為健康組、隱性乳房炎組以及臨床乳房炎組（表1）[3]。由于SCC 值呈偏態分布且方差不同質，通常將其轉化為接近正態分布的SCS，再進行遺傳統計分析。利用國際通用的轉換公式“SCS=log2（SCC/100 000）+3”將SCC 值轉化為SCS。為提高后期模型預測的精度，取大于0 并小于10 的SCS 作為數據集[10]。剔除后，共得781 611條DHI 記錄。同時，劃定200～400 頭泌乳牛規模的奶牛場為小型奶牛場（Small），400～800 頭泌乳牛為中型奶牛場（Medium），大于800 頭泌乳牛為大型奶牛場（Large）。

表1 中國荷斯坦牛健康狀況劃分

1.2.2 奶牛乳房炎單因素分析及入選變量 使用北京地區DHI 數據，利用SAS 9.2 軟件對影響奶牛乳房炎預測的各個因素進行單因素分析，對單個因素進行初步篩選，探究單個因素對奶牛乳房炎的影響大小，得到對奶牛乳房炎預測的影響有統計學差異的指標。對于分類變量的分析，通過計算對數比率即Logit（p）確定自變量進入模型的形式。

1.2.3 奶牛乳房炎多因素分析及Logistic 回歸模型建立根據北京地區DHI 數據單因素分析的結果，確定健康組、隱性乳房炎組和臨床乳房炎組組間有統計學意義的指標，即奶牛乳房炎風險因素（表2）。

表2 北京地區中國荷斯坦牛乳房炎風險因素

將這些指標作為入選變量代入Logistic 回歸方程中，利用SAS 9.2 軟件進行Logistic 回歸分析。公式如下：

1.2.4 ROC 曲線的繪制 ROC 曲線主要用于反映預測變量對響應變量的預測準確率，是一條由點連接而成的曲線，這里體現奶牛乳房炎風險評估模型的預測準確率。其中Y 軸為靈敏度，X 軸為特異度。靈敏度是指一項診斷試驗能將實際患病的奶牛正確地診斷為病牛的概率，特異度是指一項診斷試驗能將實際無病的奶牛正確診斷為非病牛的概率。靈敏度越高，說明真正的病牛被診斷出來的可能性越大；特異度越高，說明真正的非病牛被排除的可能性越大。曲線下面積（Area Under the Curve，AUC）即模型診斷價值，一般用來評價試驗的診斷性能，面積越大，說明該試驗的診斷性能越好。一般來說，AUC 在0.50～0.70 之間被認為其辨別力一般，AUC 在0.70～0.90 之間被認為該模型良好，AUC 高于0.90 被認為該模型是優秀的。

1.2.5 Logistic 回歸模型預測準確率的評估 針對構建的Logistic 風險評估模型，應用北京地區某荷斯坦牛場2019 年全年共11 338 頭次奶牛的DHI 數據，對Logistic風險評估模型預測奶牛乳房炎的準確率進行評估。其中，DHI 數據內容包括個體號、胎次、產犢日期、測定日期、產奶相關記錄和SCC 等。利用R 3.4.1 軟件對DHI 數據進行整理，并應用Logistic 回歸模型對牛場DHI 數據中各測定指標進行分析，比較每月預測情況與下月奶牛乳房炎實際情況之間的一致性以驗證模型預測結果的準確率。奶牛乳房炎發病判定標準為SCC＞10 萬/mL。其中，10 萬/mL～50 萬/mL SCC 作為隱性乳房炎標準，大于50 萬/mL SCC 作為臨床乳房炎標準。

2 結果

2.1 奶牛乳房炎風險因素（自變量）進入風險評估模型的形式 利用北京地區DHI 數據進行分析，根據計算奶牛隱性乳房炎組及臨床乳房炎組Logit（p），并繪制折線圖，以確定變量進入模型的形式（圖1）。通過計算場規模效應可知，小型牛場及中型牛場的Logit（p）較為接近，故將小型牛場及中型牛場合并為一個變量，即中小型牛場（圖1A）；同理，測定月份在6、7、8月份時Logit（p）與其他月份明顯不同，在臨床乳房炎組中尤為明顯，而6、7、8 月份北京地區處于夏季階段，故測定時間變量分為2 個，即夏季及非夏季[11]（圖1B）。胎次變量及泌乳階段的Logit（p）隨等級增加基本呈線性變化（圖1C、圖1D）。本月SCS 對應的隱性乳房炎Logit（p）（圖1E）和臨床乳房炎Logit（p）（圖1F）的平均值呈線性變化。可見，合并變量之后，各變量均隨等級增加呈線性變化，可后續分析。

2.2 奶牛隱性乳房炎組及臨床乳房炎組各變量解釋 使用北京奶牛中心提供的DHI 數據進行奶牛乳房炎風險評估。選用“場規模+胎次+測定季節+泌乳階段”為變量代入多因素Logistic 回歸方程。最終可用于診斷奶牛下一個測定月隱性乳房炎患病風險的解釋變量有4個，具體包括場規模、胎次、測定季節、泌乳階段（表3）。

隱性乳房炎風險預測結果顯示，與處于大型奶牛場中的奶牛相比，處于中小型場中的奶牛更可能在下一測定月中患隱性乳房炎（OR=1.184，P＜0.001）；處于2胎次的奶牛下一測定月患隱性乳房炎的可能性是初產奶牛的1.377 倍，處于3 胎次的奶牛下一測定月患隱性乳房炎的可能性是初產奶牛的1.647 倍（P＜0.001）；當前測定時間為夏季時，中國荷斯坦牛下一月患隱性乳房炎的可能性是非夏季的1.092 倍（P＜0.001）；隨著泌乳階段的增加，下一月中國荷斯坦牛患隱性乳房炎的可能性也增加，尤其是泌乳天數＞300 d 的牛只。

與隱性乳房炎組風險評估過程相同，同樣采用“場規模+胎次+測定季節+泌乳階段”為變量代入多因素Logistic 回歸方程，解釋奶牛下一個測定月具有患臨床乳房炎風險的最佳解釋變量同樣包括了場規模、胎次、測定季節、泌乳階段，并且趨勢相同（表4）。處于中小型奶牛場中的奶牛下一測定月更可能患臨床乳房炎（OR=1.483，P＜0.001）；隨著奶牛胎次的增加，下月患臨床乳房炎的可能性也不斷增加，且第3 胎次的奶牛下月患臨床乳房炎的可能性是初產奶牛的2.203倍（P＜0.001）；當前測定月為夏季時，奶牛下一測定月患臨床乳房炎的可能性比非夏季高（OR=1.266，P＜0.001）；隨著泌乳階段的增加，下一測定月患臨床乳房炎的可能性也隨之增加，泌乳天數高于300 d 的奶牛下月患臨床乳房炎的可能性是泌乳天數小于100 d 奶牛的2.121倍，且增加倍數與隱性乳房炎組基本相同（P＜0.001）。

2.3 奶牛隱性乳房炎及臨床乳房炎風險評估模型的預測價值 使用北京奶牛中心提供的DHI 數據進行風險評估。選用“場規模+胎次+測定季節+泌乳階段+本月SCS 值”為變量代入多因素Logistic 回歸方程，并繪制ROC 曲線（圖2）。根據分析可得，隱性乳房炎風險評估模型預測價值為0.721，準確率為67.6%，特異度為80.3%，靈敏度為50.7%；臨床乳房炎風險評估模型預測價值為0.825，準確率為83.6%，特異度為94.9%，靈敏度為42.5%，表明兩模型預測價值均良好。另外，從ROC 曲線中也可以得出，本月SCS 對奶牛乳房炎風險評估模型預測價值貢獻最大，其次為胎次及泌乳階段。

圖1 隱性乳房炎組及臨床乳房炎組自變量Logit（p）

表3 隱性乳房炎組風險因素

表4 臨床乳房炎組風險因素

圖2 隱性乳房炎組及臨床乳房炎組預測 ROC 曲線

2.4 奶牛乳房炎風險評估模型的預測準確率 利用Logistic回歸模型對北京地區某荷斯坦牛場2019 年全年共11 338頭次奶牛的乳房炎發病風險進行評估，根據評估結果進行分析（表5）發現，Logistic 回歸模型對該牛場奶牛的乳房炎風險評估的準確率平均為70.84%，具備應用前景。

3 討 論

奶牛場中使用各種控制和監測奶牛乳腺健康的方法主要是為了減少新發病奶牛的數量，減少已有的患乳房炎奶牛，并減少奶牛患病的持續時間。奶牛乳房炎風險因素的分析能夠幫助鑒定奶牛乳房健康情況及發病風險，對于牛場中奶牛乳房炎的控制及減少新感染的牛只具有重要作用[12-13]。

表5 北京地區某荷斯坦牛場2019 年全年奶牛乳房炎風險預測結果準確率

有效利用DHI 可進行乳房炎風險評估。每月測定的DHI 報告為奶牛乳房炎的監測提供了大量的數據，其中SCC 值是奶牛乳房炎檢測的重要工具。目前，SCC 值已經在預測牛場奶牛乳房炎中廣泛應用[14]。許多因素會影響SCC 值，包括奶牛場管理、胎次、泌乳階段以及季節和月份等[15-19]。本研究表明，中小型奶牛場奶牛乳房炎發病率高于大型奶牛場，可能是由于大型奶牛場管理規范，牛舍中環境相對較好，由病原菌引起的奶牛乳房炎發病率在逐漸降低；擠奶操作和治療情況的不斷改善，使得奶牛乳房炎損傷發生情況減少，故奶牛乳房炎發病率相對較低。另外，多項研究證明，奶牛的年齡及胎次越高，患奶牛乳房炎的風險就越高[16,20]。由于不同季節溫濕度的不同，如夏季屬于雨季，濕度較高，冬季較為干燥，奶牛乳房炎發生的概率也不同，奶牛乳房炎多發生于夏季[21]。隨著泌乳天數的增加，患奶牛乳房炎的風險也會增加，尤其是產奶天數超過200 d之后。奶牛乳房炎受環境中病原體的影響較大，并且被傳染性致病菌感染的奶牛為細菌提供了繁殖的場所，進一步傳染給其他健康奶牛[22]。SCC 值隨著泌乳天數的增加而增加可能是因為奶牛機體對抗乳房組織因擠奶造成的損傷，由此增加了細菌黏附于機體的機會，從而導致乳房炎[23-25]。

在奶業生產中，奶牛乳房炎是傳播廣泛并且花費極昂貴的疾病之一。乳房炎發病風險預測有利于防病于未患，進而減少相關的經濟損失[5]。已有研究利用深度學習、隨機森林等模型對奶牛乳房炎進行風險預測，與本研究所用Logistic 回歸的準確性及預測價值相比并無顯著差異[6]。Jadhav 等[26]對214 頭奶牛進行隱性乳腺炎預測，預測價值為0.930。該預測價值較本研究高，可能是由于該研究中對于奶牛健康狀態判斷標準與本研究不同且僅預測采集奶樣當天奶牛的狀態，并且模型中加入乳房衛生情況、臥床衛生情況及擠奶方法等DHI 報告中未包含的信息。在實際生產中，對乳房、臥床衛生情況打分會存在一定主觀因素，并且會增加奶牛場工作量。本研究利用DHI 報告中包含的信息預測下一月奶牛健康狀態，應用范圍廣，且奶牛乳房炎模型預測價值良好，為早期發現奶牛乳房炎提供一定參考。

4 結 論

本研究利用DHI 記錄分析荷斯坦牛乳房炎發病的風險因素，通過對覆蓋120 余個牛場近20 年的196 萬多條DHI 數據記錄進行分析，基于多因素Logistic 回歸分析的結果表明，場規模越小、奶牛胎次越高、夏季階段測定、泌乳階段越高，奶牛患乳房炎風險越高；使用DHI 測定記錄進行Logistic 回歸模型的構建，所得奶牛隱性乳房炎風險評估模型的預測價值及準確率均低于臨床乳房炎風險評估模型；Logistic 回歸分析后得到的奶牛隱性乳房炎風險評估模型和奶牛臨床乳房炎風險評估模型預測價值分別為0.721 和0.825，模型預測價值良好，可以應用于實際荷斯坦牛群的乳房炎風險評估。