2 function(host.D, para.D, HP, nperm=999, test.links=FALSE, seed=NULL, correction="none", silent=FALSE)
4 # Test of host-parasite coevolution
5 # host.D = host distance or patristic matrix (class dist or matrix)
6 # para.D = parasite distance or patristic matrix (class dist or matrix)
7 # HP = host-parasite link matrix (n.host, n.para)
9 # Pierre Legendre, May 2009
11 epsilon <- sqrt(.Machine$double.eps)
14 seed <- .Random.seed[trunc(runif(1,1,626))]
17 host.D <- as.matrix(host.D)
18 host.pc <- pcoa(host.D, correction=correction)
19 if(host.pc$correction[2] == 1) {
20 if(min(host.pc$values[,2]) < -epsilon) stop('Host D matrix has negative eigenvalues. Rerun with correction="lingoes" or correction="cailliez"')
21 sum.host.values.sq <- sum(host.pc$values[,1]^2)
22 host.vectors <- host.pc$vectors
24 sum.host.values.sq <- sum(host.pc$values[,2]^2)
25 host.vectors <- host.pc$vectors.cor
27 n.host <- nrow(host.D)
29 para.D <- as.matrix(para.D)
30 para.pc <- pcoa(para.D, correction=correction)
31 if(para.pc$correction[2] == 1) {
32 if(min(para.pc$values[,2]) < -epsilon) stop('Parasite D matrix has negative eigenvalues. Rerun with correction="lingoes" or correction="cailliez"')
33 sum.para.values.sq <- sum(para.pc$values[,1]^2)
34 para.vectors <- para.pc$vectors
36 sum.para.values.sq <- sum(para.pc$values[,2]^2)
37 para.vectors <- para.pc$vectors.cor
39 n.para <- nrow(para.D)
41 if(!silent) cat("n.hosts =", n.host, ", n.parasites =", n.para,'\n')
44 tracemax <- max(sum.host.values.sq, sum.para.values.sq)
46 if(n.host == n.para) {
47 if(!silent) cat("The function cannot check if matrix HP has been entered in the right way.",'\n')
48 if(!silent) cat("It will assume that the rows of HP are the hosts.",'\n')
51 if(temp[1] == n.host) {
52 if(temp[2] != n.para) stop("Matrices host.D, para.D and HP not comformable")
53 } else if(temp[2] == n.host) {
54 if(temp[1] != n.para) stop("Matrices host.D, para.D and HP not comformable")
56 if(!silent) cat("Matrix HP has been transposed for comformity with host.D and para.D.",'\n')
58 stop("Matrices host.D, para.D and HP not comformable")
61 p.per.h <- apply(HP, 1, sum)
62 h.per.p <- apply(HP, 2, sum)
64 # Compute and test the global statistics
65 mat.4 <- t(host.vectors) %*% HP %*% para.vectors
66 global <- sum(mat.4^2)
72 HP.perm <- apply(HP, 2, sample)
73 mat.4.perm <- t(host.vectors) %*% HP.perm %*% para.vectors
74 global.perm <- c(global.perm, sum(mat.4.perm^2))
75 if(global.perm[i+1] >= global) nGT <- nGT+1
77 global.perm <- global.perm[-1]
78 p.global <- nGT/(nperm+1)
79 } else { p.global <- NA }
82 # Test individual H-P links
84 # 1. Create the list of H-P pairs
85 list.hp <- which( t(cbind(HP,rep(0,n.host))) > 0)
86 HP.list <- cbind((list.hp %/% (n.para+1))+1, list.hp %% (n.para+1))
87 colnames(HP.list) <- c("Host","Parasite")
88 n.links <- length(list.hp)
96 # 2. Compute reference values of link statistics
98 HP.k[HP.list[k,1], HP.list[k,2]] <- 0
99 mat.4.k <- t(host.vectors) %*% HP.k %*% para.vectors
100 trace.k <- sum(mat.4.k^2)
101 stat1 <- c(stat1, (global-trace.k))
102 den <- tracemax-global
104 stat2 <- c(stat2, stat1[k+1]/den)
106 stat2 <- c(stat2, NA)
109 # 3. Test link statistics by permutations
117 HP.k.perm <- apply(HP.k, 2, sample)
118 mat.4.k.perm <- t(host.vectors) %*% HP.k.perm %*% para.vectors
119 trace.k.perm <- sum(mat.4.k.perm^2)
120 stat1.perm <- global.perm[i]-trace.k.perm
121 if(stat1.perm >= stat1[k+1]) nGT1 <- nGT1+1
123 if(!is.na(stat2[k+1])) {
124 den <- tracemax-global.perm[i]
126 stat2.perm <- stat1.perm/den
127 if(stat2.perm >= stat2[k+1]) nGT2 <- nGT2+1
130 # if(!silent) cat("In permutation #",i,"den < epsilon",'\n')
134 p.stat1 <- c(p.stat1, nGT1/(nperm+1))
135 if(!is.na(stat2[k+1])) {
136 p.stat2 <- c(p.stat2, nGT2/(nperm2+1))
146 link.table <- cbind(HP.list, stat1[-1], p.stat1[-1], stat2[-1], p.stat2[-1])
147 colnames(link.table) = c("Host","Parasite","F1.stat","p.F1","F2.stat","p.F2")
148 out <-list(ParaFitGlobal=global, p.global=p.global, link.table=link.table, para.per.host=p.per.h, host.per.para=h.per.p, nperm=nperm)
150 if(!silent) cat("Rerun the program with option 'test.links=TRUE' to test the individual H-P links",'\n')
151 out <-list(ParaFitGlobal=global, p.global=p.global, para.per.host=p.per.h, host.per.para=h.per.p, nperm=nperm)
155 a[3] <- sprintf("%2f",a[3])
156 if(!silent) cat("Computation time =",a[3]," sec",'\n')
158 class(out) <- "parafit"