/* 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]