/* -------------------------------------------------------------------------- */ /* Lab 9: Multi-way (or factorial) MANOVA and MANCOVA */ /* -------------------------------------------------------------------------- */ /* Set the library and path where the dataset is */ %let path='E:\Teaching\Mulivariate S09\Lab 9\'; libname lab9 &path; /* Format modifications */ proc format; value gender 0="Male" 1="Female"; value diagnosis 1="Mild" 2="Moderate" 3="Severe"; run; /* Getting the individual group means */ proc sort data=lab9.psych; by gender diagnosis; proc means mean std stderr min max; var test1 test2 test3 desirability; by gender diagnosis; format gender gender. diagnosis diagnosis.; proc freq; table diagnosis*gender; format gender gender. diagnosis diagnosis.; proc corr; var test1 test2 test3 desirability age; run; /* Multivariate Normality Assessment */ /* There are three sets of groups in a factorial design. Group membership can be determined by using either of the original categorical variables or the combination of the two */ data f_mild; set lab9.psych; if gender=1 & diagnosis=1; data m_mild; set lab9.psych; if gender=0 & diagnosis=1; data f_moderate; set lab9.psych; if gender=1 & diagnosis=2; data m_moderate; set lab9.psych; if gender=0 & diagnosis=2; data f_severe; set lab9.psych; if gender=1 & diagnosis=3; data m_severe; set lab9.psych; if gender=0 & diagnosis=3; run; /* Let us investigate the leverage values for each of our factorial groups */ proc reg data=f_mild noprint; model age = test1 test2 test3 desirability / influence; output out=lev_f_mild H=Leverage; proc reg data=m_mild noprint; model age = test1 test2 test3 desirability / influence; output out=lev_m_mild H=Leverage; proc reg data=f_moderate noprint; model age = test1 test2 test3 desirability / influence; output out=lev_f_moderate H=Leverage; proc reg data=m_moderate noprint; model age = test1 test2 test3 desirability / influence; output out=lev_m_moderate H=Leverage; proc reg data=f_severe noprint; model age = test1 test2 test3 desirability / influence; output out=lev_f_severe H=Leverage; proc reg data=m_severe noprint; model age = test1 test2 test3 desirability / influence; output out=lev_m_severe H=Leverage; run; quit; /* Critical leverage values we should consider */ proc iml; Nfmil = 73; Nmmil = 67; Nfmod = 72; Nmmod = 23; Nfsev = 44; Nmsev = 21; crit1_fmil = (2*4)/Nfmil; crit1_mmil = (2*4)/Nmmil; crit1_fmod = (2*4)/Nfmod; crit1_mmod = (2*4)/Nmmod; crit1_fsev = (2*4)/Nfsev; crit1_msev = (2*4)/Nmsev; crit1 = (crit1_fmil||crit1_mmil||crit1_fmod||crit1_mmod||crit1_fsev||crit1_msev)`; crit2_fmil = (2*(4+1))/Nfmil; crit2_mmil = (2*(4+1))/Nmmil; crit2_fmod = (2*(4+1))/Nfmod; crit2_mmod = (2*(4+1))/Nmmod; crit2_fsev = (2*(4+1))/Nfsev; crit2_msev = (2*(4+1))/Nmsev; crit2 = (crit2_fmil||crit2_mmil||crit2_fmod||crit2_mmod||crit2_fsev||crit2_msev)`; crit3_fmil = ((2*(gaminv(.99,2)))/(Nfmil-1))+(1/Nfmil); crit3_mmil = ((2*(gaminv(.99,2)))/(Nmmil-1))+(1/Nmmil); crit3_fmod = ((2*(gaminv(.99,2)))/(Nfmod-1))+(1/Nfmod); crit3_mmod = ((2*(gaminv(.99,2)))/(Nmmod-1))+(1/Nmmod); crit3_fsev = ((2*(gaminv(.99,2)))/(Nfsev-1))+(1/Nfsev); crit3_msev = ((2*(gaminv(.99,2)))/(Nmsev-1))+(1/Nmsev); crit3 = (crit3_fmil||crit3_mmil||crit3_fmod||crit3_mmod||crit3_fsev||crit3_msev)`; Group = {'Female Mild','Male Mild','Female Moderate', 'Male Moderate', 'Female Severe', 'Male Severe'}; print 'Leverage cut-off values'; print Group crit1 crit2 crit3; quit; /* Comparing the leverage cut-off values against the observed leverages*/ proc means data=lev_f_mild max; var Leverage; title 'Female Mild'; proc means data=lev_m_mild max; var Leverage; title 'Male Mild'; proc means data=lev_f_moderate max; var Leverage; title 'Female Moderate'; proc means data=lev_m_moderate max; var Leverage; title 'Male Moderate'; proc means data=lev_f_severe max; var Leverage; title 'Female Severe'; proc means data=lev_m_severe max; var Leverage; title 'Male Severe'; run; /* RUN THE MULTNORM2 MACRO */ /* Let's assess multivariate normality */ %multnormplt (data=f_mild, var= test1 test2 test3 desirability, title='Female Mild'); %multnormplt (data=m_mild, var= test1 test2 test3 desirability, title='Male Mild'); %multnormplt (data=f_moderate, var= test1 test2 test3 desirability, title='Female Moderate'); %multnormplt (data=m_moderate, var= test1 test2 test3 desirability, title='Male Mild'); %multnormplt (data=f_severe, var= test1 test2 test3 desirability, title='Female Severe'); %multnormplt (data=m_severe, var= test1 test2 test3 desirability, title='Male Severe'); quit; /* Bartlett's Test */ data forbart; set lab9.psych; if gender=1 & diagnosis=1 then group=1; if gender=0 & diagnosis=1 then group=2; if gender=1 & diagnosis=2 then group=3; if gender=0 & diagnosis=2 then group=4; if gender=1 & diagnosis=3 then group=5; if gender=0 & diagnosis=3 then group=6; proc discrim data=forbart pool=test; class group; var test1 test2 test3 desirability; run; /* Running the MANOVA and saving out residuals */ proc glm data=lab9.psych; class gender diagnosis; model test1 test2 test3 desirability = gender diagnosis gender*diagnosis; lsmeans gender diagnosis gender*diagnosis / stderr cl pdiff adjust=Bon; manova h=gender diagnosis gender*diagnosis/ printe printh; output out=lab9.resids r=rt1 rt2 rt3 rdes; run; quit; /* Let us create the corresponding profile plot */ title "Profile Plot for Test Data"; data psychflat; set forbart; admin="test1"; score=test1; output; admin="test2"; score=test2; output; admin="test3"; score=test3; output; admin="desir"; score=desirability; output; keep group admin score; run; proc sort data=psychflat; by group admin; run; proc means data=psychflat noprint; by group admin; var score; output out=psychmeans mean=mean; run; legend1 label=(height=1 position=top justify=center 'Groups') value= ('F_Mild' 'M_Mild' 'F_Moderate' 'M_Moderate' 'F_Severe' 'M_Severe') across=3 down=2; proc gplot data=psychmeans; axis1 length=4.5 in; axis2 length=7.5 in; plot mean*admin=group / vaxis=axis1 haxis=axis2 legend=legend1; symbol1 v=K f=special w=2 h=2 l=join i=l color=black; symbol2 v=K f=special w=2 h=2 l=join i=l color=red; symbol3 v=K f=special w=2 h=2 l=join i=l color=green; symbol4 v=K f=special w=2 h=2 l=join i=l color=blue; symbol5 v=K f=special w=2 h=2 l=join i=l color=orange; symbol6 v=K f=special w=2 h=2 l=join i=l color=brown; run; quit; /* Investigating residuals */ proc univariate data = lab9.resids plot ; by gender diagnosis; var rt1 rt2 rt3 rdes; run; quit; /* Run the multnorm2 SAS macro */ data f_mild_r; set lab9.resids; if gender=1 & diagnosis=1; keep rt1 rt2 rt3 rdes; data m_mild_r; set lab9.resids; if gender=0 & diagnosis=1; keep rt1 rt2 rt3 rdes; data f_moderate_r; set lab9.resids; if gender=1 & diagnosis=2; keep rt1 rt2 rt3 rdes; data m_moderate_r; set lab9.resids; if gender=0 & diagnosis=2; keep rt1 rt2 rt3 rdes; data f_severe_r; set lab9.resids; if gender=1 & diagnosis=3; keep rt1 rt2 rt3 rdes; data m_severe_r; set lab9.resids; if gender=0 & diagnosis=3; keep rt1 rt2 rt3 rdes; run; %multnormplt (data=f_mild_r, var= rt1 rt2 rt3 rdes, title='Female mild residuals'); %multnormplt (data=m_mild_r, var= rt1 rt2 rt3 rdes, title='Male mild residuals'); %multnormplt (data=f_moderate_r, var= rt1 rt2 rt3 rdes, title='Female moderate residuals'); %multnormplt (data=m_moderate_r, var= rt1 rt2 rt3 rdes, title='Male moderate residuals'); %multnormplt (data=f_severe_r, var= rt1 rt2 rt3 rdes, title='Female severe residuals'); %multnormplt (data=m_severe_r, var= rt1 rt2 rt3 rdes, title='Male severe residuals'); quit;