tmpar <- par( mar = c(3.1,3.1,0,0), mfrow = c(1,2) )
on.exit( {par( tmpar ); rm( data, ylim, xlab, ylab, lsm, tmp, small, tmpar )} )

xlab <- c("(b) control","(a) tree")
ylab <- c("","pure split plot score")
xlim <- c(1,7.7)
small <- c(.7,1.6)

lsm <- lsmean( trees.spaov, factors = c("soil","date") )
lsm$type <- factor( !is.na( match( lsm$soil,
   c("dtac2","nc11004","ne-332","nm-6") ) ) )
ylim <- range( lsm$pred )

lsm <- split( lsm, lsm$type )
# following values are SE for contrasts take for SAS output
lsm[[1]]$se <- rep( .42023472 / sqrt( 2 ), nrow( lsm[[1]] ) )
lsm[[2]]$se <- rep( .35818796 / sqrt( 2 ), nrow( lsm[[2]] ) )

data <- split( trees.sp, trees.sp$type )

for ( i in 2:1 ) { # "FALSE" comes before "TRUE"
   lsd.plot( trees.spaov, data[[i]], c("date","soil"), lsm[[i]],
      xpos = 6, xlim = xlim, ylim = ylim, xaxt = "n", xlab = "", ylab = "" )
   axis(1,1:6)
   mtext( paste( xlab[i], "soils by date" ), 1, 2 )
   mtext( ylab[i], 2, 2 )
   tmp <- lsm[[i]][lsm[[i]]$date==6,]
   text( rep( 6.2, nrow( tmp ) ), uncollide( tmp$pred, small[[i]] ),
      as.character( tmp$soil ), adj = 0 )
}


###########################################################################

#lsd.plot( trees.spaov, factors = c("pyrrat","crop"),
#   xpos = .6, ypos = 3.25, xlab = "", ylab = "" )
#mtext( "(b) pyr rate by crop", 1, 2 )