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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
-------------- next part --------------
`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)
}
On Thu, Aug 7, 2008 at 7:55 AM, Gavin Simpson <gavin.simpson at ucl.ac.uk> wrote:> 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.You need to use the proper subset of rows of X: `panel.procrustes` <- function(x, y, kind, choices, rotation, X, ar.col, length = 0.05, ..., subscripts) { 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[subscripts ,1], y1 = X[subscripts, 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, ..., subscripts) { if(kind == 1) { xlim <- range(x, X[subscripts, choices[1]]) ylim <- range(y, X[subscripts, choices[2]]) } else { xlim <- range(x) ylim <- range(y) } list(ylim = ylim, xlim = xlim) } -Deepayan
On Thu, 2008-08-07 at 14:11 -0700, Deepayan Sarkar wrote:> On Thu, Aug 7, 2008 at 7:55 AM, Gavin Simpson <gavin.simpson at ucl.ac.uk> wrote: > > Dear List,<snip />> > You need to use the proper subset of rows of X:[Apologies for the delay in responding --- I have been offline for several days] Thank you very much Deepayan. I thought that was the problem, but was not aware of the correct incantation to select the relevant subsets. All the best, G> > `panel.procrustes` <- > function(x, y, kind, choices, rotation, X, > ar.col, length = 0.05, ..., subscripts) > { > 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[subscripts ,1], y1 = X[subscripts, 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, ..., subscripts) { > if(kind == 1) { > xlim <- range(x, X[subscripts, choices[1]]) > ylim <- range(y, X[subscripts, choices[2]]) > } else { > xlim <- range(x) > ylim <- range(y) > } > list(ylim = ylim, xlim = xlim) > } > > -Deepayan-- %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~% 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 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%