/*  A SAS program to perform a bootstrap analysis of
    median survival times for the data from problem
    4 on assignment 9, Spring 2002.
*/


data set1;
  infile 'c:\courses\st511\snew\heart1.dat';
  input y;
  run;

/* compute the median and other summary statistics
   and check for univariate normaliy */

proc univariate data=set1 normal;
  var y;
  run;

data set2;
  infile 'c:\courses\st511\snew\heart2.dat';
  input y;
  run;

/* compute the median and other summary statistics
   and check for univariate normaliy */

proc univariate data=set2 normal;
  var y;
  run;


/*  Create an IML program to construct a percentile
    bootstrap confidence interval for the sample
    median */

proc iml;
  start bootm;

  alpha=.05;         * Set the confidence level;
  per= (1-alpha)*100;

  use set1;          * Enter the data from set1;
  read all into x1;

  n1 = nrow(x1);     * Compute the sample size;
 
  m1=median(x1);     * Compute the sample median;
  
  use set2;          * Enter the data from set2;
  read all into x2;

  n2 = nrow(x2);     * Compute the sample size;
 
  m2=median(x2);     * Compute the sample median;
  


  print,,,,," median survival time: heart transplants = " m1;
  print,,,,," median survival time: no transplant = " m2;

  
/* Compute bootstrap confidence intervals 
   First work with the transplant cases */

  nboot = 5000;     * set number of bootstrap samples;

  m1boot = J(nboot,1);
  do i1 = 1 to nboot;
    xm1 = J(n1,1,0);
      do i2 = 1 to n1;
        krow = int(uniform(0)*n1)+1;
        xm1[i2,1] = x1[krow,1];
      end;
      m1boot[i1,1] = median(xm1);
  end;

/* Now work with non-transplant cases */

 m2boot = J(nboot,1);
  do i1 = 1 to nboot;
    xm2 = J(n2,1,0);
      do i2 = 1 to n2;
        krow = int(uniform(0)*n2)+1;
        xm2[i2,1] = x2[krow,1];
      end;
      m2boot[i1,1] = median(xm2);
  end;

/* Compute differences in bootstrapped medians
   for the two samples  */

   dboot = m1boot-m2boot;

/* Sort the bootstrapped results  */

  b1sort=m1boot;       
  b1rank=rank(m1boot);
  m1boot[b1rank, ] = b1sort;

  b2sort=m2boot;       
  b2rank=rank(m2boot);
  m2boot[b2rank, ] = b2sort;

  bdsort=dboot;       
  bdrank=rank(dboot);
  dboot[bdrank, ] = bdsort;

  /*  Compute a naive bootstrap confidence interval
      for the median of the transplant population */

  nl = int(nboot*alpha/2) +1;
  nu = int(nboot*(1-alpha/2)) +1;
  r1lower = m1boot[nl,1];
  r1upper = m1boot[nu,1];
  r2lower = m2boot[nl,1];
  r2upper = m2boot[nu,1];
  rdlower = dboot[nl,1];
  rdupper = dboot[nu,1];
  print,,,,,"Naive bootstrap: " per "percent confidence intervals: ";
  print "Median survival time with transplant:    " r1lower r1upper;
  print "Median survival time without transplant: " r2lower r2upper;
  print "Difference in median survival times:     " rdlower rdupper;


  /*  Compute bias corrected bootstrap confidence intervals */

  do i=1 to nboot while(m1boot[i,1]<m1);
     cm1 = i-1;
     end;

  zo = probit(cm1/nboot);
  za = probit(1-(alpha/2));
  n1l = int(nboot*probnorm(2*zo-za)) +1;
  n1u = int(nboot*probnorm(2*zo+za)) +1;
  r1lower = m1boot[n1l,1];
  r1upper = m1boot[n1u,1];

  do i=1 to nboot while(m2boot[i,1]<m2);
     cm2 = i-1;
     end;

  zo = probit(cm2/nboot);
  za = probit(1-(alpha/2));
  n2l = int(nboot*probnorm(2*zo-za)) +1;
  n2u = int(nboot*probnorm(2*zo+za)) +1;
  r2lower = m2boot[n2l,1];
  r2upper = m2boot[n2u,1];

  print,,,,,"Bias corrected bootstrap: " per "percent confidence intervals: ";
  print "Median survival time with transplant:    " r1lower r1upper;
  print "Median survival time without transplant: " r2lower r2upper;

/*  Output the bootstrap median values and differences
    to a file */

  rboot = m1boot||m2boot||dboot;
  create bootout from rboot [colname = {"m1" "m2" "diff"}];
  append from rboot;

  finish;

  run bootm;


  proc gchart data=bootout;
     vbar m1 / type=percent space=0 frame;
     title "Bootstraped median values";
	 title2 "Heart transplant survival times";
     run;

  proc gchart data=bootout;
     vbar m2 / type=percent space=0 frame;
     title "Bootstraped median values";
	 title2 "Non-transplant survival times";
     run;

  proc gchart data=bootout;
     vbar diff / type=percent space=0 frame;
     title "Bootstraped values";
	 title2 "Difference in median survival times";
     run;