Permutation test on correlations

[Martin_blaedel]

Hi 

I am trying (and struggling) to understand SAS syntax. I wanted to do a permutation test on the Pearsons correlation coefficient, and copied the below code(from the SAS page) into my editor. 

I thought I had it figured out, but nothing happens, or at least I don't get an output.

Can you tell me why?

the below code will in my understanding only do the permutation test on teh tobacco and cancer variable. But what if I would like to run the same test with on the correlation between alcohol and Tobacco and again on cancer and alcohol?

Can I do that fairly simple?

Thanks a lot 

OPTIONS LINESIZE=95;

DATA   test;
INFILE 'C:\Users\XXXXXXXXXXX.PRN'  FIRSTOBS=2;
INPUT  Cancer Tobacco Alcohol;

PROC PRINT;
RUN;

proc iml;
use test;
   read all var "tobacco" into X;
   read all var "cancer" into Y;
close;
 
start CorrCoeff(X, Y, method="Pearson");
   return ( corr(X||Y, method)[1,2] );   /* CORR returns a matrix; extract the correlation of the pair */
finish;
 
/* permutation test for correlation between two variables */
call randseed(12544);        /* set a random number seed */
B = 10000;                   /* total number of permutations */
PCorr = j(B, 1, .);          /* allocate vector for Pearson correlations */
PCorrObs = CorrCoeff(X, Y);  /* correlation for original data */
PCorr[1] = PCorrObs; 
do i = 2 to B;
   permY = ranperm(Y)`;      /* permute order of Y */
   PCorr[i] = CorrCoeff(X, permY); /* correlation between X and permuted Y */
end;

[PaigeMiller]

Are there errors in the log? 

What output are you expecting? You would need a PRINT statement to see your data sets written to the output window or results window; you would need to use the CREATE statement to create a SAS data set with results.

[Martin_blaedel]

No errors on the log. 

Yeah, that is a good question I think I was expecting perhaps a permutation distribution curve, and/or a table with the p value. 

But I kind of suspect what you say about adding a print statement could solve some of the problem.

I tried at some point but perhaps I did it wrong. will try again.  

[Martin_blaedel]

Hi again

Now I have corrected the code and i get my histogram. The code looks like the below now. However, I would like to control the incrementation on the X axis on the permutation distribution graph.

Do you by any chance know how to do that?

proc iml;
use test;
   read all var "alcohol" into X;
   read all var "cancer" into Y;
close;
 
start CorrCoeff(X, Y, method="Pearson");
   return ( corr(X||Y, method)[1,2] );   /* CORR returns a matrix; extract the correlation of the pair */
finish;
 
/* permutation test for correlation between two variables */
call randseed(12544);        /* set a random number seed */
B = 10000;                   /* total number of permutations */
PCorr = j(B, 1, .);          /* allocate vector for Pearson correlations */
PCorrObs = CorrCoeff(X, Y);  /* correlation for original data */
PCorr[1] = PCorrObs; 
do i = 2 to B;
   permY = ranperm(Y)`;      /* permute order of Y */
   PCorr[i] = CorrCoeff(X, permY); /* correlation between X and permuted Y */
end;

print PCorrObs;
run;


pVal = sum( abs(PCorr)>= abs(PCorrObs) ) / B;   /* proportion of permutations with extreme corr */
print pVal;
 
title "Permutation Test with H0:corr(alcohol,Cancer)=0";
title2 "Distribution of Pearson Correlation (B=1000)";
refStmt = "refline " + char(PCorrObs) + " / axis=x lineattrs=(color=red);";
call histogram(PCorr) other=refStmt label="Pearson Correlation";

run;
quit;

[PaigeMiller · Accepted solution]

I would check the documentation for PROC IML and the HISTOGRAM function.

You could also create an output data set from your IML, and then run that through PROC SGPLOT, which will give you many options for controlling the appearance.

[Martin_blaedel]

Thanks - those are great ideas. 

I will start by looking into that documentation.