Merge ae2c945 into 1e5a325

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: paco
-Version: 0.3.1
-Date: 2016-11-28
+Version: 0.4.0
+Date: 2019-05-12
 Title: Procrustes Application to Cophylogenetic Analysis
 Description: Procrustes analyses to infer co-phylogenetic
     matching between pairs of (ultrametric) phylogenetic trees.
@@ -12,11 +12,11 @@ Depends:
     R (>= 3.0.0)
 Imports:
     vegan (>= 2.2-0),
+    ape,
     plyr
 Suggests:
-    ape,
-    doMC,
     testthat
 License: MIT + file LICENSE
 URL: http://www.uv.es/cophylpaco/
-RoxygenNote: 5.0.1
+Encoding: UTF-8
+RoxygenNote: 6.1.1

Makefile

@@ -16,7 +16,7 @@ test: cran/paco doc
 
 paco.tar.gz:
 	rm $@ 2>/dev/null; true
-	cd cran; tar -zcvf paco.tar.gz paco
+	cd cran; R CMD build paco
 	mv cran/paco.tar.gz $@
 
 clean:

R/add_pcoord.r

@@ -3,8 +3,8 @@
 #' Translates the distance matrices of 'host' and 'parasite' phylogenies into Principal Coordinates, as needed for Procrustes superimposition.
 #' @param D A list with objects H, P, and HP, as returned by \code{paco::prepare_paco_data}.
 #' @param correction In some cases, phylogenetic distance matrices are non-Euclidean which generates negative eigenvalues when those matrices are translated into Principal Coordinates. There are several methods to correct negative eigenvalues. Correction options available here are "cailliez", "lingoes", and "none". The "cailliez" and "lingoes" corrections add a constant to the eigenvalues to make them non-negative. Default is "none". 
-#' @return The list that was input as the argument `D' with two new elements; the Principal Coordinates of the `host' distance matrix and the Principal Coordinates of the `parasite' distance matrix.
-#' @note To find the Principal Coordinates of each distance matrix, we internally use a function, \code{coordpcoa}, that is a modified version of \code{ape::pcoa}, using \code{vegan::eigenvals}.
+#' @return The list that was input as the argument `D' with four new elements; the Principal Coordinates of the `host' distance matrix and the Principal Coordinates of the `parasite' distance matrix, as well as, a `correction' object stating the correction used for negative eigenvalues and a `note' object stating whether or not negative eigenvalues were present and therefore corrected.
+#' @note To find the Principal Coordinates of each distance matrix, we internally a modified version of the function \code{ape::pcoa} that uses \code{vegan::eigenvals} and zapsmall
 #' @export
 #' @examples
 #' data(gopherlice)
@@ -16,10 +16,34 @@
 add_pcoord <- function(D, correction='none')
 {
 HP_bin <- which(D$HP >0, arr.ind=TRUE)
-H_PCo <- coordpcoa(D$H, correction)$vectors #Performs PCo of Host distances
-P_PCo <- coordpcoa(D$P, correction)$vectors #Performs PCo of Parasite distances
-D$H_PCo <- H_PCo[HP_bin[,1],] #Adjust Host PCo vectors
-D$P_PCo <- P_PCo[HP_bin[,2],] #Adjust Parasite PCo vectors
+if(correction == 'none'){
+	H_PCo <- coordpcoa(D$H, correction) #Performs PCo of Host distances
+	P_PCo <- coordpcoa(D$P, correction) #Performs PCo of Parasite distances
+	H_PCoVec <- H_PCo$vectors
+	P_PCoVec <- P_PCo$vectors
+	D$H_PCo <- H_PCoVec[HP_bin[,1],] #Adjust Host PCo vectors
+	D$P_PCo <- P_PCoVec[HP_bin[,2],] #Adjust Parasite PCo vectors
+}else{
+	H_PCo <- coordpcoa(D$H, correction)
+	P_PCo <- coordpcoa(D$P, correction)
+	if(H_PCo$note == "There were no negative eigenvalues. No correction was applied"){
+		H_PCoVec <- H_PCo$vectors
+		D$H_PCo <- H_PCoVec[HP_bin[,1],] #Adjust Host PCo vectors
+	}else{
+		H_PCoVec <- H_PCo$vectors.cor
+		D$H_PCo <- H_PCoVec[HP_bin[,1],] #Adjust Host PCo vectors
+	}
+
+	if(P_PCo$note == "There were no negative eigenvalues. No correction was applied"){
+		P_PCoVec <- P_PCo$vectors
+		D$P_PCo <- P_PCoVec[HP_bin[,2],] #Adjust Host PCo vectors
+	}else{
+		P_PCoVec <- P_PCo$vectors.cor
+		D$P_PCo <- P_PCoVec[HP_bin[,2],] #Adjust Host PCo vectors
+	}
+
+}
 D$correction <- correction
+D$note <- list(H_note=H_PCo$note, P_note=P_PCo$note)
 return(D)
 }

R/coordpcoa.r

@@ -2,8 +2,8 @@
 # D A list with objects H, P, and HP, returned by prepare_paco_data
 #  correction The correction to apply (none, lingoes, or cailliez)
 #  rn rownames (optional)
-# Internal function coordpcoa is a modified version of ape::pcoa, utilising vegan::eigenvals
-coordpcoa <-function (D, correction = "none", rn = NULL)
+# Internal function coordpcoa is a modified version of ape::pcoa, utilising vegan::eigenvals and zapsmall
+coordpcoa <- function (D, correction = "none", rn = NULL) 
 {
     centre <- function(D, n) {
         One <- matrix(1, n, n)
@@ -20,55 +20,57 @@ coordpcoa <-function (D, correction = "none", rn = NULL)
     epsilon <- sqrt(.Machine$double.eps)
     if (length(rn) != 0) {
         names <- rn
-    }
-    else {
+    }else {
         names <- rownames(D)
     }
     CORRECTIONS <- c("none", "lingoes", "cailliez")
     correct <- pmatch(correction, CORRECTIONS)
-    if (is.na(correct)) stop("Invalid correction method")
+    if (is.na(correct)) 
+        stop("Invalid correction method")
     delta1 <- centre((-0.5 * D^2), n)
     trace <- sum(diag(delta1))
     D.eig <- eigen(delta1)
+    D.eig$values <- as.numeric(zapsmall(vegan::eigenvals(D.eig)))
     min.eig <- min(D.eig$values)
-    D.eig$values <- vegan::eigenvals(D.eig)
-    zero.eig <- which(D.eig$values < epsilon)
+    zero.eig <- which(abs(D.eig$values) < epsilon)
+    D.eig$values[zero.eig] <- 0
     if (min.eig > -epsilon) {
         correct <- 1
         eig <- D.eig$values
         k <- length(which(eig > epsilon))
         rel.eig <- eig[1:k]/trace
         cum.eig <- cumsum(rel.eig)
-        vectors <- sweep(D.eig$vectors[, 1:k], 2, sqrt(eig[1:k]), FUN = "*")
+        vectors <- sweep(D.eig$vectors[, 1:k], 2, sqrt(eig[1:k]), 
+            FUN = "*")
         bs <- bstick.def(k)
         cum.bs <- cumsum(bs)
         res <- data.frame(eig[1:k], rel.eig, bs, cum.eig, cum.bs)
-        colnames(res) <- c("Eigenvalues", "Relative_eig", "Broken_stick",
+        colnames(res) <- c("Eigenvalues", "Relative_eig", "Broken_stick", 
             "Cumul_eig", "Cumul_br_stick")
         rownames(res) <- 1:nrow(res)
         rownames(vectors) <- names
-        colnames(vectors) <- colnames(vectors, do.NULL = FALSE,
+        colnames(vectors) <- colnames(vectors, do.NULL = FALSE, 
             prefix = "Axis.")
         note <- paste("There were no negative eigenvalues. No correction was applied")
-        out <- (list(correction = c(correction, correct), note = note,
+        out <- (list(correction = c(correction, correct), note = note, 
             values = res, vectors = vectors, trace = trace))
     }
     else {
         k <- n
         eig <- D.eig$values
         rel.eig <- eig/trace
         rel.eig.cor <- (eig - min.eig)/(trace - (n - 1) * min.eig)
-        rel.eig.cor = c(rel.eig.cor[1:(zero.eig[1] - 1)], rel.eig.cor[(zero.eig[1] +
+        rel.eig.cor = c(rel.eig.cor[1:(zero.eig[1] - 1)], rel.eig.cor[(zero.eig[1] + 
             1):n], 0)
         cum.eig.cor <- cumsum(rel.eig.cor)
         k2 <- length(which(eig > epsilon))
         k3 <- length(which(rel.eig.cor > epsilon))
-        vectors <- sweep(D.eig$vectors[, 1:k2], 2, sqrt(eig[1:k2]),
+        vectors <- sweep(D.eig$vectors[, 1:k2], 2, sqrt(eig[1:k2]), 
             FUN = "*")
         if ((correct == 2) | (correct == 3)) {
             if (correct == 2) {
                 c1 <- -min.eig
-                note <- paste("Lingoes correction applied to negative eigenvalues: D' = -0.5*D^2 -",
+                note <- paste("Lingoes correction applied to negative eigenvalues: D' = -0.5*D^2 -", 
                   c1, ", except diagonal elements")
                 D <- -0.5 * (D^2 + 2 * c1)
             }
@@ -77,69 +79,80 @@ coordpcoa <-function (D, correction = "none", rn = NULL)
                 upper <- cbind(matrix(0, n, n), 2 * delta1)
                 lower <- cbind(-diag(n), -4 * delta2)
                 sp.matrix <- rbind(upper, lower)
-                c2 <- max(Re(eigen(sp.matrix, symmetric = FALSE,
+                c2 <- max(Re(eigen(sp.matrix, symmetric = FALSE, 
                   only.values = TRUE)$values))
-                note <- paste("Cailliez correction applied to negative eigenvalues: D' = -0.5*(D +", c2, ")^2, except diagonal elements")
+                note <- paste("Cailliez correction applied to negative eigenvalues: D' = -0.5*(D +", 
+                  c2, ")^2, except diagonal elements")
                 D <- -0.5 * (D + c2)^2
             }
             diag(D) <- 0
             mat.cor <- centre(D, n)
             toto.cor <- eigen(mat.cor)
+            toto.cor$values <- as.numeric(zapsmall(vegan::eigenvals(toto.cor)))
             trace.cor <- sum(diag(mat.cor))
             min.eig.cor <- min(toto.cor$values)
-            toto.cor$values <- vegan::eigenvals(toto.cor)
-            zero.eig.cor <- which((toto.cor$values < epsilon) & (toto.cor$values > -epsilon))
-
+            zero.eig.cor <- which((toto.cor$values < epsilon) & 
+                (toto.cor$values > -epsilon))
+            toto.cor$values[zero.eig.cor] <- 0
             if (min.eig.cor > -epsilon) {
                 eig.cor <- toto.cor$values
                 rel.eig.cor <- eig.cor[1:k]/trace.cor
                 cum.eig.cor <- cumsum(rel.eig.cor)
                 k2 <- length(which(eig.cor > epsilon))
-                vectors.cor <- sweep(toto.cor$vectors[, 1:k2],
+                vectors.cor <- sweep(toto.cor$vectors[, 1:k2], 
                   2, sqrt(eig.cor[1:k2]), FUN = "*")
+                rownames(vectors.cor) <- names
+                colnames(vectors.cor) <- colnames(vectors.cor, 
+                  do.NULL = FALSE, prefix = "Axis.")
                 bs <- bstick.def(k2)
                 bs <- c(bs, rep(0, (k - k2)))
                 cum.bs <- cumsum(bs)
             }
             else {
-                if (correct == 2)
-                  cat("Problem! Negative eigenvalues are still present after Lingoes","\n")
-                if (correct == 3)
-                  cat("Problem! Negative eigenvalues are still present after Cailliez","\n")
-                rel.eig.cor <- cum.eig.cor <- bs <- cum.bs <- rep(NA, n)
+                if (correct == 2) 
+                  cat("Problem! Negative eigenvalues are still present after Lingoes", 
+                    "\n")
+                if (correct == 3) 
+                  cat("Problem! Negative eigenvalues are still present after Cailliez", 
+                    "\n")
+                rel.eig.cor <- cum.eig.cor <- bs <- cum.bs <- rep(NA, 
+                  n)
                 vectors.cor <- matrix(NA, n, 2)
+                rownames(vectors.cor) <- names
+                colnames(vectors.cor) <- colnames(vectors.cor, 
+                  do.NULL = FALSE, prefix = "Axis.")
             }
-            res <- data.frame(eig[1:k], eig.cor[1:k], rel.eig.cor,
+            res <- data.frame(eig[1:k], eig.cor[1:k], rel.eig.cor, 
                 bs, cum.eig.cor, cum.bs)
-            colnames(res) <- c("Eigenvalues", "Corr_eig", "Rel_corr_eig",
+            colnames(res) <- c("Eigenvalues", "Corr_eig", "Rel_corr_eig", 
                 "Broken_stick", "Cum_corr_eig", "Cum_br_stick")
             rownames(res) <- 1:nrow(res)
             rownames(vectors) <- names
-            colnames(vectors) <- colnames(vectors, do.NULL = FALSE,
+            colnames(vectors) <- colnames(vectors, do.NULL = FALSE, 
                 prefix = "Axis.")
-            out <- (list(correction = c(correction, correct),
-                note = note, values = res, vectors = vectors,
+            out <- (list(correction = c(correction, correct), 
+                note = note, values = res, vectors = vectors, 
                 trace = trace, vectors.cor = vectors.cor, trace.cor = trace.cor))
         }
         else {
             note <- "No correction was applied to the negative eigenvalues"
             bs <- bstick.def(k3)
             bs <- c(bs, rep(0, (k - k3)))
             cum.bs <- cumsum(bs)
-            res <- data.frame(eig[1:k], rel.eig, rel.eig.cor,
+            res <- data.frame(eig[1:k], rel.eig, rel.eig.cor, 
                 bs, cum.eig.cor, cum.bs)
-            colnames(res) <- c("Eigenvalues", "Relative_eig",
-                "Rel_corr_eig", "Broken_stick", "Cum_corr_eig",
+            colnames(res) <- c("Eigenvalues", "Relative_eig", 
+                "Rel_corr_eig", "Broken_stick", "Cum_corr_eig", 
                 "Cumul_br_stick")
             rownames(res) <- 1:nrow(res)
             rownames(vectors) <- names
-            colnames(vectors) <- colnames(vectors, do.NULL = FALSE,
+            colnames(vectors) <- colnames(vectors, do.NULL = FALSE, 
                 prefix = "Axis.")
-            out <- (list(correction = c(correction, correct),
-                note = note, values = res, vectors = vectors,
+            out <- (list(correction = c(correction, correct), 
+                note = note, values = res, vectors = vectors, 
                 trace = trace))
         }
     }
     class(out) <- "pcoa"
     out
-}
+}

R/paco_links.r

@@ -16,38 +16,41 @@
 #' D <- PACo(D, nperm=10, seed=42, method="r0")
 #' D <- paco_links(D)
 
-paco_links <- function(D, .parallel = FALSE, proc.warnings=TRUE)
-{
-   correction <- D$correction
-   HP.ones <- which(D$HP > 0, arr.ind=TRUE)
-   SQres.jackn <- matrix(rep(NA, sum(D$HP)^2), sum(D$HP))# empty matrix of jackknifed squared residuals
-   colnames(SQres.jackn) <- paste(rownames(D$proc$X),rownames(D$proc$Yrot), sep="-") #colnames identify the H-P link
-   t.critical = stats::qt(0.975,sum(D$HP)-1) #Needed to compute 95% confidence intervals.
-   nlinks <- sum(D$HP)
+paco_links <- function(D, .parallel = FALSE, proc.warnings=TRUE){
+
+  correction <- D$correction
+  HP.ones <- which(D$HP > 0, arr.ind=TRUE)
+  SQres.jackn <- matrix(rep(NA, sum(D$HP)^2), sum(D$HP))# empty matrix of jackknifed squared residuals
+  colnames(SQres.jackn) <- paste(rownames(D$proc$X),rownames(D$proc$Yrot), sep="-") #colnames identify the H-P link
+  t.critical = stats::qt(0.975,sum(D$HP)-1) #Needed to compute 95% confidence intervals.
+  nlinks <- sum(D$HP)
 
   #if .parallel is TRUE
   if(.parallel==TRUE){
-  SQres.jackn <- plyr::adply(1:nlinks, 1, function(x) single_paco_link(D, HP.ones, x, correction, proc.warnings), .parallel=.parallel)
-} else{
-  #if .parallel is FALSE
-  for(i in c(1:nlinks))
-  {
-  res.Proc.ind <- single_paco_link (D, HP.ones, i, correction, proc.warnings)
-  SQres.jackn[i, ] <- res.Proc.ind
+    SQres.jackn <- plyr::adply(1:nlinks, 1, 
+      function(i) single_paco_link(D, HP.ones, i, correction, proc.warnings), 
+      .parallel=.parallel)
+    SQres.jackn <- SQres.jackn[,-1]
+    colnames(SQres.jackn)[1] <- paste(rownames(D$proc$X),rownames(D$proc$Yrot), sep="-")[1]
+  }else{
+    #if .parallel is FALSE
+    for(i in c(1:nlinks)){
+      res.Proc.ind <- single_paco_link (D, HP.ones, i, correction, proc.warnings)
+      SQres.jackn[i, ] <- res.Proc.ind
+    }
   }
-}
 
-   SQres.jackn <- SQres.jackn^2 #Jackknifed residuals are squared
-   SQres <- (stats::residuals(D$proc))^2 # Vector of original square residuals
-   #jackknife calculations:
-   SQres.jackn <- SQres.jackn*(-(sum(D$HP)-1))
-   SQres <- SQres*sum(D$HP)
-   SQres.jackn <- t(apply(SQres.jackn, 1, "+", SQres)) #apply jackknife function to matrix
-   phi.mean <- apply(SQres.jackn, 2, mean, na.rm = TRUE) #mean jackknife estimate per link
-   phi.UCI <- apply(SQres.jackn, 2, stats::sd, na.rm = TRUE) #standard deviation of estimates
-   phi.UCI <- phi.mean + t.critical * phi.UCI/sqrt(sum(D$HP))
-   D$jackknife <- list(mean = phi.mean, upper = phi.UCI)
-   return(D)
+  SQres.jackn <- SQres.jackn^2 #Jackknifed residuals are squared
+  SQres <- (stats::residuals(D$proc))^2 # Vector of original square residuals
+  #jackknife calculations:
+  SQres.jackn <- SQres.jackn*(-(sum(D$HP)-1))
+  SQres <- SQres*sum(D$HP)
+  SQres.jackn <- t(apply(SQres.jackn, 1, "+", SQres)) #apply jackknife function to matrix
+  phi.mean <- apply(SQres.jackn, 2, mean, na.rm = TRUE) #mean jackknife estimate per link
+  phi.UCI <- apply(SQres.jackn, 2, stats::sd, na.rm = TRUE) #standard deviation of estimates
+  phi.UCI <- phi.mean + t.critical * phi.UCI/sqrt(sum(D$HP))
+  D$jackknife <- list(mean = phi.mean, upper = phi.UCI)
+  return(D)
 }
 
 #PACo setting the ith link = 0

cran/paco.Rcheck/00check.log

@@ -0,0 +1,9 @@
+* using log directory ‘/home/matthew/Desktop/paco/cran/paco.Rcheck’
+* using R version 3.4.4 (2018-03-15)
+* using platform: x86_64-pc-linux-gnu (64-bit)
+* using session charset: UTF-8
+* using option ‘--as-cran’
+* checking for file ‘paco/DESCRIPTION’ ... OK
+* this is package ‘paco’ version ‘0.4.0’
+* package encoding: UTF-8
+* checking CRAN incoming feasibility ...

cran/paco.Rcheck/00install.out

@@ -0,0 +1,9 @@
+* installing *source* package ‘paco’ ...
+** R
+** data
+** preparing package for lazy loading
+** help
+*** installing help indices
+** building package indices
+** testing if installed package can be loaded
+* DONE (paco)