plot.design <-
function(x, y = NULL, fun = mean, ..., xaxt = "n",
         xlab = "Factors", ylab = what, ylim, data = NULL, subset,
         ask = TRUE)
{
    na.range <- function(x) range(x[!is.na(x)])
    Tapply <- function(y, fac, fun, ll)
    {
        k <- length(ll)
        ans <- logical(k)
        for(i in 1:k)
            ans[i] <- fun(y[fac == ll[i]])
        names(ans) <- ll
        ans
    }
    Cl <- match.call()
    if(!inherits(x, "data.frame")) {
        if(is.null(y)) {
            if(!is.null(data))
                x <- model.frame(x, data=data)
            else {
                Cl[[1]] <- as.name("model.frame")
                names(Cl)[names(Cl)=="x"] <- "formula"
                Cl$y <- Cl$fun <- Cl$"..." <- Cl$xaxt <- Cl$xlab <-
                    Cl$ylab <- Cl$ylim <- Cl$ask <- NULL
                x <- eval(Cl, parent.frame())
            }
        } else {
            x <- model.frame(x, y)
            y <- NULL
        }
    }
    else if(inherits(y, "formula")) {
        x <- model.frame(y, data = x, na.action = function(z) z)
        y <- NULL
    }

    yname <-
        if(is.null(y)) {# `x' is a model.frame (as above)
            if(length(Terms <- attr(x, "terms")) > 0 &&
               (resp <- attr(Terms, "response")) > 0) { # get response out
                y <- x[[resp]]
                x <- x[,  - resp, drop = FALSE]
                deparse(attr(Terms, "variables")[[1 + resp]])
            }
        } else deparse(substitute(y))
    nrows <- length(attr(x, "row.names"))
    if(nrows == 0) stop("invalid `x' or `y' -- no valid resulting model.frame")
    nvars <- length(x)
    class(x) <- NULL
    allfactors <- logical(nvars)
    nas <- logical(nrows)
    for(j in seq(length = nvars)) {
        xj <- x[[j]]
        allfactors[j] <- !is.null(levels(xj))
        nas <- nas | is.na(xj)
    }
    if(length(y) == nrows) {
        y <- list(y)
        names(y) <- yname
    }
    else {
        y <- design.makey(x, allfactors, y, yname)
        ## a matrix, or a column selector, or default
        yname <- names(y)
    }
    allfactors <- (1:nvars)[allfactors]
    sfun <- substitute(fun)
    what <- paste(if(mode(sfun) == "name") as.character(sfun)
                  else if(is.character(fun)) {
                      what <- fun
                      fun <- get(fun)
                      what
                  } else if(is.function(fun)) deparse(body(fun))
                  else "Summary",
                  "of", yname)
    nf <- length(allfactors)
    statslist <- xrep <- fn <- vector("list", nf)
    nresp <- length(y)
    if(nresp > 1) {
        yr <- NULL
        if(length(ylab) < nresp)
            ylab <- rep(ylab, length = nresp)
        oldpar <- par(ask = ask)
        on.exit(par(oldpar))
    }
    for(j in 1:nresp) {
        yj <- y[[j]]
        out <- is.na(yj) | nas
        if(any(out))
            yj <- yj[!out]
        for(i in 1:nf) {
            xi <- x[[allfactors[i]]]
            if(any(out))
                xi <- xi[!out]
            ll <- fn[[i]] <- levels(xi)
            xrep[[i]] <- rep(i, length(ll))
            statslist[[i]] <- Tapply(yj, xi, fun, ll)
        }
        stats <- unlist(statslist)
        if(missing(ylim))
            ylim <- na.range(stats)
        plot(c(0, nf + 1), ylim, type = "n", xaxt = xaxt, xlab = xlab,
             ylab = ylab[j], ...)
        xr <- unlist(xrep)
        text(xr - 0.1, stats, unlist(fn), adj = 1)
        segments(xr - 0.05, stats, xr + 0.05, stats)
        frange <- sapply(statslist, na.range)
        segments(1:nf, frange[1, ],
                 1:nf, frange[2, ])
        mtext(at = 1:nf, names(x)[allfactors], side = 1, line = 1)
        gm <- fun(yj) # grand mean
        segments(0.5, gm, nf + 0.5, gm)
        yr <- if(nresp > 1) cbind(yr, stats) else stats
    }
    names(xr) <- names(stats)
    if(nresp > 1)
        dimnames(yr) <- list(names(stats), yname)
    invisible(list(x = xr, y = yr, call = Cl))
}

design.makey <-
    function(x, allfactors = sapply(x, function(xx)
                inherits(xx, "factor")), y = NULL, yname = NULL)
{
    csplit <- function(y, yname)
    {
        dy <- dim(y)
        yval <- vector("list", dy[2])
        for(i in 1:dy[2])
            yval[[i]] <- y[, i]
        if(dy[2] == 1)
            names(yval) <- yname
        else names(yval) <- paste(yname, dimnames(y)[[2]])
        yval
    }
    yval <- NULL
    if(is.null(y)) {
        if(length(tt <- attr(x, "terms")))
            y <- attr(tt, "response")
        if(!length(y)) {
            y <- seq(along = allfactors)[!allfactors]
            yname <- names(x)[y]
        }
    }
    else if(inherits(y, "data.frame")) {
        class(y) <- NULL
        yval <- y
    }
    else if(length(dim(y)) == 2)
        yval <- csplit(y, yname)
    else {
        yname <- names(x)
        if(is.character(y))
            y <- pmatch(y, yname)
        yname <- yname[y]
    }
    if(is.null(yval)) {
        if(!length(y))
            stop("no response specified or in data frame")
        if(any(is.na(y)))
            stop("can't find specified response")
        class(x) <- NULL
        ny <- 0
        yval <- list()
        for(i in seq(length = length(y))) {
            ii <- y[i]
            yi <- x[[ii]]
            if(is.matrix(yi)) {
                ny <- ny + dim(yi)[2]
                yval <- c(yval, csplit(yi, yname[i]))
            }
            else {
                ny <- ny + 1
                yval[ny] <- x[ii]
                names(yval)[ny] <- yname[i]
            }
        }
    }
    yval
}
