2 function(Dmat, nmat, n, nperm=99, make.sym=TRUE, weights=NULL, mult="holm", mantel=FALSE, silent=FALSE)
4 ### Function to carry out a posteriori tests of the contribution of individual
5 ### matrices to the congruence of a group of distance matrices.
7 ### copyleft - Pierre Legendre, December 2008
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.
14 ### Parameters of the function --
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
21 ### nmat = number of distance matrices in file Dmat.
23 ### n = number of objects in each distance matrix. All matrices have same n.
25 ### nperm = number of permutations for the tests.
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.
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.
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.
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.
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.
48 ################################################################################
50 mult <- match.arg(mult, c("sidak", "holm", "bonferroni"))
52 stop("Analysis requested for a single D matrix: CADM is useless")
56 ## Check the input file
58 stop("Error in the value of 'n' or in the D matrices themselves")
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")
64 if(is.null(weights)) {
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')
75 ## Are asymmetric D matrices present?
76 asy <- rep(FALSE, nmat)
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')
96 if(!silent) cat("\nAsymmetric matrices were transformed to be symmetric",'\n')
99 if(!silent) cat("\nAnalysis carried out on asymmetric matrices",'\n')
100 D2 <- as.list(1:nmat)
103 if(!silent) cat("Analysis of symmetric matrices",'\n')
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
113 D.temp <- as.matrix(Dmat[begin:end,])
114 vec <- as.vector(as.dist(D.temp))
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))))
122 D.temp2 <- rank(D.temp)
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)))
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)
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)
146 Y <- as.matrix(Y[,-1])
147 colnames(Y) <- colnames(Y,do.NULL = FALSE, prefix = "Dmat.")
149 ## Begin calculations: compute reference value of S
151 ## Transform the distances to ranks, by column
152 Rmat <- apply(Y,2,rank)
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)
159 ## Begin a posteriori tests of individual matrices
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
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)
174 ## Test each matrix 'k' in turn
176 ## Create a new Rmat table where the permuted column has been removed
177 Rmat.mod <- Rmat[,-k]
179 ## Permutation loop: string out permuted matrix 'k' only
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])))
188 vec <- as.vector(as.dist(D1[[k]][order,order]))
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
196 ## Calculate P-values
197 counter <- counter/(nperm+1)
199 ## Correction to P-values for multiple testing
200 if(mult == "sidak") {
202 for(i in 1:p) vec.corr = c(vec.corr, (1-(1-counter[i])^p))
203 vec.corr <- vec.corr[-1]
205 if(mult == "holm") vec.corr <- p.adjust(counter, method="holm")
206 if(mult == "bonferroni") vec.corr <- p.adjust(counter, method="bonferroni")
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.")
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)
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)
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)) {
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)
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)
241 for(k in 1:(nmat-1)) {
243 vec <- as.vector(as.dist(D1[[k]][order,order]))
244 Rmat.perm <- cbind(Rmat.perm, vec)
246 vec <- as.vector(as.dist(D1[[nmat]]))
247 Rmat.perm <- cbind(Rmat.perm, vec)
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
259 Mantel.prob <- as.matrix(as.dist(Mantel.prob/(nperm+1)))
260 diag(Mantel.prob) <- 1
264 a[3] <- sprintf("%2f",a[3])
265 if(!silent) cat("Time to compute a posteriori tests (per matrix) =",a[3]," sec",'\n')
267 out <- list(A_posteriori_tests=table, Correction.type=mult)
270 out$Mantel.cor <- Mantel.cor
271 out$Mantel.prob <- Mantel.prob
274 class(out) <- "CADM.post"