tree <- read.table( "tree.dat", header = T )
tree <- tree[ tree$acc == "notauto", ]
tree$acc <- NULL
tree$type <- factor( !is.na( match( tree$soil,
   c("dtac2","nc11004","ne-332","nm-6") ) ) )
tree$reps <- factor( tree$reps )

trees <- tree[ rep( seq( nrow( tree ) ), 6 ), c("soil","reps","type") ]
trees$date <- ordered( rep( 1:6, rep( nrow( tree ), 6 ) ) )
trees$score <- unlist( tree[ , paste( "score", 1:6, sep="" ) ] )
row.names(trees) <- interaction( trees[, c("soil","reps","date") ] )

trees.soil <- aov(score ~ soil * date + Error(reps %in% soil), trees)
trees.type <- aov(score ~ type/soil * date + Error(reps %in% soil %in% type),
   trees)

# projections (*** CAUTION: SLOW ***)
trees.proj <- proj ( trees.type )
trees.ref <- data.frame(Reference=c(trees.proj[[1]]))
for (i in names(trees.proj)[-1])
   trees.ref[[i]] <- apply( trees.proj[[i]], 1, sum )

# whole plot
trees.wp <- nested( trees.proj, trees, c("type","soil","reps"), "score",
   c(1,2), trees.ref )

trees.wpaov <- aov( score ~ soil, trees.wp )

# split plot
trees.sp <- nested( trees.proj, trees, c("type","soil","reps","date"),
   "score", c(1,3), trees.ref )

trees.spaov <- aov( score ~ reps*soil + date + date:soil, trees.sp )

# combined (sort of)
trees.all <- aov( score ~ soil * date, trees )

# split plot with time trend removed
trees.sp$res <- resid( aov( score ~ reps*soil + date, trees.sp ) ) +
   lsmean( trees.spaov, fac=NULL )$pred

trees.spres <- aov( res ~ reps*soil + date + date:soil, trees.sp )