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

xlab <- c("(b) control","(a) tree")
ylab <- c("","soil sample means")
grmean <- lsmean( trees.wpaov, fac = NULL )$pred

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

lsm <- split( lsm, lsm$type )
ylim <- range( trees.wp$score )

data <- split( trees.wp, trees.wp$type )
for ( i in 2:1 ) { # "FALSE" comes before "TRUE"
   lsm[[i]]$rank <- rank( -lsm[[i]]$pred )

   # have to be careful about levels of soil in type subsets
   data[[i]]$soil <- factor( data[[i]]$soil )
   attach( data[[i]] )
   names( lsm[[i]]$rank ) <- unique( soil )
   plot( jitter( lsm[[i]]$rank[soil] ), score, xaxt = "n",
      ylim = ylim, xlab = "", ylab = "" )
   detach()
   axis( 1, lsm[[i]]$rank, row.names( lsm[[i]] ) )
   abline( h = grmean, lty = 3 )
   mtext( paste( xlab[i], "soil sample means" ), 1, 2 )
   mtext( ylab[i], 2, 2 )
}
axis( 1, lsm[[1]]["grnhous","rank"], "grnhous" )