pcoa <- function(D, correction="none", rn=NULL)
#
# Principal coordinate analysis (PCoA) of a square distance matrix D
# with correction for negative eigenvalues.
#
# References:
# Gower, J. C. 1966. Some distance properties of latent root and vector methods
# used in multivariate analysis. Biometrika. 53: 325-338.
# Gower, J. C. and P. Legendre. 1986. Metric and Euclidean properties of
# dissimilarity coefficients. J. Classif. 3: 5-48.
# Legendre, P. and L. Legendre. 1998. Numerical ecology, 2nd English edition.
# Elsevier Science BV, Amsterdam. [PCoA: Section 9.2]
#
# Pierre Legendre, October 2007
{
centre <- function(D,n)
# Centre a square matrix D by matrix algebra
# mat.cen = (I - 11'/n) D (I - 11'/n)
{ One <- matrix(1,n,n)
mat <- diag(n) - One/n
mat.cen <- mat %*% D %*% mat
}
bstick.def <- function (n, tot.var = 1, ...) # 'bstick.default' from vegan
{
res <- rev(cumsum(tot.var/n:1)/n)
names(res) <- paste("Stick", seq(len = n), sep = "")
return(res)
}
# ===== The PCoA function begins here =====
# Preliminary actions
D <- as.matrix(D)
n <- nrow(D)
epsilon <- sqrt(.Machine$double.eps)
if(length(rn)!=0) {
names <- rn
} else {
names <- rownames(D)
}
CORRECTIONS <- c("none","lingoes","cailliez")
correct <- pmatch(correction, CORRECTIONS)
if(is.na(correct)) stop("Invalid correction method")
# cat("Correction method =",correct,'\n')
# Gower centring of matrix D
# delta1 = (I - 11'/n) [-0.5 d^2] (I - 11'/n)
delta1 <- centre((-0.5*D^2),n)
trace <- sum(diag(delta1))
# Eigenvalue decomposition
D.eig <- eigen(delta1)
# Negative eigenvalues?
min.eig <- min(D.eig$values)
zero.eig <- which(abs(D.eig$values) < epsilon)
D.eig$values[zero.eig] <- 0
# No negative eigenvalue
if(min.eig > -epsilon) { # Curly 1
correct <- 1
eig <- D.eig$values
k <- length(which(eig > epsilon))
rel.eig <- eig[1:k]/trace
cum.eig <- cumsum(rel.eig)
vectors <- sweep(D.eig$vectors[,1:k], 2, sqrt(eig[1:k]), FUN="*")
bs <- bstick.def(k)
cum.bs <- cumsum(bs)
res <- data.frame(eig[1:k], rel.eig, bs, cum.eig, cum.bs)
colnames(res) <- c("Eigenvalues","Relative_eig","Broken_stick","Cumul_eig","Cumul_br_stick")
rownames(res) <- 1:nrow(res)
rownames(vectors) <- names
colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
note <- paste("There were no negative eigenvalues. No correction was applied")
out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace))
# Negative eigenvalues present
} else { # Curly 1
k <- n
eig <- D.eig$values
rel.eig <- eig/trace
rel.eig.cor <- (eig - min.eig)/(trace - (n-1)*min.eig) # Eq. 9.27 for a single dimension
rel.eig.cor = c(rel.eig.cor[1:(zero.eig[1]-1)], rel.eig.cor[(zero.eig[1]+1):n], 0)
cum.eig.cor <- cumsum(rel.eig.cor)
k2 <- length(which(eig > epsilon))
k3 <- length(which(rel.eig.cor > epsilon))
vectors <- sweep(D.eig$vectors[,1:k2], 2, sqrt(eig[1:k2]), FUN="*")
# Only the eigenvectors with positive eigenvalues are shown
# Negative eigenvalues: three ways of handling the situation
if((correct==2) | (correct==3)) { # Curly 2
if(correct == 2) { # Curly 3
# Lingoes correction: compute c1, then the corrected D
c1 <- -min.eig
note <- paste("Lingoes correction applied to negative eigenvalues: D' = -0.5*D^2 -",c1,", except diagonal elements")
D <- -0.5*(D^2 + 2*c1)
# Cailliez correction: compute c2, then the corrected D
} else if(correct == 3) {
delta2 <- centre((-0.5*D),n)
upper <- cbind(matrix(0,n,n), 2*delta1)
lower <- cbind(-diag(n), -4*delta2)
sp.matrix <- rbind(upper, lower)
c2 <- max(Re(eigen(sp.matrix, symmetric=FALSE, only.values=TRUE)$values))
note <- paste("Cailliez correction applied to negative eigenvalues: D' = -0.5*(D +",c2,")^2, except diagonal elements")
D <- -0.5*(D + c2)^2
} # End curly 3
diag(D) <- 0
mat.cor <- centre(D,n)
toto.cor <- eigen(mat.cor)
trace.cor <- sum(diag(mat.cor))
# Negative eigenvalues present?
min.eig.cor <- min(toto.cor$values)
zero.eig.cor <- which((toto.cor$values < epsilon) & (toto.cor$values > -epsilon))
toto.cor$values[zero.eig.cor] <- 0
# No negative eigenvalue after correction: result OK
if(min.eig.cor > -epsilon) { # Curly 4
eig.cor <- toto.cor$values
rel.eig.cor <- eig.cor[1:k]/trace.cor
cum.eig.cor <- cumsum(rel.eig.cor)
k2 <- length(which(eig.cor > epsilon))
vectors.cor <- sweep(toto.cor$vectors[,1:k2], 2, sqrt(eig.cor[1:k2]), FUN="*")
# bs <- broken.stick(k2)[,2]
bs <- bstick.def(k2)
bs <- c(bs, rep(0,(k-k2)))
cum.bs <- cumsum(bs)
# Negative eigenvalues still present after correction: incorrect result
} else {
if(correct == 2) cat("Problem! Negative eigenvalues are still present after Lingoes",'\n')
if(correct == 3) cat("Problem! Negative eigenvalues are still present after Cailliez",'\n')
rel.eig.cor <- cum.eig.cor <- bs <- cum.bs <- rep(NA,n)
vectors.cor <- matrix(NA,n,2)
} # End curly 4
res <- data.frame(eig[1:k], eig.cor[1:k], rel.eig.cor, bs, cum.eig.cor, cum.bs)
colnames(res) <- c("Eigenvalues", "Corr_eig", "Rel_corr_eig", "Broken_stick", "Cum_corr_eig", "Cum_br_stick")
rownames(res) <- 1:nrow(res)
rownames(vectors) <- names
colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace, vectors.cor=vectors.cor, trace.cor=trace.cor))
} else { # Curly 2
note <- "No correction was applied to the negative eigenvalues"
bs <- bstick.def(k3)
bs <- c(bs, rep(0,(k-k3)))
cum.bs <- cumsum(bs)
res <- data.frame(eig[1:k], rel.eig, rel.eig.cor, bs, cum.eig.cor, cum.bs)
colnames(res) <- c("Eigenvalues","Relative_eig","Rel_corr_eig","Broken_stick","Cum_corr_eig","Cumul_br_stick")
rownames(res) <- 1:nrow(res)
rownames(vectors) <- names
colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace))
} # End curly 2: three ways of handling the situation
} # End curly 1
class(out) <- "pcoa"
out
} # End of PCoA
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.