Dear List, I am writing a custom panel function and xyplot method to plot the results of a procrustes analysis from the vegan package.
I am having trouble getting the call to panel.arrows to work as I wish when conditioning. The attached file contains the function definitions for the xyplot method and the custom panel and prepanel functions I am using. This example, using data and functions from the vegan package illustrates the problem. require(vegan) require(lattice) data(varespec) vare.dist <- vegdist(wisconsin(varespec)) library(MASS) ## isoMDS mds.null <- isoMDS(vare.dist, tol=1e-7) mds.alt <- isoMDS(vare.dist, initMDS(vare.dist), maxit=200, tol=1e-7) vare.proc <- procrustes(mds.alt, mds.null) vare.proc groups <- factor(c(rep(1,16), rep(2,8)), labels = c("grazed","ungrazed")) source("xyplot.procrustes.R") xyplot(vare.proc, y ~ x | groups, data = as.data.frame(groups), kind = 1) The resulting plot has too many arrows on each panel - some points have multiple arrows emanating from they. panel.procrustes() is defined as: `panel.procrustes` <- function(x, y, kind, choices, rotation, X, ar.col, length = 0.05, ...) { tp <- trellis.par.get() if(missing(ar.col)) ar.col <- tp$superpose.symbol$col[1] if(kind == 1) { panel.abline(h = 0, lty = "dashed") panel.abline(v = 0, lty = "dashed") if(ncol(rotation) == 2) { ## Sometimes rotation[1,1] is 2.2e-16 above one rotation[1,1] <- min(rotation[1,1], 1) panel.abline(0, tan(acos(rotation[1, 1])), lty = "solid") panel.abline(0, 1/tan(acos(-rotation[1, 1])), lty = "solid") } else { Y <- cbind(x,y) %*% t(rotation) for (k in seq_len(ncol(Y))) { tmp <- matrix(0, nrow = 2, ncol = ncol(Y)) tmp[, k] <- range(Y[, k]) tmp <- tmp %*% rotation panel.lines(tmp[, choices], lty = 1) panel.text(tmp[2, choices[1]], tmp[2, choices[2]], as.character(k)) } } panel.xyplot(x, y, type = "p", ...) ## Problem here panel.arrows(x0 = x, y0 = y, x1 = X[,1], y1 = X[,2], length = length, col = ar.col, ends = "last", ...) ## } else if(kind == 2) { quant <- quantile(y) panel.xyplot(x, y, type = "h", ...) panel.abline(h = quant[2:4], lty = c(2,1,2)) } } The bit I am having trouble with is the call to panel.arrows. The plotting of the points (line above the panel.arrows call) works fine with the conditioning, but I'm not getting the panel.arrows call to condition correctly. So, to my question: how does one tell panel.arrows to plot only the arrows for the relevant conditioning variable? TIA G -- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% Dr. Gavin Simpson [t] +44 (0)20 7679 0522 ECRC, UCL Geography, [f] +44 (0)20 7679 0565 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk Gower Street, London [w] http://www.ucl.ac.uk/~ucfagls/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
`xyplot.procrustes` <- function(object, formula, data = NULL, kind = 1, choices = 1:2, xlab, ylab, main, ...) { require(lattice) || stop("requires package 'lattice'") if (missing(main)) main <- "Procrustes errors" if(kind == 1) { dat <- data.frame(x = object$Yrot[, choices[1]], y = object$Yrot[, choices[2]]) if(missing(xlab)) xlab <- paste("Dimension", choices[1]) if(missing(ylab)) ylab <- paste("Dimension", choices[2]) aspect <- "iso" } else { res <- residuals(object) dat <- data.frame(x = seq_along(res), y = res) if (missing(xlab)) xlab <- "Index" if (missing(ylab)) ylab <- "Procrustes residual" aspect <- "fill" } if (!is.null(data)) dat <- cbind(dat, data) if (missing(formula)) { v <- colnames(dat) formula <- as.formula(paste(v[2], "~", v[1])) } xyplot(formula, data = dat, choices = choices, kind = kind, rotation = object$rotation, X = object$X[, choices], aspect = aspect, xlab = xlab, ylab = ylab, main = main, panel = panel.procrustes, prepanel = prepanel.procrustes, ...) } `panel.procrustes` <- function(x, y, kind, choices, rotation, X, ar.col, length = 0.05, ...) { tp <- trellis.par.get() if(missing(ar.col)) ar.col <- tp$superpose.symbol$col[1] if(kind == 1) { panel.abline(h = 0, lty = "dashed") panel.abline(v = 0, lty = "dashed") if(ncol(rotation) == 2) { ## Sometimes rotation[1,1] is 2.2e-16 above one rotation[1,1] <- min(rotation[1,1], 1) panel.abline(0, tan(acos(rotation[1, 1])), lty = "solid") panel.abline(0, 1/tan(acos(-rotation[1, 1])), lty = "solid") } else { Y <- cbind(x,y) %*% t(rotation) for (k in seq_len(ncol(Y))) { tmp <- matrix(0, nrow = 2, ncol = ncol(Y)) tmp[, k] <- range(Y[, k]) tmp <- tmp %*% rotation panel.lines(tmp[, choices], lty = 1) panel.text(tmp[2, choices[1]], tmp[2, choices[2]], as.character(k)) } } panel.xyplot(x, y, type = "p", ...) panel.arrows(x0 = x, y0 = y, x1 = X[,1], y1 = X[,2], length = length, col = ar.col, ends = "last", ...) } else if(kind == 2) { quant <- quantile(y) panel.xyplot(x, y, type = "h", ...) panel.abline(h = quant[2:4], lty = c(2,1,2)) } } `prepanel.procrustes` <- function(x, y, X, choices, kind, ...) { if(kind == 1) { xlim <- range(x, X[, choices[1]]) ylim <- range(y, X[, choices[2]]) } else { xlim <- range(x) ylim <- range(y) } list(ylim = ylim, xlim = xlim) }
______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.