裂区试验设计中统计分析的SAS线性混合模型研究

董中东 张 宁 赵 磊 任 妍 孙丛苇 陈 锋

（河南农业大学农学院，河南郑州 450046）

裂区试验设计是农业试验中常用的一种试验设计方法[1-3]。SAS GLIMM功效分析表明，当试验单元误差相同时，裂区试验设计较不完全随机区组和完全随机区组有着更高的统计功效。裂区试验设计可以分为二因素裂区和三因素裂区，后者又有其中两因素组合作主区或者作副区的不同裂区试验。但是，在裂区试验的设计和分析过程中经常会存在一些问题，尤其是其统计分析。裂区试验的统计分析难点主要有2个方面：一方面，在进行方差分析各因素主效应和互作效应的检验所用到的误差项不同；另一方面，最主要的难点是其互作效应的多重比较，由于在对不同的处理组合进行比较时用到的误差方差及其自由度的计算非常繁琐，因而裂区试验在统计分析时经常会存在一些问题[4-5]。

SAS 软件（SAS Institute，Cary，NC）是一款功能非常强大的统计分析软件，尤其是近年来其开发了线性混合模型（linear mixed model）和广义线性混合模型（generalized linear mixed model）的应用程序 PROC MIXED和PROC GLIMMIX，这些程序均使用了很多近年来的统计理论和方法，使SAS对试验数据的分析有了很大提高。由于裂区试验的分析需要估计2种误差方差以及互作项的多重比较分母误差项的矫正等问题，在SAS公司出版的著作SAS for Mixed Models中不建议使用PROC GLM程序进行裂区试验的分析。当试验数据满足方差分析的基本假定时，若涉及随机效应的估计，PROC MIXED基本不需要太多额外的程序编写，用其标准程序语句就可以实现。PROC MIXED还可以实现三因素裂区、时间裂区、时空裂区、多年多点等试验的统计分析。这里仅以二因素裂区试验设计的统计分析介绍SAS mixed程序和分析结果。

1 裂区统计分析的混合线性模型程序及结果解释

首先以一个裂区试验介绍一般情况下的SAS分析程序：有一水稻不同灌水方式（A）和种植密度（B）试验，主处理为3种灌水方式，分别以A1、A2和A3表示，副处理为3种密度，分别以B1、B2和B3表示，设置3个区组，其小区产量结果如表1，本例题和数据引自《农业试验设计与统计分析》[6]。

表1 灌水方式和密度二因素裂区试验产量结果

分析程序如下：

程序解释：首先用DATA步建立一个名为spd_1的数据集，数据集包含4个变数，其名称分别为blk、A、B 和 yld，blk代表区组，分别用 1、2、3代表 3个不同区组；A代表A因素，A1、A2和A3代表A因素的3个水平，为字符型，需加$说明其为字符型；B代表B因素，B1、B2和B3代表B因素的3个水平，也为字符型，需加以说明；yld代表小区产量。@@符号表明数据连续读取，datalines下面的数据为要读取的数据，到分号代表数据输入结束。

接下来为PROC MIXED程序步，首先PROC MIXED说明调用的程序名称，nobound选项的作用是在方差组分的估计时，避免SAS把估计的负的方差组分当作0来输出和进一步计算，此选项在主区误差组分为负值时尤其重要，它会影响到主处理的检验结果。model语句等号左边yld为要分析的试验指标或性状，右边为固定效应项，注意这里和常用的PROC GLM用法不同，mixed程序中等号右边只包含固定效应，而不能有随机效应，A|B，为A、B和A*B这3项的简写。后面的选项DDFM，是进行计算分母自由度的矫正方法，有6种方法可供选择，现在一般推荐使用KENWARDROGER（可简写为KENROG或 KR）或 SATTERTHWAITE（可简写为 SATTERTH或SAT）方法。random语句是指定随机效应，其对于裂区的统计分析尤为重要，只有把随机效应项设置正确才能得到正确的分析结果。这里把区组作为随机效应进行分析，blk*A则是定义主区误差项的随机效应。lsmeans语句定义输出A、B和A*B项的最小二乘法估计的均值，后面的选项diff是进行A、B和A*B的各个水平间最小二乘法估计均值差数的t检验（LSD 法），adjust=tukey，输出 tukey法多重比较的均值差数的结果。输出结果如表2、表3、表4所示。

表2 各随机效应项的方差估计

表3 固定效应的检验结果

表4 各处理和处理组合均值差数显著性部分结果

表2为各随机效应项的方差分量估计值，由于在程序中有nobound选项，所以区组项的方差为负值，如果没有该选项，其值将会为0。另外2项分别为主区误差组分和副区误差组分。

表3是各固定效应的F检验结果，结果表明A因素各水平间差异不显著（P=0.471 8），B因素的水平间差异极显著（P=0.004 0），2个因素存在极显著的交互作用（P=0.000 8）。由于互作项存在极显著差异，主要对互作项的结果进行解释和分析。表4为各因素不同水平和处理组合差数的比较结果，对于各因素水平间的差异显著性结果，只是看其计算结果是否正确，本例题应主要看各处理组合的差异显著性分析。先看同主区不同副区的一个比较结果，以LSD法的A1B1和A1B2为例，该比较只与副区误差有关，其自由度为副区误差自由度12。再看不同主区的比较，以A1B1和A2B2为例，由于不同主区比较的误差是主区误差和副区误差的合成，所以其既不是主区误差自由度，也不是副区误差自由度，需要根据公式矫正，其矫正公式如下：

式中，b为副处理水平，MSEa为主区误差均方，MSEb为副区误差均方，dfEa为主区误差自由度，dfEb为副区误差自由度，矫正后的自由度为12.3。从分析结果可以看出，对于互作项的多重比较，程序会自动根据要比较的2个处理组合算出其所需标准误和自由度，这种比较使用PROC GLM不能实现。

通过使用PROC MIXED对二因素裂区试验结果的分析可以看出，分析程序非常简洁，而且可以得到所要分析的结果。但是，由于处理组合的比较非常多，需要把结果进一步整理。当试验需要对一些综合的效应进行比较时，比如A1水平下B1+B2的平均效应与A1水平下B3的效应进行比较时需要使用estimate或者contrast语句才能分析出结果，本文不再介绍。

2 裂区分析的一般线性模型程序及结果

由于SAS mixed程序输出的方差分析表中没有各效应项的SS和MS，为便于验证mixed程序结果的正确性，这里给出SAS GLM程序和输出的方差分析结果。

GLM程序如下：

方差分析结果如表5、表6所示。从表5和表6可以看出，A因素的F值（0.91）和对应概率（0.471 8）、B 因素的 F 值（9.06）和对应概率（0.004 0）、A×B 互作的F值（10.22）和对应概率（0.000 8）与mixed程序输出结果完全一致。处理组合的多重比较可以在程序中加入语句“lsmeans A*B/diff；”进行，在结果中有需要误差和自由度合成的结果都是不正确的。通过GLM和mixed程序的运行结果可以看出，在裂区分析中mixed程序较glm程序更加准确、方便和快捷。

表5 方差分析结果

表6 使用Ⅲ型MS作为blk*A的误差项的假设检验