/*
Demonstration SAS Program 
Analysis of Covariance 
Taken from Snedecor and Cochran 6ed. pg. 422.
SAS code from Dr. Dickey's ST512 website 
*/

options linesize=75;
data leprosy; do drug = 'A', 'B', 'C';
  input X Y @@;
d1=(drug = 'A'); d2= (drug='B');
xd1 = x*d1; xd2=x*d2; int =1;
label X = "Before Treatment";
label Y = "After treatment";
 output; end;
cards;
11  6     6  0      16 13
 8  0     6  2      13 10
 5  2     7  3      11 18
14  8     8  1       9  5
19 11    18 18      21 23
 6  4     8  4      16 12
10 13    19 14      12  5
 6  1     8  9      12 16
11  8     5  1       7  1
 3  0    15  9      12 20
 ;
/* plot bacteria counts after (vertical axis) against before (horizontal axis)
use a different color or symbol for each drug */
proc gplot ; plot y*x=drug;
run;


/* Consider a multiple linear regression: 
Y = a linear function of X
(Y=after, X=before)
with different slopes and different intercepts for each drug.
Use PROC PRINT to see the design matrix.
*/

title "Leprosy Bacilli Abundance Scores";
title2 "Before and after 48 weeks of treatment";
title3 "At 6 body sites for each of 10 patients per drug";
footnote1 "C is a control, A and B are antibiotics";
footnote2 "source: Snedecor and Cochran 6ed. pg. 422";
proc print; var y  int d1 d2 x xd1 xd2; Title4 'Y and X matrix';
options pagesize=55;
proc reg; model Y = X d1 d2 Xd1 Xd2/ss1; Title 'REG with interaction';
footnote ' ';  * -- turn off all footnotes ;

/* Note R(Xd1|X,d1,d2) + R(Xd2|X,d1,d2,Xd1)  is needed for test of equal slopes 
   Also, see PROC GLM note below.
*/

proc reg; model Y = X d1 d2;  Title 'Reg, no interaction';

proc glm; 
   class drug; 
   model y = X drug drug*x;
   *model y = X|drug;   * shorthand for the same model;
   title 'GLM test for interaction';

/*Note that SS(drug*x) could be computed from Type I SS in PROC REG interaction model above */

/* Now consider a model with equal slopes, but different intercepts */
proc glm; class drug; model y = X drug;
   means drug;  lsmeans drug/pdiff out=out1;
   output out=out2 predicted = yhat r=resid;
    Title1 'GLM no interaction';
    Title2 '(same slopes, different intercepts)';
proc gplot data=out2;  plot resid*yhat; run;
proc gplot data=out2; plot yhat*x=drug/href=10.7;
           title  'Illustrate adjusted treatment means';
run;

proc plot data=out1; plot lsmean*drug/vpos=30;  
*vpos and hpos specify # positions available on each axis ;
           title  'Adjusted treatment means';
run;
proc gplot data=out1; plot lsmean*drug;
           title  'Adjusted treatment means';
run;