关于重复测量方差分析的SAS程序请教！

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

有这样的一组试验数据：处理组别：A1(试验药高剂量）,A2（试验药低剂量）,C（安慰剂对照）；重复观测：T0（在用药之前的测量值），T1（用药之后1个月时的测量值），T2（用药之后2个月的测量值），T3（用药之后3个月的测量值）。另外A1入组12人，A2入组9人，A3入组13人，即unbalanced，想要知道3个处理组之间的治疗效果有无差异。程序：
data a;
input group $ t0 t1 t2 t3 @@;
cards;
A1 12 13 23 29
........
A2 13 24 34 29
........
A3 14 34 32 31
........
;
run;
proc glm data=a;
class group;
model t0 t1 t2 t3=group;
repeated time 4 contrast(1)/printe summary;
means group/snk lsd;
run;
quit;
其实，对于治疗效果还要考虑到基线值的差异（即清除T0的影响），所以应该把T0作为协变量来处理。在SPSS中是通过“Analyze"-"General Linear Model"-"Repeated Measures Define Factor(s)"然后再选择3个factors,作为时间因素，并把T0作为协变量，group作为between subject变量，再得出分析的结果。
那么，我用如下的程序是否就对了？
proc glm data=a;
class group;
model t1 t2 t3=group t0;
repeated time 3 contrast(1)/printe summary;
repeated time 3 contrast(2)/printe summary;
means group/snk lsd;
run;
quit;
我对了下SPSS和SAS的运行结果，出来的值（能对应的上的）都是一样的。
也欢迎各位大虾提供建议，多多批评和指正，先谢过了！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏8 回帖

关键词：重复测量方差分析 sas程序 方差分析 重复测量 unbalanced 安慰剂 测量 程序

期待回答啊。有没有在CRO公司或者在诺华、辉瑞等制药公司的统计分析部门的朋友，这种重复测量的方差分析肯定是会经常碰到的，所以还请大家多给意见和建议。

1# zhangweibeckham

为什么不试试混合线性模型呢？(proc mixed)
1. GLM是应用最小二乘法原理计算协方差参数的，这一模型对协方差结构有极端的假定。当需要估计尽可能多的方差及协方差参数时，GLM模型有一定的局限性，而基于似然函数法原理的混合效应线性模型分析方法可满足不受协方差参数限制的要求，而且MIXED模型还可以分析观察时间点不相等的资料，也能充分利用具有缺失观察值的资料，它允许资料存在某种相关性及协方差矩阵的多样性，更适应重复测量资料的特点。
2. 一般线性模型（GLM）要求误差项e具有独立一致的正态分布，而许多实际资料不能满足这一假定。Mixed 模型允许误差项具有更灵活的结构，包括相关性和方差的不齐性。
3. 混合效应线性模型（proc mixed）采用最大似然估计法（maximum likelihood，ML）和约束最大似然估计法（restricted maximum likelihood，REML）原理计算协方差矩阵。
4. 应用混合效应线性模型的步骤：
  ① 确定固定效应和随机效应；
  ② 选择协方差结构，常见的有7种。
       a. 独立结构（又称方差分量结构） VC ，矩阵中含1个协方差参数；
       b. 复合对称结构CS，矩阵中含2个协方差参数；
       c. 空间幂相关结构 SP(POW)，含2个协方差参数；
      d. 无结构（又称不规则结构） UN，含n(n+1)/2个协方差参数；
       e. 一阶自回归结构 AR(1)，含2个协方差参数；
       f. 带状主对角结构 UN(1)，含n个协方差参数；
       g. 循环相关结构 TOEP，含n个协方差参数；
  ③上述不同的协方差矩阵中，选出似然比统计量（-2 Log Likelihood）、Akaike’s Information Criterion(AIC)、及Schwartz’s Bayesian Criterion（BIC）较小的一个。如果这些统计量很近似，则选含参数个数最少的一个，通常以AIC为主要判断指标。
  ④ 选定协方差结构后，再筛选固定效应参数，剔除无统计学意义的高阶效应。

1# zhangweibeckham

代码如下：

1. Data a;
   
2. Input group pnt t0 t1 t2 t3 @@;
   
3. Cards;
   
4. 1 1 12 13 23 29
   
5. 1 2 . . . .
   
6. …
   
7. 1 12 . . . .
   
8. 2 1 13 24 34 29
   
9. 2 2 . . . .
   
10. …
    
11. 2 9 . . . .
    
12. …
    
13. 3 13 . . .
    
14. ;
    
15. Data a_b;
    
16.    Set a;
    
17.    Array t{4} t0 t1 t2 t3;
    
18.    Do i= 1 to 4;
    
19.       If i=1 then time=0;
    
20.      If i=2 then time=1;
    
21.      if i=3 then time=2;
    
22.      If i=4 then time=3;
    
23.   Y=t{i};output;
    
24.   Timepnt=time;
    
25.   End;
    
26.   Drop i t0 t1 t2 t3 ;
    
27. Run;
    
28. Proc mixed data=a_b covtest method =REML;
    
29.    CLASS group pnt;
    
30.    Model y=group time group*time time*time time*time*group /S; /*选项S输出随机效应参数的解 */
    
31.    Repeated /type= UN subject=pnt R; /* 选项R输出R矩阵的区块 */
    
32.    Lsmeans group time / tdiff pdiff;
    
33. Run;
    
34. Proc mixed data=a_b covtest method =REML;
    
35.   CLASS group pnt;
    
36.   Model y=group time group*time time*time time*time*group /S;
    
37.   Repeated /type= UN(1) subject=pnt R;
    
38.   Lsmeans group time / tdiff pdiff;Run;
    
39. Proc mixed data=a_b covtest method =REML;
    
40.   CLASS group pnt;
    
41.   Model y=group time group*time time*time time*time*group /S;
    
42.   Repeated /type= CS subject=pnt R;
    
43.   Lsmeans group time / tdiff pdiff;
    
44. Run;
    
45. Proc mixed data=a_b covtest method =REML;
    
46.   CLASS group pnt;
    
47.   Model y=group time group*time time*time time*time*group /S;
    
48.   Repeated /type= SP(POW)(timepnt) subject=pnt R;
    
49.   Lsmeans group time / tdiff pdiff;
    
50. Run;
    
51. Proc mixed data=a_b covtest method =REML;
    
52.   CLASS group pnt;
    
53.   Model y=group time group*time time*time time*time*group /S;
    
54.   Repeated /type= AR(1) subject=pnt R;
    
55.   Lsmeans group time / tdiff pdiff;
    
56. Run;
    
57. Proc mixed data=a_b covtest method =REML;
    
58.   CLASS group pnt;
    
59.   Model y=group time group*time time*time time*time*group /S;
    
60.   Repeated /type= TOEP subject=pnt R;
    
61.   Lsmeans group time / tdiff pdiff;
    
62. Run;
    
63. /* 根据AIC准则，选择AIC较小的协方差结构，采用Type 3 tests of Fixed Effects 结果 */
    
64. /* Lsmeans 语句输出不同效应间两两比较的结果 */

复制代码

没有原始数据，楼主自己运行一下吧。

sean1945 发表于 2010-12-28 11:21
1# zhangweibeckham

为什么不试试混合线性模型呢？(proc mixed)
1. GLM是应用最小二乘法原理计算协方差参数的，这一模型对协方差结构有极端的假定。当需要估计尽可能多的方差及协方差参数时，GLM模型有一定的局限性，而基于似然函数法原理的混合效应线性模型分析方法可满足不受协方差参数限制的要求，而且MIXED模型还可以分析观察时间点不相等的资料，也能充分利用具有缺失观察值的资料，它允许资料存在某种相关性及协方差矩阵的多样性，更适应重复测量资料的特点。
2. 一般线性模型（GLM）要求误差项e具有独立一致的正态分布，而许多实际资料不能满足这一假定。Mixed 模型允许误差项具有更灵活的结构，包括相关性和方差的不齐性。
3. 混合效应线性模型（proc mixed）采用最大似然估计法（maximum likelihood，ML）和约束最大似然估计法（restricted maximum likelihood，REML）原理计算协方差矩阵。
4. 应用混合效应线性模型的步骤：
  ① 确定固定效应和随机效应；
  ② 选择协方差结构，常见的有7种。
       a. 独立结构（又称方差分量结构） VC ，矩阵中含1个协方差参数；
       b. 复合对称结构CS，矩阵中含2个协方差参数；
       c. 空间幂相关结构 SP(POW)，含2个协方差参数；
      d. 无结构（又称不规则结构） UN，含n(n+1)/2个协方差参数；
       e. 一阶自回归结构 AR(1)，含2个协方差参数；
       f. 带状主对角结构 UN(1)，含n个协方差参数；
       g. 循环相关结构 TOEP，含n个协方差参数；
  ③上述不同的协方差矩阵中，选出似然比统计量（-2 Log Likelihood）、Akaike’s Information Criterion(AIC)、及Schwartz’s Bayesian Criterion（BIC）较小的一个。如果这些统计量很近似，则选含参数个数最少的一个，通常以AIC为主要判断指标。
  ④ 选定协方差结构后，再筛选固定效应参数，剔除无统计学意义的高阶效应。

理论很深的嘛，学习了！

谢谢您的回答，非常详细，不过你的程序和思路和我想要的做的分析目的不同，我是想去出第一次观测的影响（即基线值的影响）去考察3个组别之间的治疗效果有没有差异，所以T0应该不能放在因变量中考虑，你这个模型是个更一般的方差分析，约束条件少了。我已经确定我的程序处理这种剔除基线差异的重复观测的方差分析是合适的。再一次感谢您使我受益匪浅！

呵呵 要去除基线值影响的话 在PROC mixed /glimmix 里面加个RANDOM语句就OK了~

同意“sean1945"的说法，用proc mixed 是目前最好的方法。你用MANOVA得出的东西还是要用ANOVA两两时间去比，并且要考虑多重比较的alpha.

MIXED model 利用了各时点的相关性，节省了样本。如”想去出第一次观测的影响（即基线值的影响）去考察3个组别之间的治疗效果有没有差异，“，把基线值作为协变量即可。

Proc mixed data=a_b covtest method =REML;
  CLASS group pnt;
  Model y=group time group*time t0/S;
  Repeated /type= CS subject=pnt R;
  Lsmeans group time / tdiff pdiff;
Run;

LS说得对。。才发现我错了，LZ抱歉!

sean1945 发表于 2010-12-28 11:22
1# zhangweibeckham

代码如下：

请问 pnt 是哪个变量？