Split-plot设计 SAS实践

试验资源的限制
第一种统计模型
- GLM procedure
- Mixed procedure
第二种统计模型
- GLM procedure
- Mixed procedure

Split-plot设计有两大难点，第一个难点是判断什么时候用split-plot设计；第二个难点是判断哪些效应可以省略。某造纸厂工程师试图研究如何提高纸张的抗拉强度，为此他选择了3种纸浆制备方法、完成每种方法时分别使用四种不同的温度，重复试验三次，一共需要试验36组。受制备中间程序限制，每天只能完成12组试验，他打算以每一天作为一个block，分三天完成试验。下面是试验数据：

data paper;
input block method temp resp@@;
datalines;
1 1 1 30 1 1 2 35 1 1 3 37 1 1 4 36
1 2 1 34 1 2 2 41 1 2 3 38 1 2 4 42
1 3 1 29 1 3 2 26 1 3 3 33 1 3 4 36
2 1 1 28 2 1 2 32 2 1 3 40 2 1 4 41
2 2 1 31 2 2 2 36 2 2 3 42 2 2 4 40
2 3 1 31 2 3 2 30 2 3 3 32 2 3 4 40
3 1 1 31 3 1 2 37 3 1 3 41 3 1 4 40
3 2 1 35 3 2 2 40 3 2 3 39 3 2 4 44
3 3 1 32 3 3 2 34 3 3 3 39 3 3 4 45
;
run;

试验资源的限制

在用这个例子介绍Split-plot试验的分析方法之前，我们可以总结一下Split-plot设计的适用范围。随机效应与混合效应模型是因子交叉设计的分析方法，当试验经费不足，为了节约试验资源的时候可以使用nested design，但如果是试验资源的限制，就只能使用split-plot design了，常见的试验资源限制有两种：时间限制与空间限制，造纸的试验是时间限制。

在介绍裂区设计的博客中，我介绍了另一个例子（空间限制）：要研究肥料对水稻的产量影响，考虑6种不同的水稻品种与4种不同的肥料，重复试验3次（也被叫做3个blocking）；按交叉设计的思路，一共需要6×4×3=726\times 4 \times 3 = 726×4×3=72块地，假设每平方米种25墩，每墩插4株，每一个试验单位为500株，则一共需要72×500/(4×25)=36072 \times 500 /(4 \times 25) = 36072×500/(4×25)=360平方米的试验田；如果实验室的试验田只有180平方米，那就可以考虑裂区设计，把试验田平均分为3部分（field），用来做3次重复试验，每部分面积为60平方米，然后把3部分每一部分平均分成4大块地（whole-plots），每一部分中每一块地随机施一种肥料，每块地施的肥料各不相同，将每一大块地平均分为6小块（sub-plots），每一大块地中每一小块随机种一种水稻，每块地种的水稻各不相同。尽管这样每一小块地只能种(4×25)×60/24=250(4 \times 25) \times 60/24 = 250(4×25)×60/24=250株，但如果可以对裂区设计的试验数据建立适当的统计模型，我们同样可以利用这有限的资源完成试验并得出结论。在这个例子中称肥料种类为whole-plot factor，称水稻品种为sub-plot factor。

对比这个试验，对于提高纸张的抗拉强度试验这个例子，我们也可以总结出：

试验田分为三部分对应的是分三天完成试验，这既是replicates也是blocking；
三种纸浆制备方法是whole-plot factor；
四种不同的温度是sub-plot factor；

与区组化析因设计相比，Split-plot设计的factor在每个block中是嵌套（有明显的顺序、从属关系）而不是交叉的关系；并且Split-plot设计的区组一定是随机的。

第一种统计模型

GLM procedure

proc glm data=paper;
class block method temp;
model resp=block method block*method temp block*temp method*temp block*method*temp;
random block block*method block*temp block*method*temp;
test h=method e=block*method;
test h=temp e=block*temp;
test h=method*temp e=block*method*temp;
run;

这段SAS代码中，model语句model resp=block method block*method temp blocktemp methodtemp block*method*temp;用来表示block和两种因素及其所有的交互项，一共是23−1=72^3-1=723−1=7种。random语句random block block*method block*temp block*method*temp;用来指定模型中的随机项，因为区组是随机的，包含区组的交互项都是随机的。

根据这个EMS的结论，要比较纸浆制备方法与温度及其交互项的效应，分别需要block*method、block*temp、block*method*temp的均方和作为error term，三个test语句就是为了做这些检验，结果如下：

这些结果说明在0.1的显著性水平下，纸浆制备方法与制备使用的温度对纸张抗拉强度都有显著影响。

Mixed procedure

proc mixed data=paper method=type1 covtest cl;/*因为随机项都与区组有关，没啥特别意义可以不用covtest和cl这两个命令*/
class block method temp;
model resp=method temp method*temp;
random block block*method block*temp block*method*temp;
run;

红框内的结果与GLM procedure的test语句的结果相同。

第二种统计模型

与第一种统计模型相比，第二种统计模型去掉了区组与因素的交互项以及sub-plot error，这两项被计入到random error中。根据第一种统计模型EMS的结果，去掉区组与因素的交互项以及sub-plot error就不能对subplot factor以及whole-plot factor与sub-plot factor的交互项做检验。

GLM procedure

proc glm data=paper;
class block method temp;
model resp=block method block*method temp method*temp;
random block block*method;
test h=method e=block*method;
run;

第二种统计模型只能做制备方法的效应做检验，结果是第一种统计模型相同。

Mixed procedure

proc mixed data=paper method=type1;
class block method temp;
model resp=method temp method*temp;
random block block*method;
run;

尽管Mixed procedure返回完整的ANOVA，但依然只有红框内的可以参考，temp与method*temp的结果是基于残差作为error term计算出来的，与理论结果不一致，不能参考！