options ls=78 nodate; data salsinol; input patient day1 day2 day3 day4; cards; 1 .64 .70 1.00 1.40 2 .33 .70 2.33 3.20 3 .73 1.85 3.60 2.60 4 .70 4.20 7.30 5.40 5 .40 1.60 1.40 7.10 6 2.60 1.30 .70 .70 7 7.80 1.20 2.60 1.80 8 5.30 .90 1.80 .70 9 2.50 2.10 1.12 1.01 10 1.90 1.30 4.40 2.80 11 .98 .32 3.91 .66 12 .39 .69 .73 2.45 13 .31 6.34 .63 3.86 14 .50 .40 1.10 8.10 ; run; proc corr data=salsinol cov noprob outp=corrout; var day1 day2 day3 day4; run; proc iml; /* First, lets use Hotellings T^2 to test equal means across time */ use corrout; read all var {day1 day2 day3 day4} where(_TYPE_='COV') into S; read all var {day1 day2 day3 day4} where(_TYPE_='MEAN') into xbartran; xbar=xbartran`; C={1 -1 0 0, 0 1 -1 0, 0 0 1 -1}; n=14; q=4; 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=salsinol; model day1 day2 day3 day4 = / 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 4 / 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 four */ /* response variables day1 - day4 are repeated */ /* measures across 4 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 */ /* Below we reenter the data in the configuration that we need. Alternatively,*/ /* we could avoid re-entering the data by using PROC TRANSPOSE. The code to */ /* reconfigure the data using PROC TRANSPOSE is included below for your */ /* information, but is commented out. */ data sal2; input obs patient day $ salsinol; cards; 1 1 DAY1 0.64 2 1 DAY2 0.70 3 1 DAY3 1.00 4 1 DAY4 1.40 5 2 DAY1 0.33 6 2 DAY2 0.70 7 2 DAY3 2.33 8 2 DAY4 3.20 9 3 DAY1 0.73 10 3 DAY2 1.85 11 3 DAY3 3.60 12 3 DAY4 2.60 13 4 DAY1 0.70 14 4 DAY2 4.20 15 4 DAY3 7.30 16 4 DAY4 5.40 17 5 DAY1 0.40 18 5 DAY2 1.60 19 5 DAY3 1.40 20 5 DAY4 7.10 21 6 DAY1 2.60 22 6 DAY2 1.30 23 6 DAY3 0.70 24 6 DAY4 0.70 25 7 DAY1 7.80 26 7 DAY2 1.20 27 7 DAY3 2.60 28 7 DAY4 1.80 29 8 DAY1 5.30 30 8 DAY2 0.90 31 8 DAY3 1.80 32 8 DAY4 0.70 33 9 DAY1 2.50 34 9 DAY2 2.10 35 9 DAY3 1.12 36 9 DAY4 1.01 37 10 DAY1 1.90 38 10 DAY2 1.30 39 10 DAY3 4.40 40 10 DAY4 2.80 41 11 DAY1 0.98 42 11 DAY2 0.32 43 11 DAY3 3.91 44 11 DAY4 0.66 45 12 DAY1 0.39 46 12 DAY2 0.69 47 12 DAY3 0.73 48 12 DAY4 2.45 49 13 DAY1 0.31 50 13 DAY2 6.34 51 13 DAY3 0.63 52 13 DAY4 3.86 53 14 DAY1 0.50 54 14 DAY2 0.40 55 14 DAY3 1.10 56 14 DAY4 8.10 ; run; /* proc transpose data=salsinol out=sal2(rename=(_name_=day col1=salsinol)); by patient; run; proc print data=sal2; run; */ title 'RCB Analysis using PROC GLM'; proc glm data=sal2; class patient day; model salsinol=day patient; random patient; run;