中国荷斯坦牛乳房炎风险评估模型建立及预测分析

史良玉，李文龙，唐永杰，米思远，肖 炜，刘 林，张 毅，王雅春，俞 英*

（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，模型预测价值良好，可以应用于实际荷斯坦牛群的乳房炎风险评估。