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

（1.中國農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，北京 100193；2.北京市畜牧總站，北京 100029；3.北京奶牛中心，北京 100192）

奶牛乳房炎不僅使奶牛產(chǎn)奶量顯著下降，造成嚴(yán)重經(jīng)濟(jì)損失，還可導(dǎo)致乳汁成分發(fā)生改變，使牛奶的營(yíng)養(yǎng)價(jià)值和食用價(jià)值顯著下降[1-2]。奶牛發(fā)生乳房炎時(shí)，由于炎性反應(yīng)程度不同，乳汁成分的改變程度也不盡相同。根據(jù)乳房或乳汁有無肉眼可見的變化，奶牛乳房炎分為臨床乳房炎和隱性乳房炎2 類。

目前，測(cè)定乳房炎反應(yīng)的指標(biāo)是奶牛群體改良（Dairy Herd Improvement，DHI）測(cè)定記錄中的乳汁體細(xì)胞數(shù)（Somatic Cell Count，SCC）。DHI 是一套針對(duì)奶牛泌乳性能及乳成分的完整奶牛生產(chǎn)性能記錄體系。目前國際上奶牛隱性乳房炎乳SCC 的判斷標(biāo)準(zhǔn)多為20 萬/mL～50 萬/mL。馬裴裴等[3]研究表明，北京地區(qū)中國荷斯坦牛群宜以10 萬/mL～50 萬/mL SCC 作為隱性乳房炎標(biāo)準(zhǔn)。為克服SCC 在統(tǒng)計(jì)分析中的不足，通常將其轉(zhuǎn)換成遵循正態(tài)分布的體細(xì)胞評(píng)分（Somatic Cell Score，SCS）的形式[4]。

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

1 材料與方法

1.1 實(shí)驗(yàn)材料 來自北京市奶牛中心1998—2016 年的北京地區(qū)DHI 數(shù)據(jù)，共包括121 個(gè)牛場(chǎng)78 386 頭中國荷斯坦牛的1 967 310 條測(cè)定記錄。數(shù)據(jù)內(nèi)容包括個(gè)體號(hào)、胎次、產(chǎn)犢日期、測(cè)定日期、產(chǎn)奶量、乳脂量、乳脂率、乳蛋白量、乳蛋白率和SCC 等。

1.2 實(shí)驗(yàn)方法

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

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

1.2.2 奶牛乳房炎單因素分析及入選變量 使用北京地區(qū)DHI 數(shù)據(jù)，利用SAS 9.2 軟件對(duì)影響奶牛乳房炎預(yù)測(cè)的各個(gè)因素進(jìn)行單因素分析，對(duì)單個(gè)因素進(jìn)行初步篩選，探究單個(gè)因素對(duì)奶牛乳房炎的影響大小，得到對(duì)奶牛乳房炎預(yù)測(cè)的影響有統(tǒng)計(jì)學(xué)差異的指標(biāo)。對(duì)于分類變量的分析，通過計(jì)算對(duì)數(shù)比率即Logit（p）確定自變量進(jìn)入模型的形式。

1.2.3 奶牛乳房炎多因素分析及Logistic 回歸模型建立根據(jù)北京地區(qū)DHI 數(shù)據(jù)單因素分析的結(jié)果，確定健康組、隱性乳房炎組和臨床乳房炎組組間有統(tǒng)計(jì)學(xué)意義的指標(biāo)，即奶牛乳房炎風(fēng)險(xiǎn)因素（表2）。

表2 北京地區(qū)中國荷斯坦牛乳房炎風(fēng)險(xiǎn)因素

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

1.2.4 ROC 曲線的繪制 ROC 曲線主要用于反映預(yù)測(cè)變量對(duì)響應(yīng)變量的預(yù)測(cè)準(zhǔn)確率，是一條由點(diǎn)連接而成的曲線，這里體現(xiàn)奶牛乳房炎風(fēng)險(xiǎn)評(píng)估模型的預(yù)測(cè)準(zhǔn)確率。其中Y 軸為靈敏度，X 軸為特異度。靈敏度是指一項(xiàng)診斷試驗(yàn)?zāi)軐?shí)際患病的奶牛正確地診斷為病牛的概率，特異度是指一項(xiàng)診斷試驗(yàn)?zāi)軐?shí)際無病的奶牛正確診斷為非病牛的概率。靈敏度越高，說明真正的病牛被診斷出來的可能性越大；特異度越高，說明真正的非病牛被排除的可能性越大。曲線下面積（Area Under the Curve，AUC）即模型診斷價(jià)值，一般用來評(píng)價(jià)試驗(yàn)的診斷性能，面積越大，說明該試驗(yàn)的診斷性能越好。一般來說，AUC 在0.50～0.70 之間被認(rèn)為其辨別力一般，AUC 在0.70～0.90 之間被認(rèn)為該模型良好，AUC 高于0.90 被認(rèn)為該模型是優(yōu)秀的。

1.2.5 Logistic 回歸模型預(yù)測(cè)準(zhǔn)確率的評(píng)估 針對(duì)構(gòu)建的Logistic 風(fēng)險(xiǎn)評(píng)估模型，應(yīng)用北京地區(qū)某荷斯坦牛場(chǎng)2019 年全年共11 338 頭次奶牛的DHI 數(shù)據(jù)，對(duì)Logistic風(fēng)險(xiǎn)評(píng)估模型預(yù)測(cè)奶牛乳房炎的準(zhǔn)確率進(jìn)行評(píng)估。其中，DHI 數(shù)據(jù)內(nèi)容包括個(gè)體號(hào)、胎次、產(chǎn)犢日期、測(cè)定日期、產(chǎn)奶相關(guān)記錄和SCC 等。利用R 3.4.1 軟件對(duì)DHI 數(shù)據(jù)進(jìn)行整理，并應(yīng)用Logistic 回歸模型對(duì)牛場(chǎng)DHI 數(shù)據(jù)中各測(cè)定指標(biāo)進(jìn)行分析，比較每月預(yù)測(cè)情況與下月奶牛乳房炎實(shí)際情況之間的一致性以驗(yàn)證模型預(yù)測(cè)結(jié)果的準(zhǔn)確率。奶牛乳房炎發(fā)病判定標(biāo)準(zhǔn)為SCC＞10 萬/mL。其中，10 萬/mL～50 萬/mL SCC 作為隱性乳房炎標(biāo)準(zhǔn)，大于50 萬/mL SCC 作為臨床乳房炎標(biāo)準(zhǔn)。

2 結(jié)果

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

2.2 奶牛隱性乳房炎組及臨床乳房炎組各變量解釋 使用北京奶牛中心提供的DHI 數(shù)據(jù)進(jìn)行奶牛乳房炎風(fēng)險(xiǎn)評(píng)估。選用“場(chǎng)規(guī)模+胎次+測(cè)定季節(jié)+泌乳階段”為變量代入多因素Logistic 回歸方程。最終可用于診斷奶牛下一個(gè)測(cè)定月隱性乳房炎患病風(fēng)險(xiǎn)的解釋變量有4個(gè)，具體包括場(chǎng)規(guī)模、胎次、測(cè)定季節(jié)、泌乳階段（表3）。

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

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

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

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

表3 隱性乳房炎組風(fēng)險(xiǎn)因素

表4 臨床乳房炎組風(fēng)險(xiǎn)因素

圖2 隱性乳房炎組及臨床乳房炎組預(yù)測(cè) ROC 曲線

2.4 奶牛乳房炎風(fēng)險(xiǎn)評(píng)估模型的預(yù)測(cè)準(zhǔn)確率 利用Logistic回歸模型對(duì)北京地區(qū)某荷斯坦牛場(chǎng)2019 年全年共11 338頭次奶牛的乳房炎發(fā)病風(fēng)險(xiǎn)進(jìn)行評(píng)估，根據(jù)評(píng)估結(jié)果進(jìn)行分析（表5）發(fā)現(xiàn)，Logistic 回歸模型對(duì)該牛場(chǎng)奶牛的乳房炎風(fēng)險(xiǎn)評(píng)估的準(zhǔn)確率平均為70.84%，具備應(yīng)用前景。

3 討 論

奶牛場(chǎng)中使用各種控制和監(jiān)測(cè)奶牛乳腺健康的方法主要是為了減少新發(fā)病奶牛的數(shù)量，減少已有的患乳房炎奶牛，并減少奶牛患病的持續(xù)時(shí)間。奶牛乳房炎風(fēng)險(xiǎn)因素的分析能夠幫助鑒定奶牛乳房健康情況及發(fā)病風(fēng)險(xiǎn)，對(duì)于牛場(chǎng)中奶牛乳房炎的控制及減少新感染的牛只具有重要作用[12-13]。

表5 北京地區(qū)某荷斯坦牛場(chǎng)2019 年全年奶牛乳房炎風(fēng)險(xiǎn)預(yù)測(cè)結(jié)果準(zhǔn)確率

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

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

4 結(jié) 論

本研究利用DHI 記錄分析荷斯坦牛乳房炎發(fā)病的風(fēng)險(xiǎn)因素，通過對(duì)覆蓋120 余個(gè)牛場(chǎng)近20 年的196 萬多條DHI 數(shù)據(jù)記錄進(jìn)行分析，基于多因素Logistic 回歸分析的結(jié)果表明，場(chǎng)規(guī)模越小、奶牛胎次越高、夏季階段測(cè)定、泌乳階段越高，奶?；既榉垦罪L(fēng)險(xiǎn)越高；使用DHI 測(cè)定記錄進(jìn)行Logistic 回歸模型的構(gòu)建，所得奶牛隱性乳房炎風(fēng)險(xiǎn)評(píng)估模型的預(yù)測(cè)價(jià)值及準(zhǔn)確率均低于臨床乳房炎風(fēng)險(xiǎn)評(píng)估模型；Logistic 回歸分析后得到的奶牛隱性乳房炎風(fēng)險(xiǎn)評(píng)估模型和奶牛臨床乳房炎風(fēng)險(xiǎn)評(píng)估模型預(yù)測(cè)價(jià)值分別為0.721 和0.825，模型預(yù)測(cè)價(jià)值良好，可以應(yīng)用于實(shí)際荷斯坦牛群的乳房炎風(fēng)險(xiǎn)評(píng)估。