options ls=78 ps=66 nodate;

data sparrow;
  infile 'sparrow.dat';
  input birdno x1 x2 x3 x4 x5;
run;

proc univariate data=sparrow normal plot;
  var x1 x2 x3 x4 x5;
run;

proc corr data=sparrow cov outp=corrout;
  var x1 x2 x3 x4 x5;
run;

proc print data=corrout;
run;


proc iml;
  use corrout;
  read all var {x1 x2 x3 x4 x5} where(_TYPE_='COV') into S;

  print S;

  Sinv=inv(S);
  call eigen(eigvals,eigvecs,Sinv);
  print eigvals eigvecs;

  lambda1=eigvals[1,1];
  e1=eigvecs[,1];
  print lambda1 e1;

  check=Sinv*e1-lambda1*e1;
  print check; /* should equal zero */

run;

data sparrow2;
  set sparrow;
  e11=0.017652;
  e12=-0.039642;
  e13=-0.356953;
  e14=0.9265815;
  e15=-0.110219;
  ep1=0.5365005;
  ep2=0.8290154;
  ep3=0.0964961;
  ep4=0.0743522;
  ep5=0.1003044;
  
  lincomb1=x1*e11 + x2*e12 + x3*e13 + x4*e14 + x5*e15;
  lincomb2=x1*ep1 + x2*ep2 + x3*ep3 + x4*ep4 + x5*ep5;

  keep birdno x1 x2 x3 x4 x5 lincomb1 lincomb2;

run;

proc print data=sparrow2;
run;

title 'Q-Q plot for <e_1,x>, <e_p,x> using PROC UNIVARIATE';
proc univariate data=sparrow2 plot normal;
  var lincomb1 lincomb2;
run;

/* An alternative method for obtaining the Q-Q plots: */
/* Compute the expected normal quantiles and plot the observed */
/* quantiles versus the normal quantiles using PROC PLOT */
/* This procedure allows for calculation of r_q */

proc sort data=sparrow2 out=sortout;
  by lincomb1;
run;

data l1norm;
  set sortout;
  phiinv = probit((_N_-0.5)/49); /* sample size is 49 */
run;

title 'Q-Q plot for <e_1,x> using PROC PLOT';
proc plot data=l1norm;
  plot lincomb1*phiinv='*';
run;

title 'Correlation between observed and expected quantiles (r_q) for lincomb1';
proc corr data=l1norm noprob nosimple;
  var lincomb1 phiinv;
run;

proc sort data=sparrow2 out=sortout;
  by lincomb2;
run;

data l2norm;
  set sortout;
  phiinv = probit((_N_-0.5)/49); /* sample size is 49 */
run;

title 'Q-Q plot for <e_p,x> using PROC PLOT';
proc plot data=l2norm;
  plot lincomb2*phiinv='*';
run;

title 'Correlation between observed and expected quantiles (r_q) for lincomb2';
proc corr data=l2norm noprob nosimple;
  var lincomb2 phiinv;
run;