# [R] Matrix multiplication problem

Gavin Simpson gavin.simpson at ucl.ac.uk
Fri Mar 8 15:02:42 CET 2002

```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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```