split1 <- read.table("http://stat.ufl.edu/~winner/sta6207/kuehl_prob14_1.txt",
   header=F,col.names=c("hybrid","block","density","hd.wt"))
attach(split1); names(split1)

hybrid <- factor(hybrid)
block <- factor(block)
density <- factor(density,levels=c(10,15,25,40), ordered=TRUE)

### Create contrasts that match Kuehl's on p. 488

contrasts(density) <- cbind(c(-0.546,-0.327,0.109,0.764),
   c(0.513,-0.171,-0.741,0.399),c(-0.435,0.783,-0.435,0.087))

contrasts(density)

### Fit the model with the correct error term for hybrid (densityxblock)
### Obtain polynomial contrasts for density using correct error term 
prob14.mod1 <-  aov(hd.wt ~ hybrid*density + Error(block/density))
summary(prob14.mod1)
summary(prob14.mod1,split=list(density=1:3))


### Fit the model with the incorrect error term for hybrid (Error2)
### F-test for density polynomials use incorrect error term
### Obtain polynomial contrasts for hybridxdensity using its correct error term 
prob14.mod2 <- aov(hd.wt ~ block + density + block:density + hybrid + density:hybrid)
summary(prob14.mod2)
summary(prob14.mod2,split=list(density=1:3))