广灵驴体重与体尺指标的主成分分析和回归模型的建立

牛晓艳，曹 亮，詹海杰，冯国亮，李燕平，明世清，唐 贵，梁茂文，任克良*

（1.山西农业大学动物科学学院，山西太谷 030801；2.山西省畜禽繁育工作站，山西太原 030001；3.广灵县优种驴场，山西广灵 037500）

作为一种经济型草食动物，驴具有耐粗饲、抗逆性强、适应性强、易于饲养且饲料来源广的特点。与牛羊肉相比，驴肉脂肪含量明显较低，且胆固醇含量低，不饱和脂肪酸含量丰富，可预防心脑血管疾病。驴皮是熬制阿胶的主要原料，具有补血、滋阴、养肝、益气、润肺等功效。养驴业是传统畜牧业的一个重要组成部分，随着社会的发展，驴的役用价值逐渐被其肉用、皮用、奶用等经济价值所替代。驴存栏量以每年3.5%的速度下降[1]，地方品种更是面临种群数量缩减、种质退化甚至濒临灭绝的危险境地，亟需进行保护性研究。

广灵驴主要分布在山西省广灵、灵邱两县及其周围各县的边缘地区[2]，以广灵县南村镇、壶泉镇、加斗乡分布最多，质量最好，其品种的形成具有二百多年的历史[3]。该地区属于黄土高原东缘、太行山北部，全年无霜期120 d，年降水量400 mm 左右。广灵驴俗称广灵画眉驴，具有体格高大、肌肉丰满、结构匀称、“五白一黑”等特点，是改良小型驴种的优良品种之一，也是全国24 个优秀地方驴品种之一。1983 年，广灵驴被编入全国《马驴品种志》，2006 年农业部（现为农业农村部）第662 号公告确定广灵驴为国家级遗传资源保护品种[4]。根据2019 年底的统计数字，广灵驴共存栏1 960头，其中基础母驴1 020 头，公驴52 头，青年驴和驴驹888 头[3]。

目前，对于广灵驴品种特征、生产性能等方面的描述较少，所报道的数据集中于20 世纪八九十年代[2,5]，不能适应现代社会的发展要求，此外，对于其体重体尺之间主成分分析的研究尚未见报道。基于此，本研究对广灵县优种驴场现存栏的不同年龄阶段的种驴进行了体重体尺的相关分析和主成分分析，并建立了体重与体尺的最优回归模型，为种驴的测量评估和进一步的选育研究提供了方法。

1 材料与方法

1.1 实验动物 体重、体尺数据通过对山西省广灵县优种驴场2008—2018 年现存驴群进行测定和数据采集而获得。分析过程中将驴群按照年龄划分为6 个月～2 岁的幼驴阶段、3～4 岁的青年驴阶段、5 岁及以上的成年驴3 个阶段，其中幼驴29 头、青年驴20 头、成年驴51 头，共计100 头。

驴群的饲喂采用以谷草为基础粗饲料，秋季加喂青玉米秸秆（带穗），精料补充以玉米、黑豆、麸皮为主的传统饲养方法。

1.2 体重、体尺测定 分别测定不同年龄阶段驴群的体重、体高、体长、胸围、管围5 个指标。其中体重由地秤测定活重，读取仪表值。体高为鬐甲最高点到地面的垂直距离。体长为肩端到臀端的斜直线距离。胸围是指从肩胛骨后端引一垂线，绕体躯一周的周长。管围是指左前管上1/3 处至管骨最细之处的周长。同时记录采样个体的毛色、系谱、性别等指标。

1.3 统计分析 不同年龄阶段驴群体重、体尺的平均值、标准差、标准误等描述性指标采用SPSS 24.0 进行计算。

主成分分析的基本思想是把原有的多个指标转化成少数几个代表性较好的综合指标，这些少数指标能够反映原来大部分信息（85%以上），并且各个指标间保持独立，可起到降维和简化数据结构的作用[6]。采用R 语言中cor_pearson 函数分别对不同年龄组种驴5 个体重、体尺性状进行相关性分析，采用corr 函数对相关性分析结果进行显著性检验。

采用R 语言中princomp 函数进行主成分分析，采用summary 函数提取主成分的具体信息，采用loadings函数分析主成分的载荷内容，采用screeplot 函数绘制主成分分析的碎石图。

采用R 语言中lm 进行体高、体长、胸围、管围相对于体重的多元回归方程构建；采用step 函数进行方程的逐步回归，并参考AIC（赤池信息量）准则和update函数进行模型的优化，最终得到最优的拟合模型。

2 结果与分析

2.1 年龄组和性别对广灵驴体重体尺的影响 广灵驴不同年龄组和性别对其体重体尺的影响见表1～2。由表1 可知，不同年龄阶段对广灵驴体重（F=36.67，P＜0.001）、体高（F=37.85，P＜0.001）、体长（F=37.56，P＜0.001）、胸围（F=50.41，P＜0.001）、管围（F=32.31，P＜0.001）均具有极显著影响。几个体重体尺指标都表现出幼年组极显著低于青年组和成年组，但青年组和成年组间差异不显著，表明年龄对驴的生长发育性状有极显著的效应，但4 岁之后，驴的生长趋于稳定。同时也说明在广灵驴快速生长阶段（4 岁以前）应注意饲料中的营养补充，使其快速达到成年体重。

表1 不同年龄组对广灵驴体重体尺的最小二乘分析结果（最小二乘均值±标准误）

表2 不同性别对广灵驴体重体尺的最小二乘分析结果（最小二乘均值±标准误）

2.2 广灵驴体重和体尺间的相关性分析 由表3 可知，5 周岁及以上的成年驴体重与体高、体长、胸围、管围均表现出极显著的正相关关系，其中体重与体长的相关性最高（γ=0.648）。各指标间均呈现出极显著的正相关关系，其中体长和体高的正相关关系最高（γ=0.914），胸围与管围的正相关关系最低（γ=0.368）。

表3 广灵驴成年驴体重和体尺的相关性分析

由表4 可知，3～4 岁的青年驴体重与体高（P＜0.01）、体长（P＜0.05）、胸围（P＜0.05）之间存在显著或极显著正相关关系，但体重与管围之间的相关关系不显著。青年驴的体长和体高相关性最高（γ=0.908），体重和胸围之间的相关性最低（γ=0.497），胸围与体高、体长以及体重与管围，胸围与管围之间均无显著的相关关系。

表4 广灵驴青年驴体重和体尺的相关性分析

由表5 可知，6 个月～2 岁幼驴体重与各体尺指标间均存在极显著的正相关关系，且相关系数很高。其中，体重与体高的相关性最高（γ=0.945），体重与管围的相关性最低（γ=0.817）。各指标间的相关关系中，体长和体高的相关性最高（γ=0.970），胸围与体重（γ=0.908）、胸围与体高（γ=0.947）、胸围与体长（γ=0.931）的正相关关系均较高。从上述广灵驴不同生长发育阶段体重与各体尺指标的相关分析结果来看，不同阶段的相关性关系存在差异。

表5 广灵驴幼驴体重和体尺的相关性分析

2.3 广灵驴5 个性状的主成分分析 由表6～7 可知，6 个月～2 岁阶段广灵驴的特征值为4.61，贡献率为0.92。第一主成分标准变化量所表达的关系式为：Y1=-0.45BW–0.46BH–0.46BL–0.45CM–0.42CC。由于贡献率代表每个特征值所能够解释方差的比例，特征向量代表变量预测某一主成分的回归系数。从贡献率和累积贡献率来看，协方差矩阵第一个特征根的贡献率达到92%，说明第一个主成分代表了原来5 个因素90%以上的信息，所以第一主成分又称为体重因子。3～4 岁青年驴阶段广灵驴的前3 个特征值分别为3.12、0.75、0.64，前3 个主成分的贡献率分别为0.62、0.15、0.13，前3 个主成分的累积贡献率为0.90（＞0.80），所以分别取体重、体高和体长3 个主成分代表原来5 个因素的信息，第一主成分的关系式为：Y1=-0.42BW–0.52BH–0.49BL–0.38CM–0.41CC；第二主成分的关系式为：Y2=0.32BW–0.37BH–0.43BL+0.76CM–0.04CC；第 三主成分的关系式为：Y3=-0.59BW–0.15BH–0.11BL+0.16CM+0.77CC。5 岁以上成年驴阶段广灵驴的前2 个特征值分别为3.33 和0.76，前2 个主成分的贡献率分别为0.67 和0.15，累积贡献率为0.82（＞0.80），所以分别取体重和体高代表原来5 个因素的信息。第一主成分的关系式为：Y1=-0.46BW–0.48BH–0.50BL–0.36CW–0.43CC；第二主成分的关系式为：Y2=-0.28BW+0.36BH+0.35BL–0.80CW–0.43CC。

表6 广灵驴不同生长发育阶段相关矩阵的特征值、贡献率和累计贡献率

图1 为广灵驴不同年龄阶段5 个体重体尺性状的碎石图，反映了不同年龄阶段广灵驴不同因子所涵盖的信息量。从图1 可以直观地看出，对于6 个月～2 岁的幼驴，第一主成分涵盖了绝大多数的变异，代表了90% 以上的信息；对于3～4 岁的青年驴，前3 个主成分涵盖了90%以上的变异；而对于5 岁以上的成年驴，前2 个主成分代表了80%以上的变异。

图1 不同月龄广灵驴的碎石图

表7 广灵驴不同生长发育阶段入选主成分的特征向量

2.4 广灵驴体重和体尺间最优回归模型的建立 应用R语言中lm 函数和step 函数建立广灵驴体尺相对于体重的多元线性回归方程，结果列于表8。

根据AIC 准则，从一组可供选择的模型中选择最佳模型时以AIC 最小时为最佳[6]。由表8 可知，在拟合的4 个多元线性方程中，Y2和Y4的AIC 含量较小，但Y2中所采用的2 个体尺性状（体高、体长）的回归系数检验不显著，故该方程不成立，而Y4所采用的2个体尺性状（体长、胸围）回归系数检验均极显著，说明体长和胸围对广灵驴的体重具有决定性作用，拟合的方程可用于预测广灵驴体重，最终获得的拟合方程为：BW=4.35BL+2.18CM－579.39。进一步的，随机选取30 头不同年龄阶段的广灵驴体重体尺数据分别代入上述公式和经验公式：体重=+25 进行检验，结果如表9 所示。

表8 广灵驴体重相对于不同体尺性状的多元线性回归方程

由表9 可知，通过对本研究得出的拟合方程和经验公式拟合度（R2）的比较，发现2 个公式的拟合度都随着广灵驴年龄的增加而提高，对于成年驴的估计拟合度均可达到0.99 以上，较为理想。但也不难看出，对于成年驴阶段和幼驴阶段，拟合方程的拟合度均比经验公式高，说明拟合方程的效果更好；对于青年驴阶段，2种公式的拟合度接近。另外，由表9 还可以看出，拟合公式拟合效果随着个体体重的增加而提高，对青年和成年驴阶段个体的拟合效果普遍好于幼年阶段。

表9 采用拟合方程与经验公式对广灵驴体重的估计

3 讨 论

3.1 广灵驴体重、体尺分析 广灵驴属大型驴种，据《山西省家畜家禽品种志》记载，广灵驴的体型外貌特征为：体型接近正方形、体格高大、骨骼粗壮、结构匀称、耐寒耐粗饲。全身被毛短而粗密，毛色以黑化眉为主，青化眉、灰色、纯黑次之。文献报道广灵驴成年公驴平均体高为138.4 cm，体长138.5 cm，胸围为147.2 cm，管围为17.8 cm；成年母驴平均体高为134.1 cm，体长为131.6 cm，胸围为146.9 cm，管围为15.7 cm[7]。与本研究所获得数据相比，体高和体长数据一致，但胸围和管围较小，这可能反映了不同饲养环境下广灵驴生产性能的差异，也和不同种群的选择强度有关。

与其他大型驴种相比，广灵驴体型较晋南驴大，较德州驴小。据聂鸿瑶等[2]对1979—1980 年在山西省闻喜县、夏县测定的晋南驴的体重体尺数据，成年公驴体重为127.3 cm，母驴为124.3 cm；体长公驴为123.2 cm，母驴为125.2 cm；胸围公驴为136.3 cm，母驴为135.6 cm；管围公驴为15.2 cm，母驴为14.8 cm。而冯志华[8]对100 头晋南驴（其中公驴28 头，母驴72 头）体重体尺进行测定，得到公驴平均体高为135.5 cm，母驴为134.5 cm；体长公驴为132.3 cm，母驴为128.4 cm；胸围公驴为132.3 cm，母驴为128.4 cm；管围公驴为17.2 cm，母驴为15.9 cm，均略小于广灵驴，但差异不大。农博[9]描述了德州驴的品种特征，成年公驴体高为142 cm，成年母驴体高为140 cm；体长公驴可达143.6 cm，母驴达137 cm；胸围公驴为152.8 cm，母驴为160 cm；管围公驴为18.7 cm，母驴为16.4 cm，均高于广灵驴或与之相当。

不同的年龄组因素对广灵驴体重体尺性状均具有极显著影响，表现为青年组和成年组的个体极显著高于幼年组个体，但青年组和成年组之间差异不显著。这与杨莉等[10]研究结果一致，他们采用3 种曲线对德州驴的生长发育规律进行拟合分析，得到公、母驴达到成熟时的拐点月龄分别为2.43 和1.95，说明在2 岁以后，驴的生长速度趋于平稳，接近成年水平。不同性别因素对广灵驴体重、体高、体长、胸围均无显著影响，但公驴管围极显著大于母驴，这与其他一些大型驴品种相似。

3.2 广灵驴体重、体尺的相关性分析 本研究结果表明，不同发育阶段的广灵驴体重体尺之间均存在显著的正相关，且成年驴阶段体重与体长相关性最高（γ=0.648，P＜0.001），而青年和幼年阶段体重与体高的相关性最高（γ=0.606，P＜0.001；γ=0.945，P＜0.001）。这与杨莉等[10]报道的公驴和母驴体重与胸围的相关性最高不一致（公驴：γ=0.963；母驴：γ=0.971）。Pierre 等[11]计算了西非布基纳法索等4 个地区1 001 头公驴和351 头未怀孕母驴体重和体尺的相关关系，也得到体重与胸围具有最强的相关性（公驴：γ=0.89；母驴：γ=0.90）。但陆会宁等[12]在对80 只兰州大尾羊× 小尾寒羊F1代体重和体尺的相关性研究中发现，体重和体长（γ=0.875）的相关性与体重和胸围的相关性接近（γ=0.907）。说明体重和各体尺指标间的相关性与物种和样本量有关。本研究与其他研究相比样本量较少，可能会对不同指标间的相关关系计算结果产生影响。

3.3 广灵驴体重、体尺的主成分分析 主成分分析法可以在力求数据信息丢失最少的情况下，对高维的变量空间进行降维，研究指标体系的少数几个线性组合，而且这几个线性组合所构成的综合指标能够尽可能的保留原来指标的信息。目前，国内尚未见对地方驴品种体重体尺性状进行主成分分析的报道，通过本研究结果可以看出，不同年龄阶段对广灵驴生长发育起到主要作用的性状不同。在幼年阶段，体重的影响最大；青年阶段，体重、体高和体长则代表了其主要体型特征；成年阶段，体重和体高最具代表性。这说明在广灵驴的选育过程中，应根据不同阶段性状之间的关联性，并兼顾其他体尺性状，对育种目标进行优化，保留地方品种的体型外貌特征。

3.4 广灵驴体重、体尺的回归分析 采用多元回归分析中的逐步回归分析法结合AIC 准则得出广灵驴体重与体尺性状的最优回归线性模型，发现广灵驴的体重主要与其体长和胸围密切相关。杨阳等[13]对3～9 岁178 头德州驴进行了体重和体尺的回归分析，建立了德州驴体重（Y）与其体长（X1）、胸围（X2）和尻长（X3）的多元线性回归方程：Y=2.6X1+1.8X2+3.0X3–500.2，与本研究结果相似，但本研究未对广灵驴的尻长数据进行采集，所以最终的方程仅包括了体长和胸围2 个性状。肖海霞等[14]采用R 软件包对62 头吐鲁番驴进行体重和体尺分析，并建立了最优的回归模型，得出回归方程为：Y（体重）=0.71X1（体高）+1.54X2（体长）+2.84X3（胸围）–447.36，与本研究结果相似。

在实际生产中，对驴的体重测定通常较为困难，一些参考资料给出了如下计算公式：体重=+25，此公式的运算较为繁琐难记[7]，本研究给出的体重估计公式将体重作为体长和胸围的线性函数进行拟合，校正拟合度达到0.85 以上，更加简单实用。比较经验公式与线性拟合方程的拟合度可知，线性拟合方程的拟合度在广灵驴生长发育的各个阶段均高于经验公式或与之接近，这说明线性方程的拟合和预测更加准确，可代替传统经验公式在生产中进行推广使用。

本研究提出的拟合公式是通过对广灵驴体重体尺性状进行分析得出的，通过与其他地方品种比较可知，不同品种具有不同的体型特征和生长发育特点，因此此公式仅适用于广灵驴体重的预测。另外，虽然此公式是通过对6 个月～5 岁（幼年～成年）不同生长阶段个体数据的拟合而得到，但也不难发现，公式估计值的拟合度随着个体年龄的降低而迅速下降，所以本公式用于估计青年和成年驴的体重更加准确适用。最后，与其他研究者的估计方程相比，本研究未考虑尻长、十字部高等数据，在表现形式上略有差异，今后，可进一步对广灵驴的体尺数据进行收集，并扩大采样规模，得到更为完善的公式。

4 结 论

本研究基于R 语言对不同年龄阶段广灵驴的体重体尺进行了比较，计算了各体重体尺性状间的相关性，并对各性状进行了主成分分析，最后得出了体重与体尺间的最优回归线性方程。研究结果表明：①幼年和青年阶段体重与体高的相关性最高（γ=0.945，γ=0.606），成年阶段体重和体长的相关性最高（γ=0.648）；②主成分分析结果表明，幼驴阶段第一主成分的特征值为4.61，贡献率为0.92；青年驴阶段第一、第二、第三主成分的特征值分别为3.12、0.75、0.64，累积贡献率达到0.90；成年驴阶段第一、第二主成分的特征值分别为3.33 和0.76，累积贡献率为0.82；③通过逐步回归分析得到广灵驴体重和体尺的最优回归方程为：BW=4.35BL+2.18CM–579.39。