topic Re: Permutation test on correlations in SAS Programming

Permutation test on correlations

Martin_blaedel — Tue, 20 Feb 2024 18:57:54 GMT

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;

Re: Permutation test on correlations

PaigeMiller — Tue, 20 Feb 2024 19:08:01 GMT

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.

Re: Permutation test on correlations

Martin_blaedel — Tue, 20 Feb 2024 19:13:28 GMT

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.

Re: Permutation test on correlations

Martin_blaedel — Tue, 20 Feb 2024 20:15:14 GMT

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;

Re: Permutation test on correlations

PaigeMiller — Tue, 20 Feb 2024 20:24:56 GMT

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.

Re: Permutation test on correlations

Martin_blaedel — Tue, 20 Feb 2024 20:28:47 GMT

Thanks - those are great ideas.

I will start by looking into that documentation.