Change in rotated NMDS scores as a response variable

Hello again listers,
I've had no problems rotating metaMDS ordindation with 2 axes to a vector of
interest. Unfortunately, for my particular dataset stress is > 0.20 for 2
axes.  I'm wondering if it's possible to rotate the first axis for an
ordination with 3 axes somehow, maybe by specifying a 2 by n matrix? That
doesn't seem to work, though. Does anyone have any thoughts on how to
accomplish this (or whether it's a legitimate goal)?
Thanks, Erik

Erik,

It is not documented to work, and it won't work: argument 'vec' must be a
single vector. I don't know how to do this exactly. However, you can achieve
something resembling your goal with Procrustes analysis. One problem is that
you must scale your target matrix variables adequately with respect to each
other. Another problem is that you minimize the distances between your
target scores (two variables) and ordination, and this is not orthogonal to
these target variables. You may also get weird results if your target
variables are correlated. Your responsibility is to take care of all these
aspects (and perhaps some more that I haven't listed). Technically, you can
proceed like this:

xy <- with(varechem, cbind(pH, Humdepth))
xy <- scale(xy)
# check: cor(xy)
# consider: xy <- scores(prcomp(~ pH + Humdepth, data=varechem, scale=TRUE))
ord3 <- metaMDS(varespec, k=3)
pro <- procrustes(xy, ord3) # will re-scale scores of ord3!
rot3 <- fitted(pro)

Cheers, Jari Oksanen

################### EXAMPLE
library(vegan)
data(varespec); data(varechem)
ord1 <- metaMDS(varespec, k = 2)
ord2 <- metaMDSrotate(ord1, varechem$pH) # no problems here

ord3 <- metaMDS(varespec, k = 3)
ord4 <- metaMDSrotate(ord3, varechem$pH) # not "conformable"
ord5 <- metaMDSrotate(ord3, cbind(varechem$pH, varechem$Humdepth)) # doesn't
work with matrix

################ END EXAMPLE