1 pcoa <- function(D, correction="none", rn=NULL)
3 # Principal coordinate analysis (PCoA) of a square distance matrix D
4 # with correction for negative eigenvalues.
7 # Gower, J. C. 1966. Some distance properties of latent root and vector methods
8 # used in multivariate analysis. Biometrika. 53: 325-338.
9 # Gower, J. C. and P. Legendre. 1986. Metric and Euclidean properties of
10 # dissimilarity coefficients. J. Classif. 3: 5-48.
11 # Legendre, P. and L. Legendre. 1998. Numerical ecology, 2nd English edition.
12 # Elsevier Science BV, Amsterdam. [PCoA: Section 9.2]
14 # Pierre Legendre, October 2007
16 centre <- function(D,n)
17 # Centre a square matrix D by matrix algebra
18 # mat.cen = (I - 11'/n) D (I - 11'/n)
19 { One <- matrix(1,n,n)
20 mat <- diag(n) - One/n
21 mat.cen <- mat %*% D %*% mat
24 bstick.def <- function (n, tot.var = 1, ...) # 'bstick.default' from vegan
26 res <- rev(cumsum(tot.var/n:1)/n)
27 names(res) <- paste("Stick", seq(len = n), sep = "")
31 # ===== The PCoA function begins here =====
36 epsilon <- sqrt(.Machine$double.eps)
42 CORRECTIONS <- c("none","lingoes","cailliez")
43 correct <- pmatch(correction, CORRECTIONS)
44 if(is.na(correct)) stop("Invalid correction method")
45 # cat("Correction method =",correct,'\n')
47 # Gower centring of matrix D
48 # delta1 = (I - 11'/n) [-0.5 d^2] (I - 11'/n)
49 delta1 <- centre((-0.5*D^2),n)
50 trace <- sum(diag(delta1))
52 # Eigenvalue decomposition
53 D.eig <- eigen(delta1)
55 # Negative eigenvalues?
56 min.eig <- min(D.eig$values)
57 zero.eig <- which(abs(D.eig$values) < epsilon)
58 D.eig$values[zero.eig] <- 0
60 # No negative eigenvalue
61 if(min.eig > -epsilon) { # Curly 1
64 k <- length(which(eig > epsilon))
65 rel.eig <- eig[1:k]/trace
66 cum.eig <- cumsum(rel.eig)
67 vectors <- sweep(D.eig$vectors[,1:k], 2, sqrt(eig[1:k]), FUN="*")
71 res <- data.frame(eig[1:k], rel.eig, bs, cum.eig, cum.bs)
72 colnames(res) <- c("Eigenvalues","Relative_eig","Broken_stick","Cumul_eig","Cumul_br_stick")
73 rownames(res) <- 1:nrow(res)
75 rownames(vectors) <- names
76 colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
77 note <- paste("There were no negative eigenvalues. No correction was applied")
78 out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace))
80 # Negative eigenvalues present
85 rel.eig.cor <- (eig - min.eig)/(trace - (n-1)*min.eig) # Eq. 9.27 for a single dimension
86 rel.eig.cor = c(rel.eig.cor[1:(zero.eig[1]-1)], rel.eig.cor[(zero.eig[1]+1):n], 0)
87 cum.eig.cor <- cumsum(rel.eig.cor)
88 k2 <- length(which(eig > epsilon))
89 k3 <- length(which(rel.eig.cor > epsilon))
90 vectors <- sweep(D.eig$vectors[,1:k2], 2, sqrt(eig[1:k2]), FUN="*")
91 # Only the eigenvectors with positive eigenvalues are shown
93 # Negative eigenvalues: three ways of handling the situation
94 if((correct==2) | (correct==3)) { # Curly 2
96 if(correct == 2) { # Curly 3
97 # Lingoes correction: compute c1, then the corrected D
99 note <- paste("Lingoes correction applied to negative eigenvalues: D' = -0.5*D^2 -",c1,", except diagonal elements")
100 D <- -0.5*(D^2 + 2*c1)
102 # Cailliez correction: compute c2, then the corrected D
103 } else if(correct == 3) {
104 delta2 <- centre((-0.5*D),n)
105 upper <- cbind(matrix(0,n,n), 2*delta1)
106 lower <- cbind(-diag(n), -4*delta2)
107 sp.matrix <- rbind(upper, lower)
108 c2 <- max(Re(eigen(sp.matrix, symmetric=FALSE, only.values=TRUE)$values))
109 note <- paste("Cailliez correction applied to negative eigenvalues: D' = -0.5*(D +",c2,")^2, except diagonal elements")
114 mat.cor <- centre(D,n)
115 toto.cor <- eigen(mat.cor)
116 trace.cor <- sum(diag(mat.cor))
118 # Negative eigenvalues present?
119 min.eig.cor <- min(toto.cor$values)
120 zero.eig.cor <- which((toto.cor$values < epsilon) & (toto.cor$values > -epsilon))
121 toto.cor$values[zero.eig.cor] <- 0
123 # No negative eigenvalue after correction: result OK
124 if(min.eig.cor > -epsilon) { # Curly 4
125 eig.cor <- toto.cor$values
126 rel.eig.cor <- eig.cor[1:k]/trace.cor
127 cum.eig.cor <- cumsum(rel.eig.cor)
128 k2 <- length(which(eig.cor > epsilon))
129 vectors.cor <- sweep(toto.cor$vectors[,1:k2], 2, sqrt(eig.cor[1:k2]), FUN="*")
130 # bs <- broken.stick(k2)[,2]
132 bs <- c(bs, rep(0,(k-k2)))
135 # Negative eigenvalues still present after correction: incorrect result
137 if(correct == 2) cat("Problem! Negative eigenvalues are still present after Lingoes",'\n')
138 if(correct == 3) cat("Problem! Negative eigenvalues are still present after Cailliez",'\n')
139 rel.eig.cor <- cum.eig.cor <- bs <- cum.bs <- rep(NA,n)
140 vectors.cor <- matrix(NA,n,2)
143 res <- data.frame(eig[1:k], eig.cor[1:k], rel.eig.cor, bs, cum.eig.cor, cum.bs)
144 colnames(res) <- c("Eigenvalues", "Corr_eig", "Rel_corr_eig", "Broken_stick", "Cum_corr_eig", "Cum_br_stick")
145 rownames(res) <- 1:nrow(res)
147 rownames(vectors) <- names
148 colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
149 out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace, vectors.cor=vectors.cor, trace.cor=trace.cor))
153 note <- "No correction was applied to the negative eigenvalues"
155 bs <- c(bs, rep(0,(k-k3)))
158 res <- data.frame(eig[1:k], rel.eig, rel.eig.cor, bs, cum.eig.cor, cum.bs)
159 colnames(res) <- c("Eigenvalues","Relative_eig","Rel_corr_eig","Broken_stick","Cum_corr_eig","Cumul_br_stick")
160 rownames(res) <- 1:nrow(res)
162 rownames(vectors) <- names
163 colnames(vectors) <- colnames(vectors, do.NULL = FALSE, prefix = "Axis.")
164 out <- (list(correction=c(correction,correct), note=note, values=res, vectors=vectors, trace=trace))
166 } # End curly 2: three ways of handling the situation