/*  This is a program for analyzing the 
    penicillan data from Box, Hunter, 
    and Hunter.  It is posteded in the 
    file
              penclln.sas          */

/*  First enter the data  */


data set1;
  infile 'c:\courses\st511\sas\penclln.dat';
  input batch process $ yield;
  run;


/* Compute the ANOVA table, formulas for
   expectations of mean squares, process
   means and their standard errors */

proc glm data=set1;
  class batch process;
  model yield = batch process / e e3;
  random batch / q test;
  lsmeans process / stderr pdiff tdiff;
  output out=set2 r=resid p=yhat;
  run;

/* Compute a normal probability plot for 
   the residuals and the Shapiro-Wilk 
   test for normality  */

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

proc univariate data=set2 normal plot;
  var resid;
  run;
/* UNIX users can use the following 
   graphics options.  This code is uses
   the graphics options for WINDOWS  */
/*
goptions cback=white colors=(black) 
       targetdevice=ps300 rotate=portrait;
*/

goptions cback=white colors=(black)
     device=win target=winprtc rotate=portrait;

  axis1 label=(h=2.5 r=0 a=90 f=swiss 'Residuals')
        value=(f=swiss h=2.0) w=3.0 length=5.0 in;

  axis2 label=(h=2.3 f=swiss 'Standard Normal Quantiles')
        value=(f=swiss h=2.0) w=3.0 length=5.0 in;

  axis3 label=(h=2.3 f=swiss 'Production Process')
        value=(f=swiss h=2.0) w=3.0 length=5.0 in;

  symbol1 v=circle i=none h=2 w=3 c=black;

proc gplot data=set2;
  plot resid*q / vaxis=axis1 haxis=axis2;
  title h=3.0 ls=1.0in f=swiss 
        c=black 'Normal Probability Plot';
  footnote ls=0.6in '  ';
  run;

proc gplot data=set2;
  plot resid*process / vaxis=axis1 haxis=axis3;
  title h=3.0 ls=1.0in f=swiss 
        c=black 'Residual Plot';
  footnote ls=0.6in '   ';
  run;

/* Fit the same model using PROC MIXED.  Compute
   REML estimates of variance components.  Note
   that PROC MIXED provides appropriate standard
   errors for process means.  When block effects
   are random. PROC GLM does not provide correct
   standard errors for process means  */

proc mixed data=set1;
  class process batch;
  model yield = process / ddfm=satterth solution;
  random batch / solution cl alpha=.05;
  make 'solutionR' out=setr;
  make 'CovParms' out=setv;
  lsmeans process / pdiff tdiff;
run;

proc mixed data=set1;
  class process batch;
  model yield = process / ddfm=satterth solution;
  random batch / type=vc G solution cl alpha=.05;
  lsmeans process / pdiff tdiff;
run;


proc print data=setr; run;
proc print data=setv; run;