# read data
#cloning <- read.table( "cloning.dat", header=T )
cloning <- cloning[ !is.na(cloning$titer), ]
cloning$clone <- factor(cloning$code)

# ANOVA table
cloning.fit <- aov( titer ~ clone, cloning )
summary( cloning.fit )

# estimate means
attach( cloning )
cloning.means <- data.frame( clone=unique(code), pred=tapply( titer, code, mean ) )
cloning.means$sd <- tapply( titer, code, function( x ){ sqrt( var( x ) ) } )
cloning.means$n <- tapply( titer, code, length )
cloning.means$se <- cloning.means$sd / sqrt( cloning.means$n )
cloning.means$type <- type[ !duplicated( code ) ]
detach( )

# or try the PDA routine:
cloning.lsm <- lsmean( cloning.fit )

# estimate common SD
cloning.fit
std.dev( cloning.fit )

# preselected contrasts
aov(titer ~ clone, cloning, 
   subset=is.na( match( cloning$type, c("cult","etb","odd") ) ) )

# critical values for multiple comparisons (from SAS output)
cloning.crit <- data.frame( all=c(306,555,737,530,278),
   bal=c(254,435,527,414,232),
   row.names=c("LSD","BON","SCHEFFE","HSD","WALLER") )
for (i in names(cloning.crit))
   names(cloning.crit[[i]]) <- row.names(cloning.crit)

# critical values for multi-stage testing (MST) multiple comparisons
cloning.mst <- read.table( "cloncrit.dat", header=T )
tmp <- cloning.mst$set
on.exit(rm(tmp))
cloning.mst$set <- NULL
cloning.mst <- split(cloning.mst,tmp)