options ls=78 nodate;

data sparrow;
  infile 'sparrow.dat';
  input birdno x1 x2 x3 x4 x5;
  if (birdno < 22) then died=0;
  else died=1;
run;

proc sort data=sparrow;
  by died;
run;

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

proc iml;
  use corrout;
  read all var {x1 x2 x3 x4 x5} where(_type_='COV' & died=0) into S1;
  read all var {x1 x2 x3 x4 x5} where(_type_='MEAN' & died=0) into xbart1;
  read all var {x1 x2 x3 x4 x5} where(_type_='COV' & died=1) into S2;
  read all var {x1 x2 x3 x4 x5} where(_type_='MEAN' & died=1) into xbart2;
  xbar1=xbart1`;
  xbar2=xbart2`;
  n1=21;
  n2=28;
  p=5;
  alpha=.05;
  Spooled=((n1-1)*S1+(n2-1)*S2)/(n1+n2-2);
  T2=(xbar1-xbar2)`*inv(((1/n1)+(1/n2))*Spooled)*(xbar1-xbar2);
  Fcrit=finv(1-alpha,p,n1+n2-p-1)*p*(n1+n2-2)/(n1+n2-p-1);
  scaledT2=(n1+n2-p-1)*T2/p/(n1+n2-2);
  pval=1-probf(scaledT2,p,n1+n2-p-1);
  print T2 Fcrit scaledT2 pval;

  lowlim=j(p,1,0); /* define lowlim, uplim initially as 5 x 1 matrices of 0's*/
  uplim=j(p,1,0); 

  do i=1 to p;
    lowlim[i,1]=(xbar1[i,1]-xbar2[i,1])-sqrt(Fcrit)*sqrt( ((1/n1)+(1/n2))*Spooled[i,i] );
    uplim[i,1] =(xbar1[i,1]-xbar2[i,1])+sqrt(Fcrit)*sqrt( ((1/n1)+(1/n2))*Spooled[i,i] );
  end;
  
  print lowlim uplim;
  
  ahat=inv(Spooled)*(xbar1-xbar2);
  print ahat;   

run;

proc glm data=sparrow;
  class died;
  model x1 x2 x3 x4 x5=died/nouni;
  manova h=died;
run;