[R] panel.arrows problem in custom panel function

[Gavin Simpson]

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)
}