【SAS】permutation test 語法邏輯

permutation test 置換測試-->也透過原樣本重複多次的隨機排列，來驗證原樣本在總體樣本可能的情況。

步驟：

1. 計算原樣本統計結果。

2. 將原樣本劑型隨機排列(EX. 分組的隨機排列)，並計算排列後的統計結果。

3. 重複step 2.的步驟多次(Simulation) (EX. 10000次)。

4. 將每次simulation 結果與原樣本統計結果比較，如果原結果>simulation 結果，設定新欄位變項=1。

5. 計算置換測試後的p-value=simulation 中新欄位=1的個數/simulation 次數。

SAS CODE (以ANCOVA檢定 比較LSMEANS DIFF為例):

1. /* 步驟1：首先進行原始 ANCOVA 分析，取得基準 lsmeansd ;*/

   ods output diffs=diffs ; /*匯出需要的統計結果*/

   proc mixed data=a1 ;

   class trt (ref='R') ;

   model iop_change= v1 trt ;

   lsmeans trt /pdiff cl;

   run ;

   *將原樣本跑出的統計結果設定成一個巨集，後面步驟會使用到 ;

   data null; /*統計結果: LSMANS DIFF*/

   set diffs;

   if trt = 'T' then original_meandiff = estimate ;

   call symputx('original_meandiff', original_meandiff);

   run;

2. /* 步驟2：進行 permutation test，重複多次 (例如 1000次)*/ ;

   %let num_permutations = 1000;

   proc printto log="路徑" ; /*將log輸出，避免sas log 頁面不足*/

   %macro aa;

   %do iteration = 1 %to &num_permutations;

   proc surveyselect data=a1 out=perm_data group=2 /*重排資料的 group 標籤 ; /seed=&iteration /設定種子*/ noprint ;

   run;

   data perm_data ;

   set perm_data ;

   if GroupID =1 then trt="T" ; /*根據隨機後的分組Recoding*/

   else if GroupID =2 then trt="R" ;

   run ;

   *simulation 樣本進行 ANCOVA 分析，並提取每次重排後的統計結果 ;

   ods output diffs=diffs ;

   ods results off ; /不呈現result結果，避免佔用記憶體*/

   proc mixed data=perm_data ;

   class trt (ref='R') ;

   model IOP_change= v1 trt ;

   lsmeans trt /pdiff cl;run ;

   * 把每次 permutation 的 統計結果存入一個彙整檔案 ;

   %IF &iteration=1 %THEN %DO ;

   data perm_meandiff;set diffs ;

   iteration = &iteration;

   keep iteration estimate ;

   run;

   %END ;

   %ELSE %DO ;

   data estimate2;

   set diffs;

   iteration = &iteration;

   keep iteration estimate;

   run;

   data perm_meandiff;

   set perm_meandiff estimate2; /*疊上新的simulation結果*

   run;

   %END ;

   %end;

   %mend ;

   %aa ;

   proc printto ; /*log 輸出至log頁面*/*建議permutation 的結果存入永久檔;

   data aa.Perm_meandiff ;

   set Perm_meandiff ;

   run ;

3. /* 步驟3：繪製所有simulation結果的直方圖*/

   ods select all ; /*呈現所有輸出結果*/

   ods results on ; /*result頁面呈現結果*/

   proc univariate data=aa.perm_meandiff;var estimate;histogram / normal;

   inset mean std / position=ne;

   run;  

4. /* 步驟4：比較原始 meandiff 與 permutation meandiff 分佈 */

   data aa.perm_meandiff;

   set aa.perm_meandiff;

   if estimate>&original_meandiff then p=1 ;

   else p=0 ;

   run ;

5. /*步驟5：計算 permutation test 結果*/

   proc sql ;

   select sum(p)/count(*) as pvalue

   into :permutationp

   from aa.perm_meandiff;

   run ;

   %put Original meandiff : &original_meandiff ;

   %put p-value distribution from permutation test: &&permutationp ;