/* Mantel correlation with randomisation test */ /* Author: Steve Palmer, University of Aberdeen */ /* Last updated: 23 October 2012 */ %macro manteltest(indata,Nrows,row,col,X,Y,Nperms,tline); * Parameters: * indata - the name of the dataset holding the two full square matrices to be correlated, 1 row per row*column combination (diagonal elements may be omitted); * Nrows - the number of rows (and columns) in the matrix; * row - the name of the column holding the row number; * col - the name of the column holding the column number; * X - the name of the column holding the cell values of the first matrix; * Y - the name of the column holding the cell values of the second matrix; * Nperms - the number of random permutations; * tline - a single digit indicating the title line on which macro output titles are to be written (e.g. 3 for title on the third line); * Extract the required data from the input dataset and set row and column numbers to be sequential; %if &row = row or &col = col %then %do; data r00; set &indata; _row = &row; _col = &col; Xdata = &X; Ydata = &Y; keep _row _col Xdata Ydata; run; %let row = _row; %let col = _col; %put Value of row changed to &row and col changed to &col; %end; %else %do; data r00; set &indata; Xdata = &X; Ydata = &Y; keep &row &col Xdata Ydata; run; %end; proc sort data=r00; by &row &col; run; proc means data=r00 noprint nway; class &row; output out=rows n=nrows; run; data rows; set rows; if _n_ eq 1 then row = 1; else row = row+1; retain row; keep &row row; run; data r00row; merge r00 rows; by &row; run; proc sort data=r00row; by &col &row; run; proc means data=r00 noprint nway; class &col; output out=cols n=ncols; run; data cols; set cols; if _n_ eq 1 then col = 1; else col = col+1; retain col; keep &col col; run; data r0; merge r00row cols; by &col; run; proc sort data=r0; by row col; run; *proc print data=r0; %let tt1 = Input matrix for &X - not randomised; %let qtt1 = %nrbquote('&tt1'); title&tline %unquote(&qtt1); proc transpose data=r0 out=r0x prefix=c; by row &row; id col; var Xdata; run; proc print data=r0x noobs; var row &row c1-c&Nrows; run; %let tt2 = Input matrix for &Y - to be randomised; %let qtt2 = %nrbquote('&tt2'); title&tline %unquote(&qtt2); proc transpose data=r0 out=r0y prefix=c; by row &row; id col; var Ydata; run; proc print data=r0y noobs; var row &row c1-c&Nrows; run; * Generate random permutations of the rows and columns; data r1; %do p = 1 %to &Nperms; do row = 1 to &Nrows; perm = &p; rrr = ranuni(1); output; end; %end; run; proc rank data=r1 out=rowrand; by perm; var rrr; ranks row2; run; *proc print data=rowrand; run; data colrand; set rowrand; col = row; col2 = row2; drop row row2; run; *proc print data=colrand; run; data rowrand2; set rowrand; drop rrr; %do c = 1 %to &Nrows; col&c = &c; %end; run; proc transpose data=rowrand2 out=rowrand3; by perm row row2; var col1-col&Nrows; run; data rowrand3; set rowrand3; col = col1; drop col1 _name_; run; *proc print data=rowrand3; run; proc sort data=rowrand3; by perm col; run; data rowrand4; merge rowrand3 colrand; by perm col; run; *proc print data=rowrand4; run; * Combine permutations with the X data (not permuted); proc sort data=rowrand4; by row col perm; run; data r2; merge r0 rowrand4; by row col; keep perm row col row2 col2 Xdata; run; *proc print data=r2; run; * Combine permutations with the Y data (permuted); data ydata; set r0; row2 = row; col2 = col; keep row2 col2 Ydata; run; proc sort data=ydata; by row2 col2; run; proc sort data=r2; by row2 col2 perm; run; data r3; merge ydata r2; by row2 col2; run; * Append the observed data; data r4; set r0 r3; if perm eq . then perm = 0; keep perm row col row2 col2 Xdata Ydata; run; proc sort data=r4; by perm row col; run; *proc print data=r4; run; * Correlation for each permutation; proc corr data=r4 noprint outp=mcorr01; by perm; var Xdata; with Ydata; run; data mcorr02; set mcorr01; r = Xdata; keep perm r; if _type_ eq 'CORR' then do; output; if perm eq 0 then do; if r lt 0 then rnegative = 1; else rnegative = 0; call symput('rneg',left(put(rnegative,1.0))); end; end; run; %if &rneg %then %do; %put Observed correlation coefficient is negative; proc rank data=mcorr02 out=mcorr03; var r; ranks rrank; run; %end; %else %do; %put Observed correlation coefficient is positive; proc rank data=mcorr02 out=mcorr03 descending; var r; ranks rrank; run; %end; data mcorr03; set mcorr03; P = rrank / (&Nperms + 1); run; %let ttt = Mantel correlation of &X and &Y, probability from &Nperms randomisations; %let qttt = %nrbquote('&ttt'); title&tline %unquote(&qttt); proc print data=mcorr03 noobs; where perm eq 0; var r P; run; title&tline; %mend;