随机效应与混合效应模型 SAS实践

两个处理变量都具有固定效应
两个处理变量都具有随机效应
- GLM procedure
- Mixed procedure, type 1 method
- Mixed procedure, REML method
只有一个处理变量具有随机效应
- GLM procedure
- - Restricted Model
  - Unrestricted Model
  - 两两比较
- Mixed procedure

精密仪器制造非常依赖零部件尺寸的精确度。一位工程师随机选取了三种仪器对某零部件的二十个关键部位进行了测量，每一个部位重复测量了两次。这个测量试验的目的是研究零部件部位与测量仪器对测量结果的影响，零部件部位与测量仪器是处理变量。实际上这两个处理变量都应该考虑为随机的，但为了比较固定效应与随机效应的区别，下面我们分别用固定效应、随机效应与混合效应模型分析这个试验。

两个处理变量都具有固定效应

假设两个处理变量都是固定的。

proc glm data=randr;/*数据表读取比较长，放在文末，数据表名为randr*/
class operator part;/*定义仪器与部位为类型变量*/
model resp=operator|part;/*|表示使用full model，即operator part operator*part*/
run;

这个GLM模型是显著的，下面是处理变量的Type III 平方和分解：

这个结果说明测量仪器不会对测量结果造成显著的影响，但不管用什么仪器测量，每个部位测量精确度都有显著差异。

两个处理变量都具有随机效应

GLM procedure

假设两个处理变量都是随机的。

proc glm data=randr;
class operator part;
model resp=operator|part;
random operator part operator*part/test;
run;

与固定效应模型相比，随机效应模型只需要在GLM procedure中加入random语句。在这段SAS代码中，random operator part operator*part/test; 语句表示operator part operator*part这三个因素都具有随机效应，/test表示对随机效应进行F检验，之所以要加上/test是因为随机效应模型中并非所有的均方和都与残差的均方和比较：

随机效应模型也是显著的，下面是处理变量的Type III 平方和分解：

这个结果与随机模型的结果一模一样，但是我们只能用这个结果做交互项operator*part的推断，要对另外两个因素做推断，需要random语句的结果：

这两张表是random语句的结果。在第二张表中，交互项与残差的平方和分解结果与固定效应模型和随机效应模型的都相同；第一张表是random语句中/test的结果，对于operator与part而言，error term是operator*part，这个结果同样说明测量仪器不会对测量结果造成显著的影响，但每个部位测量精确度都有显著差异。

需要注意的是random语句中的/test语句可以用test语句代替：

proc glm data=randr;
class operator part;
model resp=operator|part;
random operator part operator*part;
test H=operator E=operator*part;
test H=part E=operator*part;
run;

语句test H=operator E=operator*part;表示检验operator的效应，用operator*part的均方和作为error term的均方和。语句test H=part E=operator*part;表示检验part的效应，用operator*part的均方和作为error term的均方和。这两个语句输出的结果与上面的第一张表相同：

Mixed procedure, type 1 method

随机效应模型也可以用Mixed procedure来做，具体有两种方法，type 1 method和REML method：

proc mixed cl maxiter=20 covtest method=type1 data = randr;/*covtest表示估计每个效应的方差、cl表示列出方差的置信区间*/
class operator part;
model resp = ;/*mixed procedure后的随机因素不用写出来*/
random operator part operator*part;/*mixed procedure没有/test语句*/
run;

type 1 method输出的结果与GLM procedure的random语句输出结果一样，下面是covtest与cl语句的输出，包括每个因素的方差估计、标准差以及置信区间，type 1 method估计方差的方法是最小二乘法。

Mixed procedure, REML method

proc mixed cl maxiter=20 covtest method=reml;
class operator part;
model resp = ;
random operator part operator*part;
run;

REML method估计方差的方法是残差最大似然方法，是Mixed procedure的默认方法。

只有一个处理变量具有随机效应

假设operator是固定的因素。

GLM procedure

Restricted Model

proc glm data=randr;
class operator part;
model resp=operator|part;
random  part operator*part; /*operator是固定的，随机的只有part和交互项*/
test H=operator E=operator*part;
run;

注意restricted model在做推断时，检验如下

因此上面那段SAS代码中加了一个test语句，test H=operator E=operator*part;用以对operator的效应做推断。

这是GLM procedure的输出结果，part和operator的结果是可以参考的，operator的结果需要参考test语句的输出：

Unrestricted Model

proc glm data=randr;
class operator part;
model resp=operator|part;
random  part operator*part/test;
run;

Restricted Model与Unrestricted Model相比唯一的区别就是交互项完全随机还是交互项关于固定效应的因子的边际效应是固定的。根据Unrestricted model的EMS，

part与operator的均方和都需要与它们的交互项进行比较，所以用/test语句来实现，下面是/test的输出结果

两两比较

proc glm;
class operator part;
model resp=operator|part;
random part operator*part / test;
means operator / tukey lines E=operator*part;
lsmeans operator / adjust=tukey E=operator*part tdiff stderr;
run;

一般对固定效应的因素，不同的因素取值对测量值的效应两两之间可以进行比较。means语句和lsmeans语句都可以做比较，means语句means operator / tukey lines E=operator*part;表示对operator的效应进行比较，用Tukey方法，设定error term为operator*part，lines表示用线图表示结果：

Estimate是每一种测量仪器测量结果的均值，蓝线贯穿了三个均值表示这三个均值没有显著差异。lsmeans语句lsmeans operator / adjust=tukey E=operator*part tdiff stderr;可以做更精细的比较：

Mixed procedure

proc mixed alpha=.05 cl covtest;
class operator part;
model resp=operator / ddfm=kr;
random part operator*part;
lsmeans operator / alpha=.05 cl diff adjust=tukey;
run;

上面的SAS代码中只有ddfm=kr语句需要注意一下，ddfm或者ddf用来指定交互项自由度的计算方法，一般用kr就可以了。lsmeans语句的输出包括operator均值的估计、置信区间以及两两比较：