生态学区组试验设计的方差分析及P值探讨

吕世杰，闫宝龙,2，王忠武，李治国*，，康萨如拉，刘红梅

(1.内蒙古农业大学草原与资源环境学院/草地资源教育部重点实验室/农业农村部饲草 栽培、加工与高效利用重点实验室/内蒙古自治区草地管理与利用重点实验室，呼和浩特 010019; 2.内蒙古民族大学农学院，通辽028043；3.内蒙古自治区林业科学研究院，呼和浩特 010010)

在生态学野外或田间试验中，我们常用到区组试验设计[1～3],原因是区组试验设计操作相对简单，数据分析相对容易掌握，研究者更愿意接受并用于揭示研究对象的变化特征和变化规律[4～6]。然而,事实上试验设计的操作简单不等于数据分析的科学运用。尽管有什么样的试验设计就会有什么样的数据分析方法，但数据分析方法受试验设计种类、样本容量、取样方法和模型参数等多方面的影响[4,5,7]。区组试验设计的数据分析方法首先考虑到的就是方差分析[4,5]。然而，方差分析模型的选取并不是单因素和双因素这么简单，也不是固定模型和随机模型这么容易，而是一个综合考量过程[5,6]。方差分析结果判断出显著差异之后，需要进一步做多重比较[5,6]。可是，我们在文献中经常会看到多重比较结果等同于方差分析[1～3]。因此，方差分析模型选取以及其与多重比较的关系有必要详细阐明。特别值得注意的是，近些年对于P值问题的争论引起了统计学界高度重视，对其他学科的影响也具有深远的意义[8～10]。本研究立足于生态学单因素区组试验设计，采用实例解析的方法对方差分析模型选取、多重比较、差异显著性(P值)逐一探讨，为生态学区组试验设计的数据分析提供科学合理的解决途径，也为生态学科学问题的阐释给予比较全面的理论支撑。

1 材料和方法

试验数据来源于内蒙古农业大学草原与资源环境学院四子王旗放牧试验基地，放牧试验基地采用区组试验设计[11]。2016年8月份，在每一个试验处理区随机选择10个50cm×50cm的样方，测定短花针茅植物种群的高度、盖度、密度和地上现存量以及植物群落地上现存量，并以该试验设计下的取样数据(高度、密度)展开方差分析相关问题的讨论。数据分析采用SAS 9.2软件，其中正态性检验调用UNIVARIATE过程，方差分析调用GLM过程，多重比较选用DUNCAN关键字，方差同质性选用HOVTEST关键字。

2 方差分析相关问题探讨

2.1 样本容量如何确定

从表1看出，每一载畜率下样本容量为27或30个观测数据，这一样本指的是观测样本容量，而不是试验设计样本容量。试验设计的样本容量由4个载畜率和3个区组构成，也就是总样本容量为12个数据。这一点可以从单因素区组试验设计的模型[4]中可以看到：

xij=μ+αi+βj+εij

(1)

式中，xij为观测数据(每一载畜率每一区组内短花针茅种群高度或密度)；μ为总体(短花针茅种群高度或密度)均值；αi为载畜率导致短花针茅种群高度或密度的差异；βj为区组试验设计中区组导致短花针茅种群高度或密度的差异。因此，每一载畜率每一区组内10个观测样方的观测数据不能直接用于方差分析，原因是违背了单因素区组试验设计数据分析的统计模型。那么每一载畜率每一区组内短花针茅种群高度或密度观测数据(表1中7或10个样本容量)还有没有意义，原因是生态学野外或田间试验空间异质性大，观测数据波动性大，需要增加观测重复来弥补数据的波动性，从而使获得的观测数据均值更稳定，更具有代表每一载畜率每一区组内短花针茅植物种群高度或密度指标的集中情况。

表1 不同载畜率下短花针茅植物种群高度和密度样本数据描述

2.2 正态性检验

确定了样本对象和样本容量，我们才能进入正态性检验环节，这是方差分析的前提条件之一。在荒漠草原，短花针茅高度属于数量性状数据(连续型变量)，而密度属于质量性状数据(非连续性变量)，所以有必要先假设荒漠草原短花针茅这一总体的高度、密度服从正态分布[5]。短花针茅高度、密度数据的样本容量均为12个，高度经Shapiro-wilk检验[4,6]统计量W=0.9524，P=0.6730；Kolmogorov-smirnov检验[4,6]结果显示，D=0.1098，P>0.1500；所以高度属于正态分布数据。密度经Shapiro-wilk检验统计量W=0.9679，P=0.8872；Kolmogorov-smirnov检验结果显示，D=0.1275，P>0.1500；所以密度也属于正态分布数据。因此，高度和密度指标可以进行下一步数据分析过程。

2.3 方差同质性检验

方差同质性检验调用SAS的GLM过程，一般在MEANS关键字后的待比较变量进行指定，格式为“MEANS CHL/HOVTEST DUNCAN ALPHA=0.05”，其中CHL代表载畜率变量，其余属于SAS系统关键字[12]。在这里需要明白一点，当进行多重比较时不管因素变量存在几个，只是针对其中的一个因素变量进行比较，其余因素变量均变为重复。所以，多因素变量的方差分析在进行方差同质性检验时模型只能指定一个因素变量。在本研究中，多重比较的因素变量为载畜率，相应的区组变量成为重复变量。也就是说，当进行方差同质性检验时，涉及的待比较因素只有一个，那就是CHL，即载畜率变量。当方差同质性检验通过时(一般要求P>0.05)，说明不同载畜率下3个区组的观测数据来自同一总体，即荒漠草原不同载畜率下短花针茅高度、密度来源于同一总体，其差异仅是由载畜率和区组差异引起。经检验发现，短花针茅高度F=0.80，P=0.5294；密度F=2.26，P=0.1587；P值均大于0.05，认为不同载畜率下的高度、密度指标没有因载畜率不同而产生差异，也就是说不同载畜率下高度、密度数据具有方差同质性，来源于同一个总体。

2.4 线性可加性检验

获得的样本数据，正态性和方差同质性检验通过后还需要进行线性可加性检验[5]，即公式(1)的模型检验。在进行线性模型线性可加性检验时，研究者往往过于看重F值和Pr>F值。看重F值的原因是根据相对大小进行变量取舍，以便在多因素模型中进行模型优化；看重Pr>F值的原因是判断是否继续进行多重比较的标准。然而，这种做法是欠妥当的，因我们线性可加模型是基于试验设计对观测数据的数学表述，所以拟合率的大小也是影响线性可加模型是否成立的关键指标[6]。本研究案例线性可加模型检验结果显示，高度F=2.36，P=0.1633，R2=0.6625；密度F=6.67，P=0.0194，R2=0.8475；所以高度方差分析检验没有通过，拟合率也比较低；而密度方差分析检验通过(P<0.05)，拟合率高达84.75%。

2.5 区组效应算不算一个因素

区组效应是不是一个因子，不同研究者具有不同的理解。茆诗松等[4]在编著的《试验设计》中指出，在进行线性可加模型检验时，观察区组效应是否显著大于误差效应是进行判别的依据，如果区组效应显著大于误差效应(P<0.05)，就要考虑区组效应的价值，当区组效应与误差效应无显著差异，则可直接划归为误差效应。然而，盖钧镒[7]在主编的《试验统计方法》中指出，区组效应可以看成另一个因素。本研究认为，区组试验设计前提要求区组内不同试验处理的基本条件尽可能一致，不同区组间的差异可以大一些，同时在野外或田间生态学实验中经常涉及山坡大小、水文条件甚至植被状况等自然条件限制，在安排试验时就应该尽可能考虑区组间差异，可以采用区组试验设计进行单向控制[7]。如果没有考虑，也不能进行补救的情况下，区组效应属于随机效应，如果进行单向控制，区组效应属于固定效应，这时无论是随机效应还是固定效应，方差分析线性模型不会改变，均属于双因素方差分析模型，对于试验处理的检验不再是简单的取舍问题，而是固定模型、随机模型和混合模型的问题[5]。这样随之而来会带来另外的问题，单因素区组设计当考虑区组设计为因素时，不能考虑交互作用(没有重复，缺少自由度)，所以也就不能判断处理因素强弱；进而考虑区组因素为随机因素时，不能对处理因素进行有效的检验(缺少自由度)。因此，这也成了单因素区组试验设计的缺点。

3 方差分析结果的呈现

3.1 方差分析表呈现

在方差分析结果呈现时，研究者对方差分析结果进行了精简(只有F值和Pr>F值)，我们只能看到因素的方差效应是否大于误差效应，难以看到因素效应的方差贡献，更看不到线性可加模型对原始数据的拟合效果，所以对于数据分析结果的科学性和可信性存在质疑。同时，对于存在随机因素的多因素方差分析模型，我们也无法判断采用的是固定模型、随机模型或混合模型中的哪一类。因此，综合来说方差分析结果应该采用表2的样式。对于单因素随机区组试验设计，其方差分析模型应该是双因素无交互作用方差分析模型(表2)。对于高度数据，其更符合单因素方差分析模型，此时区组效应属于随机效应，且对因素载畜率(CHL)具有干扰作用，因此采用单因素方差分析模型更为合适(此时区组效应成为随机误差效应)。尽管模型选择是以损失拟合率为代价，但是拟合率下降幅度不大，且能够表征载畜率因素对短花针茅高度的影响，所以结合多重比较结果(图1)，单因素方差分析模型更为准确。对于密度指标，采用双因素固定效应方差分析模型更为合理，此时具有较高的拟合率(R2=84.75%，即0.8475)，且载畜率对密度的影响也能够得到真实体现。所以，单因素区组试验设计的方差分析本质是双因素方差分析模型，而且处理效应与区组效应划归为固定效应，模型为固定效应模型，即全称应该为双因素固定效应方差分析模型。然而，在分析数据时受区组效应的影响，究竟是采用单因素方差分析模型还是双因素固定效应模型，需要通过检验结果进行判定。

3.2 多重比较结果呈现

在进行方差分析时，当方差分析检验结果显著时(P<0.05)需要对样本均值进行多重比较，且在概率水平下(一般P=0.0500或P=0.0100)进行差异显著性标记[5]。从多重比较结果来看，高度指标在不同载畜率下存在显著性差异(重牧的HG处理区短花针茅种群的高度显著低于CK和LG，P<0.05)，且伴随载畜率增大具有下降的变化趋势，因此多重比较结果印证了单因素方差分析模型的合理性(表2)。从密度比较结果来看，CK和LG处理区短花针茅密度显著低于MG和HG处理区(P<0.05)，结合拟合率，所以选用双因素固定效应模型更适合。

表2 荒漠草原短花针茅高度和密度方差分析表

4 讨论

4.1 方差分析和多重比较的合理性选择

多数情况下，研究者认为方差分析最简单，其最善于利用方差分析解决问题。然而,事实上方差分析最不简单，主要体现在以下几个方面:首先，方差分析3个前提假设(正态性、同质性和线性可加性均)需要进行严格的检验[5]；其次，方差分析的线性可加模型依赖于试验设计，有什么样的试验设计就会有什么样的数据分析方法，主要针对的就是方差分析[4]；第三，方差分析模型[5]存在单因素模型、双因素模型(又分有重复模型和无重复模型)和多因素模型(也分为有重复模型和无重复模型)；第四，按照因素是否为固定效应和随机效应，方差分析模型又分为固定模型、随机模型和混合模型[5]；第五，综合前四项条件，选择合理的方差分析模型并对数据进行科学分析是很困难的，甚至有时候方差分析这一方法不能应用。综上所述，方差分析需要考察样本容量对象、前提假设检验、试验设计种类、试验因素类别判定、线性可加模型及其拟合效果等诸多因素，期待在以后的研究中能够看到基于多方面考量的方差分析运用过程，进而保证数据分析佐证科学问题的可靠性和科学性。

4.2 方差分析结果呈现的发展趋势

方差分析结果更多的是以多重比较结果呈现[13～15]，但随着数据可视化理念的提倡，表格呈现方式逐渐被图形呈现方式替代。图形呈现方式伴随着计算机技术的提高，也存在发展趋势，即柱形图→柱形图+误差线→箱线图。目前，柱形图+误差线表示方法最多，但是值得注意的是误差线究竟是用标准偏差还是用标准误差，不同的研究者具有不同标注方法[16]。在统计学中，如果误差线采用的是标准偏差，表示获得样本数据为大样本数据(样本容量n≥30)，此时代表的统计学意义为样本容量中有95%或99%的样本观测数据在此区间；如果误差线采用的是标准误差，此时代表的统计学意义为样本均值有95%或99%的概率下在此区间内波动。由于是样本均值的多重比较，且在野外或田间区组试验设计中大样本数据很难获得，所以标注标准误差更为合理[16]。除了多重比较结果，还有方差分析结果，这一结果最初研究者是采用全表放置[17～19]，但是随着国际交流与合作的发展，简表形式出现；原因是在生态学领域不是研究统计学结果，而是利用统计学表征生态学专业研究结果。这一想法或说法并没有错，但是读者很难在多重比较结果和简表中读到数据拟合信息、模型选用信息，从而导致研究者在数据分析相互借鉴中产生偏差。因此，如果能交代清楚数据分析方法，建议尽可能交代清楚，比如单因素方差分析、简单的双因素方差分析；如果交代不清楚，建议将方差分析全表放上(如表2)，以便为读者提供更为全面的数据分析信息，也方便研究者之间的相互借鉴。

4.3 关于P值的界定和使用

在进行方差分析和多重比较时，均涉及P值的界定，且近几年关于P值的争论比较激烈[8,20～23]；首先是P=0.0490和P=0.0500的问题，其次P值是否合理的问题，最终可以归结为一个问题：P值究竟该怎么用，有没有必要用。这一争论最终以美国统计协会2019年给出的建议逐渐平息[9]：其认为在涉及概率统计的时候，标注P值的具体数值，不要过于强调是P=0.0500还是P=0.0100。本研究认为，P值是概率统计下的产物，是小样本推测总体特征的保障，因此不能因P=0.0490和P=0.0500区别否定整个概率统计；其次，P值显著临界点的定义是根据小概率事件是否发生定义的[5]，所以P=0.0500和P=0.0100仍然可以沿用；第三，临界点的定义从来都不是一成不变的，比如方差分析采用的临界点是P=0.0500和P=0.0100，而回归分析引入变量的临界点一般是P=0.1500，所以临界点的使用可以根据研究的专业内容进行合理的调整，比如张金屯[24]就有过相关的尝试和实例。综合来看，P值是进行概率检验的保障，不同的分析方法可以有不同的概率水平进行界定。比如方差分析，常用的临界点P=0.0500和P=0.0100仍可沿用，视具体情况也可适当调整，比如P=0.1000；但也不可能无限制的增大临界值，否则统计学的弃真和纳伪错误发生概率就会发生变化，毕竟增大了概率临界值也就增加了弃真错误的概率区间[5]。在SAS等统计软件中，P值并不是固定的，数据分析者可以根据研究内容进行指定，也就是说统计软件或统计学家也认为P值是可以调整的，从而适合不同研究专业和不同研究方向。结合美国统计学会针对P值的建议，方差分析还是要强调临界值选用的具体数值(毕竟多重比较结果是在特定临界值下计算的结果)，其他分析方法可以标注P值，然后根据自己的研究内容具体情况具体分析。此外，以后的数据分析，贝叶斯统计越来越受到统计学家的支持，生态学领域的数据分析也应该倾向于此。

5 结语

单因素区组试验设计的方差分析模型应该是双因素固定效应方差分析模型，且不能考虑试验处理效应与区组效应的交互作用；在方差分析过程中，如果存在区组效应干扰时，可调整为单因素方差分析模型；这个过程均需要对指标数据进行正态性、方差同质性检验，结合样本容量和线性拟合率综合探讨方差分析的可靠性和科学性。在进行多重比较时，需要给定具体的P值；对应的方差分析模型、处理效应检验可以根据研究情况，对P值可做出合理的调整；建议方差分析模型、处理效应检验和多重比较检验的P值最好一致。