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



More information about the R-help mailing list