/*-------------------------------------------------------------------*/ /* Univariate and Multivariate General Linear Models: */ /* Theory and Applications Using SAS(R) Software */ /* by Neil H. Timm and Tammy A. Mieczkowski */ /* Copyright(c) 1997 by SAS Institute Inc., Cary, NC, USA */ /* SAS Publications order # 55809 */ /* ISBN 1-55544-987-5 */ /*-------------------------------------------------------------------*/ /* */ /* This material is provided "as is" by SAS Institute Inc. There */ /* are no warranties, express or implied, as to merchantability or */ /* fitness for a particular purpose regarding the materials or code */ /* contained herein. The Institute is not responsible for errors */ /* in this material as it now exists or will exist, nor does the */ /* Institute provide technical support for it. */ /* */ /*-------------------------------------------------------------------*/ /* Questions or problem reports concerning this material may be */ /* addressed to the author: */ /* */ /* Neil H. Timm, Ph.D. */ /* Professor */ /* Department of Psychology in Education */ /* 5C01 Forbes Quadrangle */ /* Pittsburgh PA 15260 */ /* */ /* or by email: */ /* TIMM@vms.cis.pitt.edu */ /* */ /*-------------------------------------------------------------------*/ This file contains example code and data sets that are used in the book, Univariate and Multivariate General Linear Models: Theory and Applications Using SAS(R) Software by Neil H. Timm and Tammy A. Mieczkowski. ************************************************************************* ************************************************************************* EXAMPLE CODE USED IN THIS BOOK *************************************************** /* The following code appears on page 8, */ /* where it is used to produce output 1.6. */ *************************************************** /* Program 1_6.sas */ /* Program to create a multivariate normal data set */ options ls=80 ps=60 nodate nonumber; filename app1 '1_6.dat'; title1 'Output 1.6: Generating a Multivariate Normal Data Set'; proc iml; seed=30195; z=normal(repeat(seed,50,4)); u={10,20,30,40}; s={3 1 0 0, 1 4 0 0, 0 0 1 4, 0 0 4 20}; a=root(s); uu=repeat(u`,50,1); y=(z*a) + uu; print y; file app1; do i=1 to nrow(y); do j=1 to ncol(y); put (y[i,j]) 10.2 +2 @; end; put; end; closefile app1; quit; *************************************************** /* The following code appears on page 11, */ /* where it is used to produce output 1.7.1. */ *************************************************** /* Program 1_7_1.sas */ /* Program to create Q-Q plot of data */ options ls=80 ps=60 nodate nonumber; filename app1 '1_6.dat'; title1 'Output 1.7.1: Q-Q plot of 1.6 Data (y1)'; data ex171; infile app1; input y1-y4; proc sort; by y1; proc univariate noprint; var y1; output out=stats n=nobs mean=mean std=std; data quantile; set ex171; if _n_=1 then set stats; i+1; p=(i - .5) / nobs; z=probit(p); normal = mean + z*std; proc print; proc corr; var y1 z; run; filename out '1_7_1.cgm'; goptions device=cgmmwwc gsfname=out gsfmode=replace colors=(black) hsize=5in vsize=4in; proc gplot data=quantile; plot y1*z normal*z /overlay frame; symbol1 v=; symbol2 i=join v=none l=1; run; *************************************************** /* The following code appears on page 13, */ /* where it is used to produce output 1.7.2. */ *************************************************** /* Program 1_7_2.sas */ /* Program to create Q-Q plot of 1/(y1**2) Data */ options ls=80 ps=60 nodate nonumber; filename app1 '1_6.dat'; title1 'Output 1.7.2: Q-Q Plot of 1/(y1**2)'; data ex172; infile app1; input y1-y4; ty1=1/(y1**2); proc sort; by ty1; proc univariate noprint; var ty1; output out=stats n=nobs mean=mean std=std; data quantile; set ex172; if _n_=1 then set stats; i+1; p=(i - .5) / nobs; z=probit(p); normal = mean + z*std; proc print; proc corr; var ty1 z; run; filename out '1_7_2.cgm'; goptions device=cgmmwwc gsfname=out gsfmode=replace colors=(black) hsize=5in vsize=4in; proc gplot data=quantile; plot ty1*z normal*z /overlay frame; symbol1 v=; symbol2 i=join v=none l=1; run; *************************************************** /* The following code appears on page 15, */ /* where it is used to produce output 1.7.3. */ *************************************************** /* Program 1_7_3.sas */ /* Program to create Q-Q plot of y1 data with Outlier */ options ls=80 ps=60 nodate nonumber; filename app1 '1_6.dat'; title1 'Output 1.7.3: Q-Q plot of y1 data with an Outlier'; data ex173; infile app1; input y1-y4; if y1 ge 13.7 then y1 = 17; proc sort; by y1; proc univariate noprint; var y1; output out=stats n=nobs mean=mean std=std; data quantile; set ex173; if _n_=1 then set stats; i+1; p=(i - .5) / nobs; z=probit(p); normal = mean + z*std; proc print; proc corr; var y1 z; run; filename out '1_7_3.cgm'; goptions device=cgmmwwc gsfname=out gsfmode=replace colors=(black) hsize=5in vsize=4in; proc gplot data=quantile; plot y1*z normal*z /overlay frame; symbol1 v=; symbol2 i=join v=none l=1; run; *************************************************** /* The following code appears on page 17, */ /* where it is used to produce output 1.7.4. */ *************************************************** /* Program 1_7_4.sas */ /* Program to create Q-Q plot of a dataset */ /* To run this program on your own dataset change */ /* the name of the file in the file=_____statement */ /* and the number of columns in the p=___ statement */ options ls=80 ps=60 nodate nonumber; %let file = ycondx.dat; %let p = 3; /* macro to expand the string of variables that are processed */ %macro expand(cols); %do j=1 %to &cols; v&j %end; %mend expand; /* macro to perform Q-Q plotting of the variables */ %macro qq(cols); %do i=1 %to &cols; proc sort data=ex174; by v&i; proc univariate noprint data=ex174; var v&i; output out=stats n=nobs mean=mean std=std; data quantile; set ex174; if _n_=1 then set stats; k+1; pr=(k - .5) / nobs; z=probit(pr); normal = mean + z*std; proc print data=quantile; title "Output 1.7.4: Q-Q plot, variable V&i, &file"; proc corr data=quantile; var v&i z; filename out "1_7_4_&i'.cgm"; goptions device=cgmmwwc gsfname=out gsfmode=replace colors=(black) hsize=5in vsize=4in; proc gplot data=quantile; title "Output 1.7.4: Q-Q plot, variable V&i, &file"; plot v&i*z normal*z /overlay frame; symbol1 v=; symbol2 i=join v=none l=1; run; %end; %mend qq; data ex174; infile "&file"; input %expand(&p); proc print data=ex174; title "Output 1.7.4: &file"; %qq(&p); *************************************************** /* The following code appears on page 22, */ /* where it is used to produce output 1.8.1. */ *************************************************** /* Program 1_8_1.sas */ /* Program to create Chi-Square Plot of Y Data */ /* Data set from 1_6.sas is used, with a column */ /* of observation numbers added to the file */ options ls=80 ps=60 nodate nonumber; filename app1 '1_6.da2'; title1 'Output 1.8.1: Chi-Square Plot of the 1.6 Dataset'; data ex181; infile app1; input tag $ y1 - y4; label tag = 'id' y1 = 'var1' y2 = 'var2' y3 = 'var3' y4 = 'var4'; %let id=tag; %let var=y1 y2 y3 y4; proc iml; reset; start dsquare; use _last_; read all var {&var} into y [colname=vars rowname=&id]; n=nrow(y); p=ncol(y); r1=&id; print y; m=y[ :,]; d=y - j(n,1) * m; s=d` * d / (n-1); dsq=vecdiag(d * inv(s) * d`); r=rank(dsq); val=dsq; dsq [r, ] = val; val=r1; &id [r] = val; z=((1:n)` - .5) / n; chisq = 2 * gaminv(z,p/2); result = dsq||chisq; cl={'dsq' 'chisq'}; create dsquare from result [colname=cl rowname=&id]; append from result [rowname=&id]; finish; run dsquare; quit; proc print data=dsquare; var tag dsq chisq; run; filename out '1_8_1.cgm'; goptions device=cgmmwwc gsfname=out gsfmode=replace colors=(black) hsize=5in vsize=4in; proc gplot data=dsquare; plot dsq*chisq /frame; symbol1 v=; run; *************************************************** /* The following code appears on page 24, */ /* where it is used to produce output 1.8.2. */ *************************************************** /* Program 1_8_2.sas */ /* Program to create Chi-Square Plot */ /* of Residuals from Rhower data */ options ls=80 ps=60 nodate nonumber; filename rhower 'ycondx.da2'; title1 'Output 1.8.2: Chi-Square Plot of Residuals'; data ex182; infile rhower; input tag $ yc1-yc3; label tag = 'id' yc1 = 'var1' yc2 = 'var2' yc3 = 'var3'; %let id=tag; %let var=yc1 yc2 yc3; proc iml; reset; start dsquare; use _last_; read all var {&var} into y [colname=vars rowname=&id]; n=nrow(y); p=ncol(y); r1=&id; print y; m=y[ :,]; d=y - j(n,1) * m; s=d` * d / (n-1); dsq=vecdiag(d * inv(s) * d`); r=rank(dsq); val=dsq; dsq [r, ] = val; val=r1; &id [r] = val; z=((1:n)` - .5) / n; chisq = 2 * gaminv(z,p/2); result = dsq||chisq; cl={'dsq' 'chisq'}; create dsquare from result [colname=cl rowname=&id]; append from result [rowname=&id]; finish; run dsquare; quit; proc print data=dsquare; var tag dsq chisq; run; filename out '1_8_2.cgm'; goptions device=cgmmwwc gsfname=out gsfmode=replace colors=(black) hsize=5in vsize=4in; proc gplot data=dsquare; plot dsq*chisq /frame; symbol1 v=; run; *************************************************** /* The following code appears on page 26, */ /* where it is used to produce output 1.9.1. */ *************************************************** /* Program 1_9_1.sas */ /* Program to produce bivariate scatter plot of Rhower Data */ options ls=80 ps=60 nodate nonumber; filename rhower '5_1.dat'; title 'Output 1.9.1: Bivariate Scatter Plots of Rhower Data'; data ex191; infile rhower; input y1-y3 x0-x5; proc plot; plot y1*x1; plot y1*x2; plot y1*x3; plot y1*x4; plot y1*x5; run; *************************************************** /* The following code appears on page 29, */ /* where it is used to produce output 1.9.2. */ *************************************************** /* Program 1_9_2.sas */ /* Program to create 3-D Plots of bivariate normal distributions*/ options ls=80 ps=60 nodate nonumber; title 'Output 1.9.2: Bivariate Normal Distribution'; title2 'with u=(0, 0), var(y1)=3, var(y2)=4, cov(y1,y2)=1, r=.289'; data bivar; vy1=3; vy2=4; r=.289; keep y1 y2 z; cons=1/(2*3.14159265*sqrt(vy1*vy2*(1-r*r))); do y1=-10 to 10 by .2; do y2=-10 to 10 by .2; zy1=y1/sqrt(vy1); zy2=y2/sqrt(vy2); d=((zy1**2)+(zy2**2)-2*r*zy1*zy2)/(1-r**2); z=cons*exp(-d/2); if z > .001 then output; end; end; run; filename out1 '1_9_2_1.cgm'; goptions device=cgmmwwc gsfname=out1 gsfmode=replace colors=(black) hsize=6in vsize=5in; proc g3d data=bivar; plot y1*y2=z; run; /* A Second Plot*/ title2 'with u=(0, 0), var(y1)=3, var(y2)=20, cov=0, r=0'; data bivar2; vy1=3; vy2=20; r=0; keep y1 y2 z; cons=1/(2*3.14159265*sqrt(vy1*vy2*(1-r*r))); do y1=-10 to 10 by .2; do y2=-10 to 10 by .2; zy1=y1/sqrt(vy1); zy2=y2/sqrt(vy2); d=((zy1**2)+(zy2**2)-2*r*zy1*zy2)/(1-r**2); z=cons*exp(-d/2); if z > .001 then output; end; end; run; filename out2 '1_9_2_2.cgm'; goptions device=cgmmwwc gsfname=out2 gsfmode=replace colors=(black) hsize=6in vsize=5in; proc g3d data=bivar2; plot y1*y2=z; run; *************************************************** /* The following code appears on page 32, */ /* where it is used to produce output 1.10. */ *************************************************** /* Program 1_10.sas */ /* Program to calculate Mardia's measure of multivariate */ /* skewness and kurtosis */ options ls=80 ps=60 nodate nonumber; title 'Output 1.10: Mardias Multivariate Skewness & Kurtosis'; data Rhower; infile '5_1.dat'; input y1-y3 x0-x5; proc print data=Rhower; proc iml; use Rhower; v={y1 y2 y3}; w={x0 x1 x2 x3 x4 x5}; read all var v into y; read all var w into x; beta=inv(x`*x)*x`*y; n=nrow(y); p=ncol(y); k=ncol(x); s=(y`*y-y`*x*beta)/(n-k); s_inv=inv(s); q=i(n)-x*(inv(x`*x)*x`); d=q*y*s_inv*y`*q; b1=(sum(d#d))/(n*n); b2=trace(d#d)/n; kappa1= n*b1/6; kappa2=(b2-p*(p+2))/sqrt(8*p*(p+2)/n); dfchi=p*(p+1)*(p+2)/6; pvalskew=1-probchi(kappa1,dfchi); pvalkurt=2*(1-probnorm(abs(kappa2))); print s; print s_inv; print b1; print kappa1; print pvalskew; print b2; print kappa2; print pvalkurt; quit; *************************************************** /* The following code appears on page 34, */ /* where it is used to produce output 1.11. */ *************************************************** /* Program 1_11.sas */ /* Program to compute Box Cox Transformations */ /* To run this program on your own dataset change */ /* the name of the file in the file=____ statement */ /* the number of rows in the n=____ statement and */ /* the number of columns (variables) in the p=___ statement */ options ls=80 ps=60 nodate nonumber mprint; %let file=c:\exp1_6.dat; %let n=50; %let p=1; /*macro to expand the string of variables that are processed */ %macro expand(cols); %do j=1 %to &cols; x&j %end; %mend expand; /*macro to perform the Box-Cox transformation on the data matrix */ %macro loop(cols); %do i=1 %to &cols; proc iml; use matrix; read all var {x&i}; in=i(&n); allh={-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3}; one=j(&n,1,1); c=in-(one*(inv(one`*one))*one`); do k=1 to 23 by 1; h=allh[k,1]; xh=x&i##h; hinv=1/h; vhinv=j(&n,1,hinv); y=(xh-one)#vhinv; my=(one`*y)/&n; ycy=y`*c*y; lnx=log(x&i); slnx=one`*lnx; ycyn=ycy/&n; if ycyn > 0 then lhp1=-(&n/2)*log(ycyn); else lhp1=.; lhp2=(h-1)*slnx; if ycyn > 0 then lh=lhp1+lhp2; else lh=.; lhs=lhs//lh; end; Lambda=allh||lhs; print, "Lambda and corresponding likelihood for variables x&i", lambda; call pgraf(lambda, '*', 'lambda', 'likelihood', "plot of lambda vs likelihood for variable x&i"); quit; %end; %mend loop; /*input the data and process the macro */ data matrix; infile "&file"; input (%expand(&p)) (&p*:25.) ; title "Output 1.11: Box-Cox Transformation plots of &file"; proc print; %loop(&p) run; *************************************************** /* The following code appears on page 51, */ /* where it is used to produce output 2.6. */ *************************************************** /* Program 2_6.sas */ /* Program to perform Multiple Linear Regression Analysis of homicide data */ /* Data from exercise 4, p 117, Kleinbaum, Kupper, Muller */ options ls=80 ps=60 nodate nonumber; title 'Output 2.6: Multiple Linear Regression Analysis of Homicide data'; data reg; infile 'c:\2_6.dat'; input city y x1 x2 x3; label y='homicide rate' x1='population size' x2='percent low income' x3='unemployment rate'; proc print; proc univariate; var y; proc corr; var y x1 x2 x3; proc reg; model y = x1 x2 x3 /vif; proc reg; model y = x1 x2 x3 /partial selection = backward; model y = x2 x3 /vif r collin influence; paint student. ge 2 or student. le -2 /symbol = 'o'; plot r.*x2 r.*x3 /hplots=2; plot r.*p. student.*p.; run; *************************************************** /* The following code appears on page 59, */ /* where it is used to produce output 2.7. */ *************************************************** /* Program 2_7.sas */ /* Program to perform analysis of a One-Way Design with */ /* Unequal n, using Full Rank model */ options ls=80 ps=60 nodate nonumber; filename treat '2_7.dat'; title 'Output 2.7: One-Way Design with Unequal n, Full Rank Model'; data treat; infile treat; input treat y; proc print; proc glm; class treat; model y=treat /noint e xpx; contrast 'treatments' treat 1 0 0 -1, treat 0 1 0 -1, treat 0 0 1 -1; estimate 'treat1-treat4' treat 1 0 0 -1; estimate 'treat2-treat4' treat 0 1 0 -1; estimate 'treat3-treat4' treat 0 0 1 -1; run; proc glm; class treat; model y=treat /noint e xpx; contrast 'treats123 vs 4' treat -1 -1 -1 3; estimate 'treats123 vs 4' treat -1 -1 -1 3; run; *************************************************** /* The following code appears on page 63, */ /* where it is used to produce output 2.8. */ *************************************************** /* Program 2_8.sas */ /* Program to perform analysis of a nested design with */ /* Unequal n, using Full Rank model */ options ls=80 ps=60 nodate nonumber; filename nested 'c:\2_8.dat'; title 'Output 2.8: Nested Design with Unequal n, Full Rank Model'; data nest; infile nested; input a b y; label a='course' b='section'; proc print; proc glm; class a b; model y = a b(a) /noint e; contrast 'ha_1' a 1 -1; estimate 'ha_1' a 1 -1 /e; contrast 'hb_1' b(a) 1 -1 0 0 0; estimate 'hb_1vs2' b(a) 1 -1 0 0 0 /e; contrast 'hb_2' b(a) 0 0 1 0 -1, b(a) 0 0 0 1 -1; estimate 'hb_1vs3' b(a) 0 0 1 0 -1 /e; estimate 'hb_2vs3' b(a) 0 0 0 1 -1 /e; contrast 'hb_1and2' b(a) 1 -1 0 0 0, b(a) 0 0 1 0 -1, b(a) 0 0 0 1 -1; proc glm; class a b; model y = a b(a); run; *************************************************** /* The following code appears on page 69, */ /* where it is used to produce output 2.9. */ *************************************************** /* Program 2_9.sas */ /* Program to perform analysis of covariance with different slopes and */ /* Unequal n, using Full Rank model */ /* the one-way Unbalanced Intraclass covariance model */ options ls=80 ps=60 nodate nonumber; filename ex29 'c:\2_9.dat'; title 'Output 2.9: One-way Unbalanced Intraclass covariance model'; data cov; infile ex29; input y a z a1 a2 a3 z1 z2 z3; proc print; /* Using the full rank model */ proc reg; model y = a1 a2 a3 z1 z2 z3 /noint; parallel: mtest z1-z2=0, z1-z3=0 /print details; run; /* Using the less than full rank model */ proc glm; class a; model y=a z a*z /e xpx solution; lsmeans a /stderr pdiff cov; run; *************************************************** /* The following code appears on page 79, */ /* where it is used to produce output 3.3. */ *************************************************** /* Program 3_3.sas */ /* Two-Way Factorial Design, Unequal Cell Freq */ /* with and without interaction */ /* data from Timm and Carlson (1975, page 58) */ options ls=80 ps=60 nodate nonumber; title 'Output 3.3: 2-Way Factorial Design'; data two; infile 'c:\3_3.dat'; input a b y1 x11 x12 x13 x14 x21 x22 x23 x24; proc print data=two; run; proc reg; title2 'Unrestricted: 2-Way w/Interaction (univar)--Proc Reg'; model y1 = x11 x12 x13 x14 x21 x22 x23 x24 /noint; AB: mtest x11-x12-x21+x22=0, x12-x13-x22+x23=0, x13-x14-x23+x24=0 /print details; A: mtest x11+x12+x13+x14-x21-x22-x23-x24=0 /print details; B: mtest x11+x21-x14-x24=0, x12+x22-x14-x24=0, x13+x23-x14-x24=0 /print details; wtA: mtest .1875*x11+.3125*x12+.1875*x13+.3125*x14- .2500*x21-.2500*x22-.2500*x23-.2500*x24=0 /print details; wtB: mtest .4286*x11+.5714*x21-.5556*x14-.4444*x24=0, .5556*x12+.4444*x22-.5556*x14-.4444*x24=0, .4286*x13+.5714*x23-.5556*x14-.4444*x24=0 /print details; run; proc glm; title2 'Unrestricted: 2-Way w/Interaction (univar)--Proc GLM'; class a b; model y1 = a b a*b /ss1 ss2 ss3; proc glm; class a b; model y1 = b a a*b /ss1 ss2 ss3; run; proc reg; title2 'Restricted: 2-Way w/o Interaction (univar)--Proc Reg'; model y1 = x11 x12 x13 x14 x21 x22 x23 x24 /noint; restrict x11 - x12 - x21 + x22 = 0, x12 - x13 - x22 + x23 = 0, x13 - x14 - x23 + x24 = 0; A: mtest x11+x12+x13+x14-x21-x22-x23-x24=0 /print details; B: mtest x11+x21-x14-x24=0, x12+x22-x14-x24=0, x13+x23-x14-x24=0 /print details; wtA: mtest .1875*x11+.3125*x12+.1875*x13+.3125*x14- .2500*x21-.2500*x22-.2500*x23-.2500*x24=0 /print details; wtB: mtest .4286*x11+.5714*x21-.5556*x14-.4444*x24=0, .5556*x12+.4444*x22-.5556*x14-.4444*x24=0, .4286*x13+.5714*x23-.5556*x14-.4444*x24=0 /print details; run; proc glm; title2 'Restricted: 2-Way w/o Interaction (univar)--Proc GLM'; class a b; model y1 = a b /ss1; proc glm; class a b; model y1 = b a /ss1; run; *************************************************** /* The following code appears on page 94, */ /* where it is used to produce output 3.4. */ *************************************************** /* Program 3_4.sas */ /* Latin Square Design, Full Rank Model */ options ls=80 ps=60 nodate nonumber; title 'Output 3.4: Latin Square Design, Full Rank Model'; data Latin; infile 'c:\3_4.dat'; input a b c y1 x112 x123 x131 x213 x221 x232 x311 x322 x333; proc print data=latin; run; proc reg; model y1 = x112 x123 x131 x213 x221 x232 x311 x322 x333 /noint; restrict x112 - x131 - x213 + x221 - x322 + x333 = 0, -x112 + x123 - x221 + x232 + x311 - x333 = 0; A: mtest x112 + x123 + x131 - x311 - x322 - x333 = 0, x213 + x221 + x232 - x311 - x322 - x333 = 0 /print details; B: mtest x112 - x131 + x213 - x232 + x311 - x333 = 0, x123 - x131 + x221 - x232 + x322 - x333 = 0 /print details; C: mtest -x123 + x131 - x213 + x221 + x311 - x333 = 0, x112 - x123 - x213 + x232 + x322 - x333 = 0 /print details; run; *************************************************** /* The following code appears on page 102 */ /* where it is used to produce output 3.5.1. */ *************************************************** /* Program 3_5_1.sas */ /* Split-Plot (Two-Group Profile) Design */ options ls=80 ps=60 nodate nonumber; title 'Output 3.5.1: Split-Plot Repeated Measures Design'; data split; infile 'c:\3_5_1.dat'; input y x111 x121 x131 x112 x122 x132 x211 x221 x231 x212 x222 x232 x311 x321 x331 x312 x322 x332; proc print data=split; run; /* Using full rank model */ proc reg; title2 'Full Rank Model'; model y = x111 x121 x131 x112 x122 x132 x211 x221 x231 x212 x222 x232 x311 x321 x331 x312 x322 x332 /noint; restrict x111 - x121 - x112 + x122 = 0, x121 - x131 - x122 + x132 = 0, x211 - x221 - x212 + x222 = 0, x221 - x231 - x222 + x232 = 0, x311 - x321 - x312 + x322 = 0, x321 - x331 - x322 + x332 = 0; A: mtest x111+x121+x131+x112+x122+x132-x211-x221-x231-x212-x222-x232=0, x211+x221+x231+x212+x222+x232-x311-x321-x331-x312-x322-x332=0 /print details; B: mtest x111-x131+x112-x132+x211-x231+x212-x232+x311-x331+x312-x332=0, x121-x131+x122-x132+x221-x231+x222-x232+x321-x331+x322-x332=0 /print details; AB: mtest x111-x121+x112-x122-x221+x231-x222+x232=0, x121-x131+x122-x132-x221+x231-x222+x232=0, x211-x221+x212-x222-x311+x321-x312+x322=0, x221-x231+x222-x232-x321+x331-x322+x332=0 /print details; AS: mtest x111+x121+x131-x112-x122-x132=0, x211+x221+x231-x212-x222-x232=0, x311+x321+x331-x312-x322-x332=0 /print details; run; /* Using less than full rank model */ data split2; input a b1 b2 b3; cards; 1 3 4 3 1 2 2 1 2 3 7 7 2 5 4 6 3 3 4 6 3 2 3 5 ; run; proc glm; title2 'Less than Full Rank Model'; class a; model b1 b2 b3 = a; repeated b 3 profile /nom summary; run; *************************************************** /* The following code appears on page 109, */ /* where it is used to produce output 3_5_2. */ *************************************************** /* Program 3_5_2.sas */ /* Split-Plot (Two-Group Profile) Design, Full Rank Model */ /* Timm (1975, p244) data */ options ls=80 ps=60 nodate nonumber; title 'Output 3.5.2: Two Group Profile Analysis'; data splitb; infile 'c:\3_5_2.dat'; input grp p1 p2 p3 p4 p5; proc print data=splitb; run; proc glm; class grp; model p1 p2 p3 p4 p5 = grp; repeated position 5 profile /nom summary; run; *************************************************** /* The following code appears on page 111, */ /* where it is used to produce output 3.5.3. */ *************************************************** /* Program 3_5_3.sas */ /* Tests of Covariance Structures */ options ls=80 ps=60 nodate nonumber mprint; title 'Output 3.5.3: Tests of Covariance Structures'; data struct; infile 'c:\3_5_2.dat'; input grp p1 p2 p3 p4 p5; proc print data=struct; run; proc iml; use struct; read all var {p1 p2 p3 p4 p5} where (grp=1) into y1; read all var {p1 p2 p3 p4 p5} where (grp=2) into y2; g=2; n1=nrow(y1); n2=nrow(y2); n=n1+n2; p=ncol(y1); df1=n1-1; df2=n2-1; s1=(y1`*(i(n1)-(1/n1)*j(n1,n1))*y1)/df1; s2=(y2`*(i(n2)-(1/n2)*j(n2,n2))*y2)/df2; s=(1/(n-g))*(df1*s1+df2*s2); /* the hypothesis test matrix a` */ ap={1 -1 0 0 0, 0 1 -1 0 0, 0 0 1 -1 0, 0 0 0 1 -1}; a=ap`; q=nrow(ap); call gsorth(as,t,lindep,a); asp=as`; as1a=as`*s1*as; as2a=as`*s2*as; asa=(1/(n-g))*(df1*as1a+df2*as2a); print "A` and A*`", ap asp; print s1 s2; print "Reduced Covariance Matrices", as1a as2a; print "Reduced Pooled Covaraince Matrix", asa; /* Step 1 : Test of Equality of the Reduced Covariance Matrices */ m1=(n-g)*(log(det(asa))); m2=df1*(log(det(as1a)))+df2*(log(det(as2a))); m=m1-m2; print m; e1=(2*q*q+3*q-1)/(6*(q+1)*(g-1)); e2=(1/df1)+(1/df2)-(1/(n-g)); e=e1*e2; print e; xb=(1-e)*m; v1=(q*(q+1)*(g-1))/2; probxb=1-probchi(xb,v1); print "Test of Equality of Reduced Covariance Matrices", xb v1; print "Probability of Xb with df=v1", probxb; print "The above test works well for ni>20 and q<6 and g<6"; eo1=(1/(df1*df1))+(1/(df2*df2))-(1/((n-g)*(n-g))); eo2=(q-1)*(q+2)/(6*g-6); eo=eo2*eo1; vo=(v1+2)/(eo-(e*e)); fb=((1-e-(v1/vo))/v1)*m; probfb=1-probf(fb,v1,vo); print "Test of Equality of Reduced Covariance Matrices", fb v1 vo; print "Probability of Fb with df=v1,vo", probfb; print "The above test can be used for small ni and/or q and/or g >6"; /* Step 2 : Test of Circularity */ const=-( (n-1)-((2*q*q+q+2)/(6*q)) ); lds=log(det(asa)); qlt=q*(log((trace(asa))/q)); print const lds qlt; xs=const*(lds-qlt); dfxs=((q*q+q)/2)-1; probxs=1-probchi(xs,dfxs); print "Tests of Circularity (Mauchly)", xs dfxs; print "Probability of Xs with df=dfxs", probxs; sasa=asa*asa; v=trace(sasa)/((trace(asa))*(trace(asa))); vs=((q*q*n)/2)*(v-(1/q)); probvs=1-probchi(vs,dfxs); print "Test of Circularity (L.B.I.)", vs, dfxs; print "Probability of Vs with df=dfxs", probvs; quit; run; *************************************************** /* The following code appears on page 119, */ /* where it is used to produce output 3.6.1. */ *************************************************** /* Program 3_6_1.sas */ /* Program to perform ANCOVA with one covariate and */ /* Unequal n, using Full Rank model */ options ls=80 ps=60 nodate nonumber; filename examp4 'c:\2_9.dat'; title 'Output 3.6.1: ANCOVA, One Covariate, Unequal n'; data onecov; infile examp4; input y a z a1 a2 a3 z1 z2 z3; proc print; run; /* Using the full rank model */ proc reg; title2 'Full Rank Model'; model y = a1 a2 a3 z1 z2 z3 /noint; restrict z1-z2=0, z1-z3=0; A: mtest a1-a3=0, a2-a3=0 /print details; REG: mtest z3=0 /print details; run; /* Using the less than full rank model */ proc glm; title2 'Less than Full Rank Model'; class a; model y=a z /e xpx solution; lsmeans a /stderr pdiff cov; run; *************************************************** /* The following code appears on page 124, */ /* where it is used to produce output 3.6.2. */ *************************************************** /* Program 3_6_2.sas */ /* Program to perform ANCOVA with two covariates and */ /* Unequal n, using Full Rank model */ options ls=80 ps=60 nodate nonumber; title 'Output 3.6.2: ANCOVA, Two Covariates, Unequal n'; data twocov; infile 'c:\3_6_2.dat'; input y a z1 z2 a1 a2 a3 z11 z12 z13 z21 z22 z23; proc print; run; /* Using the full rank model */ proc reg; title2 'Full Rank Model'; model y = a1 a2 a3 z11 z12 z13 z21 z22 z23 /noint; parallel: mtest z11-z13=0, z12-z13=0, z21-z23=0, z22-z23=0 /print details; run; proc reg; model y = a1 a2 a3 z11 z12 z13 z21 z22 z23 /noint; restrict z11-z13=0, z12-z13=0, z21-z23=0, z22-z23=0; A: mtest a1-a2=0, a2-a3=0 /print details; REG: mtest z13=0, z23=0 /print details; run; /* Using the less than full rank model */ proc glm; title2 'Less than full rank model'; class a; model y=a z1 z2 a*z1 a*z2 /e solution; run; proc glm; class a; model y=a z1 z2 /e xpx solution; lsmeans a /stderr pdiff cov; run; *************************************************** /* The following code appears on page 143, */ /* where it is used to produce output 4.5. */ *************************************************** /* Program 4_5.sas */ /* WLSE for data with heteroscedasticity */ /* Data are from Neter, Wasserman, and Kutner, 1990, p. 421 */ options ls=78 ps=60 nodate nonumber; title1 'Output 4.5: Data with Heteroscedasticity'; data heter; infile 'c:\4_5.dat'; input age blpress; if age lt 30 then agegrp=1; else if age lt 40 then agegrp=2; else if age lt 50 then agegrp=3; else agegrp=4; run; proc sort; by agegrp; run; proc means mean std var n; var blpress; by agegrp; run; /* Ordinary Least Square Regression */ proc reg data=heter; title2 'Ordinary Least Squares Regression '; model blpress=age/ p cli; output out=resid1 r=olsresid; run; /* Plot of OLS residuals */ filename out1 'c:\4_5_1.cgm'; goptions device=cgmmwwc gsfname=out1 gsfmode=replace colors=(black) vsize=5in hsize=4in htitle=1; proc gplot data=resid1; title1 'OLS Residuals'; title2; plot olsresid*age /frame; symbol1 v=; run; /* Add weight vector to the SAS dataset */ data wlse; set heter; if agegrp=1 then w=(1/sqrt(24.9231)); else if agegrp=2 then w=(1/sqrt(50.5256)); else if agegrp=3 then w=(1/sqrt(96.3524)); else if agegrp=4 then w=(1/sqrt(133.1410)); run; /* Weighted Least Squares Regression */ proc reg data=wlse; title1 'Output 4.5: Data with Heteroscedasticity'; title2 'Weighted Least Squares Regression '; model blpress=age/ p cli; weight w; output out=resid2 r=wlsresid; run; data resid3; set resid2; wwlsresi=w*wlsresid; label wwlsresi = 'weighted residual'; run; /* Plot of weighted WLS residuals */ filename out2 'c:\4_5_2.cgm'; goptions device=cgmmwwc gsfname=out2 gsfmode=replace colors=(black) vsize=5in hsize=4in htitle=1; proc gplot data=resid3; title1 'Weighted WLS Residuals'; title2; plot wwlsresi*age /frame; symbol1 v=; run; *************************************************** /* The following code appears on page 156, */ /* where it is used to produce output 4.7.2. */ *************************************************** /* Program 4_7_2.sas */ /* Data are from Marascuilo and McSweeney, 1977, p.174 */ /* Using PROC CATMOD to evaluate Marginal Homogeniety */ options ls=78 ps=60 nodate nonumber; title1 'Output 4.7.2: Test of Marginal Homogeniety'; /* Using PROC IML */ title2 'Using PROC IML'; proc iml; start cat1; /* Input variables and data */ n=124; s=3; r=3; freq={50 9 10, 2 6 29, 0 0 18}; print freq; /* Form proportions and shape p */ prop=freq/n; p=shape(prop,r*s,1); print p; /* Construct Covariance matrix and test matrix */ omega=(diag(p)-p*p`)/n; qr=2; A={0 1 1 -1 0 0 -1 0 0, 0 -1 0 1 0 1 0 -1 0}; S=A*Omega*A`; print S; InvS=Inv(S); print InvS; F=A*p; print F; /* Calculate Wald statistic and p-value */ Q=F`*InvS*F; pv=1-probchi(Q,qr); print 'Chi-square Test of Marginal Homogeniety',Q, 'p-value',pv; finish cat1; run cat1; quit; run; /* Using PROC CATMOD */ data school; input t55 t65 count @@; cards; 1 1 50 1 2 9 1 3 10 2 1 2 2 2 6 2 3 29 3 1 0 3 2 0 3 3 18 ; run; proc catmod; title2 'Using PROC CATMOD'; weight count; response marginals; model t55*t65=_response_; repeated year 2; run; *************************************************** /* The following code appears on page 162, */ /* where it is used to produce output 4.7.3. */ *************************************************** /* Program 4_7_3.sas */ /* Data are from Grizzle et al.,1969 */ /* Testing for Homogeniety of Proportions */ options ls=78 ps=60 nodate nonumber; title1 'Output 4.7.3: Testing Homogeniety of Proportions, One-Way ANOVA'; /* Using PROC IML */ title2 'Using PROC IML'; proc iml; start cat2; freq={23 7 2, 23 10 5, 20 13 5, 24 10 6, 18 6 1, 18 6 2, 13 13 2, 9 15 2, 8 6 3, 12 4 4, 11 6 2, 7 7 4, 12 9 1, 15 3 2, 14 8 3, 13 6 4}; s=nrow(freq); r=ncol(freq); wt={1 2 3}; rowsum=freq[,+]; prob=freq/(rowsum*repeat(1,1,r)); p=shape(prob[,1:r],0,1); A=I(16)@(repeat(1,1,r)#wt); D=diag(rowsum); V=(diag(p)-p*p`)#(Inv(D)@repeat(1,r,r)); S=A*V*A`; InvS=Inv(S); F=A*p; print F; X={1 1 0 0 1 0 0, 1 1 0 0 0 1 0, 1 1 0 0 0 0 1, 1 1 0 0 -1 -1 -1, 1 0 1 0 1 0 0, 1 0 1 0 0 1 0, 1 0 1 0 0 0 1, 1 0 1 0 -1 -1 -1, 1 0 0 1 1 0 0, 1 0 0 1 0 1 0, 1 0 0 1 0 0 1, 1 0 0 1 -1 -1 -1, 1 -1 -1 -1 1 0 0, 1 -1 -1 -1 0 1 0, 1 -1 -1 -1 0 0 1, 1 -1 -1 -1 -1 -1 -1}; Beta=Inv(X`*InvS*X)*X`*InvS*F; print Beta; Q=F`*InvS*F-Beta`*(X`*InvS*X)*Beta; dfq=9; pvq=1-probchi(Q,dfq); print 'Chi-square Test of Fit',Q,'p-value for fit',pvq; CT={ 0 0 0 0 1 0 0, 0 0 0 0 0 1 0, 0 0 0 0 0 0 1}; dfct=3; CH={ 0 1 0 0 0 0 0, 0 0 1 0 0 0 0, 0 0 0 1 0 0 0}; dfch=3; L={0 0 0 0 3 2 1}; dfL=1; Wh=(CH*Beta)`*Inv(CH*Inv(X`*InvS*X)*CH`)*(CH*Beta); pvWh=1-probchi(Wh,dfch); Wt=(CT*Beta)`*Inv(CT*Inv(X`*InvS*X)*CT`)*(CT*Beta); pvWt=1-probchi(Wt,dfct); W=(L*Beta)`*Inv(L*Inv(X`*InvS*X)*L`)*(L*Beta); pvL=1-probchi(W,dfL); print 'Hosp Chi-Square Test',WH, 'p-value Hosp',pvWh, 'Treat Chi-Square Test',WT,'p-value Treat',pvWt, 'Linear Trend',W,'p-value Linear',pvL; finish cat2; run cat2; quit; run; /* Using PROC CATMOD */ data operate; input treat hosp $ severity $ wt @@; cards; 1 a N 23 1 a S 7 1 a M 2 1 b N 18 1 b S 6 1 b M 1 1 c N 8 1 c S 6 1 c M 3 1 d N 12 1 d S 9 1 d M 1 2 a N 23 2 a S 10 2 a M 5 2 b N 18 2 b S 6 2 b M 2 2 c N 12 2 c S 4 2 c M 4 2 d N 15 2 d S 3 2 d M 2 3 a N 20 3 a S 13 3 a M 5 3 b N 13 3 b S 13 3 b M 2 3 c N 11 3 c S 6 3 c M 2 3 d N 14 3 d S 8 3 d M 3 4 a N 24 4 a S 10 4 a M 6 4 b N 9 4 b S 15 4 b M 2 4 c N 7 4 c S 7 4 c M 4 4 d N 13 4 d S 6 4 d M 4 ; run; proc catmod order=data; title2 'Using PROC CATMOD' ; weight wt; response 1 2 3; model severity = hosp treat/ freq oneway; contrast 'Linear Trend - Treatment' treat 3 2 1; run; *************************************************** /* The following code appears on page 169, */ /* where it is used to produce output 4.7.4. */ *************************************************** /* Program 4_7_4.sas */ /* Testing for Independence */ options ls=78 ps=60 nodate nonumber; title 'Output 4.7.4: Test of Independence ' ; data independ; input office $ perform $ count @@; cards; Y Success 50 Y Usuccess 20 N Success 10 N Usuccess 20 ; run; proc catmod; weight count; model office*perform=_response_/ freq wls; loglin office perform; run; *************************************************** /* The following code appears on page 191, */ /* where it is used to produce output 5.6. */ *************************************************** /* Program 5_6.sas */ /* Rhower data from Timm(1975), pp.281, 345 */ options ls = 80 ps = 60 nodate nonumber; title1 'Output 5.6: Multivariate Regression'; data rhower; infile 'c:\5_6.dat'; input y1-y3 x0-x5; proc print data=rhower; run; /* Calculations for Regression and Multivariate Cooks Distance */ proc iml; title2 'Using PROC IML, including Cooks Distance'; use rhower; v={y1 y2 y3}; w={x0 x1 x2 x3 x4 x5}; read all var v into y; read all var w into x; beta=inv(x`*x)*x`*y; print 'Regression Coefficients' beta; n=nrow(y); p=ncol(y); k=ncol(x); S=(y`*y-y`*x*beta)/(n-k); print 'Estimated Covarianc Matrix' S; v=inv(k*s)@(x`*x); b=shape(beta`,p*k,1); m=n-1; x1=x[1:m,1:k]; y1=y[1:m,1:p]; beta1=inv(x1`*x1)*x1`*y1; b1=shape(beta1`,p*k,1); d1=(b-b1)`*v*(b-b1); obs=32; index=obs; do i=2 to 31; j=n-i; x2=x[1:j,1:k]; x1=x[j+2:n,1:k]; x2=x2//x1; y2=y[1:j,1:p]; y1=y[j+2:n,1:p]; y2=y2//y1; beta2=inv(x2`*x2)*x2`*y2; b2=shape(beta2`,p*k,1); d2=(b-b2)`*v*(b-b2); d1=d1//d2; obs=obs-1; index=index//obs; end; x1=x[2:n,1:k]; y1=y[2:n,1:p]; beta1=inv(x1`*x1)*x1`*y1; b1=shape (beta1`,p*k,1); d2=(b-b1)`*v*(b-b1); d1=d1//d2; obs=1; index=index//obs; D=index||d1; print 'Influence Obs # and Multivariate Cooks Distance', D; quit; /* Multivariate Regression using PROC REG */ proc reg data=rhower; title2 'Using PROC REG'; /* Full Model Analysis */ model y1-y3 = x1-x5; Gammaa:mtest x1-x5/ print; /* test that all coefficients equal zero */ Gamma1:mtest x1; /* test that x1 equals zero */ Gamma2:mtest x2; /* test that x2 equals zero */ Gamma3:mtest x3; /* test that x3 equals zero */ Gamma4:mtest x4; /* test that x4 equals zero */ Gamma5:mtest x5; /* test that x5 equals zero */ Beta:mtest intercept,x1-x5; /* test that int and all coeff are zero */ /* Reduced Model Analysis*/ model y1-y3 = x2-x4; Gamma:mtest x2-x4; /* test that gamma=0 for reduced model */ Gamma2:mtest x2; Gamma3:mtest x3; Gamma4:mtest x4; run; /* Calculation of Multivariate Eta Squared for Full and Reduced Models */ proc iml; title2 'Multivariate Eta Squared for Full and Reduced Models'; n=32; g1=6; g2=4; p=3; LmdaF=.81206193; LmdaR=.70761904; Full=1-n*LmdaF/(n-g1+LmdaF); Reduced=1-n*LmdaR/(n-g2+LmdaR); print 'Eta Squared full model' Full, 'Eta Squared reduced model' Reduced; quit; /* Multivariate Regression using PROC GLM and Full Rank Model */ proc glm data=rhower; title2 'Using PROC GLM'; model y1-y3=x2-x4/ nouni; manova h=x2-x4/ printe printh; run; *************************************************** /* The following code appears on page 208, */ /* where it is used to produce output 5.7. */ *************************************************** /* Program 5_7.sas */ /* Nonorthogonal MANOVA Design */ options ls=80 ps=60 nodate nonumber; title1 'Output 5.7: Nonorthogonal Three Factor MANOVA'; data three; infile 'c:\5_7.dat'; input a b c y1 y2; proc print data=three; run; /* Unweighted MANOVA Analysis */ proc glm; title2 'Unweighted Analysis'; class a b c; model y1 y2 = a b c a*b a*c b*c a*b*c/ ss3 e; contrast 'a1 vs a3' a 1 0 -1; estimate 'a1 vs a3' a 1 0 -1/ divisor=3 e; contrast 'c1 vs c2' c .5 -.5; estimate 'c1 vs c2' c .5 -.5/ e; lsmeans a b c a*b a*c b*c a*b*c; manova h=a|b|c; run; /* Full Rank Model for MANOVA Design */ proc glm; title2 'Full Rank Model with Unweighted Contrast'; class a b c; model y1 y2 = a*b*c/ noint; contrast 'a*b at c1' a*b*c 1 0 -1 0 0 0 -1 0 1 0 0 0 0 0 0 0 0 0, a*b*c 0 0 1 0 -1 0 0 0 -1 0 1 0 0 0 0 0 0 0, a*b*c 0 0 0 0 0 0 1 0 -1 0 0 0 -1 0 1 0 0 0, a*b*c 0 0 0 0 0 0 0 0 1 0 -1 0 0 0 -1 0 1 0/ e; manova h=a*b*c/ printe printh; run; /* Weighted MANOVA Analysis */ proc glm data=three; title2 'Weighted Analysis'; class a b c; model y1 y2 = a b a*b c a*c b*c a*b*c/ ss1; means a b c a*b a*c b*c a*b*c; manova h=a a*b a*b*c; run; proc glm; class a b c; model y1 y2 = b a a*b c c*b c*a a*b*c/ ss1; manova h=b; run; proc glm; class a b c; model y1 y2 = c b b*c a a*c a*b a*b*c/ ss1; manova h= c c*b; run; proc glm; class a b c; model y1 y2 = a c a*c b a*b b*c a*b*c/ ss1; manova h=a*c; run; /* Cell Means Model for MANOVA Design */ proc glm; title2 'Full Rank Model with Weighted Contrast'; class a b c; model y1 y2 = a*b*c/ noint nouni; contrast 'a*b at c1'a*b*c 1 0 -1 0 0 0 -1 0 1 0 0 0 0 0 0 0 0 0, a*b*c 0 0 1 0 -1 0 0 0 -1 0 1 0 0 0 0 0 0 0, a*b*c 0 0 0 0 0 0 1 0 -1 0 0 0 -1 0 1 0 0 0, a*b*c 0 0 0 0 0 0 0 0 1 0 -1 0 0 0 -1 0 1 0/ e; contrast 'a1 vs a3'a*b*c .2143 .1429 .1429 .2143 .1429 .1429 0 0 0 0 0 0 -.1333 -.0667 -.200 -.200 -.200 -.200/ e; contrast 'c1 vs c2'a*b*c .1429 -.1053 .0952 -.1579 .0952 -.1053 .0952 -.0526 .0476 -.1053 .1429 -.1053 .0952 -.0526 .1429 -.1579 .1429 -.1579/ e; manova h=a*b*c/printe printh; run; *************************************************** /* The following code appears on page 219, */ /* where it is used to produce output 5.8. */ *************************************************** /* Program 5_8.sas */ /* MANCOVA Designs */ options ls=80 ps=60 nodate nonumber; title 'Output 5.8: Multivariate MANCOVA'; data mancova; infile 'c:\5_8.dat'; input drugs $ apgar1 apgar2 x; proc print; run; proc glm; title2 'Test of Parallelism'; class drugs; model Apgar1 Apgar2 = drugs x drugs*x; manova h= drugs*x; run; proc glm; title2 'MANCOVA Tests'; class drugs; model apgar1 apgar2 = drugs x/ solution; contrast '1 vs. 4' drugs 1 0 0 -1; manova h=drugs x; means drugs; lsmeans drugs/ stderr pdiff tdiff cov out=adjmeans; run; proc print data=adjmeans; title2 'Adjusted Means'; run; *************************************************** /* The following code appears on page 225, */ /* where it is used to produce output 5.9. */ *************************************************** /* Program 5_9.sas */ /* Data from Smith, Gnanadesikan and Hughes (1962) */ /* Stepdown Analysis--MANCOVA */ options ls=80 ps=60 nodate nonumber; title1 'Output 5.9: Multivariate MANCOVA with Stepdown Analysis'; data mancova; infile 'c:\5_9.dat'; input group $ y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 x1 x2; proc print data=mancova; run; /* Test of Parallelism */ proc glm data=mancova; title2 'Test of Parallelism'; class group; model y1-y11 = group x1 x2 group*x1 group*x2/ nouni ss3; manova h = group *x1 group*x2; run; /* MANCOVA Tests */ proc glm data=mancova; title2 'MANCOVA Tests'; class group; model y1-y11 = x1 x2 group/ e nouni ss3; contrast 'Ov Mean' intercept 1/ e; manova h=group x1 x2/ printh; means group; lsmeans group/ stderr pdiff tdiff cov out=adjmeans; run; proc print data=adjmeans; title2 'Adjusted Means'; run; /* Mancova stepdown analysis for additional information */ proc glm data=mancova; class group; title2 'Test of Gamma 2.1'; model y2 y3 y4 y6 y9-y11 = group y1 y5 y7 y8 x1 x2/ nouni ss3; manova h = group y1 y5 y7 y8 x1 x2; run; proc glm data=mancova; class group; title2 'Test of Gamma1'; model y1 y5 y7 y8 = group x1 x2/ nouni ss3; manova h = group x1 x2; run; *************************************************** /* The following code appears on page 231, */ /* where it is used to produce output 5.10. */ *************************************************** /* Program 5_10.sas */ /* Data from Timm(1975, page 244) */ /* Repeated Measures Analysis */ options ls=80 ps=60 nodate nonumber; title1 'Output 5.10: Repeated Measurement Analysis'; data timm; infile 'c:\5_10.dat'; input group y1 y2 y3 y4 y5; proc print data=timm; run; proc glm; title2 'Multivariate Tests'; class group; model y1-y5 = group/ intercept nouni; means group; manova h = group/ printe printh; /*test of group diffs for means */ manova h =_all_ m=(1 -1 0 0 0, 0 1 -1 0 0, 0 0 1 -1 0, 0 0 0 1 -1) prefix = diff/ printe printh; /* test for parallel profiles */ proc glm; class group; model y1-y5= group/noint nouni; contrast 'Mult Cond' group 1 0, group 0 1; manova m=(1 -1 0 0 0, 0 1 -1 0 0, 0 0 1 -1 0, 0 0 0 1 -1) prefix = diff/ printe printh; /* test of condtions as vectors */ run; proc glm; title2 'Tests given Parallelism of Profiles'; class group; model y1 - y5 = group/ noint nouni; contrast 'Univ Gr' group 1 -1; manova m=(.2 .2 .2 .2 .2) prefix=Gr/ printe printh; /* test of group means given parallel profiles */ run; /* Univariate Test given Parallelism and Sphericity */ proc glm; title2 'Univariate Test given Parallelism and Sphericity'; class group; model y1 - y5 = group/ nouni; repeated cond 5 (1 2 3 4 5) profile/printm summary; manova h=group m=(1 -1 0 0 0, 0 1 -1 0 0, 0 0 1 -1 0, 0 0 0 1 -1) prefix = diff / printe printh; run; *************************************************** /* The following code appears on page 243, */ /* where it is used to produce output 5.11. */ *************************************************** /* Program 5_11.sas */ /* Data from Timm(1975, page 454) */ /* Extended Linear Hypotheses */ options ls=80 ps=60 nodate nonumber; title1 'Output 5.11: Extended Linear Hypotheses'; data timm; infile 'c:\5_11.dat'; input group y1 y2 y3 x1 x2 x3; proc print data=timm; run; proc glm data=timm; title2 'Multivariate Test of Group--Using PROC GLM'; class group; model y1-y3 = group/nouni; means group; manova h=group/printh printe; run; /* IML procedure for Extended Linear Hypothesis */ title2 'Multivariate Test of Group--Using PROC IML'; proc iml; use timm; a={x1 x2 x3}; b={y1 y2 y3}; read all var a into x; read all var b into y; beta=inv(x`*x)*x`*y; print beta; n=nrow(y); p=ncol(y); k=ncol(x); nu_h=2; u=3; nu_e=n-k; s0=min(nu_h,u); r=max(nu_h,u); alpha=.05; denr=(nu_e-r+nu_h); roy_2=(r/denr)*finv(1-alpha,r,denr); rvalue=sqrt(roy_2); m0=(abs(nu_h-u)-1)/2; n0=(nu_e-u-1)/2; num=s0**2*(2*m0+s0+1); dent=2*(s0*n0+1);df=s0*(2*m0+s0+1); t0_2=(num/dent)*finv(1-alpha,df,dent); tovalue=sqrt(t0_2); print s0 m0 n0; e=(y`*y-y`*x*beta); co={1 -1 0, 0 1 -1}; ao=i(3); eo=ao`*e*ao; bo=co*beta*ao; wo=co*inv(x`*x)*co`; ho=bo`*inv(wo)*bo; print,"Overall Error Matrix", eo , "Overall Hypothesis Test Matrix",ho; /* c`c=eo where c is upper triangle Cholesky matrix */ c=root(eo); f=inv(c`)*ho*inv(c); eig=Eigval(round(f,.0001)); vec=inv(c)*eigvec(round(f,.0001)); print,"Eigenvalues & Eigenvectors of Overall Test of Ho (Groups)", eig vec; /* Extended Linear Hypothesis following Overall Group Test */ m={1 .5, -.5 .5, -.5 -1}; g=ao*m*co; print, "Extended Linear Hypothesis Test Matrix",m g; psi=m*bo; psi_hat=trace(psi); tr_psi=abs(psi_hat); h=m*wo*m`; einv=inv(eo); c=root(einv); f=inv(c`)*h*inv(c); xeig=Eigval(round(f,.0001)); print, "Eigenvalues of Extended Linear Hypothesis", xeig; denrt=sum(sqrt(xeig)); dentr=sqrt(sum(xeig)); to_2=tr_psi/dent; print, "Extended To**2 Statistic", to_2; print, "Extended To**2 Critical Value", tovalue; root=tr_psi/denr; print, "Extended Largest Root Statistic", root; print, "Extended Largest Root Critical Value", rvalue; print psi_hat alpha; ru=psi_hat+rvalue*denrt; rl=psi_hat-rvalue*denrt; vu=psi_hat+tovalue*dentr; vl=psi_hat-tovalue*dentr; print 'Approximate Simultaneous Confidence Intervals'; print 'Contrast Significant if interval does not contain zero'; print 'Extended Root interval: ('rl ',' ru ')'; print 'Extended Trace interval: ('vl ',' vu ')'; /* Multiple Extended Linear Hypothesis using To**2 */ m1={1 0,0 0,0 0}; m2={0 0,1 0,0 0}; m3={0 0,0 1,0 0}; m4={0 0,0 0,0 1}; print,"Multiple Extended Linear Hypothesis Test Matrices", m1,m2,m3,m4; g1=ao*m1*co; g2=ao*m2*co; g3=ao*m3*co; g4=ao*m4*co; t1=trace(m1*bo); t2=trace(m2*bo); t3=trace(m3*bo); t4=trace(m4*bo); tau=t1//t2//t3//t4; t11=trace(m1*wo*m1`*eo); t21=trace(m2*wo*m1`*eo); t22=trace(m2*wo*m2`*eo); t31=trace(m3*wo*m1`*eo); t32=trace(m3*wo*m2`*eo); t33=trace(m3*wo*m3`*eo); t41=trace(m4*wo*m1`*eo); t42=trace(m4*wo*m2`*eo); t43=trace(m4*wo*m3`*eo); t44=trace(m4*wo*m4`*eo); r1=t11||t21||t31||t41; r2=t21||t22||t32||t42; r3=t31||t32||t33||t43; r4=t41||t42||t43||t44; t=r1//r2//r3//r4; print tau,t; to_4=tau`*inv(t)*tau; print, "Extended Linear Hypothesis Criterion To**2 Squared", to_4; print, "Extended To**2 Critical Value", t0_2; quit; /* Multivariate test of Parallelism */ proc glm data=timm; title2 'Multivariate Test of Parallelism--Using PROC GLM'; class group; model y1-y3 = group/nouni; manova h = group m = ( 1 -1 0, 0 1 -1) prefix = diff/ printe printh; run; title2 'Multivariate Test of Parallelism--Using PROC IML'; proc iml; use timm; a={x1 x2 x3}; b={y1 y2 y3}; read all var a into x; read all var b into y; beta=inv(x`*x)*x`*y; n=nrow(y); p=ncol(y); k=ncol(x); nu_h=2; u=2; nu_e=n-k; s0=min(nu_h,u); r=max(nu_h,u); alpha=.05; denr=(nu_e-r+nu_h); roy_2=(r/denr)*finv(1-alpha,r,denr); rvalue=sqrt(roy_2); m0=(abs(nu_h-u)-1)/2; n0=(nu_e-u-1)/2; num=s0**2*(2*m0+s0+1); dent=2*(s0*n0+1); df=s0*(2*m0+s0+1); t0_2=(num/dent)*finv(1-alpha,df,dent); tovalue=sqrt(t0_2); print s0 m0 n0; e=(y`*y-y`*x*beta); co={1 -1 0, 0 1 -1}; ao={1 0, -1 1, 0 -1}; eo=ao`*e*ao; bo=co*beta *ao; wo=co*inv(x`*x)*co`; ho=bo`*inv(wo)*bo; c=root(eo); f=inv(c`)*ho*inv(c); eig=eigval(round(f,.0001)); vec=inv(c)*eigvec(round(f,.0001)); print,"Eigenvalues & Eigenvectors of Overall test of Ho (Parallelism)", eig vec; /* Extended Linear Hypothesis following overall Parallelism test */ m={0 1,1 0}; g=ao*m*co; print, "Extended Linear Hypothesis Test Matrix", m g; psi=m*bo; psi_hat=trace(psi); tr_psi=abs(psi_hat); h=m*wo*m`; einv=inv(eo); c=root(einv); f=inv(c`)*h*inv(c); xeig=eigval(round(f,.0001)); print, "Eigenvalues of Extended Linear Hypothesis", xeig; to_2=tr_psi/sqrt(sum(xeig)); print, "Extended To**2 Statistic", to_2; print,"Extended To**2 Critical Value", tovalue; root=tr_psi/sum(sqrt(xeig)); print, "Extended Largest Root Statistic", root; print, "Extended Largest Root Critical Value", rvalue; print psi_hat alpha; ru=psi_hat+rvalue*sum(sqrt(xeig)); rl=psi_hat-rvalue*sum(sqrt(xeig)); vu=psi_hat+tovalue*sqrt(sum(xeig)); vl=psi_hat-tovalue*sqrt(sum(xeig)); print 'Approximate Simultaneous Confidence Intervals'; print 'Contrast Significant if interval does not contain zero'; print 'Extended Root Interval: ('rl ',' ru ')'; print 'Extended Trace Interval: ('vl '.' vu ')'; quit; /* Multivariate test of Conditions as vectors */ proc glm data=timm; title2 'Multivariate Test of Conditions--Using PROC GLM'; class group; model y1-y3 = group/noint nouni; contrast 'Mult Cond' group 1 0 0, group 0 1 0, group 0 0 1; manova m=(1 -1 0, 0 1 -1) prefix = diff/ printe printh; run; title2 'Multivariate Test of Conditions--Using PROC IML'; proc iml; use timm; a={x1 x2 x3}; b={y1 y2 y3}; read all var a into x; read all var b into y; beta=inv(x`*x)*x`*y; n=nrow(y); p=ncol(y); k=ncol(x); nu_h=3; u=2; nu_e=n-k; s0=min(nu_h,u); r=max(nu_h,u); alpha=.05; denr=(nu_e-r+nu_h); roy_2=(r/denr)*finv(1-alpha,r,denr); rvalue=sqrt(roy_2); m0=(abs(nu_h-u)-1)/2; n0=(nu_e-u-1)/2; num=s0**2*(2*m0+s0+1); dent=2*(s0*n0+1); df=s0*(2*m0+s0+1); t0_2=(num/dent)*finv(1-alpha,df,dent); tovalue=sqrt(t0_2); print s0 m0 n0; e=(y`*y-y`*x*beta); co=i(3); ao={1 0, -1 1, 0 -1}; eo=ao`*e*ao; bo=co*beta*ao; wo=co*inv(x`*x)*co`; ho=bo`*inv(wo)*bo; c=root(eo); f=inv(c`)*ho*inv(c); eig=eigval(round(f,.0001)); vec=inv(c)*eigvec(round(f,.0001)); print, "Eigenvalues & Eigenvectors of Overall test of Ho (Conditions)", eig vec; m={1 0 1, 0 1 1}; g=ao*m*co; print, "Extended Linear Hypothesis Test Matrix", m g; psi=m*bo; psi_hat=trace(psi); tr_psi=abs(psi_hat); h=m*wo*m`; einv=inv(eo); c=root(eo); f=inv(c`)*h*inv(c); xeig=eigval(round(f,.0001)); print, "Eigenvalues of Extended Linear Hypothesis", xeig; to_2=tr_psi/sqrt(sum(xeig)); print, "Extended To**2 Statistic", to_2; print, "Extended T0**2 Critical Value", tovalue; root= tr_psi/sum(sqrt(xeig)); print, "Extended Largest Root Statistic", root; print, "Extended Largest Root Critical Value", rvalue; print psi_hat alpha; ru=psi_hat+rvalue*sum(sqrt(xeig)); rl=psi_hat-rvalue*sum(sqrt(xeig)); vu=psi_hat+tovalue*sqrt(sum(xeig)); vl=psi_hat-tovalue*sqrt(sum(xeig)); print 'Approximate Simultaneous Confidence Intervals'; print 'Contrast Significant if interval does not contain zero'; print 'Extended Root Interval: ('rl ',' ru ')'; print 'Extended Trace Interval: ('vl ',' vu ')'; quit; *************************************************** /* The following code appears on page 253, */ /* where it is used to produce output 5.12. */ *************************************************** /* Program 5_12.sas */ /* Data are from Timm(1975, p.264) */ /* The SAS routine power.pro was supplied by Muller, La Vange */ /* Ramey, and Ramey (1992) as shown in their JASA (1992) article */ /* "Power calculations for general linear models including repeated */ /* measures applications", pp 1209-1226 */ options ls=80 ps=60 nodate nonumber; title 'Output 5.12: Power Analysis for MANOVA design'; proc iml symsize=1000; %include 'c:\power.pro'; alpha = .01; c={1 -1}; u=i(2); beta={20.38 -38.33, 0 0}; sigma={307.08 280.83, 280.83 421.67}; essencex={1 0, 1 0, 1 0, 1 0, 1 0, 0 1, 0 1, 0 1, 0 1, 0 1}; run power; quit; *************************************************** /* The following code appears on page 263, */ /* where it is used to produce output 6.5. */ *************************************************** /* Program 6_5.sas */ /* Example from Timm(1980b) and Boik(1988,1991) - Zullo Dental data */ options ls=80 ps=60 nodate nonumber; title 'Output 6.5: Double Multivariate Linear Model'; data dmlm; infile 'c:\mmm.dat'; input group y1 - y9; proc print data=dmlm; run; proc glm; title2 ' Double Multivariate Model Analysis'; class group; model y1 - y9 = group/ nouni; means group; manova h=group/ printh printe; run; /* Multivariate test of Parallelism */ proc glm; class group; model y1 - y9 = group / nouni; contrast 'Parallel' group 1 -1; manova m=(.7071 0 -.7071 0 0 0 0 0 0, 0 0 0 -.408 .816 -.408 0 0 0, 0 0 0 0 0 0 .7071 0 -.7071, -.408 .816 -.408 0 0 0 0 0 0, 0 0 0 .7071 0 -.7071 0 0 0, 0 0 0 0 0 0 -.408 .816 -.408) prefix = parl/ printh printe; run; /* Multivariate test of Conditions as vectors */ proc glm; class group; model y1 - y9 = group/ noint nouni; contrast 'Mult Cond' group 1 0, group 0 1; manova m=(1 -1 0 0 0 0 0 0 0, 0 1 -1 0 0 0 0 0 0, 0 0 0 1 -1 0 0 0 0, 0 0 0 0 1 -1 0 0 0, 0 0 0 0 0 0 1 -1 0, 0 0 0 0 0 0 0 1 -1) prefix = diff/ printh printe; run; /* Multivariate test of Conditions given Parallelism */ proc glm; class group; model y1- y9 = group/noint nouni; contrast 'Cond|Parl' group .5 .5; manova m=(.7071 0 -.7071 0 0 0 0 0 0, 0 0 0 -.408 .816 -.408 0 0 0, 0 0 0 0 0 0 .7071 0 -.7071, -.408 .816 -.408 0 0 0 0 0 0, 0 0 0 .7071 0 -.7071 0 0 0, 0 0 0 0 0 0 -.408 .816 -.408) prefix=cond/ printh printe; run; proc glm; class group; model y1 - y9 = group/noint nouni; contrast 'Group|Parl' group 1 -1; manova m=(.577 .577 .577 0 0 0 0 0 0, 0 0 0 .577 .577 .577 0 0 0, 0 0 0 0 0 0 .577 .577 .577) prefix=ovall/ printh printe; /* Multivariate Mixed Model Analysis */ data mix; infile 'c:\mixed.dat'; input group subj cond y1 y2 y3; proc print data=mix; run; proc glm; title2 'Multivaraiate Mixed Model Analysis'; class group subj cond; model y1 - y3 = group subj(group) cond cond*group; random subj(group); contrast 'Group' group 1 -1/ e=subj(group); manova h = cond group*cond/ printh printe; run; /* Test for Multivariate Spericity and calculation of Epsilon for MMM*/ proc iml; title2 'Multivariate Sphericity Test'; print 'Test of Multivariate Sphericity UsingChi-Square and Adjusted Chi- Square Statistics'; e={ 9.6944 7.3056 -6.7972 -4.4264 -0.6736 3.7255, 7.3056 8.8889 -4.4583 -3.1915 -3.2396 2.9268, -6.7972 -4.4583 18.6156 2.5772 0.8837 -10.1363, -4.4264 -3.1915 2.5772 5.3981 1.4259 -1.8546, -0.6736 -3.2396 0.8837 1.4259 18.3704 -.7769, 3.7255 2.9268 -10.1363 -1.8546 -0.7769 6.1274}; print e; n=18; p=3; t=3; k=2; u=6; q=u/p; nu_e=n-k; nu_h=1; e11=e[1:3,1:3];print e11; e22=e[4:6,4:6];print e22; dn=(e11+e22)/2; b=eigval(dn); print b; a=eigval(e); print a; b=log(b); a=log(a); chi_2=n#(q#sum(b)-sum(a)); df=p#(q-1)#(p#q+p+1)/2; pvalue=1-probchi(chi_2,df); print chi_2 df pvalue; c1=p/(12#q#nu_e#df); rho= 1-c1#(2#p##2#(q##4-1)+3#p#(q##3-1)-(q##2-1)); ro_chi_2=(rho#nu_e/n)#chi_2; print rho; c2=1/(2#rho##2); c3=((p#q-1)#p#q#(p#q+1)#(p#q+2))/(24#nu_e##2); c4=((p-1)#p#(p+1)#(p+2))/(24#q##2#nu_e##2); c5=df#(1-rho)##2/2; omega=c2#(c3-c4-c5); print omega; p1=1-probchi(ro_chi_2,df); p2=1-probchi(ro_chi_2,df+4); cpvalue=(1-omega)#p1+omega#p2; print ro_chi_2 cpvalue; s=e/nu_e; s11=s[1:3,1:3]; s12=s[1:3,4:6]; s21=s[4:6,1:3]; s22=s[4:6,4:6]; enum=trace((s11+s22)*(s11+s22))+(trace(s11+s22))##2; eden=q#( trace(s11)##2+trace(s11*s11)+trace(s12)##2+trace(s12*s12)+ trace(s21)##2+trace(s21*s21)+trace(s22)##2+trace(s22*s22)); epsilon=enum/eden; nu_h=nu_h#q; nu_e=nu_e#q; s0=min(nu_h,p); Mnu_h=nu_h#epsilon; Mnu_e=nu_e#epsilon; ms0=min(mnu_h,p); m0=(abs(mnu_h-p)-1)/2; n0=(mnu_e-p-1)/2; denom=s0##2#(2#m0+ms0+1); numer=2#(ms0#n0+1); df1=ms0#(2#m0+ms0+1); df2=2#(ms0#n0+1); print 'Epsison adjusted F-Statistics for cond and group X cond MMM tests'; print epsilon; f_cond = df2#13.75139851/(ms0#df1); f_gXc = df2#0.19070696/(ms0#df1); print f_cond f_gXc df1 df2; p_cond=1-probf(f_cond,df1,df2); p_gXc=1-probf(f_gXc,df1,df2); print 'Epsilon adjusted pvalues for MMM tests using T0**2 Criterion'; print p_cond p_gXc; quit; *************************************************** /* The following code appears on page 305, */ /* where it is used to produce output 7.5. */ *************************************************** /* Program 7_5.sas */ options ls=80 ps=60 nodate nonumber; title1 'Output 7.5: Restricted Nonorthogonal 3 Factor MANOVA Design'; data three; infile 'c:\7_5.dat'; input a b c y1 y2 u1 - u18; proc print data=three; run; /* Using PROC REG */ proc reg; title2 'Analysis using PROC REG'; model y1 y2 = u1 - u18/noint; restrict .3333*u1+.3333*u2+.3333*u3-.3333*u4- .3333*u5-.3333*u6-.3333*u7-.3333*u8- .3333*u9+.3333*u10+.3333*u11+.3333*u12=0, .3333*u7+.3333*u8+.3333*u9-.3333*u10- .3333*u11-.3333*u12-.3333*u13-.3333*u14- .3333*u15+.3333*u16+.3333*u17+.3333*u18=0, .3333*u1-.3333*u2-.3333*u4+.3333*u5+ .3333*u7-.3333*u8-.3333*u10+.3333*u11+ .3333*u13-.3333*u14-.3333*u16+.3333*u17=0, .3333*u2-.3333*u3-.3333*u5+.3333*u6+ .3333*u8-.3333*u9-.3333*u11+.3333*u12+ .3333*u14-.3333*u15-.3333*u17+.3333*u18, u1-u2-u4+u5-u7+u8+u10-u11=0, u2-u3-u5+u6-u8+u9+u11-u12=0, u7-u8-u10+u11-u13+u14+u16-u17=0, u8-u9-u11+u12-u14+u15+u17-u18=0; A: mtest u1+u2+u3+u4+u5+u6-u13-u14-u15-u16-u17-u18=0, u7+u8+u9+u10+u11+u12-u13-u14-u15-u16-u17-u18=0/print details; B: mtest u1-u3+u4-u6+u7-u9+u10-u12+u13-u15+u16-u18=0, u2-u3+u5-u6+u8-u9+u11-u12+u14-u15+u17-u18=0/print details; C: mtest u1+u2+u3-u4-u5-u6+u7+u8+u9-u10-u11-u12+u13+u14+u15-u16-u17-u18=0 /print details; AB: mtest u1-u2+u4-u5-u7+u8-u10+u11=0, u2-u3+u5-u6-u8+u9-u11+u12=0, u7-u8+u10-u11-u13+u14-u16+u17=0, u8-u9+u11-u12-u14+u15-u17+u18=0/print details; run; /* Using PROC GLM */ proc glm; title2 'Analysis using PROC GLM'; class a b c; model y1 y2 = a b c a*b/ss3; lsmeans a b c a*b; manova h=a b c a*b/printe printh; run; *************************************************** /* The following code appears on page 311, */ /* where it is used to produce output 7.6. */ *************************************************** /*Program 7_6.sas */ options ls=80 ps=60 nodate nonumber; title1 'Output 7.6: Restricted Intra-class Covariance Design'; data intra; infile 'c:\7_6.dat'; input a b y1 y2 u1-u6 z1-z6; proc print data=intra; run; /* Testing Equality of intra-class regression coefficients */ proc reg; title2 'Testing Equality of Inta-class Regression Coefficients'; model y1 y2 = u1 - u6 z1 - z6/noint; Parallel: mtest z1-z6=0, z2-z6=0, z3-z6=0, z4-z6=0, z5-z6=0; Coin: mtest z1-z6=0,z2-z6=0,z3-z6=0,z4-z6=0,z5-z6=0, u1-u6=0,u2-u6=0,u3-u6=0,u4-z6=0,u5-u6=0; run; /* MANCOVA model with unweighted and weighted tests*/ proc reg; title2 'MANCOVA Analyis using PROC REG'; model y1 y2= u1 - u6 z1 - z6/noint; restrict z1-z6=0, z2-z6=0, z3-z6=0, z4-z6=0, z5-z6=0; A: mtest u1+u2+u3-u4-u5-u6=0; B: mtest u1-u3+u4-u6=0, u2-u3+u5-u6=0; AB: mtest u1-u3-u4+u6=0, u2-u3-u5+u6=0; Reg: mtest z6=0; Awt: mtest .4167*u1+.3333*u2+.25*u3-.3*u4-.4*u5-.3*u6=0; Bwt: mtest .625*u1-.5*u3+.375*u4-.5*u6=0, .500*u2-.5*u3+.500*u5-.5*u6=0; run; /* MANOVA model with unweighted tests */ proc reg; title2 ' Restricted Unweighted MANOVA tests using PROC REG'; model y1 y2= u1 - u6 z1 -z6/noint; restrict z1=z2=z3=z4=z5=z6=0; A: mtest u1+u2+u3-u4-u5-u6=0; B: mtest u1-u3+u4-u6=0, u2-u3+u5-u6=0; AB: mtest u1-u3-u4+u6=0, u2-u3-u5+u6=0; run; data mancova; infile 'c:\7_6a.dat'; input a b y1 y2 z; proc print data=mancova; run; /* Unweighted MANOVA Analysis using GLM */ proc glm; title2 'Unweighted MANOVA tests using PROC GLM'; class a b; model y1 y2=a b a*b/ss3; lsmeans a b a*b; manova h= a b a*b/htype=3 etype=3; run; /* Unweighted MANCOVA Analysis using GLM */ proc glm; title2 'Unweighted MANCOVA tests using PROC GLM'; class a b; model y1 y2=a b a*b z/ss3; manova h=a b a*b z/htype=3 etype=3; run; *************************************************** /* The following code appears on page 320, */ /* where it is used to produce output 7.7. */ *************************************************** /* Program 7_7.sas */ /* Growth Curve Analysis (Grizzle and Allen, 1969) */ options ls=80 ps=60 nodate nonumber; title1 'Output 7.7: Growth Curve Analysis - Grizzle and Allen Data'; data grizzle; infile'c:\7_7.dat'; input group y1 y2 y3 y4 y5 y6 y7; proc print data=grizzle; run; proc summary data=grizzle; class group; var y1 -y7; output out=new mean=mr1-mr7; proc print data=new; run; data plot; set new; array mr(7) mr1-mr7; do time = 1 to 7; response = mr(time); output; end; drop mr1-mr7; run; proc plot; title2 'Growth Curve Plots of Group Means'; plot response*time=group; run; proc glm data=grizzle; title2 'Transformed Data Polynomials'; class group; model y1 - y7 = group/nouni; repeated time polynomial/summary nom nou; run; proc iml; use grizzle; read all var {y1 y2 y3 y4 y5 y6 y7} into y; read all var {Group} into gr; /* Create Orthogonal Polynomials of degree q-1=6 */ x={1 2 3 4 5 6 7}; P_prime=orpol(x,6); Y0=Y*P_prime; print 'P prime matrix',P_prime; t=y0||gr; /* Create New Transformed data set */ varnames={yt1 yt2 yt3 yt4 yt5 yt6 yt7 group}; create trans from t (|colname=varnames|); append from t; quit; run; proc print data=trans; run; /* Test of model fit */ proc glm data=trans; title2 'Test of Cubic Fit'; model yt5 - yt7=/nouni; manova h=intercept; run; proc glm data=trans; title2 'Test of Quadratic Fit'; model yt4 - yt7=/nouni; manova h=intercept; run; /* Using a Rao-Khatri MANCOVA model */ /* Test for group differences (coincidence) */ proc glm data=trans; title2 'Test of Group Differences'; class group; model yt1 - yt4 = group yt5 - yt7/nouni; manova h=group/printh printe; run; /*Test for Parallelism of Profiles */ proc glm data=Trans; title2 ' Test of Parallelism'; class group; model yt1 - yt4 = group yt5 - yt7/noint nouni; contrast 'parallel' group 1 -1 0 0, group 1 0 -1 0, group 1 0 0 -1; manova m=(0 1 0 0, 0 0 1 0, 0 0 0 1) prefix = parll/printe printh; run; /* Estimate of matrix B and contrasts */ proc glm data=trans; title2 ' Estimate of Matrix B'; class group; model yt1 - yt4 = group yt5 - yt7; estimate 'beta1' intercept 1 group 1 0 0 0; estimate 'beta2' intercept 1 group 0 1 0 0; estimate 'beta3' intercept 1 group 0 0 1 0; estimate 'beta4' intercept 1 group 0 0 0 1; estimate '1vs2' group 1 -1 0 0; estimate '1vs3' group 1 0 -1 0; estimate '1vs4' group 1 0 0 -1; run; *************************************************** /* The following code appears on page 332, */ /* where it is used to produce output 7.8. */ *************************************************** /* Program 7_8.sas */ /* Growth Curve Analysis (Timm, 1980a) */ options ls=80 ps=60 nodate nonumber; title1 'Output 7.8: Multiple Response Growth Curve Analysis- Zullo data'; data zullo; infile 'c:\mmm.dat'; input group y1 - y9; proc print data=zullo; run; data test; set zullo; array z{9} y1 - y9; k=1; do rvar = 1,2,3; do time = 1,2,3; grmeans = z{k}; k=k+1; output; end; end; keep rvar group time grmeans; run; proc summary nway; class group rvar time; var grmeans; output out=new mean=grmeans; proc print data = new; run; proc plot; title2 'Plot of Means'; plot grmeans*time=group; run; proc iml; use zullo; read all var {y1 y2 y3 y4 y5 y6 y7 y8 y9} into y; read all var {Group} into gr; /* Create Orthogonal Polynomials of degree q-1=2 */ x={1 2 3}; P=orpol(x,2); print P; P_prime=block(P,P,P); Y0=Y*P_prime; t=y0||gr; /* Create New Transformed data set */ varnames={yt1 yt2 yt3 yt4 yt5 yt6 yt7 yt8 yt9 group}; create trans from t (|colname=varnames|); append from t; quit; proc print data=trans; run; /* Test model adequacy or fit */ proc glm data=trans; title2 'Test of Fit - Linear'; model yt3 yt6 yt9=; manova h=intercept; run; proc glm data=trans; title2 'Test of Fit - Quadratic(1)/Linear(2)'; model yt6 yt9=/nouni; manova h=intercept; run; /* Test of significance of Beta weights */ proc glm data=trans; title2 'Test of Quadratic Beta Weights'; class group; model yt1 - yt9 = group/noint nouni; contrast 'order' group 1 0, group 0 1; manova m=(0 0 1 0 0 0 0 0 0, 0 0 0 0 0 1 0 0 0, 0 0 0 0 0 0 0 0 1); run; /* Using a Rao-Khatri MANCOVA model */ /* Test for group differences (coincidence) p=q and q