SAS Estimate 或 Contrast

看到一個簡明例子，分享給大家。

 

 

 

 

 

1.

lsmeans vacgrp ;
estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;
estimate '用estimate計算lsmean PBO' intercept 1 vacgrp 0 1;

 

 

 

兩種結果計算結果相同。一定記得在estimate中寫 intercept 1.這是截距。第一個estimate是估計ACT的最小二乘均值，第二個是估計PBO的。

estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;

負號表示被減，H0: PBO - ACT  = 0;

 

 

estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0 visit .25 .25 .25 .25 visit*vacgrp .25 .25 .25 .25 0 0 0 0;
estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1 visit*vacgrp -.25 -.25 -.25 -.25 .25 .25 .25 .25;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;

 

estimate '用estimate計算lsmean ACT' intercept 4 vacgrp 4 0 visit 1 1 1 1 visit*vacgrp 1 1 1 1 0 0 0 0 / divisor = 4;
estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -4 4 visit*vacgrp -1 -1 -1 -1 1 1 1 1 / divisor = 4;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;
estimate 'HI vs PB at Site 1' vacgrp vacgrp*siteid 1 0 -1
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 1 0 0 0     *HI at site 1 2 3 4;
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 0 0 0 0     *LO at site 1 2 3 4;
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 -1 0 0 0     *PB at site 1 2 3 4;

 

 

上述兩段code作用完全一樣。visit*vacgrp是說：vacgrp = ACT，對應的VISIT*VACGRP標上需要得數字。

 Code中 0.25那的意義是每個visit占的權重是一樣的，當計算均值時。每個visit樣本量占總體的比例。

當使用estimate估計均值時，計算得出的統計量和means中得出的統計量一致。參考方差分析 - GLM.

 

 

 

 

class visit vacgrp pat ;
lsmeans visit*vacgrp;
estimate '用estimate計算lsmean ACT' intercept 4 vacgrp 4 0 visit 1 1 1 1 visit*vacgrp 1 0 1 0 1 0 1 0 / divisor = 4;
estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -4 4 visit*vacgrp -1 1 -1 1 -1 1 -1 1 / divisor = 4;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;

class中變量出現得順序，影響交互作用相乘時得順序。

 

 

proc iml;
      a = {0 1};
      b = {0.2 0.2 0.2 0.2 0.2};
      ab = b@a;
      print ab;
quit;

如果變量比較多，不太好區分交互作用乘積，用上述代碼算。

data arthr;
 input vacgrp $ pat mo1 mo2 mo3 mo4;
 datalines;
ACT 101 6 3 0 3
ACT 103 7 3 1 2
ACT 104 4 1 2 8
ACT 107 8 4 3 1
PBO 102 6 5 5 7
PBO 105 9 4 6 6
PBO 106 5 3 4 3
PBO 108 6 2 3 9
; 

data discom; set arthr;
 keep vacgrp pat visit score;
 score = mo1; visit = 1; output;
 score = mo2; visit = 2; output;
 score = mo3; visit = 3; output;
 score = mo4; visit = 4; output;
run;

ods html;
proc mixed data = discom;
class vacgrp pat visit ;
model score = vacgrp visit visit*vacgrp;
repeated visit / subject=pat(vacgrp) type=un r rcorr;
/*lsmeans vacgrp ;*/
/*estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;*/
/*estimate '用estimate計算lsmean PBO' intercept 1 vacgrp 0 1;*/
/*estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;*/
lsmeans visit*vacgrp;
estimate '用estimate計算lsmean ACT' intercept 4 vacgrp 4 0 visit 1 1 1 1 visit*vacgrp 1 1 1 1 0 0 0 0 / divisor = 4;
estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -4 4 visit*vacgrp -1 -1 -1 -1 1 1 1 1 / divisor = 4;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;
estimate 'Month 3 Change from Baseline: ACT v. PBO' 
  vacgrp * visit -1 0 0 1 1 0 0 -1;
run;

12proc mixed data = discom;
class visit vacgrp pat ;
model score = vacgrp visit visit*vacgrp;
repeated visit / subject=pat(vacgrp) type=un r rcorr;
class visit vacgrp pat ;
lsmeans visit*vacgrp;
estimate '用estimate計算lsmean ACT' intercept 4 vacgrp 4 0 visit 1 1 1 1 visit*vacgrp 1 0 1 0 1 0 1 0 / divisor = 4;
estimate '用estimate計算lsmean ACT' intercept 1 vacgrp 1 0;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -4 4 visit*vacgrp -1 1 -1 1 -1 1 -1 1 / divisor = 4;
estimate '用estimate計算t檢驗 ACT PBO' vacgrp -1 1;
run;

proc iml;
      a = {1 -1};
      b = {-1 0 0 1};
      ab = a@b;
      print ab;
quit;
*負號表示被減;

負號表示被減, b是visit 3 - baseline得。