options ls=78 nodate;
data hwk4data;
input patient time1 time2 time3 time4 time5;
cards;
1 11 18 15 18 15
2 33 27 31 21 17
3 20 28 27 23 29
4 28 26 18 18 18
5 22 23 22 16 10
6 20 22 16 15 12
7 24 27 22 21 24
8 30 26 30 24 20
;
run;
proc corr data=hwk4data cov noprob outp=corrout;
var time1-time5;
run;
proc iml;
/* First, lets use Hotellings T^2 to test equal means across time */
use corrout;
read all var {time1 time2 time3 time4 time5} where(_TYPE_='COV') into S;
read all var {time1 time2 time3 time4 time5} where(_TYPE_='MEAN') into xbartran;
xbar=xbartran`;
C={1 -1 0 0 0, 0 1 -1 0 0, 0 0 1 -1 0, 0 0 0 1 -1};
n=8;
q=5;
alpha=.05; /* For test of equal means across time */
T2=n*(C*xbar)`*inv(C*S*C`)*C*xbar;
Fcrit=(q-1)*(n-1)*finv(1-alpha,q-1,n-q+1)/(n-q+1);
scaledT2=T2*(n-q+1)/((q-1)*(n-1));
pval=1-probf(scaledT2,q-1,n-q+1);
print T2 Fcrit pval;
/* Now compute test of H-F conditions */
QQ=orpol(1:q,q-1); /* creates a matrix with q columns, that are orthogonal*/
/* polynomials of degree 0, 1, ...., q-1 */
PP=QQ(|1:q,2:q|); /* I only need columns 2 through q, so let PP equal the */
/* matrix consisting of all rows and columns 2-q of QQ */
print QQ PP;
Sigstar=PP`*S*PP; /* This is an estimate of PP`*Sigma*PP where Sigma is the */
/* pop var-cov matrix. Sigma satisfies H-F conditions if */
/* and only if PP`*Sigma*PP satisfies a condition known */
/* as sphericity. Note that Sigstar is (q-1) x (q-1) */
p=q-1; /* Lets call the dimension of Sigstar p */
Lambda= det(sigstar)/(( (1/p)*trace(sigstar) )**p);
/* Lambda is the likelihood ratio test statistic */
/* The p-value for this test statistic can be approximated by a formula */
/* given in Srivistava and Carter, 1983 p.327, which is implemented below */
m=n-1-(2*p*p+p+2)/(6*p);
a=(p+1)*(p-1)*(p+2)*(2*p*p*p+6*p*p+3*p+2)/(288*p*p);
f=p*(p+1)/2-1;
z=-m*log(Lambda);
pval=(1-probchi(z,f))+(a/(m*m))*(probchi(z,f)-probchi(z,f+4));
print, "LR test statistic =" lambda, "p-value for H-F conditions = " pval;
quit;
/* The same LR test of the H-F conditions can be obtained with SAS's PROC GLM */
title 'Hotellings T2 and test of H-F conditions using PROC GLM';
proc glm data=hwk4data;
model time1-time5 = / nouni;/* Remember we used PROC GLM to compute */
/* Hotellings T2 for the kite data (see */
/* kite2.sas). PROC GLM doesnt compute */
/* T2, but instead computes Wilks Lambda*/
/* which is equivalent and gives the */
/* same p-value. Notice that the p-value*/
/* for Wilks Lambda given by PROC GLM */
/* for the hypothesis of no TIME effect */
/* is the same as obtained above using */
/* PROC IML */
repeated time 5 / printe;/* the printe requests the test of H-F conditions */
/* which is called Mauchlys criterion applied to */
/* orthogonal components in the output. */
/* The repeated statement tells SAS that the five */
/* response variables time1-time5 are repeated */
/* measures across 5 time points */
run;
/* Now do randomized complete block analysis of these data as though H-F */
/* conditions were satisfied. For this we need the data in a different form */
proc transpose data=hwk4data out=hwk4d2(rename=(_name_=time col1=response));
by patient;
run;
title 'RCB Analysis using PROC GLM';
proc glm data=hwk4d2;
class patient time;
model response=time patient;
random patient;
run;