Practicalities of Using ESTIMATE and CONTRAST Statements

原文載點：http://support.sas.com/resources/papers/proceedings10/269-2010.pdf

在 PROC GLM 和 PROC MIXED 裡面，最令一般 SAS 使用者感到頭痛的不是模型的檢定而是事後分析(post hoc test)。這牽涉到 contrast 和 estimate 兩個語法。在我接觸過所有非統計科班出身但需要用到線性模型以及事後分析的人，沒有人不對這兩個語法感到聞風喪膽。多年人我曾經試圖用各種方法來教人如何使用這兩種語法，但效果都不怎麼顯著。簡單來講，Contrast 和 estimate 無論是在 PROC GLM 還是 PROC MIXED 裡面，功能完全一樣，語法內的選項也一樣，而兩者之間唯一的差別是，contrast 可以只能算出 p-value，而 estimate 可以把指定的點估計量及其 p-value 一次輸出。David J. Pasta 在 SAS Global Forum 2010 發表了一篇技術文件將這兩種語法做一個詳盡的介紹，堪稱一大貢獻。在本文中，我會介紹幾個比較實用的範例。
原文作者套用一個醫學資料並詳述其背景，但為了簡化初學者的時間，我將會以程式來主。第一個範例資料只有三個變數，分別是作為應變數的離散型變數 sex (兩層) 和 race (五層)，而依變數則為一個連續型變數 y1。有了對資料的基本概念，首先來看一個最簡單的例子：

proc glm data=anal;
      class sex;
      model y1 = sex / solution;
      lsmeans sex / stderr tdiff e;
run;

這個例子只單獨分析 sex 對 y1 的影響，而 lsmeans 是用來算出每個性別對 y1 的最小平方估計量。後面的 tdiff 則是去求兩個性別的差異。這段程式會跟用 PROC TTEST 做出來的結果一模一樣。如果不用 lsmeans 而改用 estimate 語法，則可以這樣寫：

proc glm data=anal;
      class sex;
      model y1 = sex / solution;
      estimate 'male' intercept 1 sex 0 1;
      estimate 'female' intercept 1 sex 1 0;
      estimate 'female-male' sex 1 -1;
run;

頭兩個 estimate 語法可算出不同性別在 y1 上的平均值，後面接著 intercept 和 sex 是為了要告訴 SAS 估計值是從哪些變數來的。一般使用者最常忘記的就是把 intercept 1 加上去，因為除非在 model 的選項上面加上 noint 抑制掉截距項的計算，否則一切用 estimate 求算的點估計量都要加上 intercept 1。不過你也不用擔心不加會怎麼樣，因為 log 視窗裡面自然會出現紅色的錯誤訊號提醒你加上。而 sex 後面跟著「0  1」和「1  0」則是要準確地將要計算的性別定位出來，當然你也要對自己如何將 sex 編碼要非常了解，而不是亂安插一通。在本例中，女性的編碼順位是排在男生前面（若真的不清楚編碼順位，看一下輸出報表一開始的 Class Level Information 裡面就會清清楚楚地把所有類別變數的編碼順位列出來），所以當 sex 後面寫「1  0」時，自然是去算女生的值，反之則是算男生的值。

而當要比較女生和男生的差異時，則有兩點要注意。一是 intercept 可以不用寫了，因為兩兩相減的情況下，截距項自然就被減掉了。二是 sex 後面要改成「1  -1」，很明顯地這就是指用女生的估計值去減男生的估計值。如果要算男生的估計值減女生的估計值，把「1  -1」改成「-1  1」即可。偷懶的方法是不改，但把 output 輸出的估計值變號即可，至於 p-value 是不會改變的。

當考慮兩個應變數時，程式如下：

proc glm data=anal;
      class race sex;
      model y1 = race sex / solution e;
      lsmeans race sex / stderr tdiff e;
      lsmeans race sex / stderr tdiff e om;
run;

第一個 lsmeans 語法會一次把兩個變數的最小平方估計量通通算出來，並且會把成對比較的結果也列印出來。第二個 lsmeans 語法和第一個的不同是後面多了 om 這個選項。它的功用是去計算不同類別的權重，並在計算平均數的時候考量到權重大小。簡單來講就是去計算加權的估計量。以 sex 為例，在沒有考量加權的情況下，每個 race 的權重都會是 0.2。而考量權重的話，每個 race 的權重則等於每個種族在全部樣本裡面所佔的比例。

若用 estimate 語法則程式如下：

proc glm data=anal;
      class race sex;
      model y1 = race sex / solution e;
      estimate 'male' intercept 1 sex 0 1;
      estimate 'female' intercept 1 sex 1 0;
      estimate 'female-male' sex 1 -1;
run;

寫法跟第一個程式範例一樣，但我會習慣加上「race 0.2 0.2 0.2 0.2 0.2」，只是說不加的話也沒關係，因為程式會自動用等同的權重去平均。寫到這邊只有兩點要強調：一是類別變數有幾層，estimate 語法內變數後面要跟的數字就有幾個。有人習慣將 1 後面可能還會出現的 0 都省略掉，但我強烈建議補上，因為多寫幾個 0 不會花太多時間，而且有可好處是可以把所有的 estimate 語法排列的工工整整，要 debug 也比較方便。二是腦筋一定要清楚知道變數的編碼順位，在寫之前最好再去對一下 Class Level Information 表格內的訊息。

當考慮交互作用項時，estimate 語法就沒有辦法簡化了：

proc glm data=anal;
      class sex race;
      model y1 = race sex race*sex / solution;
      estimate 'male'   intercept 1 sex 0 1 race .2 .2 .2 .2 .2  race*sex 0 0 0 0 0  .2 .2 .2 .2 .2;
      estimate 'female' intercept 1 sex 1 0 race .2 .2 .2 .2 .2  race*sex .2 .2 .2 .2 .2  0 0 0 0 0;
      estimate 'female-male' sex 1 -1 race*sex .2 .2 .2 .2 .2  -.2 -.2 -.2 -.2 -.2;
run;

在 sex 的點估計中，除了「race 0.2 0.2 0.2 0.2 0.2」一定要加上外，交互作用項 race*sex 後面要寫的數字向量變將整個程式複雜程度提昇了許多。由於 race 有五層，sex 有兩層，所以交互作用項就會是 5X2 = 10 層。因此開始寫這段數字向量時，先要知道到底該寫幾個數字進去，這是第一個要注意的事情。其次，將 sex 的向量和 race 的向量用 Kronecker product 的方式去算。如果不清楚 Kronecker product，我在這邊提供一個小程式碼：

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

自行替換 a 和 b 裡面的數字就可以算出兩個變數的交互作用項的向量。然後就直接從 output 視窗裡面複製貼上即可。在成對比較上面，交互作用的數字向量寫法也是一樣可以用這個小程式去算，只要把{0 1}改成{1 -1}即可。但在寫程式時，「race 0.2 0.2 0.2 0.2 0.2」因為在相減的過程中被扣掉了，所以不用寫上去。一個比較白痴的方法是把第一個 estimate 和 第二個 estimate 排好，然後下面減上面，自然就得到最後一個的數字向量。只是當你有大量變數要比較時，這種作法就有點慢了。

可是當考慮權重時，程式如下：

proc glm data=anal;
      class sex race;
      model y1 = race sex race*sex / solution e;
      estimate 'male' intercept 1  sex 0 1  race  .104 .142 .128 .047 .579 race*sex 0 0 0 0 0  .104 .142 .128 .047 .579;
      estimate 'female' intercept 1  sex 1 0 race  .104 .142 .128 .047 .579 race*sex .104 .142 .128 .047 .579  0 0 0 0 0;
      estimate 'female-male' sex 1 -1 race*sex .104 .142 .128 .047 .579 -.104 -.142 -.128 -.047 -.579;
      estimate 'female-male (off)' sex 1 -1 race*sex .105 .142 .128 .047 .579 -.105 -.142 -.128 -.047 -.579;
run;

紅色部分就是 race 在五個不同總族下的權重，這可以用 PROC FREQ 額外算出來。只要這個部分解決，後面交互作用項就簡單了，因為可以直接套用上面我提供的小程式，之後便是一連串的剪貼。兩兩比較的技巧也如同之前所述。最後一個標記為 'female-male (off)' 的程式，只是要提醒，在計算權重時，我們可能會使用到進位的權重，但進位後可能會造成無法估計的結果。本例中，由於進位的關係，race*sex 後面的數字總和雖然還是零，但是 race 權重的總和變成 1.001 而不是 1.000 了。因此你就會看到一段紅色錯誤訊息出現在 log 視窗告訴你這個事後比較無法估計。解決的方法沒有特別，就是手動去微調裡面的數據，東加一點西減一點，調到總合是 1 即可。
有時候事後比較的目的是要看某離散型的變數是否有線性趨勢，則程式如下：

proc glm data=anal;
      class sex educat;
      model y1 = sex educat / solution;
      lsmeans sex educat / stderr e;
      lsmeans sex educat / stderr e om;
      estimate 'educat linear (no divisor)' educat -2 -1 0 1 2;
      estimate 'educat linear (divisor=10)' educat -2 -1 0 1 2 / divisor=10;
      estimate 'educat linear (divisor=40)' educat -4 -2 0 2 4 / divisor=40;
      estimate 'educat by year (not centered)' educat 8 12 14 16 20;
      estimate 'educat by year (centered)' educat -6 -2 0 2 6 / divisor=80;
run;

此程式把 race 換成 educat，但一樣是五層。線性趨勢的事後檢定所使用的數字向量為「-2 -1 0 1 2」，這個比例大致上是不會改變的，當然不同數目的層次有不同的寫法，請自行孤狗。有時會有比較特殊的寫法，就比方'educat by year (not centered)'的寫法，不過情況比較少。大部分還是會寫成像'educat by year (centered)'所示，也就是讓數字向量內的數字是對秤的，而不是單調遞升或遞減的。至於 divisor 的用法，在此例沒有特殊意義，只是當做是分母把要用來估計的數字向量一次全部除掉，換句話說，以第二個'educat linear (divisor=10)'為例，你也可以寫成

estimate 'educat linear (divisor=10)' educat -0.2 -0.1 0 0.1 0.2;

但是 divisor 的真正威力，是在下面這個範例：

proc glm data=anal;
      class sex race;
      model y1 = race sex race*sex / solution e;
      estimate 'male (333)' intercept 1  sex 0 1  race  .333 .333 .333 race*sex 0 0 0  .333 .333 .333; 
      estimate 'female (333)' intercept 1  sex 1 0 race  .333 .333 .333 race*sex  .333 .333 .333  0 0 0; 
      estimate 'female-male (333)' sex 1 -1 race*sex  .333 .333 .333  -.333 -.333 -.333; 
      estimate 'male (334)' intercept 1  sex 0 1  race  .333 .333 .334 race*sex 0 0 0  .333 .333 .334; 
      estimate 'female (334)' intercept 1  sex 1 0 race  .333 .333 .334 race*sex  .333 .333 .334  0 0 0; 
      estimate 'female-male (334)' sex 1 -1 race*sex  .333 .333 .334  -.333 -.333 -.334; 
      estimate 'male (d)' intercept 3  sex 0 3  race  1 1 1  race*sex 0 0 0  1 1 1 / divisor=3; 
      estimate 'female (d)' intercept 3  sex 3 0 race  1 1 1  race*sex 1 1 1  0 0 0 / divisor=3; 
      estimate 'female-male (d)' sex 3 -3 race*sex 1 1 1   -1 -1 -1 / divisor=3; 
      estimate 'male (o)' intercept 849  sex 0 849  race  142 128 579 race*sex 0 0 0  142 128 579 / divisor=849; 
      estimate 'female (o)' intercept 849  sex 849 0 race  142 128 579 race*sex 142 128 579  0 0 0 / divisor=849; 
      estimate 'female-male (o)' sex 849 -849 race*sex 142 128 579 -142 -128 -579 / divisor=849; 
run;

此例中，race 不再是五層，而是三層！所以當要計算 sex 內男女分別的點估計量，SAS 裡面是不允許寫分數而是得直接寫小數上去，造成 race 後面的數字向量不能是{1/3  1/3  1/3}而是{0.333  0.333  0.333}。這種情況使得三個數字的總和變成 0.999 而不是 1.000，然後就會看到 SAS 的 log 視窗給你錯誤的訊息了。有一種解決方法是改成{0.333  0.333  0.334}，但一旦層數變多了，如七層、九層、十一層、二十一層....，則不但整條 estimate 又臭又長，每次都要去微調裡面的數字也很麻煩，所以此時就要靠 divisor 這個選項幫我們一次除盡所有數。首先，先把原本每個變數後面的數字都乘上 3，接著在最後面加上 / divisor = 3 即可。如果考量到 race 的權重，請回憶國小國中學過的最小公倍數。你可能會問，一對小數點怎麼去求最小公倍數呢？在此額外提供一個 SAS 小程式方便你的計算：

data _null_;
    x = lcm(a1,a2,....,aN);
    put x;
run;

把你算出來的權重全部丟到 lcm function 裡面，這個函數顧名思義就是最小公倍數 least common multiple 的縮寫，馬上就可以算出 divisor 該填什麼數字進去了。

以上就是很基本但卻很實用的 estimate 寫法。contrast 的用法一模一樣，所以沒有多做介紹。原文後面還有一些更進階的語法，但我個人的經驗是極少數的情況會使用到，因此大家熟練以上的程式就夠了。不過有時候離散變數內的層數很多，在建立 estimate 語法的時候，後面會跟著一長串的數字向量，非常驚人（恐怖）。原文作者建議將一些會一直重複用到的數字向量用 %let 先定義好，等到要用的時候再呼叫即可，比方說：

%let S_1    = 1 0 0 0 0 0 0 0 0 0;
%let S_2    = 0 1 0 0 0 0 0 0 0 0;
%let S_3    = 0 0 1 0 0 0 0 0 0 0;
%let S_4    = 0 0 0 1 0 0 0 0 0 0;
%let S_5    = 0 0 0 0 1 0 0 0 0 0;
%let S_6    = 0 0 0 0 0 1 0 0 0 0;
%let S_7    = 0 0 0 0 0 0 1 0 0 0;
%let S_8    = 0 0 0 0 0 0 0 1 0 0;
%let S_9    = 0 0 0 0 0 0 0 0 1 0;
%let S_10   = 0 0 0 0 0 0 0 0 0 1;
%let S_none = 0 0 0 0 0 0 0 0 0 0;
%let P_U    = .1 .1 .1 .1 .1 .1 .1 .1 .1 .1;
%let P_O    = 0.0666048524 0.0776077886 0.0888270746 
              0.0955339206 0.1017462525 0.1059496214 
              0.1124710246 0.1164580436 0.1168907433 
              0.1179106784 ;
  *** Intercept at the end of the BEFORE and at the beginning of the AFTER ***;
  estimate 'LS:1_Before'   intercept             1
                           decile                &S_1.
                           decile*tafter         &S_none.;      
  estimate 'LS:1_After'    intercept             1
                           decile                &S_1.
                           decile*tafter         &S_1.;   
  estimate 'LS:2_Before'   intercept             1
                           decile                &S_2.
                           decile*tafter         &S_none.;      
  estimate 'LS:2_After'    intercept             1
                           decile                &S_2.
                           decile*tafter         &S_2.;      
 . . . 
  estimate 'LS:10_Before'  intercept             1
                           decile                &S_10. 
                           decile*tafter         &S_none.;      
  estimate 'LS:10_After'   intercept             1
                           decile                &S_10.
                           decile*tafter         &S_10.;      
  estimate 'LS:Before_U'   intercept             1
                           decile                &P_U.
                           decile*tafter         &S_none.; 
  estimate 'LS:After_U'    intercept             1
                           decile                &P_U.
                           decile*tafter         &P_U.;  
  estimate 'LS:Before_O'   intercept             1
                           decile                &P_O.
                           decile*tafter         &S_none.;   
  estimate 'LS:After_O'    intercept             1
                           decile                &P_O.
                           decile*tafter         &P_O.;

試想，如果不用 %let 去額外定義巨集變數，而把那些 0 和 1 或 -1 等數字全部填回到 estimate 後面，則整串程式碼絕對會讓人暈眩。因此，我強烈建議使用 %let 來定義好 estimate 後面該用到的數字向量。此外，當要考量權重時，雖然可以自行算好寫上，但當層數多的時候（如本例），一個個複製貼上也是很麻煩，所以原文作者提供了一段程式碼讓 SAS 自動幫你把把這些權重貼到巨集變數裡面。

proc sort data=temp01(keep=patid_age decile) out=temp02 nodupkey;
  by patid_age decile;
run;
data _null_;
  set temp02;
  by patid_age decile;
  if not(first.patid_age and last.patid_age) 
then error 'ERROR: decile varies within patid_age';
run;
ods listing close;
ods output onewayfreqs = temp03(keep=decile frequency cumfrequency);
proc freq data = temp02;
  tables decile;
run;
ods listing;
data _null_;
  set temp03;
  call symput('N_'||left(trim(put(decile,2.))),frequency);
run;

data _null_;
  set temp03;
  if (decile eq 10) then call symput('Ntot',cumfrequency);
run;
data _null_;
  call symput('P_1', &N_1. / &Ntot.));
  call symput('P_2', &N_2. / &Ntot.);
  call symput('P_3', &N_3. / &Ntot.);
  call symput('P_4', &N_4. / &Ntot.);
  call symput('P_5', &N_5. / &Ntot.);
  call symput('P_6', &N_6. / &Ntot.);
  call symput('P_7', &N_7. / &Ntot.);
  call symput('P_8', &N_8. / &Ntot.);
  call symput('P_9', &N_9. / &Ntot.);
  call symput('P_10', &N_10. / &Ntot.);
run;
*** Macro variable for model specifications ***;
%let P_O = &P_1. &P_2. &P_3. &P_4. &P_5. &P_6. &P_7. &P_8. &P_9. &P_10.;

這段程式碼簡而言之，就是先用 PROC FREQ 把次數分配表的結果算出來，然後把這必要的變數（也就是 frequency 和 cumfrequency）存出來計算每一層的比例，再把這比例結果用 symput 函數存成不同的巨集變數，然後一次丟進 %let 裡面整合成一個巨集變數供 estimate 語法使用。不過我覺得原文作者寫的這段程式稍嫌複雜些，可能因為他用的資料結構比較複雜的因素。在一般的情況下，用 PROC FREQ 可以立刻求出類別變數內各層的比例：

proc freq data=xxx;
      tables varcat;
      ods output onewayfreqs = temp(keep=percent);
run;

用 ods output 另存的資料裡面有個 percent 變數就是每一層的比例。然後用 proc transpose 轉置：

proc transpose data=temp out=temp;
      var percent;
run;

假設有十層的話，轉置後的 percent 會被存在 COL1~COL10 裡面。然後用 symput 函數分別將他們存成不同的巨集變數：

data _null_;
      set temp;
      call symput('P_1', COL1);
      call symput('P_2', COL2);
      call symput('P_3', COL3);
      call symput('P_4', COL4);
      call symput('P_5', COL5);
      call symput('P_6', COL6);
      call symput('P_7', COL7);
      call symput('P_8', COL8);
      call symput('P_9', COL9);
      call symput('P_10', COL10); 
run;

最後再全部丟到 %let 後面用一個巨集變數整合他們：

%let P = &P_1 &P_2 &P_3;

如此一來，在寫 estimate 時需要用到這串數字向量時，輕鬆填上&P 即可！

CONTACT INFORMATION 
Your comments and questions are valued and encouraged.  Contact the author at:
 David J. Pasta
 Vice President, Statistics & Data Operations
ICON Clinical Research
 188 Embarcadero, Suite 200
 San Francisco, CA 94105
 +1.415.371.2111  
 david.pasta@iconplc.com
 www.iconplc.com