基于因子分析和聚类分析的30份山黄皮种质资源评价

张涛 贺鹏 宋海云 许鹏 莫庆道 覃振师 韦媛荣 汤秀华 王文林 郑树芳

摘  要：选取单果重、果实纵径、果实横径、果形指数、单种重、种子纵径、种子横径、可食率、可溶性固形物含量9个性状指标对30份山黄皮种质资源进行评价。结果表明，山黄皮果实的单果重变异幅度最大，果形指数的变异幅度最小;除可溶性固形物含量与其他8个性状相关性不显著外，其余各指标之间大都存在显著或极显著的相关性;采用因子分析和聚类分析将30份山黄皮种质聚为3大类，山黄皮种质资源的果实性状与地理起源有一定的相关性。其中，来自广西龙州县和宁明县的13份山黄皮种质的果实性状表现最好。

关键词：山黄皮;种质资源;果实性状;因子分析;聚类分析中图分类号：S663. 9      文献标识码：A

Evaluation of 30 Species of Clausena anisum-olens （Blanco） Merr. Germplasm Resources Based on Factor Analysis and Cluster Analysis

ZHANG Tao， HE Peng， SONG Haiyun， XU Peng， MO Qingdao， QIN Zhenshi， WEI Yuanrong，TANG Xiuhua， WANG Wenlin， ZHENG Shufang^*

South Subtropic Agricultural Scientific Research Institute of Guangxi， Longzhou， Guanxi 532415， China

Abstract： In order to evaluate the germplasm resources of 30 species ofClausena anisum-olens（Blanco） Merr.， nine fruit traits including weight of single fruit， fruit equatorial diameter， fruit transverse diameter， fruit shape index， weight of single seed， seed equatorial diameter， seed transverse diameter， edible rate， soluble solid content were measured and analyzed by statistical analysis. The variation of weight of single fruit ofC. anisum-olens（Blanco） Merr. was the most significant， while the variation of fruit shape index was the most non-significant. The soluble solid content was not correlated with the other eight indices significantly， and there was significant and extremely significant correlation between the other indexes. Thirty species ofC. anisum-olens （Blanco） Merr. germplasm resources were clustered into three categories by factor analysis and cluster analysis. There was a certain correlation between fruit trait and geographical origin. Fruit traits of 13 germplasm from Longzhou County and Ningming County were the best.

Keywords： Clausena anisum-olens（Blanco） Merr.; germplasm resources; fruit trait; factor analysis; cluster analysis

DOI： 10.3969/j.issn.1000-2561.2020.07.007

山黄皮[Clausena anisum-olens（Blanco） Merr.]为芸香科（Rutaceae）黄皮属（Clausena）多年生常绿大灌木或小乔木，俗称“鸡皮果”，又名细叶黄皮，主要分布于北回归线以南的南亚热带地区，主产于我国广西西南部的百色和龙州、云南南部的蒙自和河口、广东的新会和鹤山及越南北部和菲律宾^[1-2]。山黄皮果实为浆果，果皮、果肉可食，果实中糖类、蛋白质、维生素等含量丰富，还含有铁、钙等各种营养元素，营养价值高。山黄皮以其可口的酸甜味及香气在芸香科水果中独具一格，可加工成山黄皮干、果脯、果酱、饮料等制品，其产品消腻开胃，风味独特，具有化痰止咳、祛风除湿的功效，深受消费者的欢迎^[2-3]。

目前对山黄皮的研究主要集中在种质资源收集评价及其遗传多样性分析方面。覃振师等^[4-5]从广西崇左市5个县（区）收集引进了50份野生山黄皮种质资源，对单穗挂果数、平均单果重、可食率、产量、外观品质、成熟期等方面的因素简单综合分析后，初步选出6个优良单株，后对收集的山黄皮种质资源进行ISSR遗传多样性分析，表明大部分山黄皮种质亲缘关系较接近，且与地理起源有一定的相关性，部分遗传差异较大的种质可作为山黄皮优良品种和单株选育及良种繁育等研究的重要材料。而多元统计分析中的主成分分析、因子分析和聚类分析方法早已经在苹果^[6-8]、大蒜^[9]、冬枣^[10]、砂梨^[11]、核桃^[12]、桃^[13]、歐李^[14]、香梨^[15]、葡萄^[16]、猕猴桃^[17-18]、西瓜^[19]、茄子^[20]、普洱茶^[21]等作物品质性状分析和综合评价方面大量应用，但应用多元统计分析方法对山黄皮种质果实性状进行评价还未见报道。

本研究以本所资源圃内从广西崇左市各地引进的30份山黄皮种质资源为研究对象，在前人研究的基础上，对其果实性状进行观察和测定，应用因子分析对各种质山黄皮的果实性状进行综合评价，在综合评价的基础上再进行聚类分析，以期为本地区山黄皮种质的资源保存和利用、品种品质改良以及优质新品种选育提供一定的参考依据。

1  材料与方法

1.1 材料

供试材料选自广西南亚热带农业科学研究所山黄皮种质资源圃的30份山黄皮种质资源，种质均采用嫁接保存，树龄3～5年。资源编号：SHP1～SHP2、SHP4～SHP15、SHP17～SHP18、SHP21～SHP22、SHP26、SHP28～SHP32、SHP34、SHP36～SHP40。各种质生长期间管理水平一致。试验在自然开花和自然成熟时采收选择能代表种质特性、成熟度一致的果实。供试山黄皮种质的采集地信息如表1所示。

仪器与设备：电子数显（快闪）卡尺，桂林量具刃具有限责任公司;WYT-IO-80%手持折射仪，成都华西光学电子仪器厂;BP211D电子分析天平，广州广一科学仪器有限公司。

1.2 方法

参考相关水果的果实性状研究的测定指标^[19，22-26]，以及《热带、南亚热带果树种质资源描述规范》和《热带、南亚热带果树种质资源数据质量控制规范》^[27-28]，选取9个性状指标进行测定;每个种质随机选取3棵树，于果实成熟期进行采样，分别采摘树冠外围中上部的果实，每个种质选取30～40个，样品果存入4 ℃冰箱中，之后每次随机选取10个样品果进行各指标的测定，重复测试3次，记录测定数据。各指标测定方法如下：

单果重（g）：用电子分析天平称10个样品果的总重量，计算平均单果重。

果实纵径（mm）：用游标卡尺分别测量10个样品果纵向最长处的直径，计算平均果实纵径。

果实横径（mm）：用游标卡尺分别测量10个样品果横向最宽处的直径，计算平均果实横径。

果形指数：果实纵径与果实横径的比值即果形指数。

单种重（g）：用电子分析天平称20粒种子的总重量，计算平均单种重。

种子长（mm）：用游标卡尺分别测量10粒种子纵向最长处的直径，计算平均种子长。

种子宽（mm）：用游标卡尺分别测量10粒种子横向最宽处的直径，计算平均种子宽。

可食率：取样品果10个，洗净，拭干，分别称量果重（W）和种子重量（W₁），计算平均可食率：可食率=（W-W₁）/W×100%

可溶性固形物含量（%）^[29]：用WYT-IO-80%手持折射仪测量果汁的可溶性固形物含量。

1.3数据处理

采用Excel 2007对9个果实性状指标数据进行多样性分析和相关性分析;采用Excel 2007对这些数据进行标准化处理，然后采用SPSS 20.0软件进行因子分析和聚类分析。

2  结果与分析

2.1 山黄皮种质果实品质性状的变异情况分析

30份山黄皮种质果实的9个品质性状的变异情况分析结果如表2所示。30份山黄皮种质果实的9个品质性状间表现出不同程度的变异。其中，单果重的变异幅度最大，变异系数为23.45%;其次为单种重、可溶性固形物含量、种子宽、可食率、种子长、果实横径、果实纵径，变异系数在7%～19%;果形指数的变异系数最小，仅为4.73%。

2.2 山黄皮种质果实性状的相关性分析

山黄皮果实性状之间的相关性分析结果如表3所示，单果重与果实纵径、果实横径、种子长宽以及可食率呈极显著的正相关，与单种重呈显著的正相关，与果形指数呈显著的负相关;果实纵径与果实横径、单种重、种子长以及可食率呈极显著的正相关，与种子宽呈显著的正相关;果

注：^*表示相關性显著（P<0.05），^**表示相关性极显著（P<0.01）。

Note：^*means significant correlation （P<0.05），^**means extremely significant correlation （P<0.01）.

实横径与种子长宽、可食率呈极显著的正相关，与单种重呈显著的正相关，与果形指数呈显著的负相关;果形指数与可食率呈显著的负相关;单种重与种子长呈极显著的正相关，与种子宽呈显著的正相关;种子长与种子宽呈极显著的正相关;可溶性固形物与其他8个性状相关不显著。

2.3 因子分析

2.3.1  数据标准化处理  对30个种质的9个性状指标数据按照公式1进行无量纲标准化处理^[30-34]，公式如下：

（1）

式中，为标准化值;表示第j个样品（种质）第i个指标的原始测定值;表示第i个指标的平均值;表示第i个指标的标准差（SD值）。标准化后的数据如表4所示。

2.3.2  因子分析的适宜性检验  对经过无量纲标准化处理后的数据进行KMO检验和Bartlett球度检验^[33-38]，结果如表5所示。KMO值为0.561，根据统计学家Kaiser给出的KMO度量标准可知，原始变量适合进行因子分析。Bartlett球度检验得出的相伴概率为0.000，小于显著性水平0.05，因此拒绝Bartlett球度检验的零假设，认为变量间存在较强的相关性，适合做因子分析。综合以上检验结果，本研究的试验数据可以进行因子分析。

2.3.3  公因子的选择  对标准化后的数据进行因子分析，各个主因子的特征值、累积贡献率如表6所示。从表6中可以看出，前3个因子变量的特征值均大于1，并且经过最大方差旋转后它们的方差贡献率分别为42.052%、23.907%、15.887%，累计方差贡献率已达到81.846%，即这3个因子已

注：表中加括弧的数据表示为负值。

Note： Data in brackets was negative values.

经把山黄皮种质的全部9项果实性状指标的绝大部分信息反映出来，因此可以选取这3个因子作为山黄皮种质果实性状的综合评价指标。

为了更加突出主要公因子的作用，将因子进行方差最大正交旋转，旋转后各个性状指标在空间的分布如图1所示。

公因子所包含的指标因子载荷若为正，表明该因子对公因子产生正向影响;因子载荷若为负，表明该因子对公因子产生负向影响。图1可直观地解释各性状指标对3个公因子的贡献。第1公因子F₁主要由A1（单果重）、A2（果实纵径）、A3（果实横径）、A8（可食率）4个因子综合作用的指标，主要反映山黄皮果实的大小、性状以

及果实的食用品质;第2公因子F₂主要由A6（种子长）、A7（种子宽）2个因子综合作用的指标，主要反映山黄皮果实种子的大小;第3公因子F₃主要由A9（可溶性固形物含量）1个因子作用的指标，主要反映山黄皮果实可溶性糖类的含量。

2.3.4  山黄皮果实性状的公因子得分  旋转后的各公因子得分系数如表7所示，30份山黄皮种质果实性状各公因子的得分根据公式（2）进行计算;30份山黄皮种质果实性状指标的综合评价得分，根据3个公因子旋转后的方差贡献率及对应的公因子得分进行计算，详见公式（3）。

（2）

（3）

公式中為公因子得分，为综合评价得分，，…，为第n份种质果实样品的9个性状指标标准化后的值，，…，为旋转后9个果实性状指标因子得分系数。

2.3.5  山黄皮果实性状的分布状况  以各主因子得分绘制二维排序图，可直观地反映种质与性状指标之间关系及山黄皮果实性状的分布状况。以第1主因子得分做横坐标，第2主因子得分做纵坐标，绘制二维排序图（图2）：1、2、4、6、8、10、16、19、25、26（SHP1、SHP2、SHP5、SHP7、SHP9、SHP11、SHP18、SHP26、SHP34、SHP36）10个种质（系）坐落在公因子F₁得分的正向区间、公因子F₂得分的负向区间，表明这10个种质果实大、种子小、可食率高;12、15、23（SHP13、SHP17、SHP31）3个种质坐落在公因子F₁、公因子F₂得分的負向区间，表明这些种质果实的种子虽然小，但果实整体也偏小、可食率也低。其中15（SHP17）种质果实的种子相较于其他种质明显偏小。

编号对应的山黄皮种质名称见表1。

Number corresponding to the name ofClausena anisum-olens（Blanco） Merr. germplasm resources was shown in Table 1.

图2F₁主因子与F₂主因子得分二维排序

Fig. 2  Scatter plot based on principal factor F₁and F₂

以第1主因子得分做横坐标，第3主因子得分做纵坐标，绘制二维排序图（图3）：1、2、3、4、5、6、8、10、11、16、19、21、30（SHP1、SHP2、SHP4、SHP5、SHP6、SHP7、SHP9、SHP11、SHP12、SHP18、SHP26、SHP29、SHP40）13个种质坐落在公因子F₁得分的正向区间、公因子F₃得分的负向区间，表明这13个种质果实大、可食率高，可溶性糖类含量也高;12、20、23、28、29（SHP13、SHP28、SHP31、SHP38、SHP39）5个种质坐落

编号对应的山黄皮种质名称见表1。

Number corresponding to the name ofClausena anisum-olens（Blanco） Merr. germplasm resources was shown in Table 1.

图3F₁因子与F₃因子得分二维排序

Fig. 3  Scatter plot based on principal factor F₁and F₃

在公因子F₁、公因子F₃得分的负向区间，表明这些种质果实偏小、可食率低，但可溶性糖类含量高。

2.4 山黄皮种质的聚类分析

选用因子分析中筛选出对山黄皮果实性状影响较大的7个性状，采用欧氏距离类平均法对30份山黄皮种质进行系统聚类分析，聚类结果如图4所示，各类群山黄皮种质果实性状指标的平均值如表8所示。

聚类结果（图4）显示，以欧氏距离10为临界值，可将30份山黄皮资源分为3大类群。在欧氏距离10～15之间，编号为15的种质首先从30份种质中分离出来，其果实偏小，且可食率低，单独为一类，即第Ⅰ类群;剩余29份种质可分为2大类群，即有20份种质组成的第Ⅱ大类群，这一类果实质量平均为2.838g，果实大，可食率高，可溶性固形物含量也高;其他9份种质组成第Ⅲ大类群，这一类平均果粒质量为1.873 g，可溶性固形物含量低。在欧氏距离5～10之间，第Ⅱ大类群可分为2个亚类，第1个亚类由13份种质组成，包括编号分别为1、2、3、4、6、8、10、11、13、16、19、20、26，这一类果实平均质量为2.658 g，果实大，可食率较高，可溶性固形物含量高;第2个亚类有7份种质，包括编号分别为5、9、17、18、21、25、30，这一类表现为果实平均质量为

3.171 g，可食率高，可溶性固形物含量较高。同样，可将第Ⅲ大类群分为2个亚类，第1个亚类有3份种质，包括编号分别为7、14、22，这一类种质果实平均质量为2.287 g，果实较大，但可食率、可溶性固形物含量较低;剩余6份种质为第2亚类，包括编号分别为12、23、24、27、28、29，这一类种质果实平均质量为1.667 g，果实小，可食率、可溶性固形物含量较低。

3  讨论

本研究通过对所选取的单果重、单种重、果实横径、果实纵径、可溶性固形物含量等9个果实性状进行描述性统计发现，单果重的变异幅度最大，变异系数在20%以上，说明单果重这个性状具有丰富的遗传信息和选择潜力，在育种改良上还有很大空间;而果实横径、果实纵径及果形指数变异范围较小，变异系数均在10%以下，表明这些性状遗传较为稳定，在现有种质的基础上进行改良的空间较小。

覃振师等^[5]对从广西崇左市5个县（区）收集引进的50份野生山黄皮种质资源进行ISSR遗传多样性分析，将50份山黄皮种质分成A、B两大组，A组包含40份种质，可分2个亚组：第Ⅰ亚组包含28份种质，其中16份种质来自龙州县，7份种质来自宁明县，5份种质来自凭祥市;第Ⅱ亚组包含12份种质，其中6份种质来自江州区，4份种质来自扶绥县，2份种质来自宁明县;B组包含10份种质，可分2个亚组：第Ⅰ亚组包含8份种质，其中4份种质来自天等县，4份种质来自大新县;第Ⅱ亚组包含2份种质，均来自大新县;来自宁明的NM9、来自天等的TD1、来自大新的DX2和来自江州区的JZ5在形态和遗传上差异较大。

本研究从这50份中选取了30份，通过对9个果实性状的因子分析和聚类分析，将这30份山黄皮种质聚为3大类，聚类分析结果得出的各类群性状特征与各主因子得分绘制二维分布图分析结果基本一致。编号为15，原始编码为DX2的种质被归到第Ⅰ大类，该种质来自大新县硕龙镇，该类种质的果实偏小，且可食率低，与其他种质差异较大，这与覃振师等^[5]研究结果DX2在形态和遗传上差异较大相吻合。20份种质被归到第Ⅱ大类群，第Ⅱ大类群又可分为2个亚类，其中13份种质被归到第一个亚类，来自龙州县的有9份，来自宁明县的有4份，该类种质的果实大，可食率较高，可溶性固形物含量高;另7份种质被归到第二亚类，来自宁明县有4份，来自凭祥市有3份，该类种质的可食率高，可溶性固形物含量较高。第Ⅱ大类群与覃振师等^[5]研究结果A组的第Ⅰ亚组相吻合，第Ⅱ大类群分为2个亚类相当于对覃振师等^[5]研究结果A组的第Ⅰ亚组进行了再次分类。剩下的9份种质被归到第Ⅲ大类群，第Ⅲ大类群也分为2个亚类，其中3份种质被归到第一个亚类，来自崇左市江州区，该类种质的果实较大，但可食率、可溶性固形物含量较低;另6份种质被归到第二亚类，来自崇左市江州区的有3份，来自扶绥县渠旧镇的有3份，该类种质的果实小，可食率、可溶性固形物含量较低。第Ⅲ大类群与覃振师等^[5]研究结果A组的第Ⅱ亚组相吻合，第Ⅲ大类群分为2个亚类相当于对覃振师等^[5]研究结果A组的第Ⅱ亚组进行了再次分类。表明山黄皮种质资源的果实性状与地理起源有一定的相关性。

综合分析来看，第Ⅱ大类群的第一个亚类的13份种质，果实大，可食率较高，可溶性固形物含量高，果实性状综合表现最优，因子分析和聚类分析的两种方法的综合利用有助于了解各种质间的亲缘关系和优良种质的分布情况，为本地区山黄皮品种品质改良、优质新品种选育以及种质资源保护利用提供了科学、客观的依据。

4  结论

通过因子分析和聚类分析将30份山黄皮种质聚为3大类，山黄皮种质资源的果实性状与地理起源有一定的相关性。来自龙州县和宁明县的13份种质果实大，可食率较高，可溶性固形物含量高，果实性状表现最优。

果实品质的组成因子较多，由外在和内在的复合因子构成，本研究主要集中在山黄皮外在品质指标的分析，下一步将增加对内在品质指标及内外在品质指标间的相关性分析，探索合理地简化果实品质评价指标，构建科学的山黄皮果实综合评价体系。

参考文献

• 江苏新医学院. 中药大词典[M]. 上海： 上海科学技术出版社， 1985： 196.
• 中国科学院中国植物志编辑委员会. 中国植物志： 43卷2册[M]. 北京： 科学出版社， 2004： 138.
• 广西中医药研究所. 广西药用植物名录[M]. 南宁： 广西人民出版社， 1986： 308.
• 覃振师， 韦持章， 何铣扬， 等. 广西崇左市野生山黄皮种质资源调查与开发利用[J]. 广东农业科学， 2012， 39（5）： 138-139.
• 覃振师， 邓立宝， 王文林， 等. 广西崇左市山黄皮种质资源调查及ISSR遗传多样性分析[J]. 南方农业学报， 2016， 47（7）： 1071-1076.
• 公丽艳， 孟宪军， 刘乃侨， 等. 基于主成分与聚类分析的苹果加工品质评价[J]. 农业工程学报， 2014， 30（13）： 276-285.
• 孟宪军， 公丽艳， 毕金峰， 等. 中早熟苹果果实品质因子的选择[J]. 食品科学， 2013， 34（23）： 72-76.
• 冯  娟， 任小林， 田建文， 等. 不同产地富士苹果品质分析与比较[J]. 食品工业科技， 2013， 34（14）： 108-112.
• 陈书霞， 周  静， 申晓青， 等. 大蒜种质产量和品质性状主成分聚类分析与综合评价[J]. 植物遗传资源学报， 2012， 13（3）： 429-434.
• 马庆华， 李永红， 梁丽松， 等. 冬枣优良单株果实品质的因子分析与综合评价[J]. 中国农业科学， 2010， 43（12）： 2491-2499.
• 董星光， 田路明， 曹玉芬， 等. 我国南方砂梨主产区主栽品种果实品质因子分析及综合评价[J]. 果树学报， 2014， 31（5）： 815-822.
• 潘学军， 张文娥， 李琴琴， 等. 核桃感官和营养品质的主成分及聚类分析[J]. 食品科学， 2013， 34（8）： 195-198.
• 徐臣善， 高东升. 基于主成分分析的设施桃果实品质综合评价[J]. 食品工业科技， 2014， 35（23）： 84-88， 94.
• 李卫东， 李  欧， 和银霞， 等. 基于TXRF法的欧李果肉中营养元素特征分析[J]. 食品科学， 2015， 36（4）： 164-167.
• 李龙飞， 林彩霞， 吐尔逊阿依·达吾提， 等. 库尔勒香梨杂交品种（系）果实品质测定与综合评价[J]. 新疆农业大学学报， 2014， 37（2）： 153-158.
• 法洁琼， 张振文. 甘肃河西走廊产区主栽酿酒葡萄品质比较研究[J]. 北方园艺， 2013， 37（3）： 26-30.
• 秦红艳， 许培磊， 艾  军， 等. 软枣猕猴桃种质资源果实品质、表型性状多样性及主成分分析[J]. 中国农学通报， 2015， 31（1）： 160-165.
• 曲雪艳， 郎彬彬， 钟  敏， 等. 野生毛花猕猴桃果实品质主成分分析及综合评价[J]. 中国农学通报， 2016， 32（1）： 92-96.
• 赵卫星， 李晓慧， 徐小利， 等. 西瓜品种果实性状的变异性及主成分和聚类分析[J]. 河南农业科学， 2017， 46（1）： 105-109.
• Arivalagan M， Gangopadhyay K K， Gunjeet K， et al. Variability in mineral composition of India eggplant （Solanum melongena L.） genotypes[J]. Journal of Food Composition and Analysis， 2012（26）： 173-176.
• 刘兴勇， 师  江， 邵金良， 等. 普洱茶晒青样游离氨基酸和色泽差异研究[J]. 食品科学， 2015， 36（1）： 46-50.
• 张振境， 蔡小林， 潘介春， 等. 4种番荔枝树形树体生长和果实性状比较[J]. 安徽农业科学， 2017， 45（27）： 55-56， 67.
• 钟金仙， 罗  英， 曾仁杰， 等. 黄瓜种质资源果实性状的主成分分析与综合评价[J]. 中国农学通报， 2017， 33（16）： 46-52.
• 解德宏， 张翠仙， 张永超， 等. 澜沧江流域芒果种质资源果实性状多样性分析[J].南方农业学报， 2018， 49（2）： 214-221.
• 唐  克， 单金友， 吴雨蹊， 等. 引进的沙棘品种实生苗和果实性状研究[J]. 黑龙江农业科学， 2018（2）： 60-62.
• 张阳梅， 曾建生， 赵志平， 等. 菠萝种质资源果实主要数量性状的主成分分析[J]. 热带农业科技， 2016， 39（2）： 10-13， 28.
• 陈业渊. 热带、南亚热带果树种质资源描述规范[M]. 北京： 中国农业出版社， 2005.
• 陈业渊， 贺军虎. 热带、南亚热带果树种质资源数据质量控制规范[M]. 北京： 中国农业出版社， 2006.
• 全国果品标准化技术委员会. 水果和蔬菜可溶性固形物含量的测定——折射仪法： NY/T 2637-2014[S]. 北京： 中国标准出版社， 2014.
• 白世践， 李  超， 蔡军社， 等. 吐鲁番地区新征集葡萄资源果实主要品质性状的因子分析和聚类分析[J]. 西北农业学报， 2016， 25（7）： 1006-1016.
• 古丽尼沙·卡斯木， 木合塔尔·扎热， 张东亚， 等. 基于因子分析的无花果引进品种果实品质性状综合评价[J]. 食品科学， 2018， 39（1）： 99-104.
• 孙亚强， 吴翠云， 王  德， 等. 野生酸枣资源果实品质因子分析及评价指标选择[J]. 食品科学， 2016， 37（9）： 29-34.
• 杜  强， 贾丽艳. SPSS统计分析从入门到精通[M]. 北京： 人民邮电出版社， 2009.
• 王益民， 张  珂， 许飞华， 等. 不同品种枸杞子营养成分分析及评价[J]. 食品科学， 2014， 35（1）： 34-38.
• 薛  毅， 陳立萍. 统计建模与R软件[M]. 北京： 清华大学出版社， 2007.
• 李静萍. 多元统计分析——原理与基于SPASS的应用[M]. 2版. 北京： 中国人民大学出版社， 2015.
• 李静萍， 谢邦昌. 多元统计分析方法与应用[M]. 北京： 中国人民大学出版社， 2008.
• 袁志发， 周静芋. 多元统计分析[M]. 北京： 科学出版社， 2002.