【SAS】permutation test 語法邏輯

Ping
·
(修改过)
·
IPFS
·
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 ; 

CC BY-NC-ND 4.0 授权

喜欢我的作品吗?别忘了给予支持与赞赏,让我知道在创作的路上有你陪伴,一起延续这份热忱!

Ping最近... 一直寫SAS程式 | 持續增肌減脂、桌球訓練 | 剛考到健身C級證照 | 將要看《相聲瓦舍-狐說》《表演工作坊-那一年,我們下凡》
  • 来自作者
日常堆積成立體的 「Ping」
2 篇作品

2025.04.11-意外且小有成就感的點

日常堆積成立體的 Ping