diet <- read.table( "diet.dat", header = T, na.strings = "." )
diet$trt <- factor(diet$trt)
diet$block <- factor(diet$block)#
# ordinary anova
diet.fit <- aov(dmi ~ trt, diet)
# summaries of anova fit
print( diet.fit )
print( summary( diet.fit ))
# splus routine calculates group means
# but chokes on standard errors because of imbalance
print( model.tables( diet.fit, type = "means", se = T ))
# another way to get means, as predicted values of observations
print( predict( diet.fit, se = T ))
# or use Yandell's PDA library to get means and SEs
library(pda)
print( lsmean( diet.fit ))
# histogram (plotter) or stem-and-leaf (printer)
# be sure to turn on plotter first, such as motif()
print( stem( diet$dmi ))
hist( diet$dmi, nclass = 20, prob = T )
# plot of group means against individual values
# with some jitter to make points more readable
diet.jit <- jitter( predict( diet.fit ))
plot( diet.jit, diet$dmi, type = "n",
xlab = "group means", ylab = "DMI" )
text( diet.jit, diet$dmi, as.character( diet$trt ))
abline( a = 0, b = 1, lty = 2 )
# plot of predicted with jitter against residuals
plot( diet.jit, resid( diet.fit ), type = "n",
xlab = "group means", ylab = "residuals" )
text( diet.jit, resid( diet.fit ), as.character( diet$trt ))
abline( h = 0, lty = 2 )
abline( h = std.dev( diet.fit ) * c(-1,1), lty = 3 )
# or consider the default plots of Splus
plot( diet.fit )
# polynomial contrasts
contrasts(diet$trt)<-contr.poly( 6 )
diet.poly <- aov( dmi ~ trt, diet )
print( summary( diet.poly ))
print( summary.lm( diet.poly ))
# coefficients
print( coef( diet.poly ))
# reconstruct treatment means
print( cbind(1,contrasts(diet$trt)) %*% matrix(coef(diet.poly),6) )
# orthogonal contrasts for unbalanced data
# note that these do not give the same coefficients as above: why?
diet$trtlin <- unfactor( diet$trt )
diet$trtquad <- diet$trtlin ^ 2
diet.orth <- aov( dmi ~ trtlin + trtquad + trt, diet )
print( diet.orth )
print( coef( diet.orth ))
print( cbind( 1, 1:6, (1:6)^2, contrasts(diet$trt)[,3:5] ) %*%
matrix( coef( diet.orth ), 6 ))