understanding better argument scale=T in rda

[aguilart]

For really long I thought that the argument of scale=T in rda was doing z scaling: (xi-mean(x))/sd(x) prior to pca calcuations . However, recently (after reading about how pca is calculating uising rda) I found that z scaling is always calcuated (even when scale=F). For example, using varechem:

example<-rda(varechem)

The z-scaled values are already calculated:

all(round(example$Ybar[,1],3)==  
round(((varechem$N-mean(varechem$N))/sd(varechem$N))/0.86754,3))

that is TRUE

(Side note: I stil do not know why the number are not identical and there is some scaling factor, here being 0.8675)

Now, if I run:

example2<-rda(varechem,scale=T)

The scaled values are different as in example1. They were just "re-scaled" (in this case by 0.1808945)

example2$Ybar[,1]/example$Ybar[,1]

In the documentation it says that the argument scale "Scale species to unit variance (like correlations)". So my main questions are: what variance exactly? the total variance of all the data? how is this calculated? And why is this necessary since all the data is z-scaled already to units of standard deviation?

Thanks for all the support!

[jarioksa]

> m <- rda(varechem)
> zapsmall(colSums(m$Ybar^2))
       N        P        K       Ca       Mg        S       Al       Fe 
   30.56   223.39  4204.53 59332.17  1681.93   136.14 14962.28  3654.07 
      Mn       Zn       Mo Baresoil Humdepth       pH 
 1149.98     8.90     0.06   270.26     0.44     0.05 
> zapsmall(apply(varechem, 2, var))
       N        P        K       Ca       Mg        S       Al       Fe 
   30.56   223.39  4204.53 59332.17  1681.93   136.14 14962.28  3654.07 
      Mn       Zn       Mo Baresoil Humdepth       pH 
 1149.98     8.90     0.06   270.26     0.44     0.05 
> m1 <- rda(varechem, scale=TRUE)
> colSums(m1$Ybar^2)
       N        P        K       Ca       Mg        S       Al       Fe 
       1        1        1        1        1        1        1        1 
      Mn       Zn       Mo Baresoil Humdepth       pH 
       1        1        1        1        1        1 
> head(m$Ybar/m1$Ybar, 3)
          N       P        K      Ca       Mg        S       Al      Fe
18 5.528084 14.9464 64.84236 243.582 41.01134 11.66783 122.3204 60.4489
15 5.528084 14.9464 64.84236 243.582 41.01134 11.66783 122.3204 60.4489
24 5.528084 14.9464 64.84236 243.582 41.01134 11.66783 122.3204 60.4489
         Mn       Zn        Mo Baresoil  Humdepth        pH
18 33.91134 2.983939 0.2399804 16.43944       NaN 0.2140026
15 33.91134 2.983939 0.2399804 16.43944       NaN 0.2140026
24 33.91134 2.983939 0.2399804 16.43944 0.6646346 0.2140026
> all.equal(crossprod(m$Ybar), cov(varechem))
[1] TRUE
> all.equal(crossprod(m1$Ybar), cor(varechem))
[1] TRUE

[aguilart]

Thanks a lot! It is now clearer: Using scale=T I am re-scaling each deviation to get a proportion/percentage of the standard deviation. That is in your example:

head(m$Ybar/m1$Ybar, 3)
          N       P        K      Ca       Mg        S       Al      Fe
18 5.528084 14.9464 64.84236 243.582 41.01134 11.66783 122.3204 60.4489
15 5.528084 14.9464 64.84236 243.582 41.01134 11.66783 122.3204 60.4489
24 5.528084 14.9464 64.84236 243.582 41.01134 11.66783 122.3204 60.4489
         Mn       Zn        Mo Baresoil  Humdepth        pH
18 33.91134 2.983939 0.2399804 16.43944       NaN 0.2140026
15 33.91134 2.983939 0.2399804 16.43944       NaN 0.2140026
24 33.91134 2.983939 0.2399804 16.43944 0.6646346 0.2140026

all those values are the square root the variances of each variables (so,the sd).

I have a follow up question. As I understand, the next step in rda is to calculate the eigenvectors and eigenvalues from the covariance matrix of m$Ybar. In which cases is better to have, then, m$Ybar scaled or not? I still cannot understand what is the advantage of the scaling.

[jarioksa]

You need to use scaling if the variables are not measured in same units. The varechem is such a case. If you change some of the chemical variables to use, say, millimoles instead of milligrams, or to use micrograms instead of milligrams the variance will change and so will change the analysis. Changing the unit of measurement will change the variance, and the higher the variance, the higher the weight of the variable in the analysis. RDA/PCA explain variance, and if you use unscaled data, all that they need to do is to model Ca and Al, but you do not need to care about pH since its contribution to the total variance is so small. With scale=TRUE all variables have equal weight and will be treated equally. Community data is a special case: all abundances are measure in same units. However, even with community data, the most abundant species usually have the highest variance and the analysis will focus on them. With scale=TRUE all species have equal importance in the analysis. Outside community ecology, RDA/PCA is almost always made on scaled and equalized data because the data usually are measured in arbitrary and different units.

Compare

> biplot(rda(varechem), dis="sp")
> biplot(rda(varechem, scale=TRUE), dis="sp")

Or for community data:

> biplot(rda(varespec), dis="sp")
> biplot(rda(varespec, scale=TRUE), dis="sp")