/* Read the bone density data from a text file */


  data set1;
    infile 'c:\courses\st511\snew\bone2.density.txt';
    input id group bmd weight age caffine calcium
         growmilk teenmilk twenmilk;        
   run;


/* Examine the contents of the file */

proc contents data=set1; run;


/* Fit a model  */

proc glm data=set1;
  class group;
  model bmd = age group / solution ss1 ss3;
  output out=set1 residual=r predicted=p;
  lsmeans group /stderr pdiff;
  run;

/* Compute plotting positions for a 
   normal probability plot  */

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



/* Windows users can use the following */	
goptions device=WIN target=WINPRTC rotate=landscape;


axis1 label=(f=swiss h=1.8 a=90 r=0 "Predicted BMD (g/cm^2)")  
      value=(f=swiss h=1.8)  w=3.0  
      length= 6.5 in;

axis2 label=(f=swiss h=2.0 "Residuals")  
      value=(f=swiss h=1.8) w= 3.0  
      length = 5.0 in;

axis1 label=(f=swiss h=1.8 a=90 r=0 "q")
      order -3 to 3 by 1 
      value=(f=swiss h=1.8)  w=3.0  
      length= 6.5 in;

 SYMBOL1 V=CIRCLE H=1.7 w=3 l=1 i=none ;
 SYMBOL2 V=DIAMOND H=1.7 w=3 l=3 i=join ;
 SYMBOL3 V=square H=1.7 w=3 l=9 i=join ;
 

 PROC GPLOT DATA=set1; 
   PLOT r*p / vaxis=axis2 haxis=axis1;
   TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Prdeicted Values";
	   footnote ls=0.01in;
   RUN;

 PROC GPLOT DATA=set1; 
   PLOT r*q / vaxis=axis2 haxis=axis1;
   TITLE1 ls=0.01in H=2.0 F=swiss "Normal Probability Plot";
	   footnote ls=0.01in;
   RUN;