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;