data Sheffield;
  input fat method lab rep ;
cards;
5.19  1  1  1
5.09  1  1  2
4.09  1  2  1
3.99  1  2  2
3.75  1  2  3
4.04  1  2  4
4.06  1  2  5
4.62  1  3  1
4.32  1  3  2
4.35  1  3  3
4.59  1  3  4
3.71  1  4  1
3.86  1  4  2
3.79  1  4  3
3.63  1  4  4
3.26  2  1  1
3.48  2  1  2
3.24  2  1  3
3.41  2  1  4
3.35  2  1  5
3.04  2  1  6
3.02  2  2  1
3.32  2  2  2
2.83  2  2  3
2.96  2  2  4
3.23  2  2  5
3.07  2  2  6
3.08  2  3  1
2.95  2  3  2
2.98  2  3  3
2.74  2  3  4
3.07  2  3  5
2.70  2  3  6
2.98  2  4  1
2.89  2  4  2
2.75  2  4  3
3.04  2  4  4
2.88  2  4  5
3.20  2  4  6
;
run;

/* Some plots to visually compare fat content from the two methods: */

goptions reset=all;
symbol v=dot h=.8 c=blue;
axis1 order=(0 to 5 by 1) label=(a=90 'Laboratory');
axis2 order=(2.5 to 5.5 by 1);
proc gplot data=sheffield;
  by method;
  plot lab*fat/ vaxis=axis1 haxis=axis2;
run;
quit;

/* Plotting the means for each method-laboratory combination */

proc sql;
  create table plot as
  select method, lab, mean(fat) as mfat
  from sheffield
  where method=1
  group by lab;
quit;
proc sql;
  create table plot1 as
  select method, lab, mean(fat) as mfat1
  from sheffield
  where method=2
  group by lab;
quit;
data combo;
  set plot plot1;
run;
symbol1 c=blue v=dot h=.8 i=join;
symbol2 c=red v=square h=.8 i=join;
axis1 label=(a=90 'Mean');
legend1 label=none value=(height=1 font=swiss 'Method 1' 'Method 2' ) 
        position=(bottom right inside) mode=share cborder=black;
proc gplot data=combo;
  plot (mfat mfat1)*lab / vaxis=axis1 overlay legend=legend1;
run;
quit; 

/* PROC MIXED gives ML estimates of the model parameters.                */
/* Note the constraints used are a bit different from those in the book. */

PROC MIXED DATA=Sheffield method=ml covtest;
CLASS method lab;
MODEL fat = method / s;
LSMEANS method/adj=tukey pdiff;
RANDOM lab method*lab;
run;


/* Omnibus test of random effects is possible with PROC GLIMMIX.                             */
/* The default estimation method in PROC GLIMMIX (and PROC MIXED, actually) is REML, not ML. */
PROC GLIMMIX data=Sheffield;
CLASS method lab;
MODEL fat = method / s;
LSMEANS method/ pdiff adjust=tukey;
RANDOM lab method*lab;
COVTEST zerog; /* Tests whether all random effect variances are zero */ 
run;