# read data,  fit models
feed <- read.table( "feed.dat", header=T )
feed$trt <- factor( feed$trt )

# adjusted treatment -- analysis of covariance
feed.ancova <- aov( bwg ~ trt + fi, feed )
# regression type fit
feed.reg <- aov( bwg ~ fi + la + cla, feed )
# simple regression
feed.fi <- aov( bwg ~ fi, feed)

# removing factor or covariate
feed$bwgfi <- resid( feed.fi ) + lsmean( feed.fi )$pred
feed.fi2 <- aov( bwgfi ~ fi + trt, feed )
feed.fi1 <- aov( bwgfi ~ trt, feed )
feed.trt <- aov( bwg ~ trt, feed )
feed$bwgtrt <- resid( feed.trt ) + lsmean( feed.trt, fac = NULL )$pred
feed.trt2 <- aov( bwgtrt ~ fi + trt, feed )
feed.trt1 <- aov( bwgtrt ~ fi, feed )

# checking for any interaction
feed.int <- aov( bwg ~ fi * trt, feed )
feed$res <- resid( feed.ancova )
feed.res <- aov( res ~ fi * trt, feed )
# checking for high interaction
feed$hi <- as.numeric((feed$cla + feed$la) == .5)
feed.hi <- aov( bwg ~ fi + trt + fi:hi, feed )