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

data set1;
  infile 'c:\courses\st511\sas\dogs.clinics.dat';
  input drug disease y clinic;
  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;

/*  UNIX users should delete the next line and replace
    it with the lines in the next comment  */

goptions device=WIN target=WINPRTC rotate=landscape;

/*
 filename graffile pipe 'lpr -Dpostscript';
 goptions gsfmode=replace gsfname=graffile
    targetdevice=ps300 rotate=landscape;
*/

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

axis2 label=(f=swiss h=2.2)
      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=' ';
       label variety = ' ';
 run;


/*  Perform analysis of variance where the disease factor 
    is entered into the model before the drug factor.  
    Compute an ANOVA table and expectations of mean squares  */
 

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


 /*  Repeat the analysis with the MIXED procedure  */
     


  proc mixed data=set1;
    class clinic drug disease;
    model y = disease drug disease*drug / solution
                e1 e3 outpm=set1 outp=set2;;
    random clinic / G solution;
    ods output solutionr=setr;
	ods output covparms=setv;
    lsmeans disease*drug / slice=disease;
    lsmeans disease*drug / pdiff tdiff;
    lsmeans disease*drug / pdiff adjust=Tukey; 
    run;

	proc print data=setr; run;

	proc print data=setv; run;


proc rank data=setr normal=blom out=setr ; 
  var clinic; ranks q;  run;

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

axis2 label=(f=swiss h=2.0 r=90 a=270)
      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=setr;
   plot clinic*q / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Normal Plot: Clinic Effects";
       label q='Standard Normal Quantiles ';
       label clinic = 'Effects ';
 run;

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

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

axis2 label=(f=swiss h=2.0 r=90 a=270)
      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=set1;
   plot resid*q / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Normal Plot: Residuals";
       label q='Standard Normal Quantiles ';
       label resid = 'Residuals';
 run;

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