summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
require(plyr) #New version of length which can handle NA's: if na.rm==T,
don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else       length(x)
}

# This is does the summary; it's not easy to understand...
datac <- ddply(data, groupvars, .drop=.drop,
           .fun= function(xx, col, na.rm) {
                   c( N    = length2(xx[,col], na.rm=na.rm),
                      mean = mean   (xx[,col], na.rm=na.rm),
                      sd   = sd     (xx[,col], na.rm=na.rm)
                      )
                  },
            measurevar,
            na.rm
     )

qplot(temp,wl,data=cleant, facets=.~habitat, geom=c("point","smooth"), method="lm")

model <- lme(shell.size ~habitat * temp *city, random = ~1|colony, data = FL.snails)

anova(model, type = 'marginal')

28L, 28L,....., 28L, 26L,...... 26L, 26L, 30L,...... 30L, ..........0.683, 1.283)), .Names = c("colony", "individual", "city", "habitat", "temp", "shell.size"), class = "data.frame", row.names = c(NA, -5471L))