/* SAS program to accompany the article 
   "Fisher's transformation of the correlation coefficient"
   by Rick Wicklin, published 20SEP2017 on The DO Loop blog
   http://blogs.sas.com/content/iml/2017/09/20/fishers-transformation-correlation.html

   This program:
   1) Graphs the Fisher transformation and the inverse transforamtion
   2) Simulates bivariate normal samples of size N from MVN(0,Sigma) where 
      Sigma={1 rho, rho 1}. Compute each sample correlation and plot
      a histogram of the sample correlations to create an approximate 
      sampling distribution.
   3) Transform each sample correlation and plot the transformed sampling distributions.
   4) Show how to use the FISHER option in PROC CORR in SAS
*/

/* Fisher's normalizing transformation 
   z = arctanh(r)
*/
title;footnote;
data artanh;
do r = -1 to 1 by 0.01;
   z = artanh(r);
   output;
end;
run;
ods graphics / width=240px height=400px;
title "Fisher's Z Transformation";
title2 "z = arctanh( r )";
proc sgplot data=artanh noautolegend;
   refline 0 / axis=x; refline 0 / axis=y;
   xaxis grid; yaxis grid min=-2 max=2 values=(-2 to 2 by 0.5);
   *lineparm x=0 y=0 slope=1 / lineattrs=(color=gray pattern=dashed);
   series x=r y=z / lineattrs=(thickness=2);
run;

data tanh;
do z = -2 to 2 by 0.01;
   r = tanh(z);
   output;
end;
run;
ods graphics / width=360px height=240px;
title "INVERSE Z Transformation";
title2 "r = tanh(z)";
proc sgplot data=tanh noautolegend;
   refline 0 / axis=x; refline 0 / axis=y;
   yaxis grid; xaxis grid min=-2 max=2 values=(-2 to 2 by 0.5);
   series x=z y=r / lineattrs=(thickness=2);
   *lineparm x=0 y=0 slope=1 / lineattrs=(color=gray pattern=dashed);
run;

/**********************************************/

/* Simulate bivariate normal samples of size N with correlation rho */
%let N = 20;             /* sample size */
%let NumSamples = 2500;  /* number of simulated samples */
proc iml;
N = &N;              /* sample size */
NumSamples = &NumSamples;
rho = .; ID = .; x = .; y = .;
create SimCorr var {"rho" ID X Y};  /* define 4 numerical output vars */

call randseed(123456);
ID = colvec( repeat( T(1:NumSamples), 1, N) );
Sigma = I(2);
rhoList = do(0.2, 0.8, 0.2);
do i = 1 to ncol(rhoList);
   Sigma[{2 3}] = rhoList[i];
   Z = randnormal(N*NumSamples, {0 0}, Sigma);
   x = Z[,1]; y = Z[,2]; 
   rho = j(N*NumSamples, 1, rhoList[i]);
   append;
end;
close;
quit;

proc corr data=SimCorr noprint
          outP=OutCorr(where=(_TYPE_="CORR") rename=(y=Corr));
   by rho ID;
   var x y;
run;

data DistribCorr;      /* Retain R[1,2] cell only */
   set OutCorr;
   if mod(_N_,2)=1;   /* Keep odd-numbered rows; delete even rows */
   /* apply Fisher's z transformation */
   z = artanh(Corr);
   z0 = artanh(rho);
   drop _NAME_ x;
   label Corr = "Sample Correlation (r)"   z = "z = arctanh(r)";
   if _N_=1 then call symputx("sigma", 1/sqrt(&N-3));
run;

%put &=sigma;


/* http://blogs.sas.com/content/iml/2016/03/09/comparative-panel-overlay-histograms-sas.html
*/
ods graphics / width=640px height=480px ANTIALIASMAX=20000;
title "Distribution of Sample Correlation";
title2 "Random Bivariate Normal Samples (N=&N)";
footnote j=L "Reference lines at values of the population correlation";
proc sgpanel data=DistribCorr noautolegend;
  panelby rho / onepanel layout=rowlattice rowheaderpos=right;
  histogram Corr;
  density Corr / type=kernel;
  refline rho / axis=x lineattrs=GraphData2(thickness=2);
  colaxis grid values=(-0.5 to 1 by 0.1);
  rowaxis offsetmin=0 values=(0 to 20 by 10) valueshint;
run;


title "Distribution of Fisher Transformation of Sample Correlation";
title2 "Random Bivariate Normal Samples (N=&N)";
footnote j=L "Reference lines at transformed values of the population correlation";
proc sgpanel data=DistribCorr noautolegend;
  format z0 4.2;
  panelby rho / onepanel layout=rowlattice rowheaderpos=right;
  histogram z;
  density z / type=normal(sigma=&sigma);
  inset z0 / border opaque;
  refline z0 / axis=x lineattrs=GraphData2(thickness=2);
  colaxis grid values=(-0.5 to 2 by 0.5);
  rowaxis offsetmin=0 values=(0 to 20 by 5) valueshint;
run;


/*******************/
/* Use the FISHER option in PROC CORR to compute CIs and test hypotheses */
title; footnote;
proc corr data=sashelp.iris fisher(rho0=0.75 biasadj=no) nosimple;
   where Species='Versicolor';
   var PetalLength PetalWidth;
run;