Color-coding 95% confidence ellipses for centroids

Ideally, the function could return the x and y coordinates used in the ellipses. Here's a simple modification of ordiellipse() that adds xy to each list element of the function output:

ordiellipse2 <- function (ord, groups, display = "sites", kind = c("sd", "se"),        conf, draw = c("lines", "polygon", "none"),        w = weights(ord, display), col = NULL, alpha = 127,        show.groups, label = FALSE,  ...){    weights.default <- function(object, ...) NULL    kind <- match.arg(kind)    draw <- match.arg(draw)    pts <- scores(ord, display = display, ...)    ## ordiellipse only works with 2D data (2 columns)    pts <- as.matrix(pts)    if (ncol(pts) > 2)        pts <- pts[ , 1:2, drop = FALSE]    if (ncol(pts) < 2)        stop("ordiellipse needs two dimensions")    w <- eval(w)    if (length(w) == 1)        w <- rep(1, nrow(pts))    if (is.null(w))        w <- rep(1, nrow(pts))    ## make semitransparent fill    if (draw == "polygon" && !is.null(col))        col <- rgb(t(col2rgb(col)), alpha = alpha, maxColorValue = 255)    if (!missing(show.groups)) {        take <- groups %in% show.groups        pts <- pts[take, , drop = FALSE]        groups <- groups[take]        w <- w[take]    }    out <- seq(along = groups)    inds <- names(table(groups))    res <- list()    if (label)        cntrs <- names <- NULL    ## Remove NA scores    kk <- complete.cases(pts)    for (is in inds) {        gr <- out[groups == is & kk]        if (length(gr) > 2) {            X <- pts[gr, ]            W <- w[gr]            mat <- cov.wt(X, W)            if (kind == "se")                mat$cov <- mat$cov/mat$n.obs            if (missing(conf))                t <- 1            else t <- sqrt(qchisq(conf, 2))            xy <- vegan:::veganCovEllipse(mat$cov, mat$center, t)            if (draw == "lines")                vegan:::ordiArgAbsorber(xy, FUN = lines,                        col = if(is.null(col)) par("fg") else col,                        ...)            else if (draw == "polygon")                 vegan:::ordiArgAbsorber(xy[, 1], xy[, 2], col = col, FUN = polygon,                        ...)            if (label && draw != "none") {                cntrs <- rbind(cntrs, mat$center)                names <- c(names, is)            }            mat$scale <- t            mat$xy <- xy   # this is the only substantive change I made            res[[is]] <- mat        }    }    if (label && draw != "none") {        if (draw == "lines")            vegan:::ordiArgAbsorber(cntrs[,1], cntrs[,2], labels=names, col = col,                      FUN = text, ...)    else         vegan:::ordiArgAbsorber(cntrs, labels = names, col = NULL,                FUN = ordilabel, ...)    }    class(res) <- "ordiellipse"    invisible(res)}

Now you can create your stats object, per your example, but specify draw = "none", because you'll want to play withstats` on your own to have full control of the ellipse plotting (no need to specify graphical arguments):

stats<-with(dat.fac, ordiellipse2(ord, Site_TP, kind = "se", conf = 0.95, draw = "none"))

notice now that stats has the same form, but there is an element xy for each ellipse:

names(stats[[1]])[1] "cov"    "center" "n.obs"  "wt"     "scale"  "xy"    

You can plot these in a loop using the polygon() function, but better would be to place all the xy coordinates in a single matrix or data.frame, so that you can just make a single call to polygon(), but with vectorized arguments, such as col, lwd. I'll extract the xy matrices using lapply(), then collapse them together as a data.frame. Key is the separate the individual ellipses with NAs, so the polygon drawing knows where one stops and the next starts:

xy.mat <- t(as.data.frame(lapply(stats, function(x){            t(rbind(x$xy, NA))        })))

Plotting becomes easy:

plot(ord, display = "sites", type = "n")polygon(xy.mat[, 1], xy.mat[, 2], col = c(rep("#0000ff20", 5), rep("#ff000030", 5)), border = c(rep("blue", 5), rep("red", 5)))

Now it's just a matter of ordering the graphical arguments to polygon() correctly.