I'm running a principal component analysis and I was told that the vector of the scaling 1 are supposed to be of length 1. Here they are enormously bigger than 1. In scaling 2, it's suppose to be less than 1.
Am I doing something wrong?
Also, is it possible, if the data are standardized to draw the equilibrium circle of descriptors in scaling 2?
The data are beak measurements on different species of birds.
the dataset looks like this (named pca.bird.clean):
MBL MBW MBD
1 9.13 7.21 7.55
2 13.97 12.73 14.33
3 8.15 6.43 7.01
4 10.80 9.35 10.02
5 ... ... ...
The PCA:
pca.all.sp=vegan::rda(pca.bird.clean[,c("MBL","MBW", "MBD")], scale=FALSE)
This is the code to produce the image (see here)
cleanplot.pca(pca.all.sp)
Here is a full reproducible example:
library(vegan)
nb.ind.per.sp = 100
nb.sp = rep(nb.ind.per.sp,4)
nind = sum(nb.sp)
sp.id = rep(c(1,2,3,4),each=nb.ind.per.sp)
nsp=4
nb.trait = 3
mu.trait = structure(list(bird1 = c(9.83, 11.64, 11.09),
bird2 = c(6.84, 8.51, 7.09),
bird3 = c(13.17, 14.48, 15.27),
bird4 = c(8.26, 14.21, 8.59)),
.Names = c("bird1", "bird2", "bird3", "bird4"),
row.names = c("MBW_mean", "MBL_mean", "MBD_mean"), class = "data.frame")
sigma.trait= structure(list(bird1 = c(0.99, 0.94, 1.75),
bird2 = c(0.10, 0.22, 0.16), #c(0.10, 0.22, 0.16)
bird3 = c(0.98, 0.76, 1.99),
bird4 = c(0.19, 0.89, 0.29)), # c(0.19, 0.89, 0.29)
.Names = c("bird1", "bird2", "bird3", "bird4"),
row.names = c("MBW_var", "MBL_var", "MBD_var"), class = "data.frame")
ind.trait <- matrix(0, nrow=nind, ncol=dim(mu.trait)[1])
nb.id.sp.last <- 1
for(nb.id.sp in 1:nsp) {
nb.id.sp.seq <- as.numeric(summary(as.factor(sp.id))[nb.id.sp])
for(nbtraits in 1:(dim(mu.trait)[1])) {
traits <- rnorm(n=nb.id.sp.seq,
mean = mu.trait[nbtraits,nb.id.sp],
sd = sigma.trait[nbtraits,nb.id.sp])
ind.trait[nb.id.sp.last:(nb.id.sp.last+nb.id.sp.seq-1),nbtraits] = traits
}
nb.id.sp.last <- nb.id.sp.last + nb.id.sp.seq
}
colnames(ind.trait) = c("MBW","MBL","MBD")
pca.all.sp=vegan::rda(ind.trait[,], scale=FALSE)
"cleanplot.pca" <- function(res.pca, ax1=1, ax2=2, point=FALSE,
ahead=0.07, cex=0.7)
{
# A function to draw two biplots (scaling 1 and scaling 2) from an object
# of class "rda" (PCA or RDA result from vegan's rda() function)
#
# License: GPL-2
# Authors: Francois Gillet & Daniel Borcard, 24 August 2012
require("vegan")
par(mfrow=c(1,2))
p <- length(res.pca$CA$eig)
# Scaling 1: "species" scores scaled to relative eigenvalues
sit.sc1 <- scores(res.pca, display="wa", scaling=1, choices=c(1:p))
spe.sc1 <- scores(res.pca, display="sp", scaling=1, choices=c(1:p))
plot(res.pca, choices=c(ax1, ax2), display=c("wa", "sp"), type="n",
main="PCA - scaling 1", scaling=1)
if (point)
{
points(sit.sc1[,ax1], sit.sc1[,ax2], pch=20)
text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex, pos=3, scaling=1)
}
else
{
text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex, scaling=1)
}
text(res.pca, display="sp", choices=c(ax1, ax2), cex=cex, pos=4,
col="red", scaling=1)
arrows(0, 0, spe.sc1[,ax1], spe.sc1[,ax2], length=ahead, angle=20, col="red")
pcacircle(res.pca)
# Scaling 2: site scores scaled to relative eigenvalues
sit.sc2 <- scores(res.pca, display="wa", choices=c(1:p))
spe.sc2 <- scores(res.pca, display="sp", choices=c(1:p))
plot(res.pca, choices=c(ax1,ax2), display=c("wa","sp"), type="n",
main="PCA - scaling 2")
if (point) {
points(sit.sc2[,ax1], sit.sc2[,ax2], pch=20)
text(res.pca, display="wa", choices=c(ax1 ,ax2), cex=cex, pos=3)
}
else
{
text(res.pca, display="wa", choices=c(ax1, ax2), cex=cex)
}
text(res.pca, display="sp", choices=c(ax1, ax2), cex=cex, pos=4, col="red")
arrows(0, 0, spe.sc2[,ax1], spe.sc2[,ax2], length=ahead, angle=20, col="red")
}
"pcacircle" <- function (pca)
{
# Draws a circle of equilibrium contribution on a PCA plot
# generated from a vegan analysis.
# vegan uses special constants for its outputs, hence
# the 'const' value below.
eigenv <- pca$CA$eig
p <- length(eigenv)
n <- nrow(pca$CA$u)
tot <- sum(eigenv)
const <- ((n - 1) * tot)^0.25
radius <- (2/p)^0.5
radius <- radius * const
symbols(0, 0, circles=radius, inches=FALSE, add=TRUE, fg=2)
}
cleanplot.pca(pca.all.sp)
