# Full set of cow measurements (3 periods per cow )
forage <- read.table( "forage.dat", header = T )
forage$code <- factor( forage$code )
forage$trt <- ordered( forage$trt )
forage$per <- ordered( forage$per, c("Early","Middle","Late") )

forage.fit <- aov( log10( dmi )~hc*trt*per+Error( cow/hc*trt ), forage, qr = T )
print( summary( forage.fit ) )

# projections
forage.proj <- proj( forage.fit )
forage.wpint <- apply( forage.proj$cow[, c("hc","trt","hc:trt") ], 1, sum )
forage.spint <- apply( forage.proj$Within[,
	c("per","hc:per","trt:per","hc:trt:per") ], 1, sum )
forage.mean <- forage.proj[[1]][1]