DM 'LOG;CLEAR;OUT;CLEAR;';
/* options ls=79 pagesize=1000 nodate sasautos = '~/gauss/sas_macros/'; */
options ls=79 pagesize=1000 nodate sasautos = 'C:/rjb/stat_macros/';
data new;
/*  infile 'glucose.dat'; */
  infile 'C:/rjb/stat.537/glucose.dat';
  input x1-x3 y1-y3;
proc glm data = new noprint;
    model y1-y3=x1-x3;
    output out=res_out rstudent=resy1-resy3;
run;



 %transform(data=new,response=y1-y3, 
    explanatory = x1-x3, time1=1 2 3, levels = 3, Polydeg=2);
run;


proc univariate data = res_out normal;
    var resy1-resy3;
    title 'Smoothed Histogram of RStudent';
    histogram resy1-resy3 /kernel(l=1 color = black);
    inset mean std skewness kurtosis;
run;
proc univariate data = res_out normal noprint;
    var resy1-resy3;
    qqplot resy1-resy3 /normal(mu=est sigma=est) cframe=ligr
                 pctlaxis(grid lgrid=35 label='Normal Percentiles');
    inset mean std / cfill=white format=3.0 header='Normal Parameters'
                    position=(95,10) refpoint=br;
    title1 'Normal Quantile-Quantile Plot';
run;
proc iml;
 small=1.e-10;

/*-----------------------------------------*/
/* Module to compute the rank of a matrix  */
/*-----------------------------------------*/
start rankM(A) global(small);
  m = nrow(A); n = ncol(A);
  call svd(U,S,V,A);
  if m = 1
    then S = S[1];
    else if m = 0 then S = {0}; 
  tol = max(S) * small;
  r = max(loc(S>tol)); /*rank of A = number of nonzero singular values */
  return(r);
finish rankM;

 use new;
 read all var {y1 y2 y3} into Y;
 read all var {x1 x2 x3} into X1;
 reset noprint;
 n=nrow(Y); d=ncol(Y);
 X=J(n,1) || X1;
 r=rankM(X);
 H=X*ginv(X`*X)*X`;
 E=Y-H*Y;
 S=(E`*E)/(n-r);
 A=E*inv(S)*E`;
 DD=vecdiag(A);
 h1=vecdiag(I(n)-H);
 u=(DD#H1##(-1))/(n-r);
/*
   Under multivariate normality, the entries of u are each;
   distributed as a beta random variable with parameters;
   d/2 and (n-r-d)/2, where r = rank(X);
   For this data set d/2=1.5 and (n-r-d)/2 = 21.5;
*/

 D2_max=max(DD);
 B_max=max(u);
 F_max=((n-r-d)/d)*B_max/(1-B_max);
 prob_F =n*(1- probf(F_max,d,n-r-d));
 print D2_max B_max F_max prob_F;
 create tdata from u [colname = 'u'];
 append from u;
 close tdata;


data total;
 merge res_out tdata;

proc univariate data = total;
 var u;
*;
*  If d = 3, n = 50, and r=4 then alpha = 1.5 and beta = 21.5;
*  For other values of (n,d,r), it is necessary to change the values;
*  of alpha and beta to equal alpha = d/2 and beta =(n-d-r)/2;
*;
qqplot u/beta(alpha=1.5 beta=21.5 threshold=0 scale=1);
histogram u / beta(alpha=1.5 beta=21.5);
 inset mean (5.3) std = 'Std Dev' (5.3) Skewness (5.3)
   Kurtosis (5.3) /header = 'Summary Statistics' pos = ne;
   title1 'Plot of Mahalanobis Distances';
run;


proc model;
    parms b01 b11 b21 b31 b02 b12 b22 b32 b03 b13 b23 b33;
    instrument x1 x2 x3;
    Y1 = b01 +b11*x1 + b21*x2 + b31*x3;
    Y2 = b02 +b12*x1 + b22*x2 + b32*x3;
    Y3 = b03 +b13*x1 + b23*x2 + b33*x3;
    fit Y1-Y3 /normal;
run;