/*******************************************************************************************************/
/*                                                                                                     */
/*  This macro implements a test described in                                                          */
/*                                                                                                     */
/*  Piepho, H.P. (2005): Permutation tests for the correlation among genetic distances and measures    */
/*  of heterosis. Theoretical and Applied Genetics (to appear)                                         */
/*                                                                                                     */
/*  The macro performs Mantel's permutation test of the correlation of two association matrices.       */
/*  The association matrices may refer to distances, measures of heterosis, specific combining ability */
/*  or other. The association matrices must be provided as SAS datasets and organised in upper         */
/*  triangular form. Variables/columns must be labelled col1-coln, where n is the number of columns    */
/*  of the association matrix. No other variables are allowed. A simple example for such a dataset     */
/*  is as follows:                                                                                     */
/*                                                                                                     */
/*  data a1;                                                                                           */
/*  input col1-col4;                                                                                   */
/*  datalines;                                                                                         */
/*  . 7 2 3                                                                                            */
/*  . . 3 1                                                                                            */
/*  . . . 4                                                                                            */
/*  ;                                                                                                  */
/*                                                                                                     */
/*  The call of the macro requires spefication of the two datasets containing the two association      */
/*  matrices (d1=, d2=). In addition, one may specify the number of permutations (default: p_max=1000) */
/*  and the seed of the random number generator (default: seed=-1, uses clock to initialize seed).     */
/*                                                                                                     */
/*******************************************************************************************************/

%macro permut(d1=,d2=,p_max=1000, seed=-1);

data &d1;
set &d1 nobs=n;
call symput('n',trim(left(put(n,8.))));
run;

data append;
array col col1-col&n;
do i=1 to &n;
  col[i]=.;
end;
drop i;

data d1a;
set &d1 append;
run;

data d2a;
set &d2 append;
run;

proc iml;
n=&n+1;
order=j(1,n,0);
use d1a;
read all var _num_ into d1;
close d1a;
use d2a;
read all var _num_ into d2;
close d2a;


/*compute r*/
d1_mean=0;
d2_mean=0;
d1_ssq=0;
d2_ssq=0;
r=0;
m=0;
do i=1 to n-1;
  do ii=i+1 to n;
    d1[ii,i]=d1[i,ii];
    d2[ii,i]=d2[i,ii];
    prod=d2[i,ii]*d1[i,ii];
	if prod = . then do;
	end;
	else do;
	   m=m+1;
       r=r+prod;
	   d1_mean=d1_mean+d1[i,ii];
	   d2_mean=d2_mean+d2[i,ii];
	   d1_ssq=d1_ssq+(d1[i,ii])**2;
	   d2_ssq=d2_ssq+(d2[i,ii])**2;
	end;
  end;
end;

d1_mean=d1_mean/m;
d2_mean=d2_mean/m;

s1=sqrt( (d1_ssq-d1_mean**2*m)/(m-1));
s2=sqrt( (d2_ssq-d2_mean**2*m)/(m-1));
r=(r-d1_mean*d2_mean*m)/s1/s2/(m-1);

/*do permutations*/
d2_perm=d2;
pval=0; pval2=0;
do p=1 to &p_max;
  do i=1 to n;
    order[i]=ranuni(&seed);
  end;
  order=rank(order);
  do i=1 to n-1;
    do ii=i+1 to n;
      d2_perm[i,ii]=d2[order[i],order[ii]];
      d2_perm[ii,i]=d2[order[ii],order[i]];
	end;
  end;

  d1_mean=0;
  d2_mean=0;
  d1_ssq=0;
  d2_ssq=0;
  r_perm=0;
  m=0;
  do i=1 to n-1;
    do ii=i+1 to n;
	  prod=d2_perm[i,ii]*d1[i,ii];
	  if prod = . then do;
	  end;
      else do; 
        r_perm=r_perm+prod;
		m=m+1;
 	    d1_mean=d1_mean+d1[i,ii];
	    d2_mean=d2_mean+d2_perm[i,ii];
	    d1_ssq=d1_ssq+(d1[i,ii])**2;
	    d2_ssq=d2_ssq+(d2_perm[i,ii])**2;
	  end;
	end;
  end;
  d1_mean=d1_mean/m;
  d2_mean=d2_mean/m;
  s1=sqrt( (d1_ssq-d1_mean**2*m)/(m-1));
  s2=sqrt( (d2_ssq-d2_mean**2*m)/(m-1));
  r_perm=(r_perm-d1_mean*d2_mean*m)/s1/s2/(m-1);
  if r_perm>=r then pval=pval+1; 
  if abs(r_perm)>=abs(r) then pval2=pval2+1; 

end;

pval=(pval+1)/(&p_max+1);
pval2=(pval2+1)/(&p_max+1);

print 'one-sided p-value' pval;
print 'two-sided p-value' pval2;

run;
quit;

%mend;


/*Example*/
data a1;
input col1-col7;
datalines;
. 0.29 0.13 0.21 0.72 0.22 0.44
. .    0.64 0.42 0.38 0.24 0.43
. . .       -0.02 0.54 -0.05 0.30
. . . .           0.45 -0.18 0.27
. . . . .               0.39 0.55
. . . . . .                  0.31
;

data a2;
input col1-col7;
datalines;
. 0.294 0.234 0.226 0.248 0.247 0.239
. .     0.301 0.257 0.224 0.247 0.256
. . .         0.220 0.269 0.253 0.239
. . . .             0.217 0.214 0.220
. . . . .                 0.213 0.232
. . . . . .                     0.222
;

%permut(d1=a1,d2=a2, p_max=10000);