R/pcoa.R

   1 pcoa <- function(D, correction="none", rn=NULL)
   2 #
   3 # Principal coordinate analysis (PCoA) of a square distance matrix D
   4 # with correction for negative eigenvalues.
   5 #
   6 # References:
   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]
  13 #
  14 #      Pierre Legendre, October 2007
  15 {
  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
  22 }
  23
  24 bstick.def <- function (n, tot.var = 1, ...)   # 'bstick.default' from vegan
  25 {
  26     res <- rev(cumsum(tot.var/n:1)/n)
  27     names(res) <- paste("Stick", seq(len = n), sep = "")
  28     return(res)
  29 }
  30
  31 # ===== The PCoA function begins here =====
  32
  33 # Preliminary actions
  34         D <- as.matrix(D)
  35         n <- nrow(D)
  36         epsilon <- sqrt(.Machine$double.eps)
  37         if(length(rn)!=0) {
  38                 names <- rn
  39                 } else {
  40                 names <- rownames(D)
  41                 }
  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')
  46
  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))
  51
  52 # Eigenvalue decomposition
  53         D.eig <- eigen(delta1)
  54
  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
  59
  60 # No negative eigenvalue
  61         if(min.eig > -epsilon) {   # Curly 1
  62                 correct <- 1
  63                 eig <- D.eig$values
  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="*")
  68                 bs <- bstick.def(k)
  69                 cum.bs <- cumsum(bs)
  70
  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)
  74
  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))
  79
  80 # Negative eigenvalues present
  81         } else {   # Curly 1
  82                 k <- n
  83                 eig <- D.eig$values
  84                 rel.eig <- eig/trace
  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
  92
  93 # Negative eigenvalues: three ways of handling the situation
  94         if((correct==2) | (correct==3)) {   # Curly 2
  95
  96                 if(correct == 2) {   # Curly 3
  97                         # Lingoes correction: compute c1, then the corrected D
  98                         c1 <- -min.eig
  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)
 101
 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")
 110                                 D <- -0.5*(D + c2)^2
 111                         }   # End curly 3
 112
 113         diag(D) <- 0
 114         mat.cor <- centre(D,n)
 115         toto.cor <- eigen(mat.cor)
 116         trace.cor <- sum(diag(mat.cor))
 117
 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
 122
 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]
 131                 bs <- bstick.def(k2)
 132                 bs <- c(bs, rep(0,(k-k2)))
 133                 cum.bs <- cumsum(bs)
 134
 135                 # Negative eigenvalues still present after correction: incorrect result
 136                 } else {
 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)
 141                 }   # End curly 4
 142
 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)
 146
 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))
 150
 151         } else {   # Curly 2
 152
 153         note <- "No correction was applied to the negative eigenvalues"
 154         bs <- bstick.def(k3)
 155         bs <- c(bs, rep(0,(k-k3)))
 156         cum.bs <- cumsum(bs)
 157
 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)
 161
 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))
 165
 166                         }   # End curly 2: three ways of handling the situation
 167                 }       # End curly 1
 168         class(out) <- "pcoa"
 169         out
 170 }       # End of PCoA