/*********************************************************
* Chapter 6 Exercise 6.19 in Johnson and Wichern 5th ed    *
* Comparison of Multivariate Mean Vectors                       *
**********************************************************/
DM "output;clear;log;clear";
Options pagesize=50 linesize=80 PageNo=1 NoDate;
OPTIONS FORMCHAR="|----|+|---+=|-/\<>*";
Title1 "Example 5.5 (Data Table 5.2)";
ODS HTML body=    "p6_19-body.html"
         contents="p6_19-contents.html"
         frame=   "p6_19-frame.html"
         page=    "p6_19-page.html"
         headtext="<title>Comparison of Multivariate Mean Vectors</title>"
         anchor="p6_19";
Data p6_19;
 infile 'C:\Teaching\Math593\Data\T6-10.dat';
 individual+1;
 Input Y1-Y3 truck $;
 label  Y1='Fuel' Y2='Repair' Y3='Capital' truck='Truck Type';


title "Comparison of Multivariate Mean Vectors for Milk Transportation Data";
Proc IML;
    Use p6_19;
    Read ALL var{Y1 Y2 Y3} where(truck='gasoline') into X1; 
	Read ALL var{Y1 Y2 Y3} where(truck='diesel') into X2; 
	N1=NROW(X1); N2=NROW(X2); P=NCOL(X1);
    close p6_19;
    One1=SHAPE(1,N1,1); One2=SHAPE(1,N2,1);
    MeanVec1=(One1`*X1)/N1; MeanVec2=(One2`*X2)/N2;  
    M1=REPEAT(MeanVec1,N1,1); M2=REPEAT(MeanVec2,N2,1);
    S1=(X1-M1)`*(X1-M1)/(N1-1); S2=(X2-M2)`*(X2-M2)/(N2-1); 
	Spool=((N1-1)*S1+(N2-1)*S2)/(N1+N2-2);
	X1bar=MeanVec1`;X2bar=MeanVec2`;
	print X1bar '   ' X2bar, S1, S2, Spool;

	*Test under equal covariance matrices;
	T2=N1*N2/(N1+N2)*(X1bar-X2bar)`*inv(Spool)*(X1bar-X2bar);
    C=(((N1+N2-2)*P)/(N1+N2-P-1))*FINV(0.99, P, N1+N2-P-1);
	F=(N1+N2-P-1)/((N1+N2-2)*P)*T2;
	Pval=1-probF(F,P,(N1+N2-P-1));
	print 'Test under equal covariance matrices';
    print T2 C Pval;
	CV=inv(Spool)*(X1bar-X2bar);
	print 'Coeff vector for the linear combination most responsible for rejection';
	print CV;

	*Test under unequal covariance matrices;
	T2=(X1bar-X2bar)`*inv(S1/N1+S2/N2)*(X1bar-X2bar);
    ChiSq=CINV(0.99, P);
	Pval=1-probChi(T2,P);
	print 'Test under unequal covariance matrices';
    print T2 ChiSq Pval;

    Start SCI(a,Mean1,Mean2,S,level,P,N1,N2);  *begin module named SCI;
    C=(((N1+N2-2)*P)/(N1+N2-P-1))*FINV(level, P, N1+N2-P-1);
    Margin=sqrt(C)*sqrt((a`*S*a)*(1/N1+1/N2));      
	D=Mean1-Mean2;
	LB=a`*D-Margin;
    UB=a`*D+Margin;
    contr=a`;

    Print contr '   ' LB UB;
    Finish;

    reset noprint;          

	print 'Simultaneous Confidence Intervals';
    Mean1=X1bar; Mean2=X2bar; S=Spool;
    a={1 0 0}`;
    Run SCI(a,Mean1,Mean2,S,.99,P,N1,N2);; *execute module SCI for various coefficient vectors;
    a={0 1 0}`;
    Run SCI(a,Mean1,Mean2,S,.99,P,N1,N2);;
    a={0 0 1}`;
    Run SCI(a,Mean1,Mean2,S,.99,P,N1,N2);;
    
run;

title 'Comparing two mean vectors using MANOVA';
PROC GLM data=p6_19;
class truck;
model Y1 Y2 Y3 = truck / nouni;
contrast 'gasoline vs. diesel' truck 1 -1;
manova h=truck / printe printh;

run;