CADM.post.R

   1 `CADM.post` <-
   2         function(Dmat, nmat, n, nperm=99, make.sym=TRUE, weights=NULL, mult="holm", mantel=FALSE, silent=FALSE)
   3 {
   4 ### Function to carry out a posteriori tests of the contribution of individual
   5 ### matrices to the congruence of a group of distance matrices.
   6 ###
   7 ### copyleft - Pierre Legendre, December 2008
   8 ###
   9 ### Reference -
  10 ### Legendre, P. and F.-J. Lapointe. 2004. Assessing congruence among distance
  11 ### matrices: single malt Scotch whiskies revisited. Australian and New Zealand
  12 ### Journal of Statistics 46: 615-629.
  13 ###
  14 ### Parameters of the function --
  15 ###
  16 ### Dmat = A text file listing the distance matrices one after the other, with
  17 ###        or without blank lines.
  18 ###        Each matrix is in the form of a square distance matrix with 0's
  19 ###        on the diagonal.
  20 ###
  21 ### nmat = number of distance matrices in file Dmat.
  22 ###
  23 ### n = number of objects in each distance matrix. All matrices have same n.
  24 ###
  25 ### nperm = number of permutations for the tests.
  26 ###
  27 ### make.sym = TRUE: turn asymmetric matrices into symmetric matrices by
  28 ###            averaging the two triangular portions.
  29 ###          = FALSE: analyse asymmetric matrices as they are.
  30 ###
  31 ### weights = a vector of positive weights for the distance matrices.
  32 ###           Example: weights = c(1,2,3)
  33 ###         = NULL (default): all matrices have same weight in calculation of W.
  34 ###
  35 ### mult = method for correcting P-values due to multiple testing. The methods
  36 ###        are "holm" (default), "sidak", and "bonferroni". The Bonferroni
  37 ###        correction is overly conservative; it is not recommended. It is
  38 ###        included to allow comparisons with the other methods.
  39 ###
  40 ### mantel = TRUE: Mantel statistics are computed from ranked distances,
  41 ###          as well as permutational P-values.
  42 ###        = FALSE (default): Mantel statistics and tests are not computed.
  43 ###
  44 ### silent = TRUE: informative messages will not be printed, except stopping
  45 ###          messages. Option useful for simulation work.
  46 ###        = FALSE: informative messages will be printed.
  47 ###
  48 ################################################################################
  49
  50         mult <- match.arg(mult, c("sidak", "holm", "bonferroni"))
  51         if(nmat < 2)
  52                 stop("Analysis requested for a single D matrix: CADM is useless")
  53
  54         a <- system.time({
  55
  56     ## Check the input file
  57     if(ncol(Dmat) != n)
  58         stop("Error in the value of 'n' or in the D matrices themselves")
  59     nmat2 <- nrow(Dmat)/n
  60     if(nmat2 < nmat)  # OK if 'nmat' < number of matrices in the input file
  61         stop("Number of input D matrices = ",nmat2,"; this value is < nmat")
  62
  63     nd <- n*(n-1)/2
  64     if(is.null(weights)) {
  65         w <- rep(1,nmat)
  66         } else {
  67         if(length(weights) != nmat)
  68                 stop("Incorrect number of values in vector 'weights'")
  69         if(length(which(weights < 0)) > 0)
  70                 stop("Negative weights are not permitted")
  71         w <- weights*nmat/sum(weights)
  72         if(!silent) cat("Normalized weights =",w,'\n')
  73         }
  74
  75     ## Are asymmetric D matrices present?
  76     asy <- rep(FALSE, nmat)
  77         asymm <- FALSE
  78     end <- 0
  79     for(k in 1:nmat) {
  80         begin <- end+1
  81         end <- end+n
  82         D.temp <- Dmat[begin:end,]
  83         if(sum(abs(diag(as.matrix(D.temp)))) > 0)
  84                 stop("Diagonal not 0: matrix #",k," is not a distance matrix")
  85         vec1 <- as.vector(as.dist(D.temp))
  86         vec2 <- as.vector(as.dist(t(D.temp)))
  87         if(sum(abs((vec1-vec2))) > 0) {
  88                 if(!silent) cat("Matrix #",k," is asymmetric",'\n')
  89                 asy[k] <- TRUE
  90                 asymm <- TRUE
  91                 }
  92         }
  93     D1 <- as.list(1:nmat)
  94     if(asymm) {
  95         if(make.sym) {
  96                 if(!silent) cat("\nAsymmetric matrices were transformed to be symmetric",'\n')
  97                 } else {
  98                 nd <- nd*2
  99                 if(!silent) cat("\nAnalysis carried out on asymmetric matrices",'\n')
 100                 D2 <- as.list(1:nmat)
 101                 }
 102         } else {
 103         if(!silent) cat("Analysis of symmetric matrices",'\n')
 104         }
 105     Y <- rep(NA,nd)
 106
 107     ## String out the distance matrices (vec) and assemble them as columns into matrix 'Y'
 108     ## Construct also matrices of ranked distances D1[[k]] and D2[[k]] for permutation test
 109     end <- 0
 110     for(k in 1:nmat) {
 111         begin <- end+1
 112         end <- end+n
 113         D.temp <- as.matrix(Dmat[begin:end,])
 114         vec <- as.vector(as.dist(D.temp))
 115         if(asymm) {
 116                 if(!make.sym) {
 117                         ## Analysis carried out on asymmetric matrices:
 118                         ## The ranks are computed on the whole matrix except the diagonal values.
 119                         ## The two halves are stored as symmetric matrices in D1[[k]] and D2[[k]]
 120                         vec <- c(vec, as.vector(as.dist(t(D.temp))))
 121                         diag(D.temp) <- NA
 122                         D.temp2 <- rank(D.temp)
 123                         diag(D.temp2) <- 0
 124                         # cat("nrow =",nrow(D.temp2)," ncol =",ncol(D.temp2),'\n')
 125                                 # cat("Matrix ",k," min =",min(D.temp2)," max =",max(D.temp2),'\n')
 126                                 # cat("Matrix ",k," max values #",which(D.temp2 == max(D.temp2)),'\n')
 127                         D1[[k]] <- as.matrix(as.dist(D.temp2))
 128                         D2[[k]] <- as.matrix(as.dist(t(D.temp2)))
 129                         } else {
 130                         ## Asymmetric matrices transformed to be symmetric, stored in D1[[k]]
 131                         vec <- (vec + as.vector(as.dist(t(D.temp)))) / 2
 132                                 D.temp2 <- (D.temp + t(D.temp)) / 2
 133                                 D.temp2 <- as.dist(D.temp2)
 134                         D.temp2[] <- rank(D.temp2)
 135                                 D.temp2 <- as.matrix(D.temp2)
 136                         D1[[k]] <- D.temp2
 137                         }
 138                 } else {
 139                 ## Symmetric matrices are stored in D1[[k]]
 140                 D.temp2 <- as.dist(D.temp)
 141                 D.temp2[] <- rank(D.temp2)
 142                 D1[[k]] <- as.matrix(D.temp2)
 143                 }
 144         Y <- cbind(Y, vec)
 145         }
 146     Y <- as.matrix(Y[,-1])
 147     colnames(Y) <- colnames(Y,do.NULL = FALSE, prefix = "Dmat.")
 148
 149     ## Begin calculations: compute reference value of S
 150
 151                 ## Transform the distances to ranks, by column
 152                 Rmat <- apply(Y,2,rank)
 153
 154                 ## Compute the S = Sum-of-Squares of the row-marginal sums of ranks (eq. 1a)
 155                 ## The ranks are weighted during the sum by the vector of matrix weights 'w'
 156                 sumRanks <- as.vector(Rmat%*%w)
 157                 S <- (nd-1)*var(sumRanks)
 158
 159     ## Begin a posteriori tests of individual matrices
 160
 161     ## Statistics displayed for each matrix: "Mantel.mean" and "W.per.matrix"
 162         ## Calculate the mean of the Mantel correlations on ranks for each matrix
 163         Mantel.cor <- cor(Rmat)
 164         diag(Mantel.cor) <- 0
 165         spear.mean <- as.vector(Mantel.cor%*%w)/(nmat-1)
 166         ## Calculate Kendall's W for each variable
 167         ## W.var <- ((nmat-1)*spear.mean+1)/nmat
 168
 169         ## P-value for each matrix: test of S, permuting values in matrix[[k]] only
 170         ## as in program CADM.f (2004)
 171         ## Initialize the counters
 172         counter <- rep(1,nmat)
 173
 174         ## Test each matrix 'k' in turn
 175         for(k in 1:nmat) {
 176                 ## Create a new Rmat table where the permuted column has been removed
 177                 Rmat.mod <- Rmat[,-k]
 178
 179                 ## Permutation loop: string out permuted matrix 'k' only
 180                 for(j in 1:nperm) {
 181                         order <- sample(n)
 182                         if(asymm & !make.sym) {
 183                                 ## For asymmetric matrices: permute the values within each triangular
 184                                 ## portion, stored as square matrices in D1[[]] and D2[[]]
 185                                 vec <- as.vector(as.dist(D1[[k]][order,order]))
 186                             vec <- c(vec, as.vector(as.dist(D2[[k]][order,order])))
 187                                 } else {
 188                                 vec <- as.vector(as.dist(D1[[k]][order,order]))
 189                                 }
 190                         Rmat.perm <- cbind(Rmat.mod, vec)
 191                         S.perm <- (nd-1)*var(as.vector(Rmat.perm%*%w))
 192                         if(S.perm >= S) counter[k] <- counter[k]+1
 193                 }
 194         }
 195
 196         ## Calculate P-values
 197         counter <- counter/(nperm+1)
 198
 199         ## Correction to P-values for multiple testing
 200                 if(mult == "sidak") {
 201                         vec.corr = NA
 202                         for(i in 1:p) vec.corr = c(vec.corr, (1-(1-counter[i])^p))
 203                         vec.corr <- vec.corr[-1]
 204                         }
 205                 if(mult == "holm") vec.corr <- p.adjust(counter, method="holm")
 206                 if(mult == "bonferroni") vec.corr <- p.adjust(counter, method="bonferroni")
 207
 208         ## Create a data frame containing the results
 209                 # table <- rbind(spear.mean, W.var, counter, vec.corr)
 210                 # rownames(table) <- c("Mantel.mean", "W.per.matrix", "Prob", "Corrected prob")
 211                 table <- rbind(spear.mean, counter, vec.corr)
 212                 rownames(table) <- c("Mantel.mean", "Prob", "Corrected.prob")
 213                 colnames(table) <- colnames(table,do.NULL = FALSE, prefix = "Dmat.")
 214
 215         ## Mantel tests
 216         if(mantel) {
 217                 diag(Mantel.cor) <- 1
 218                 rownames(Mantel.cor) <- colnames(table)
 219                 colnames(Mantel.cor) <- colnames(table)
 220                 Mantel.prob <- matrix(1,nmat,nmat)
 221                 rownames(Mantel.prob) <- colnames(table)
 222                 colnames(Mantel.prob) <- colnames(table)
 223
 224                 for(j in 1:nperm) {   # Each matrix is permuted independently
 225                                   # There is no need to permute the last matrix
 226                         Rmat.perm <- rep(NA,nd)
 227                         ##
 228                         if(asymm & !make.sym) {
 229                                 ## For asymmetric matrices: permute the values within each triangular
 230                                 ## portion, stored as square matrices in D1[[]] and D2[[]]
 231                                 for(k in 1:(nmat-1)) {
 232                                         order <- sample(n)
 233                                         vec <- as.vector(as.dist(D1[[k]][order,order]))
 234                                 vec <- c(vec, as.vector(as.dist(D2[[k]][order,order])))
 235                                         Rmat.perm <- cbind(Rmat.perm, vec)
 236                                         }
 237                                         vec <- as.vector(as.dist(D1[[nmat]]))
 238                                 vec <- c(vec, as.vector(as.dist(D2[[nmat]])))
 239                                         Rmat.perm <- cbind(Rmat.perm, vec)
 240                                 } else {
 241                                 for(k in 1:(nmat-1)) {
 242                                         order <- sample(n)
 243                                         vec <- as.vector(as.dist(D1[[k]][order,order]))
 244                                         Rmat.perm <- cbind(Rmat.perm, vec)
 245                                         }
 246                                         vec <- as.vector(as.dist(D1[[nmat]]))
 247                                         Rmat.perm <- cbind(Rmat.perm, vec)
 248                                 }
 249                         # Remove the first column of Rmat.perm containing NA
 250                         Rmat.perm <- as.matrix(Rmat.perm[,-1])
 251                         # Compute Mantel correlations on ranks under permutation
 252                         Mantel.cor.perm <- cor(Rmat.perm)
 253                         for(j2 in 1:(nmat-1)) { # Compute prob in the upper tail
 254                                 for(j1 in (j2+1):nmat) {
 255                                         if(Mantel.cor.perm[j1,j2] >= Mantel.cor[j1,j2]) Mantel.prob[j1,j2] <- Mantel.prob[j1,j2]+1
 256                                         }
 257                                 }
 258                         }
 259                 Mantel.prob <- as.matrix(as.dist(Mantel.prob/(nperm+1)))
 260                 diag(Mantel.prob) <- NA # Corrected 08feb13
 261                 }
 262
 263         })
 264         a[3] <- sprintf("%2f",a[3])
 265         if(!silent) cat("Time to compute a posteriori tests (per matrix) =",a[3]," sec",'\n')
 266
 267         out <- list(A_posteriori_tests=table, Correction.type=mult)
 268
 269         if(mantel) {
 270                 out$Mantel.cor  <- Mantel.cor
 271                 out$Mantel.prob <- Mantel.prob
 272                 }
 273         out$nperm <- nperm
 274         class(out) <- "CADM.post"
 275         out
 276 }