/*****************************************************
Wei-Johnson (1985, Biometrika)
data format need to be one recored per patient
please see the contents of sample dataset

Arguments for %weijohnson() :
1) indata -- sas data set name (one recored per patients)
2) id  ----- subject identification
3) arm ----- treatment group to compare (active = 1 and control=2, direction: code 1 - code 2)
4) vars ---- repeated measures (eg., var1 var2 var3 var4)
5) strata -- stratification variable (option)
6) kernel -- specify either TTEST or Wilcoxon
             TTEST performs test based on T-test
             Wilcoxon performs test based on Wilcoxon test
7) iter ---- Number of permutaions to get the null distribution (1000 would be recommended)
8) seed_1 -- seed for permutation
9) delta_h0 --- will be used for interval estimation. please do not change the default=0
******************************************************/

*==================================================== ;
%macro weijohnson(indata, id=, arm=, vars=, strata=., kernel=, iter=1000, seed_1=12345, delta_h0=0) ;

%if &strata eq .  %then %do ; data inds ; set &indata ; st=1 ;                      run ; %end ;
%if &strata ne . and &strata ne st
                  %then %do ; data inds ; set &indata ; st=&strata ; drop &strata ; run ; %end ;
%if &strata eq st %then %do ; data inds ; set &indata ;                             run ; %end ;
proc sort data=inds ; by st ; run ;

*---- observed statitics ----- ;
%wjstat(inds , id=&id, arm=&arm, vars=&vars, strata=st, kernel=&kernel, delta_h0=&delta_h0, arm_org=&arm) ;
data observed ; set out2 (rename=(stat=observed)) ; diff=0 ; run ;

ods listing close ;

*---- null distribution (via perumunation) -----;
data armdat ; set inds (keep=&arm) ; run ;

data permuted; set _null_ ; stat=.; diff=. ; run ;

%let rep=0;
%do %while (&rep<&iter);
    %let rep=%eval(&rep+1);

    *---- permuation ---;
    data wk1;
*     set inds (drop=&arm);
      set inds ; arm_org=&arm ; drop=&arm ;
      retain seed_1 %eval(&seed_1+&rep);
      call ranuni(seed_1, temp);
    run;
    proc sort data=wk1 ; by st temp ; run ;
    data wk2 ; merge wk1(drop=seed_1 temp) armdat ; run ;

    *---- statistics with permuted data ---;
    %wjstat(wk2, id=&id, arm=&arm, vars=&vars, strata=st, kernel=&kernel, delta_h0=&delta_h0, arm_org=arm_org) ;
        data out2 ; set out2 ; diff=0 ; run ;
    proc append base=permuted data=out2 ; run ;

%end ;

data tmp ; merge observed permuted ; by diff ;
   if stat ge observed then diff=1 ; run ;

ods output Summary=Summary ;
proc means data=tmp ; var diff ; run ;
ods output close ;

data pval ; set Summary ;
 pval1_upp=diff_Mean ;
 pval1_low=1-diff_Mean ;
 if diff_Mean ge 0.5 then pval2= (1-diff_Mean)*2 ;
 if diff_Mean lt 0.5 then pval2=  diff_Mean*2 ;
 label pval1_low="One-sided p-value (H1:Delta<0)" ;
 label pval1_upp="One-sided p-value (H1:Delta>0)" ;
 label pval2="Two-sided p-value" ;
 keep pval1_low pval1_upp pval2 ;
run ;

ods listing ;
title "Wei-Johnson (&kernel)" ; run ;
proc print label data=pval ; run ; title; run ;

%mend ;
*==================================================== ;

*==================================================== ;
%macro wjstat(indata, id=, arm=, vars=, strata=., kernel=, delta_h0=, arm_org=) ;

%if &strata eq .  %then %do ; data tmp0 ; set &indata ; st=1 ;                      run ; %end ;
%if &strata ne . and &strata ne st
                  %then %do ; data tmp0 ; set &indata ; st=&strata ; drop &strata ; run ; %end ;
%if &strata eq st %then %do ; data tmp0 ; set &indata ;                             run ; %end ;

proc sort data= tmp0 out=tmp1    ; by &id &arm st ; run ;
data tmp1 ; set tmp1 ; keep &id &arm &vars st ; run ;

proc transpose data=tmp1 out=tmp2  ; by &id &arm st ; run ;
data long ; set tmp2 (rename=(col1=res)) ; run ;

proc sort  data=long ; by _name_ st ; run ;

*--- null ---;
data long ; set long ; if &arm_org eq 1 then res = res - &delta_h0 ; run ;  

*----------;
* T- TEST  ;
*----------;
%if &kernel=TTEST %then %do ;
ods output Statistics=Statistics ;
proc ttest data=long ; class &arm ; var res ; by _name_ st ; run ;
ods output close ;
data out1 ; set Statistics ; by _name_ st ;
 variance=stderr**2 ;
 if last.st or last._name_ then output ;
 keep _name_ mean stderr variance ;
run ;
%end ;

*----------;
* Wilcoxon ;
*----------;
%if &kernel=Wilcoxon %then %do ;

 proc sort data=long ; by _name_ st ; run ;
 proc rank data=long out=long_rank percent ; var res ; by _name_ st ; run ;

 ods output Statistics=Statistics ;
 proc ttest data=long_rank ; class &arm ; var res ; by _name_ st ; run ;
 ods output close ;
 data out1 ; set Statistics ; by _name_ st ;
  variance=stderr**2 ;
  if last.st or last._name_ then output ;
  keep _name_ mean stderr variance ;
 run ;
%end ;

*--------------------------------------;

data out2 ; set _null_ ; run ;
data uno1 ; set out1 ;
  retain tmp 0 ;
  junk=1;
  tmp=tmp + 1/variance ; run ;
data uno2 ; set uno1 ; by junk ;
  totalweight=tmp ;
  if last.junk then output ;
  keep junk totalweight ; run ;
data uno3 ; merge uno1 uno2 ; by junk ; retain stat 0 ;
  weight = (1/variance)/totalweight ;
  stat = stat + mean*weight ;
data out2; set uno3 ; by junk ;
  if last.junk then output;
  keep stat ;
run ;

%mend ;
*==================================================== ;