options ls=78;
title "Profile Analysis - Women's Nutrition Data";

 /* After reading in the data, each of the original
  * variables is divided by its null value to convert 
  * to a common scale without specific units.
  * The differences are then defined for each successive
  * pair of variables.
  */

data nutrient;
  infile "D:\Statistics\STAT 505\data\nutrient.csv" firstobs=2 delimiter=',';
  input id calcium iron protein a c;
  calcium=calcium/1000;
  iron=iron/15;
  protein=protein/60;
  a=a/800;
  c=c/75;
  diff1=iron-calcium;
  diff2=protein-iron;
  diff3=a-protein;
  diff4=c-a;
  run;

 /* The iml code to compute the hotelling t2 statistic
  * is similar to that for the one-sample t2 statistic
  * except that by working with differences of variable,
  * the null values are all 0s, and the corresponding
  * null hypothesis is that all variable means are equal
  * to each other.
  */

proc iml;
  start hotel;
    mu0={0,0,0,0};
    one=j(nrow(x),1,1);
    ident=i(nrow(x));
    ybar=x`*one/nrow(x);
    s=x`*(ident-one*one`/nrow(x))*x/(nrow(x)-1.0);
    print mu0 ybar;
    print s;
    t2=nrow(x)*(ybar-mu0)`*inv(s)*(ybar-mu0);
    f=(nrow(x)-ncol(x))*t2/ncol(x)/(nrow(x)-1);
    df1=ncol(x);
    df2=nrow(x)-ncol(x);
    p=1-probf(f,df1,df2);
    print t2 f df1 df2 p;
  finish;
  use nutrient;
  read all var{diff1 diff2 diff3 diff4} into x;
  run hotel;