Dear List,
I am having trouble with some R code I have written to perform
Redundancy Analysis (RDA) on a matrix of species abundance data (Y) and
a matrix of environmental data (X).
RDA is a constrained form of PCA and can be thought of as a PCA of the
fitted values of a regression of each variable in Y on all variables in
X.
For info, the first use of RDA is in:
Rao, C.R, 1964. The use and interpretation of principal component
analysis in applied research. Sankhyaa, Ser. A 56: 329-358.
Rao proposed the problem of RDA as an exercise at the end of Chapter 8
in his book on multivariate analysis (Rao, C.R., 1973. Linear
statistical inference and its applications. 2nd edition. Wiley, New
York).
The algorithm I am following is given in Chapter 11 pages 580-587 of:
Legendre, P. and Legendre, L. 1998. Numerical Ecology. 2nd English
edition. Elsevier.
The problem I am having is with some matrix multiplication within the
code. You can see this problem I have by running the following example
using my code at the end of the email:
spp.dat <- matrix(c(4,2,0,2,0,5,1,3,2,4,0,2,2,0,3,1),4,4)
env.dat <- matrix(c(1.5,2.3,2,1.6,0.9,0.8,1.2,1.5),4,2)
test.rda <- rda(spp.dat, env.dat)
Results in the following error:
Error in df %*% t(Yhat) : non-conformable arguments
This is referring to the part of the code where I wish to do the
following:
S<sub>Yhat'Yhat</sub> = [1/(n-1)] Yhat' Yhat,
where [1/(n-1)] = df = degrees of freedom. is effectively dividing my
t(Yhat) %*% Yhat by the degrees of freedom, n being the number of
observations or sites.
df is a column matrix of the form:
[,1]
[1,] 0.3333333
[2,] 0.3333333
[3,] 0.3333333
[4,] 0.3333333
With Yhat' and Yhat both being matrices of the form:
Yhat'
[,1] [,2] [,3] [,4]
[1,] 2.1462687 -0.5641791 -0.9761194 -0.6059701
[2,] -1.9487562 1.5880597 0.8587065 -0.4980100
[3,] 0.2646766 0.9791045 -0.1472637 -1.0965174
[4,] 0.2701493 -0.6134328 -0.1089552 0.4522388
And Yhat:
[,1] [,2] [,3] [,4]
[1,] 2.1462687 -1.9487562 0.2646766 0.2701493
[2,] -0.5641791 1.5880597 0.9791045 -0.6134328
[3,] -0.9761194 0.8587065 -0.1472637 -0.1089552
[4,] -0.6059701 -0.4980100 -1.0965174 0.4522388
I am stumped as to what I am doing wrong and how to go about solving it.
If anyone could indicate how to go about solving this problem I would be
most grateful.
Gavin Simpson
Windows 2000
Version:
_
platform i386-pc-mingw32
arch x86
os Win32
system x86, Win32
status
major 1
minor 4.1
year 2002
month 01
day 30
language R
code is below:
rda <- function(y, x, sqroot=FALSE, stdy=FALSE, stdx=FALSE)
{
y <- as.matrix(y) #species data matrix
x <- as.matrix(x) #environmental or constraining variables
n <- length(y[,1]) #number of observations
if (stdy==FALSE){
y <- scale(y, center=TRUE, scale=FALSE)
} else {
y <- scale(y, center=TRUE, scale=TRUE)
}
if (stdx==FALSE){
x <- scale(x, center=TRUE, scale=FALSE)
} else {
x <- scale(x, center=TRUE, scale=TRUE)
}
if (sqroot==TRUE){
y <- sqrt(y)
}
B <- solve(t(x) %*% x) %*% t(x) %*% y
Yhat <- x %*% B
df <- as.matrix(rep((1/(n-1)), times=n))
S <- df %*% t(Yhat) %*% Yhat
svd.fitted <- svd(S)
qr.rank <- qr(S)$rank
eig <- svd.fitted$d[1:qr.rank]
U1 <- svd.fitted$u[,1:qr.rank]
U2 <- U1 %*% diag(svd.fitted$d[1:qr.rank]^0.5)
F1 <- y %*% svd.fitted$u[,1:qr.rank]
F2 <- F1 %*% diag(1/svd.fitted$d[1:qr.rank]^0.5)
Z1 <- Yhat %*% svd.fitted$u[,1:qr.rank]
Z2 <- Z1 %*% diag(1/svd.fitted$d[1:qr.rank]^0.5)
r <- cor(F1, Z1)
C <- B %*% svd.fitted$u[,1:qr.rank]
bip.xvars <- cor(x, Z1)
tmp <- list(U1 = U1, U2 = U2, F1 = F1, F2 = F2, Z1 = Z1, Z2 = Z2,
eig = eig, r = r, C = C, bip.xvars = bip.xvars)
class(tmp) <- "rda"
return(tmp)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._