/*  A SAS program to perform a bootstrap analysis of
    the correlation between average LSAT scores and GPA's
    for students in a sample of 15 law schools.
    The data are posted on the course web page as

                lawschl.dat

    and the program code is posted as

                lawschl.sas
*/


data set1;
  infile 'lawschl.dat';
  input school lsat gpa;
  run;

/*  Compute the correlation coefficient  */

proc corr data=set1;
  var lsat gpa;
  run;

/* Check for univariate normaliy */

proc univariate data=set1 normal;
  var lsat gpa;
  run;


/*  Create an IML program to perform a bootstrap
    analysis   */


proc iml;
  start bootcorr;
  use set1;          * Enter the data from set1;
  read all into zz;

  n = nrow(zz);     * Compute the number of rows and;
  pp1 = ncol(zz);    * columns in the data set;

  x = zz[ ,2:pp1];   * Select only the relavent columns;
  p = ncol(x);       * from the data file;

  a = t(x)*x-t(x[+, ])*x[+, ]/n;   * Compute the correlation matrix;
  scale = inv(sqrt(diag(a)));
  rr = scale*a*scale;

  r = rr[1,2];       * Compute a t-test for zero correlation;
  df = n-2;
  tr = r*sqrt(df)/sqrt(1-r**2);
  pval = (1-probt(abs(tr),df))*2;

  print,,,,," correlation = " r;
  print     "      t-test = " tr;
  print     "          df = " df;
  print     "     p-value = " pval;

  /* Construct a normal based confidence interval */

  alpha = .10;                       * Set the level of confidence;
  z = 0.5*log((1+r)/(1-r));
  zl = z - probit(1-alpha/2)/sqrt(n-3);
  zu = z + probit(1-alpha/2)/sqrt(n-3);
  rl = (exp(2*zl)-1)/(exp(2*zl)+1);
  ru = (exp(2*zu)-1)/(exp(2*zu)+1);

  per = (1-alpha)*100;
  print, per "% confidence interval: " rl " to " ru;

* Compute bootstrap confidence intervals ;

  nboot = 5000;            * set the number of bootstrap samples;

  rboot = J(nboot,1);
  do i1 = 1 to nboot;
     xc = J(p,p,0);
     xm = J(1,p,0);
      do i2 = 1 to n;
        krow = int(uniform(0)*n)+1;
        xc = xc + t(x[krow, ])*x[krow, ];
        xm = xm + x[krow, ];
      end;
      cov = xc - t(xm)*xm/n;
      rboot[i1,1] = cov[1,2]/sqrt(cov[1,1]*cov[2,2]);
  end;

  bsort=rboot;                     * Sort the correlations stored in rboot;
  brank=rank(rboot);
  rboot[brank, ] = bsort;

  /*  Compute a naive bootstrap confidence interval for the correlation */

  nl = int(nboot*alpha/2) +1;
  nu = int(nboot*(1-alpha/2)) +1;
  rlower = rboot[nl,1];
  rupper = rboot[nu,1];
  print,,,,,"Naive bootstrap " per "percent confidence interval: " rlower rupper;

  /*  Compute a bias corrected bootstrap confidence interval */

  do i=1 to nboot while(rboot[i,1]<r);
     cr = i-1;
     end;

  zo = probit(cr/nboot);
  za = probit(1-(alpha/2));
  nl = int(nboot*probnorm(2*zo-za)) +1;
  nu = int(nboot*probnorm(2*zo+za)) +1;
  rlower = rboot[nl,1];
  rupper = rboot[nu,1];

  print,,,,,"Bias corrected bootstrap " per "percent confidence interval: " rlower rupper;

/*  Output the bootstrap correlation values to a file */

  create bootout from rboot [colname = "bcorr"];
  append from rboot;

  finish;

  run bootcorr;


  proc gchart data=bootout;
     vbar bcorr / type=percent space=0 frame;
     title "Bootstraped correlation values";
     run;