options nocenter ls=80 ps=10000 FORMCHAR='|----|+|---+=|-/\<>*'; ****************************************************************************************** ; Title "demo 1"; data demo1; input group time1 time2 time3 ; cards; 1 10 10 10 1 10 10 10 1 10 10 10 1 10 10 10 2 15 15 15 2 15 15 15 2 16 15 15 2 15 15 15 run; proc glm data=demo1; class group; model time1 time2 time3 = group; repeated time 3 ; lsmeans group; run; quit; ****************************************************************************************** ; Title "demo 2"; data demo2; input group time1 time2 time3 ; cards; 1 14 19 29 1 15 25 26 1 16 16 31 1 12 24 32 2 10 21 24 2 17 26 35 2 19 22 32 2 15 23 34 run; proc glm data=demo2; class group; model time1 time2 time3 = group; repeated time 3 ; lsmeans group; run; quit; ****************************************************************************************** ; Title "demo 3"; data demo3; input group time1 time2 time3 ; cards; 1 35 25 16 1 32 23 12 1 36 22 14 1 34 21 13 2 57 43 22 2 54 46 26 2 55 46 23 2 60 47 25 run; proc glm data=demo3; class group; model time1 time2 time3 = group; repeated time 3 ; lsmeans group; run; quit; ****************************************************************************************** ; Title "demo 4"; data demo4; input group time1 time2 time3 ; cards; 1 35 25 12 1 34 22 13 1 36 21 18 1 35 23 15 2 31 43 57 2 35 46 58 2 37 48 51 2 32 45 53 run; proc glm data=demo4; class group; model time1 time2 time3 = group; repeated time 3 ; lsmeans group; run; quit; ***************************************************************************************** ; * Now show getting plots in SAS ; /* We use the out option in the lsmeans statement to get the data set with all the means */ proc glm data=demo4; class group; model time1 time2 time3 = group; repeated time 3 ; lsmeans group / out=means; run; quit; proc print data=means; run; /* We reset all the graph options, create different symbols for each level of the variable group and create a nice label for the vertical axis. */ goptions reset=all; symbol1 c=blue v=star h=.8 i=j; symbol2 c=red v=dot h=.8 i=j; axis1 label=(a=90 'Means'); axis2 value=('1' '2' '3') label=('Time'); proc gplot data=means; plot lsmean*_name_ = group/ vaxis=axis1 haxis=axis2; run; quit; ********************************************************************************* ; * Excersise example ; data exercise; input id exertype diet time1 time2 time3; cards; 1 1 1 85 85 88 2 1 1 90 92 93 3 1 1 97 97 94 4 1 1 80 82 83 5 1 1 91 92 91 6 1 2 83 83 84 7 1 2 87 88 90 8 1 2 92 94 95 9 1 2 97 99 96 10 1 2 100 97 100 11 2 1 86 86 84 12 2 1 93 103 104 13 2 1 90 92 93 14 2 1 95 96 100 15 2 1 89 96 95 16 2 2 84 86 89 17 2 2 103 109 90 18 2 2 92 96 101 19 2 2 97 98 100 20 2 2 102 104 103 21 3 1 93 98 110 22 3 1 98 104 112 23 3 1 98 105 99 24 3 1 87 132 120 25 3 1 94 110 116 26 3 2 95 126 143 27 3 2 100 126 140 28 3 2 103 124 140 29 3 2 94 135 130 30 3 2 99 111 150 ; run; ****************************************************************************************** ; title "Exercise model 1: time and diet" ; proc glm data=exercise; class diet; model time1 time2 time3 = diet ; repeated time 3 / printe; run; quit; ****************************************************************************************** ; title "Exercise model 2: time and exertype" ; proc glm data=exercise; class exertype; model time1 time2 time3 = exertype; repeated time 3 ; run; quit; * Showing the graph below; proc glm data=exercise; class exertype; model time1 time2 time3 = exertype; repeated time 3 ; lsmeans exertype / out=means; run; quit; proc print data=means; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j; symbol2 c=red v=dot h=.8 i=j; symbol3 c=green v=square h=.8 i=j; axis1 order=(60 to 150 by 30) label=(a=90 'Means'); axis2 value=('1' '2' '3') label=('Time'); proc gplot data=means; plot lsmean*_name_= exertype / vaxis=axis1 haxis=axis2; run; quit; ****************************************************************************************** ; * Look at covariances ; proc corr data=exercise cov; var time1 time2 time3; run; ****************************************************************************************** ; * Reshaping exercise data; proc transpose data=exercise out=long; by id diet exertype; run; data long; set long (rename=(col1=pulse) ); time = substr(_NAME_, 5, 1 )+0; drop _name_; run; proc print data=long (obs=20); var id diet exertype time pulse; run; ****************************************************************************************** ; Title "model 2 with cs structure"; proc mixed data=long; class exertype time; model pulse = exertype time exertype*time; repeated time / subject=id type=cs; run; ****************************************************************************************** ; Title "model 2 with unstructure"; proc mixed data=long; class exertype time; model pulse = exertype time exertype*time; repeated time / subject=id type=un; run; ****************************************************************************************** ; Title "model 2 with ar(1) structure"; proc mixed data=long; class exertype time; model pulse = exertype time exertype*time; repeated time / subject=id type=ar(1); run; ****************************************************************************************** ; Title "model 2 with arh(1) structure"; proc mixed data=long; class exertype time; model pulse = exertype time exertype*time; repeated time / subject=id type=arh(1) ; run; title " "; ****************************************************************************************** ; * How to get a graph using only proc mixed. Note: proc mixed does not have an out option ; * in the lsmeans statement so instead we use ODS to create the dataset of the means. ; ods output LSMeans=means1; proc mixed data=long; class exertype time; model pulse = exertype time exertype*time; repeated time / subject=id type=ar(1); lsmeans time*exertype; run; proc print data=means1; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j; symbol2 c=red v=dot h=.8 i=j; symbol3 c=green v=square h=.8 i=j; axis1 order=(60 to 150 by 30) label=(a=90 'Means'); axis2 value=('1' '2' '3') label=('Time'); proc gplot data=means1; plot estimate*time=exertype / vaxis=axis1 haxis=axis2; run; quit; ****************************************************************************************** ; * model 3: time, diet and exertype ; proc glm data=exercise; class diet exertype; model time1 time2 time3 = diet exertype diet*exertype; repeated time 3 ; run; quit; ****************************************************************************************** ; * show graph ; proc glm data=exercise; class diet exertype; model time1 time2 time3 = diet|exertype; repeated time 3 ; lsmeans diet*exertype / out=means; run; quit; proc print data=means; run; proc sort data=means out=sortdiet; by diet; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j; symbol2 c=red v=dot h=.8 i=j; symbol3 c=green v=square h=.8 i=j; axis1 order=(60 to 150 by 30) label=(a=90 'Means'); axis2 value=('1' '2' '3') label=('Time'); proc gplot data=sortdiet; by diet; plot lsmean*_name_=exertype / vaxis=axis1 haxis=axis2; run; quit; ****************************************************************************************** ; * model 3: time, diet and exertype using proc mixed and arh(1) var-cov structure ; proc mixed data=long; class exertype diet time; model pulse = exertype|diet|time; repeated time / subject=id type=arh(1) ; run; ****************************************************************************************** ; * graph of exertype over time by diet using proc mixed and arh(1) var-cov structure ; ods output LSMeans = means; proc mixed data=long; class exertype diet time; model pulse = exertype|diet|time; repeated time / subject=id type=arh(1) ; lsmeans time*diet*exertype; run; proc print data=means; run; proc sort data=means; by diet; run; goptions reset=all; symbol1 c=black v=dot i=j; symbol2 c=blue v=circle i=j; symbol3 c=red v=square i=j; axis1 order=(60 to 150 by 30) label=(a=90 'Means'); proc gplot data=means; by diet; format estimate 8.; plot estimate*time=exertype / vaxis=axis1; run; quit; ****************************************************************************************** ; *Contrasts and contrast interactions in proc mixed.; proc mixed data=long; class diet exertype time; model pulse = exertype|diet|time; repeated time / subject=id type=arh(1) ; estimate 'exer 12 v 3' exertype -.5 -.5 1; /* across time and across diet groups */ estimate 'exer 1 v 2' exertype -1 1 0; /* across time and across diet groups */ estimate 'diet' diet -1 1; /* across time and across exercise types */ estimate 'diet 1v2 & exertype 12v3' diet*exertype -.5 -.5 1 .5 .5 -1; /* across time only */ estimate 'runners only, diet 1 v 2' diet 1 -1 diet*exertype 0 0 1 0 0 -1; /* across time only */ lsmeans diet*exertype / slice=diet; /*testing for differences among exertype for each level of diet */ lsmeans exertype*time / slice=time; /*differences in exertype for each time point*/ lsmeans exertype*diet*time / slice=time*diet; run; quit; ****************************************************************************************** ; * Unequally Spaced Time Points ; * Modeling Time as a Linear Predictor of Pulse ; data study2; input id timec exertype diet pulse time; cards; 1 1 1 1 90 0 1 2 1 1 92 228 1 3 1 1 93 296 1 4 1 1 93 639 2 1 1 1 90 0 2 2 1 1 92 56 2 3 1 1 93 434 2 4 1 1 93 538 3 1 1 1 97 0 3 2 1 1 97 150 3 3 1 1 94 295 3 4 1 1 94 541 4 1 1 1 80 0 4 2 1 1 82 121 4 3 1 1 83 256 4 4 1 1 83 575 5 1 1 1 91 0 5 2 1 1 92 161 5 3 1 1 91 252 5 4 1 1 91 526 6 1 1 2 83 0 6 2 1 2 83 73 6 3 1 2 84 320 6 4 1 2 84 570 7 1 1 2 87 0 7 2 1 2 88 40 7 3 1 2 90 325 7 4 1 2 90 730 8 1 1 2 92 0 8 2 1 2 94 205 8 3 1 2 95 276 8 4 1 2 95 761 9 1 1 2 97 0 9 2 1 2 99 57 9 3 1 2 96 244 9 4 1 2 96 695 10 1 1 2 100 0 10 2 1 2 97 143 10 3 1 2 100 296 10 4 1 2 100 722 11 1 2 1 86 0 11 2 2 1 86 83 11 3 2 1 84 262 11 4 2 1 84 566 12 1 2 1 93 0 12 2 2 1 103 116 12 3 2 1 104 357 12 4 2 1 104 479 13 1 2 1 90 0 13 2 2 1 92 191 13 3 2 1 93 280 13 4 2 1 93 709 14 1 2 1 95 0 14 2 2 1 96 112 14 3 2 1 100 219 14 4 2 1 100 367 15 1 2 1 89 0 15 2 2 1 96 96 15 3 2 1 95 339 15 4 2 1 95 639 16 1 2 2 84 0 16 2 2 2 86 92 16 3 2 2 89 351 16 4 2 2 89 508 17 1 2 2 103 0 17 2 2 2 109 196 17 3 2 2 114 213 17 4 2 2 120 634 18 1 2 2 92 0 18 2 2 2 96 117 18 3 2 2 101 227 18 4 2 2 101 614 19 1 2 2 97 0 19 2 2 2 98 70 19 3 2 2 100 295 19 4 2 2 100 515 20 1 2 2 102 0 20 2 2 2 104 165 20 3 2 2 103 302 20 4 2 2 103 792 21 1 3 1 93 0 21 2 3 1 98 100 21 3 3 1 110 396 21 4 3 1 115 498 22 1 3 1 98 0 22 2 3 1 104 104 22 3 3 1 112 310 22 4 3 1 117 518 23 1 3 1 98 0 23 2 3 1 105 148 23 3 3 1 118 208 23 4 3 1 121 677 24 1 3 1 87 0 24 2 3 1 122 171 24 3 3 1 127 320 24 4 3 1 133 633 25 1 3 1 94 0 25 2 3 1 110 57 25 3 3 1 116 268 25 4 3 1 119 657 26 1 3 2 95 0 26 2 3 2 126 163 26 3 3 2 143 382 26 4 3 2 147 501 27 1 3 2 100 0 27 2 3 2 126 70 27 3 3 2 140 347 27 4 3 2 148 737 28 1 3 2 103 0 28 2 3 2 124 61 28 3 3 2 140 263 28 4 3 2 143 588 29 1 3 2 94 0 29 2 3 2 135 164 29 3 3 2 130 353 29 4 3 2 137 560 30 1 3 2 99 0 30 2 3 2 111 114 30 3 3 2 140 362 30 4 3 2 148 501 ; run; * 2. overall scatterplot of everyone ; proc sort data=study2; by id time; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j r=10; symbol2 c=red v=dot h=.8 i=j r=10; symbol3 c=green v=square h=.8 i=j r=10; axis1 order=(60 to 150 by 30) label=(a=90 'Pulse'); proc gplot data=study2; plot pulse*time=id / vaxis=axis1; run; * 3. Fit data with linear model ; proc mixed data=study2 covtest noclprint; class id exertype ; model pulse = time exertype time*exertype / solution outp=pred1r outpm = pred1f; random intercept time / subject = id; run; * The output file pred1f contains the predicted values based on the ; * fixed part of the model. We can illustrate what the predicted values ; * of pulse look like using this model below. ; goptions reset=all; symbol1 c=blue v=star h=.8 i=j; symbol2 c=red v=dot h=.8 i=j; symbol3 c=green v=square h=.8 i=j; axis1 order=(60 to 150 by 30) label=(a=90 'Predicted Pulse'); proc gplot data=pred1f; plot pred*time=exertype /vaxis=axis1; run; quit; * We can include the observed pulse as well and see that ; * this model is not fitting very well at all. ; * The green line is fitting curved data with a straight line. ; proc sort data=pred1f; by time; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j w=10; symbol2 c=red v=dot h=.8 i=j w=10; symbol3 c=green v=square h=.8 i=j w=10; symbol4 c=blue v=star h=.8 i=j r=10; symbol5 c=red v=dot h=.8 i=j r=10; symbol6 c=green v=square h=.8 i=j r=10; axis1 order=(60 to 150 by 30) label=(a=90 'Predicted and Observed Pulse'); proc gplot data=pred1f; plot pred*time=exertype / vaxis=axis1 ; plot2 pulse*time = id / vaxis=axis1 ;; run; quit; * Modeling Time as a Quadratic Predictor of Pulse * quadratic model ; proc mixed data=study2 covtest noclprint; class id exertype; model pulse = time exertype time*exertype time*time / solution outp=pred2r outpm=pred2f ; random intercept time / subject = id; run; * Below we see the predicted values from this model with the quadratic effect of time.; proc sort data=pred2f; by time; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j ; symbol2 c=red v=dot h=.8 i=j ; symbol3 c=green v=square h=.8 i=j ; axis1 order=(60 to 150 by 30) label=(a=90 'Predicted Pulse'); proc gplot data=pred2f; plot pred*time=exertype /vaxis=axis1 ; run; quit; * Again, we can plot the predicted values against the actual values of pulse.; * We see that this model fits better, but it appears that the predicted values ; * for the green group have too little curvature and the red and blue group have * too much curvature.; proc sort data=pred2f; by time; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j w=10; symbol2 c=red v=dot h=.8 i=j w=10; symbol3 c=green v=square h=.8 i=j w=10; symbol4 c=blue v=star h=.8 i=j r=10; symbol5 c=red v=dot h=.8 i=j r=10; symbol6 c=green v=square h=.8 i=j r=10; axis1 order=(60 to 150 by 30) label=(a=90 'Predicted and Observed Pulse'); proc gplot data=pred2f; plot pred*time=exertype / vaxis=axis1 ; plot2 pulse*time = id / vaxis=axis1 ;; run; quit; * Modeling Time as a Quadratic Predictor of Pulse, Interacting by Exertype ; * quadratic model , model 3 ; proc mixed data=study2 covtest noclprint; class id exertype; model pulse = time exertype time*exertype time*time time*time*exertype / solution outp=pred3r outpm=pred3f ; random intercept time / subject = id; run; * predicted vs. actual ; proc sort data=pred3f; by time; run; goptions reset=all; symbol1 c=blue v=star h=.8 i=j w=10; symbol2 c=red v=dot h=.8 i=j w=10; symbol3 c=green v=square h=.8 i=j w=10; symbol4 c=blue v=star h=.8 i=j r=10; symbol5 c=red v=dot h=.8 i=j r=10; symbol6 c=green v=square h=.8 i=j r=10; axis1 order=(60 to 150 by 30) label=(a=90 'Predicted and Observed Pulse'); proc gplot data=pred3f; plot pred*time=exertype / vaxis=axis1 ; plot2 pulse*time = id / vaxis=axis1 ;; run; quit;