anova.mean = function(depvar,   ..., fun=mean, sep="XXXX")
{
# Function applies a function (default)
# mean to each unique combination of factors
# Returns a data frame
# depvar: the dependent variable, a numeric vector
# ... at least 2 factors separated by commas
# sep: the separator used in pasting together the factors
# inside the function

cmat = NULL
for(j in list(...)){
cmat = cbind(cmat, as.character(j))
}

p = cmat[,1]
if(ncol(cmat) < 2)
stop("you must supply at least two factors")
for(j in 2:ncol(cmat)){
p = paste(p, cmat[,j],  sep=sep)
}

pu = unique(p)
m = tapply(depvar, p, fun)
pusplit = unlist(strsplit(names(m), sep))
labs = matrix(pusplit, ncol=ncol(cmat), byrow=T)
names(m) = 1:length(m)
data.frame(m, labs)
}