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