Two Independent Samples T-square Analysis
Milk Transportation Data
Options PS=55 LS=90 PageNo=1 NoDate
FORMCHAR='|----|+|---+=|-/\<>*';
title1 'Table 6.6 -- Milk Transportation-Cost Data';
/* Page 263, Johnson and Wichern, 3rd Ed. */
data milk;
input truck $ x1 x2 x3 @@;
list;
truck=upcase(truck);
label truck='Truck Type' x1='Fuel' x2='Repair' x3='Capital';
cards;
G 16.44 12.43 11.23 G 7.19 2.70 3.92 G 9.92 1.35 9.75
G 4.24 5.78 7.78 G 11.20 5.05 10.67 G 14.25 5.78 9.88
G 13.50 10.98 10.60 G 13.32 14.27 9.45 G 29.11 15.09 3.28
G 12.68 7.61 10.23 G 7.51 5.80 8.13 G 9.90 3.63 9.13
G 10.25 5.07 10.17 G 11.11 6.15 7.61 G 12.17 14.26 14.39
G 10.24 2.59 6.09 G 10.18 6.05 12.14 G 8.88 2.70 12.23
G 12.34 7.73 11.68 G 8.51 14.02 12.01 G 26.16 17.44 16.89
G 12.95 8.24 7.18 G 16.93 13.37 17.59 G 14.70 10.78 14.58
G 10.32 5.16 17.00 G 8.98 4.49 4.26 G 9.70 11.59 6.83
G 12.72 8.63 5.59 G 9.49 2.16 6.23 G 8.22 7.95 6.72
G 13.70 11.22 4.91 G 8.21 9.85 8.17 G 15.86 11.42 13.06
G 9.18 9.18 9.49 G 12.49 4.67 11.94 G 17.32 6.86 4.44
D 8.50 12.26 9.11 D 7.42 5.13 17.15 D 10.28 3.32 11.23
D 10.16 14.72 5.99 D 12.79 4.17 29.28 D 9.60 12.72 11.00
D 6.47 8.89 19.00 D 11.35 9.95 14.53 D 9.15 2.94 13.68
D 9.70 5.06 20.84 D 9.77 17.86 35.18 D 11.61 11.75 17.00
D 9.09 13.25 20.66 D 8.53 10.14 17.45 D 8.29 6.22 16.38
D 15.90 12.90 19.09 D 11.94 5.69 14.77 D 9.54 16.77 22.66
D 10.43 17.65 10.66 D 10.87 21.52 28.47 D 7.13 13.22 19.44
D 11.88 12.18 21.20 D 12.03 9.22 23.09
;
proc print data=milk label;
run;
proc plot data=milk;
plot x1*x2=truck x1*x3=truck x2*x3=truck;
run;
proc iml;
reset nolog;
use milk;
read all var{x1 x2 x3} where(truck='G') into X1;
read all var{x1 x2 x3} where(truck='D') into X2;
close milk;
n1=nrow(X1); n2=nrow(X2); p=ncol(X1);
one1=J(n1,1,1); one2=J(n2,1,1);
xbar1=X1`*one1/n1;
xbar2=X2`*one2/n2;
xstar=X1-one1*xbar1`;
s1=xstar`*xstar/(n1-1);
xstar=x2-one2*xbar2`;
s2=xstar`*xstar/(n2-1);
spool=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2);
delta=xbar1-xbar2;
*-----------------------------------------*
| Treat as Homogeneous Variances |
*-----------------------------------------*;
T2=n1*n2/(n1+n2)*delta`*inv(Spool)*delta;
F=(n1+n2-p-1)/((n1+n2-2)*p)*T2;
pvalue=1-probf(f,p,(n1+n2-p-1));
print n1 n2 p,;
print xbar1 xbar2,;
print s1,s2,Spool,;
print 'Assume Homogeneous Variances: ',;
print delta ' ' T2 ' ' F ' ' pvalue,;
/* Compute Vector quantifying largest population difference */
l=inv(Spool)*delta; l=l/sqrt(l`*l);
print ,,'Vector Proportional to Largest Population Difference',;
print l,;
/* Compute 99% Confidence Region Ellipsoid Semiaxes */
alpha=0.01; /* set the significance level */
F=Finv((1-alpha),p,(n1+n2-p-1));
c=sqrt((n1+n2-2)*p/(n1+n2-p-1)*F);
call eigen(lambda,evecs,Spool);
semiaxis=J(p,p,0);
do j=1 to p;
semiaxis[,j]=c*sqrt(lambda[j]*(1/n2+1/n2))*evecs[,j];
end;
print 'Semiaxes for 99% Confidence Ellipse';
print c ' ' lambda, evecs, semiaxis,;
/* Compute Simultaneous T2 Intervals */
var=vecdiag(Spool)*(1/n1+1/n2); /* include multiplier of N's here */
t2inter=J(p,2,0);
t2inter[,1]=delta-c*sqrt(var);
t2inter[,2]=delta+c*sqrt(var);
print,,'T2 Simultaneous Confidence Intervals',;
print f alpha c var,t2inter,;
/* Bonferroni Simultaneous Intervals for each mean */
c=Tinv((1-alpha/(2*p)),(n1+n2-2));
boninter=J(p,2,0);
boninter[,1]=delta-sqrt(var)*c; /* lower limit */
boninter[,2]=delta+sqrt(var)*c; /* upper limit */
print 'Bonferroni Intervals:',delta var alpha c,boninter,;
/* Compute Usual Univariate Intervals for each mean */
c=Tinv((1-alpha/2),(n1+n2-2));
stdinter=J(p,2,0);
stdinter[,1]=delta-sqrt(var)*c; /* lower limit */
stdinter[,2]=delta+sqrt(var)*c; /* upper limit */
print 'Univariate Intervals:',delta var alpha c,stdinter,;
*-----------------------------------------*
| Examine Equality of Variances Assumption|
*-----------------------------------------*;
/* Compute Bartlett test for E1=E2 */
df1=n1-1; df2=n2-1; df12=df1+df2;
V1=df1*S1; V2=df2*S2; V12=V1+V2;
num=(df1/2)*log(det(V1))+(df2/2)*log(det(V2))+(p*df12/2)*log(df12);
denom=(df12/2)*log(det(V12))+(p*df1/2)*log(df1)+(p*df2/2)*log(df2);
lambda=exp(num-denom);
a=(df12**2-df1*df2)*(2*p**2+3*p-1)/(12*(p+1)*df1*df2);
m=2/df12*(df12-2*a);
X2=-m*log(lambda);
df=p*(p+1)/2;
pvalue=1-probchi(X2,df);
print 'Test for Equality of Covariance Matrices: ',;
print lambda a m X2 df pvalue,;
free lambda a m x2 df pvalue num denom df1 df2 df12 v1 v2 v12 test;
*-----------------------------------------*
| Treat as Different Variances |
*-----------------------------------------*;
/* Compute the Two Sample, E1 ne E2 test */
print ,'Two Sample Test with Unequal Variances',;
S=S1/n1+S2/n2;
T2=delta`*inv(S)*delta;
n=min(n1,n2); /* Conservative d.f. */
F=(n-p)/((n-1)*p)*T2;
Fpvalue=1-probF(f,p,(n-p));
/* compute effective d.f. using generalized Welch procedure */
df= 1/(n1**3-n1**2)*(delta`*inv(S)*S1*inv(S)*delta)**2*1/T2**2
+1/(n2**3-n2**2)*(delta`*inv(S)*S2*inv(S)*delta)**2*1/T2**2;
df=1/df;
Fadj=(df-p+1)/(df*p)*T2;
AFpvalue=1-probF(Fadj,p,(df-p+1));
Xpvalue=1-probchi(T2,p);
print n1 n2 p xbar1 xbar2;
print s1,s2,S,;
print 'Assume Unequal Variances and Use Conservative F test: ',;
print ' ' t2 ' ' F ' ' Fpvalue,;
print 'Using Adjusted Degrees of Freedom Method',;
print t2 ' ' df ' ' Fadj ' ' AFpvalue,;
print 'Large Sample Chi-square Results: ',t2 ' ' Xpvalue,;
reset log;
quit;
title2 'MANOVA Approach';
proc glm data=milk;
class truck;
model x1 x2 x3 = truck / nouni;
output out=errors residual=e1 e2 e3;
manova h=truck / printe printh;
run; quit;
Table 6.6 -- Milk Transportation-Cost Data 1
Truck
OBS Type Fuel Repair Capital
1 G 16.44 12.43 11.23
2 G 7.19 2.70 3.92
3 G 9.92 1.35 9.75
4 G 4.24 5.78 7.78
5 G 11.20 5.05 10.67
6 G 14.25 5.78 9.88
7 G 13.50 10.98 10.60
8 G 13.32 14.27 9.45
9 G 29.11 15.09 3.28
10 G 12.68 7.61 10.23
11 G 7.51 5.80 8.13
12 G 9.90 3.63 9.13
13 G 10.25 5.07 10.17
14 G 11.11 6.15 7.61
15 G 12.17 14.26 14.39
16 G 10.24 2.59 6.09
17 G 10.18 6.05 12.14
18 G 8.88 2.70 12.23
19 G 12.34 7.73 11.68
20 G 8.51 14.02 12.01
21 G 26.16 17.44 16.89
22 G 12.95 8.24 7.18
23 G 16.93 13.37 17.59
24 G 14.70 10.78 14.58
25 G 10.32 5.16 17.00
26 G 8.98 4.49 4.26
27 G 9.70 11.59 6.83
28 G 12.72 8.63 5.59
29 G 9.49 2.16 6.23
30 G 8.22 7.95 6.72
31 G 13.70 11.22 4.91
32 G 8.21 9.85 8.17
33 G 15.86 11.42 13.06
34 G 9.18 9.18 9.49
35 G 12.49 4.67 11.94
36 G 17.32 6.86 4.44
37 D 8.50 12.26 9.11
38 D 7.42 5.13 17.15
39 D 10.28 3.32 11.23
40 D 10.16 14.72 5.99
41 D 12.79 4.17 29.28
42 D 9.60 12.72 11.00
43 D 6.47 8.89 19.00
44 D 11.35 9.95 14.53
45 D 9.15 2.94 13.68
46 D 9.70 5.06 20.84
47 D 9.77 17.86 35.18
48 D 11.61 11.75 17.00
49 D 9.09 13.25 20.66
50 D 8.53 10.14 17.45
Table 6.6 -- Milk Transportation-Cost Data 2
Truck
OBS Type Fuel Repair Capital
51 D 8.29 6.22 16.38
52 D 15.90 12.90 19.09
53 D 11.94 5.69 14.77
54 D 9.54 16.77 22.66
55 D 10.43 17.65 10.66
56 D 10.87 21.52 28.47
57 D 7.13 13.22 19.44
58 D 11.88 12.18 21.20
59 D 12.03 9.22 23.09
Table 6.6 -- Milk Transportation-Cost Data 3
Plot of X1*X2. Symbol is value of TRUCK.
|
30 +
| G
|
|
|
| G
|
25 +
|
|
|
|
|
|
20 +
|
|
|
| G G
F | G
u | G D
e 15 + G
l | G
| G G
| D G GG
| G D G D D G
| G G D D
| D D
10 + G G DG GG G G D D
| G D G G D D D
| G D G D D G
| G G
| G D D
| D
|
5 +
| G
|
|
|
|
|
0 +
|
--+--------+--------+--------+--------+--------+--------+--------+--------+--------+-
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5
Repair
NOTE: 2 obs hidden.
Table 6.6 -- Milk Transportation-Cost Data 4
Plot of X1*X3. Symbol is value of TRUCK.
|
30 +
| G
|
|
|
| G
|
25 +
|
|
|
|
|
|
20 +
|
|
|
| G G
F | G
u | G D
e 15 + G
l | G
| G G G
| G G G D
| G G D D
| G G D D
| D D
10 + G G GGG D G G D D
| G G G D D D D
| G D GG D D
| G
| G D D
| D
|
5 +
| G
|
|
|
|
|
0 +
|
---+----------+----------+----------+----------+----------+----------+----------+--
0 5 10 15 20 25 30 35
Capital
NOTE: 4 obs hidden.
Table 6.6 -- Milk Transportation-Cost Data 5
Plot of X2*X3. Symbol is value of TRUCK.
|
|
22.5 +
|
| D
|
20.0 +
|
|
| D
17.5 + D G
| D
|
|
15.0 + G D
| G G
| G
R | G DD D
e 12.5 + D DG
p | G D D
a | G G G
i | G
r 10.0 + G D D
| G D
| G D
| GG
7.5 + G G
| G
| G G D
| GG G D
5.0 + GG GD D
| G G D
| G
| D D
2.5 + G G G
| G
| G
|
0.0 +
|
---+----------+----------+----------+----------+----------+----------+----------+--
0 5 10 15 20 25 30 35
Capital
Table 6.6 -- Milk Transportation-Cost Data 6
N1 N2 P
36 23 3
XBAR1 XBAR2
12.218611 10.105652
8.1125 10.762174
9.5902778 18.167826
S1
23.013361 12.366395 2.906609
12.366395 17.544111 4.7730821
2.906609 4.7730821 13.963334
S2
4.3623166 0.7598872 2.3620992
0.7598872 25.851236 7.6857322
2.3620992 7.6857322 46.6544
SPOOL
15.814712 7.8866902 2.6964473
7.8866902 20.75037 5.8972629
2.6964473 5.8972629 26.580938
Assume Homogeneous Variances:
DELTA T2 F PVALUE
2.1129589 50.912787 16.375458 1.0005E-7
-2.649674
-8.577548
Vector Proportional to Largest Population Difference
L
0.5931085
-0.31176
-0.742313
Semiaxes for 99% Confidence Ellipse
Table 6.6 -- Milk Transportation-Cost Data 7
C LAMBDA
3.5959644 32.950652
20.29867
9.8966986
EVECS
0.3836645 0.4926409 0.7810931
0.5910074 0.518933 -0.617591
0.7095857 -0.69858 0.0920582
SEMIAXIS
2.3353401 2.3535894 2.6056403
3.5974218 2.4792 -2.060216
4.3191998 -3.337461 0.3070959
T2 Simultaneous Confidence Intervals
F ALPHA C VAR
4.1590807 0.01 3.5959644 1.1268937
1.4785891
1.8940524
T2INTER
-1.704346 5.930264
-7.022268 1.7229199
-13.52648 -3.628618
Bonferroni Intervals:
DELTA VAR ALPHA C
2.1129589 1.1268937 0.01 3.0639489
-2.649674 1.4785891
-8.577548 1.8940524
BONINTER
-1.139584 5.3655016
-6.375352 1.0760037
-12.7943 -4.360802
Univariate Intervals:
DELTA VAR ALPHA C
Table 6.6 -- Milk Transportation-Cost Data 8
2.1129589 1.1268937 0.01 2.6648705
-2.649674 1.4785891
-8.577548 1.8940524
STDINTER
-0.715941 4.9418589
-5.890083 0.5907353
-12.24506 -4.910032
Test for Equality of Covariance Matrices:
LAMBDA A M X2 DF PVALUE
8.6119E-8 1.7438853 1.8776221 30.544284 6 0.000031
Two Sample Test with Unequal Variances
N1 N2 P XBAR1 XBAR2
36 23 3 12.218611 10.105652
8.1125 10.762174
9.5902778 18.167826
S1
23.013361 12.366395 2.906609
12.366395 17.544111 4.7730821
2.906609 4.7730821 13.963334
S2
4.3623166 0.7598872 2.3620992
0.7598872 25.851236 7.6857322
2.3620992 7.6857322 46.6544
S
0.828926 0.3765495 0.1834391
0.3765495 1.6113032 0.4667479
0.1834391 0.4667479 2.4163226
Assume Unequal Variances and Use Conservative F test:
T2 F FPVALUE
Table 6.6 -- Milk Transportation-Cost Data 9
43.176395 13.083756 0.0000591
Using Adjusted Degrees of Freedom Method
T2 DF FADJ AFPVALUE
43.176395 37.509307 13.624742 4.5067E-6
Large Sample Chi-square Results:
T2 XPVALUE
43.176395 2.2577E-9
Table 6.6 -- Milk Transportation-Cost Data 10 MANOVA Approach
General Linear Models Procedure
Class Level Information
Class Levels Values
TRUCK 2 D G
Number of observations in data set = 59
E = Error SS&CP Matrix
X1 X2 X3
X1 901.43859577 449.54134239 153.6974965
X2 449.54134239 1182.7710663 336.1439837
X3 153.6974965 336.1439837 1515.1134885
Table 6.6 -- Milk Transportation-Cost Data 11 MANOVA Approach
General Linear Models Procedure
Multivariate Analysis of Variance
Partial Correlation Coefficients from the Error SS&CP Matrix / Prob > |r|
DF = 57 X1 X2 X3
X1 1.000000 0.435363 0.131515
0.0001 0.0006 0.3251
X2 0.435363 1.000000 0.251103
0.0006 0.0001 0.0573
X3 0.131515 0.251103 1.000000
0.3251 0.0573 0.0001
Table 6.6 -- Milk Transportation-Cost Data 12 MANOVA Approach
General Linear Models Procedure
Multivariate Analysis of Variance
H = Type III SS&CP Matrix for TRUCK
X1 X2 X3
X1 62.655678803 -78.57091527 -254.3504762
X2 -78.57091527 98.528798102 318.95831461
X3 -254.3504762 318.95831461 1032.5347352
Characteristic Roots and Vectors of: E Inverse * H, where
H = Type III SS&CP Matrix for TRUCK E = Error SS&CP Matrix
Characteristic Percent Characteristic Vector V'EV=1
Root
X1 X2 X3
0.89320679 100.00 -0.01771513 0.00931173 0.02217160
0.00000000 0.00 0.02953750 0.00270494 0.00644057
0.00000000 0.00 -0.01354161 0.03163593 -0.01310838
Manova Test Criteria and Exact F Statistics for the Hypothesis of no Overall TRUCK Effect
H = Type III SS&CP Matrix for TRUCK E = Error SS&CP Matrix
S=1 M=0.5 N=26.5
Statistic Value F Num DF Den DF Pr > F
Wilks' Lambda 0.52820432 16.3755 3 55 0.0001
Pillai's Trace 0.47179568 16.3755 3 55 0.0001
Hotelling-Lawley Trace 0.89320679 16.3755 3 55 0.0001
Roy's Greatest Root 0.89320679 16.3755 3 55 0.0001