brass <- read.table( "brassica.dat", header = T, na.strings = "." )
brass <- data.frame(brass[ brass$crop != 2, ])
brass$yldkga <- unfactor( brass$yldkga )
brass$trt <- factor( as.character( brass$trt ) )
brass$crop <- factor( brass$crop )
brass$pyrrat <- ordered( brass$pyrrat )
brass$block <- factor( brass$block )
brass$year <- ordered( brass$year %% 10 )

brass.full <- aov ( yldkga ~ crop * pyrrat * year +
   Error( year:block + crop:year:block ), brass )

brass.red <- aov ( yldkga ~ crop * year + pyrrat + crop:pyrrat +
   Error( year:block + crop:year:block ), brass, qr = T )

brass.proj <- proj ( brass.red )
brass.ref <- data.frame(Reference=c(brass.proj[[1]]))
for (i in names(brass.proj)[-1])
   brass.ref[[i]] <- apply( brass.proj[[i]], 1, sum )

# whole plot
brass.wp <- nested( brass.proj, brass, c("block","year","crop"), "yldkga",
   c(1,3), brass.ref )

brass.wpaov <- aov( yldkga ~ year*block + crop + crop:year, brass.wp,
   weight = brass.wp$weight)

# split plot
brass.sp <- nested( brass.proj, brass, c("block","year","crop","pyrrat"),
   "yldkga", c(1,4), brass.ref )

brass.spaov <- aov( yldkga ~ year*block*crop + pyrrat + crop:pyrrat, brass.sp )