Permutation Tests in Statistics and Clinical Trials Development: Sample Code

 Description:

 1. Validate the permutation test table using the measured values. There are some missing values. You can adjust the dpl dataset for other imputed data 

  2. Use the planned allocation ratio of placebo to Active treatment as 2:1. Although the actual data may not perfectly match this ratio, treatment groups in simulations are assigned based on this planned allocation using total number of subjects with available data for change from baseline.

  3. Permutation strategy: 

• If the sample size is = 14, perform full permutations.

• If the sample size is > 14, perform 10,000 iterations.

 4. For the two-sample t-test, use the pooled variance method, which assumes equal variances between the two treatment groups.

 5. To run the analysis, adjust the dataset named dpl accordingly.

/*****Sample SAS code  ******/

** One sample test and inference:

** 0. Calculated Observed t value for change

** 1. Get the random sample between 0 and 1. if >=0.5 then let rand=1, if <0.5 then rand=-1.

** 2. Get change value to be abs(change)* rand (from #1) for one sample test.

** 3. For 2-ample test, if rand=-1 then set arm to be placebo, if rand=1 then set arm=casi.  in addition, the seed should

**    change for each test. I am using 1 for casi one sample test, 2 for placebo one sample test, and 3 for two arm tests.

** 4. P-value for re-randomization test:

** (1) Get the t statistic for each of the 10000 permutations

** (2) Count how many t values >= observed t value and then divided by 10000=p1; How many t values <= observed t value and then divided by 10000=p2

** (3) Pvalue = min of (p1, p2)*2

**************************************************************************************************************************/

libname test "xxx";

data admi;

set test.admi;run;

proc contents data=admi;run;

data dpl;

set admi;

run;

/*sample size for each group and total */

proc sql noprint;

select count(distinct usubjid) into : n1  from dpl(where=(TR01PG1="Active treatment"));

select count(distinct usubjid) into : n2  from dpl(where=(TR01PG1="Placebo"));

select count(distinct usubjid) into : n3  from dpl(where=(TRT02P="Total"));

quit;

%let n1=&n1;

%let n2=&n2;

%let n3=&n3;

%put n1=&n1 n2=&n2 n3=&n3;

* subset three groups;

data s p t;

set dpl;

output t;

if TR01PG1="Active treatment" then output s;

else if TR01PG1="Placebo" then output p;

run; 

/*One-sample test: within group */

/*macros for permutation test;*/

/*sample size n<14;*/

%macro perm(data, TRT01P, order);

data d2;

set &data;

chgr=-(chg); *baseline-week xx;

n=_n_;

subj="subj"||strip(put(n, best.));

run;

/*creating the structure of permutations;*/

%macro comb(j);

data comb;

%do i=1 %to %cmpres(&&n&j); *n1, n2 or n3;

do subj&i= 0, 1;

%end;

output;

%do i=1 %to %cmpres(&&n&j);

end;

%end;

run;

%mend;

%comb(&order); 

data comb;

set comb;

replicate=_n_;

run;

proc transpose data=comb out=comb2(rename=(_name_=subj col1=c));

by replicate;

run;

/* merge the permutation structure with original data;*/

proc sql;

    create table D3 as

    select 

       b.replicate, b.c, a.*

    from COMB2 as b

    left join D2 as a

    on b.subj = a.subj

    order by replicate,n;

quit;

data d4;

set &data(in=a) d3(in=b);

if a then do; replicate=0; value=chg; end;

if c=0 then value=chg; *original change;

else if c=1 then value=chgr; *baseline and week 48 swapped;

run;

%mend perm;

/*sample size n>=14;*/

%macro rand(data, seed, repeats=10000);

data d2;

retain seed &seed; drop seed;

set &data;

do replicate=1 to &repeats;

call ranuni(seed, randu);

output;

end;

run;

proc sort data=d2; by replicate randu; run;

data d4;

set &data(in=a) d2(in=b);

if a then do; replicate=0; value=chg; end;

if not missing(randu) then do;

if randu < 0.5 then c=-1;/*change*/

if randu >= 0.5 then c=1;

value=abs(chg)*c;

end;

run;

%mend rand;

%macro results(data, TR01PG1, order, seed);

%if &&n&order >=14 %then %do; /*if sample sie >=14, derive d4 using rand macro */

%rand(&data, &seed); 

%end;

%else %do;

%perm(&data, &TR01PG1, &order);/*if sample sie <14, derive d4 using perm macro*/

%end; 

ods exclude all;

proc ttest data=d4 sides=2 h0=0;

by replicate;

    var value;

  ods output ttests=PROBF;

run;

ods exclude none;

/*count total and calculate p-value;*/

data pf&order.(rename=(col1=&TR01PG1));

set probf;

by variable replicate;

format tvalue level 13.10;

retain level count tot counte;

if first.variable and replicate=0 then do ;

level=tvalue;

tot=0;

count=0;

counte=0;

end;

else do ;

tot+1;

if round(tvalue,.0000000001)>=round(level,.0000000001) then count+1;

if round(tvalue,.0000000001)=round(level,.0000000001) then counte+1;

end;

if last.variable;

  length hd $ 80 col1 $20;

  countuse=min(count, (tot-count+counte));

col1=put(round(min(1,countuse/(tot)*2), 0.0001),7.4); *2-sided;

hd="P-value";

keep col1 hd ;

run;

%mend results;

%perm(p,Placebo, 2)

%results(s, Active treatment, 1, 10001);

%results(p, Placebo, 2, 20001);

/*Two-sample test: between group*/

/*missing values, still use allocation ratio 1:2 */

macro rand(data,  seed,repeats=10000);

    /* Step 1: Count total number of subjects */

    proc sql noprint;

        select count(*) into :n_total from &data;

    quit;

    %let n1 = %sysevalf(&n_total / 3, floor);  /* Group 1 size (1/3) ,  the result is rounded down to the nearest integer.*/

 

    /* Step 2: Create re-randomized assignments for n permutations */

    data d2;

        retain seed &seed; drop seed;

    set &data;

        do replicate = 1 to &repeats;

            call ranuni(seed, randu);

            output;

        end;

    run;

    proc sort data=d2;by replicate randu;run;

    /* Step 3: relabel the treatment group based on sorted order */

data d3;  

    length dosegroup $50;

        set d2;

        by replicate;

        retain count;

        if first.replicate then count = 0;

        count + 1;

        if count <= &n1. then do; dosegroup = "Placebo";value=chg;end;/*need to change here */

        else do; dosegroup = "Active treatment";value=chg;end;

    run;

data d4;

set &data(in=a) d3(in=b);

if a then do; replicate=0; value=chg;dosegroup=TRT01P; end;run;

%mend rand;

%macro results_grp(data, TR01PG1, order, seed);

%if &&n&order >=14 %then %do; /*if sample sie >=14, derive d4 using rand macro */

%rand(&data, &seed); 

%end;

ods exclude all;

proc ttest data=d4 sides=2 h0=0;

by replicate;

class dosegroup;

    var value;

  ods output ttests=PROBF;

run;

ods exclude none;

*count total and calculate p-value;

data pf&order.;

set probf;

by variable replicate;

format tvalue level 13.10;

retain level count tot counte;

if Method="Pooled";/*Assume equal variance*//**************caution**************************/

if first.variable and replicate=0 then do ;

level=tvalue;

tot=0;

count=0;

counte=0;

end;

else do ;

tot+1;

if round(tvalue,.0000000001)>=round(level,.0000000001) then count+1;

if round(tvalue,.0000000001)=round(level,.0000000001) then counte+1;

end;

run;

data pf&order.(rename=(col1=&TR01PG1));

set pf&order.;

by variable replicate;

if last.variable;

  length hd $ 80 col1 $20;

  countuse=min(count, (tot-count+counte));

col1=put(round(min(1,countuse/(tot)*2), 0.0001),7.4); *2-sided;

hd="P-value";

keep col1 hd ;

*keep col1 hd countuse tot;

run;

%mend results_grp;

%results_grp (t, Active treatment, 3, 30001);