tukey <- read.table("tukey.dat",header=T)
tukey$temp <- ordered(tukey$temp)
tukey$humid <- ordered(tukey$humid)

# Fit Additive Model
tukey.fit <- aov(y~temp+humid,tukey)
print(tukey.fit)

# Fit Tukey Interaction Model
tukey$inter <- predict(tukey.fit)^2
# only examine the test for "inter" covariate
tukey$pred <- predict(aov(y~temp+humid+inter,tukey))

# set up 2x2 plots
tmpar <- par(mfrow=c(2,2))

# Tukey Interaction Plots
attach(tukey)
# Interaction Plot
tmp1 <- margin.plot(humid,temp,y,method="add",orient="both",
   ylim=range(y),type="b",
   ylab=c("","yield"),xlab=c("(a) humid by temp","(b) temp by humid"),
   trace.label=rep("",2))
# Margin Plots
tmp2 <- tukey.plot(humid,temp,y,orient="both",ylim=range(y),type="b",
   ylab=c("","expected yield"),xlab=c("(c) humid by temp","(d) temp by humid"),
   trace.label=rep("",2))
detach()

# Mandel Interaction Plots
attach(tukey)
# Margin Plots
tmp1 <- mandel.plot(humid,temp,y,orient="both",type="b",
   xlab=c("(a) humid by temp","(b) temp by humid"),
   ylab=c("","expected yield"),
   trace.label=c("",""))
tmp2 <- mandel.plot(temp,humid,y,orient="both",type="b",
   xlab=c("(c) temp by humid","(d) humid by temp"),
   ylab=c("","expected yield"),
   trace.label=c("",""))
detach()

# reset plot parameters
par(tmpar)