# this is a comment
# NOTE: lsmean, mplot and std.dev are from library(pda)
library( pda )
# must turn on graphics device before plotting
motif( ) # on Splus or x11( ) on S

# read data
forcow <- read.table( "forcow.dat", header=T )
# trt is ordered factor
forcow$trt <- ordered( forcow$trt )

# I do heifer, you do cow
forcowh <- forcow[forcow$hc=="H", ]

# stem-and-leaf plot of raw data
stem( forcowh$dmi, 5 )
# box-plots by treatment
boxplot( split( forcowh$dmi, forcowh$trt ) )

# analysis of variance for heifer
forcowh.aov <- aov( dmi ~ trt, forcowh )
forcowh.aov
summary( forcowh.aov )

# linear model for heifer -- same idea
forcowh.lm <- lm( dmi ~ trt, forcowh )
forcowh.lm
summary( forcowh.lm )		# linear models summary
summary.aov( forcowh.lm )		# anova-type summary

# ordinary means
tapply( forcowh$dmi, forcowh$trt, mean )
# predicted values for all responses
data.frame( trt=forcowh$trt, mean=predict( forcowh.aov ) )
# predicted values and their standard errors by trt
forcowh.mean <- predict( forcowh.aov, se=T )
forcowh.mean <- data.frame( trt=forcowh$trt, mean=forcowh.mean$fit,
   se=forcowh.mean$se.fit )
forcowh.mean[ !duplicated( forcowh$trt ), ]
# least squares means
lsmean( forcowh.aov )

# polynomial contrasts
coef( forcowh.aov ) / sqrt(c(1,10,14,10,70))	# or coef( forcowh.lm )
# use summary.lm to get standard errors
summary.lm( forcowh.aov )		# or summary( forcowh.lm )

# S-Plus plot of anova object
plot( forcowh.aov )

# residual plot with symbols
mplot( predict( forcowh.aov ), resid( forcowh.aov ),
   group=forcowh$trt,
   xlab="mean by trt", ylab="residual" )
# residual plot with jittered symbols
mplot( jitter( predict( forcowh.aov ) ), resid( forcowh.aov ),
   group=forcowh$trt,
   xlab="mean by trt", ylab="residual" )
# add horizontal lines for 0 and +/- 1 SD
abline( h=0, lty=2 )
abline( h=c(-1,1) * std.dev( forcowh.aov ), lty=3 )