options ps=40 ls=80;

data a;
   infile 'preferb.dat' firstobs=2;
   input ID reps type $ shape $ y;

/* balanced analysis of variance (fast, but limited) */
proc anova;
   class reps type shape;
   model y = reps type shape type * shape;
   means type shape / lsd;
   means type * shape;

/* general linear model (slower, but quite versatile) */
proc glm;
   class reps type shape;
   model y = reps type shape type * shape;
   lsmeans type shape type * shape / stderr pdiff out=lsmean;
   output out=diag p=py r=ry;

/* diagnostic plots of residuals vs. predicted by type and shape */
proc plot data=diag;
   plot ry*py=type;
   plot ry*py=shape;

/* interaction plot using least squares means */
proc plot data=lsmean;
   plot lsmean*shape=type;