baking <- read.table("baking.sdat",header=T)
baking$fat <- factor(baking$fat)
baking$surf <- factor(baking$surf)
baking$flour <- factor(baking$flour)

baking.fit <- aov(spvol~flour+fat*surf,baking)
title chapter 12 - fat*surfactant;
data; set;
   drop flour1-flour4;
   flour=1; spvol=flour1; output;
   flour=2; spvol=flour2; output;
   flour=3; spvol=flour3; output;
   flour=4; spvol=flour4; output;
proc glm;
   classes flour fat surf;
   model spvol=flour fat|surf;
   lsmeans fat|surf / stderr pdiff;

effect.plot <- function(fit,data=eval(fit$call$data),
	factors = names(predictors[predictors]))
{
        predictors <- unlist(lapply(data, is.factor))[all.names(formula(fit), 
                unique = T, functions = F)[-1]]
        int.factors <- interaction(get.list(data, factors), drop = T)
        design <- !duplicated(int.factors)
        newdata <- data[design,  ]      # average covariates to mean value
	dimnames(design) <- list(as.character(unique(int.factors)),dimnames(newdata)[[2]])
	x <- predict(fit,newdata,type="terms")#[,factors]
#	x <- predict(fit,newdata)
	x
}