廣靈驢體重與體尺指標的主成分分析和回歸模型的建立

牛曉艷，曹 亮，詹海杰，馮國亮，李燕平，明世清，唐 貴，梁茂文，任克良*

（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。