/*  SAS code for analyzing data from the dog blood
    pressure experiment on assignment 7 in 2002.  */

data set1;
  infile 'c:\courses\st511\sas\dogs.dat';
  input drug disease y;
  run;


/*  Print the data  */

proc print data=set1;
run;



/* Compute sample means for all factor combinations
   with data.  Make a profile plot. */


proc sort data=set1; by drug disease;
proc means data=set1 noprint; by drug disease;
  var Y;
  output out=means mean=my;
  run;



goptions  cback=white colors=black device=WIN 
  target=WINPRTC rotate=landscape;

/* UNIX users can use the following options:

goptions cback=white colors=(black)
    targetdevice=ps300 rotate=landscape;
*/

axis1 label=(f=swiss h=2.0)
      value=(f=swiss h=2.0)  w=3.0  length= 6.0 in;

axis2 label=(f=swiss h=1.8 a=90 r=0)
      value=(f=swiss h=2.0) w= 3.0  length = 4.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=1 i=join ;
       SYMBOL2 V=diamond H=2.0 w=3 l=3 i=join ;
       SYMBOL3 V=square H=2.0 w=3 l=9 i=join ;

 proc gplot data=means;
   plot my*drug=disease / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Sample means";
       label my='Mean Response ';
       label drug = 'Drug ';
 run;


/*  Perform analysis of variance where the disease 
    factor is entered into the model before the drug 
    factor.  Use the LSMEANS statement to compare means 
    for different combinations of drug and disease factors  */

  proc glm data=set1;
    class drug disease;
    model y =  disease drug disease*drug / solution
               ss1 ss2 ss3 e e1 e2 e3 p;
    lsmeans disease*drug / pdiff tdiff stderr;
    run;

/*  Perform analysis of variance where the drug factor 
    is entered into the model before the disease
    factor.  Use the LSMEANS statement to compare means 
    for different combinations of drug and disease factors  */

  proc glm data=set1;
    class drug disease;
    model y = drug disease drug*disease / solution
               ss1 ss2 ss3 e e1 e2 e3;
    output out=set2 predicted=yhat residual=resid;
    run;

    
 /* Compute a normal probability plot for
	the resdiuals  */

 proc rank data=set2 normal=blom out=set2;
   var resid; ranks q;  run;


axis1 label=(f=swiss h=2.0)
      value=(f=swiss h=2.0)  w=3.0  length= 6.0 in;

axis2 label=(f=swiss h=1.8 a=90 r=0)
      value=(f=swiss h=2.0) w= 3.0  length = 4.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=1 i=none ;
     

 proc gplot data=set2;
   plot resid*q / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Normal Probability Plot";
       label resid='Ordered Residuals ';
       label q = 'Standard Normal Quantiles ';
 run;


 /* Check variance homogeneity  */

 
axis1 label=(f=swiss h=2.0)
      value=(f=swiss h=2.0)  w=3.0  length= 6.0 in;

axis2 label=(f=swiss h=1.8 a=90 r=0)
      value=(f=swiss h=2.0) w= 3.0  length = 4.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=1 i=none ;
    

 proc gplot data=set2;
   plot resid*yhat / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Residual Plot";
       label resid='Residuals ';
       label q = 'Estimated Means ';
 run;