wasp <- read.table( "wasp.dat", header = T )

# analysis of covariance
wasp.cov <- aov( HW ~ G1Wa * caste, wasp )
wasp.par <- aov( HW ~ G1Wa + caste, wasp )

wasp.cova <- aov( HW ~ TW + G1Wa + caste, wasp)

wasp$Hres <- resid( aov( HW ~ TW, wasp ) )
wasp$Gres <- resid( aov( G1Wa ~ TW, wasp ) )

wasp.covr <- aov( Hres ~ TW + G1Wa + caste, wasp )
wasp.covwr <- aov( Hres ~ G1Wa + caste, wasp )

# Multivariate Stuff!

library(MASS)

# discriminant analysis
wasp.lda <- lda(wasp[,-1],wasp$caste)
wasp$DA <- predict(wasp.lda,wasp[,-1])$x

# principal component analysis
wasp.pr  <- prcomp(as.matrix(wasp[,-1]))
wasp$PCA <- wasp.pr$x[,1]

## t statistics
tmp <- wasp$caste == "Q"
wasp.t <- apply( wasp[,-1], 2, function( x ) {
   ( mean( x[tmp] ) - mean( x[!tmp] ) ) /
   sqrt( ( var( x[tmp] ) + var( x[!tmp] ) ) / sum( tmp ) ) } )
rm( tmp )
wasp.t <- sort( wasp.t )