[Adephylo-commits] r78 - in pkg: . R man

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2008-12-03 11:40:14 +0100 (Wed, 03 Dec 2008)
New Revision: 78

Modified:
   pkg/R/ppca.R
   pkg/TODO
   pkg/man/orthobasis.Rd
   pkg/man/ppca.Rd
Log:
ppca is messed up. Results differ from original results that used multispati. Eigen analysis is wrong (small differences in eigenvalues/eigenvectors) 


Modified: pkg/R/ppca.R
===================================================================
--- pkg/R/ppca.R	2008-12-02 19:49:34 UTC (rev 77)
+++ pkg/R/ppca.R	2008-12-03 10:40:14 UTC (rev 78)
@@ -6,7 +6,7 @@
 ################
 # Function ppca
 ################
-ppca <- function(x, prox=NULL, method=c("patristic","nNodes","Abouheif","sumDD"), a=1,
+ppca <- function(x, prox=NULL, method=c("patristic","nNodes","oriAbouheif","Abouheif","sumDD"), a=1,
                  center=TRUE, scale=TRUE, scannf=TRUE, nfposi=1, nfnega=0){
 
     ## handle arguments
@@ -16,6 +16,7 @@
 
     tre <- as(x, "phylo4")
     method <- match.arg(method)
+    NEARZERO <- 1e-10
 
     ## proximity matrix
     if(is.null(prox)){ # have to compute prox
@@ -44,21 +45,19 @@
 
     ## replace NAs
     f1 <- function(vec){
-        if(any(is.na(vec))) {
-            warning("Replacing missing values (NA) by mean values")
-            m <- mean(vec,na.rm=TRUE)
-            vec[is.na(vec)] <- m
-        }
-
+        m <- mean(vec,na.rm=TRUE)
+        vec[is.na(vec)] <- m
         return(vec)
     }
 
-    X <- as.data.frame(lapply(X, f1))
+    if(any(is.na(X))) {
+        warning("Replacing missing values (NA) by mean values")
+        X <- as.data.frame(apply(X, 2, f1))
+    }
+
     X <- scalewt(X, center=center, scale=scale) # centring/scaling of traits
 
-    ##
 
-
     ## main computation ##
 
     ## make a skeleton of dudi
@@ -67,13 +66,18 @@
 
     ## computations of the ppca
     X <- as.matrix(X)
-    decomp <- eigen((t(X) %*% W %*% X)/N, sym=TRUE)
+    decomp <- eigen( ((t(X) %*% W %*% X)/N), sym=TRUE)
     U <- decomp$vectors # U: principal axes
-    p <- ncol(U)
     lambda <- decomp$values
 
+    ## remove null eigenvalues and corresponding vectors
+    toKeep <- (abs(lambda) > NEARZERO)
+    lambda <- lambda[toKeep]
+    U <- U[, toKeep]
+    p <- ncol(U)
+
     if(scannf){ # interactive part
-        barplot(lambda[1:res$rank])
+        barplot(lambda)
         cat("Select the number of global axes: ")
         nfposi <- as.integer(readLines(n = 1))
         cat("Select the number of local axes: ")

Modified: pkg/TODO
===================================================================
--- pkg/TODO	2008-12-02 19:49:34 UTC (rev 77)
+++ pkg/TODO	2008-12-03 10:40:14 UTC (rev 78)
@@ -55,13 +55,14 @@
 - procella.Rd -- done
 - tithonia.Rd -- done
 - ungulates.Rd -- done
+* build an orthobasis from treePart dummy variables
 * add a Moran's index and a test (is it needed, as we have morangeary ?)
 
 
 
 # TESTING:
 ==========
-* test distTips a little more (see Abouheif with an outgroup, otherwise seems ok)
+* test distTips a little more (see Abouheif with an outgroup, otherwise seems ok) -- done (TJ)
 * test ppca: are results identical to those with previous implementation?
 
 

Modified: pkg/man/orthobasis.Rd
===================================================================
--- pkg/man/orthobasis.Rd	2008-12-02 19:49:34 UTC (rev 77)
+++ pkg/man/orthobasis.Rd	2008-12-03 10:40:14 UTC (rev 78)
@@ -34,8 +34,9 @@
     - \code{patristic}: (inversed sum of) branch lengths \cr
     - \code{nNodes}: (inversed) number of nodes on the path between the
     nodes \cr
-    - \code{oriAbouheif}: original Abouheif's proximity, with diagonal (see details) \cr
-    - \code{Abouheif}: Abouheif's proximity (see details) \cr
+    - \code{oriAbouheif}: original Abouheif's proximity, with diagonal
+  (see details in \code{\link{proxTips}}) \cr
+    - \code{Abouheif}: Abouheif's proximity (see details in \code{\link{proxTips}}) \cr
     - \code{sumDD}: (inversed) sum of direct descendants of all nodes on the path
     (see details in \code{\link{proxTips}}).}
   \item{a}{the exponent used to compute the proximity (see ?\code{\link{proxTips}}).}

Modified: pkg/man/ppca.Rd
===================================================================
--- pkg/man/ppca.Rd	2008-12-02 19:49:34 UTC (rev 77)
+++ pkg/man/ppca.Rd	2008-12-03 10:40:14 UTC (rev 78)
@@ -36,8 +36,10 @@
     ambiguity) specifying the method used to compute proximities;
     possible values are:\cr
     - \code{patristic}: (inversed sum of) branch lengths \cr
-    - \code{nNodes}: (inversed) number of nodes on the path between the nodes \cr
-    - \code{Abouheif}: Abouheif's proximity (see details) \cr
+    - \code{nNodes}: (inversed) number of nodes on the path between the
+    nodes \cr
+    - \code{oriAbouheif}: original Abouheif's proximity, with diagonal (see details in \code{\link{proxTips}}) \cr
+    - \code{Abouheif}: Abouheif's proximity (see details in \code{\link{proxTips}}) \cr
     - \code{sumDD}: (inversed) sum of direct descendants of all nodes on the path
     (see details in \code{\link{proxTips}}).}
   \item{a}{the exponent used to compute the proximity (see ?\code{\link{proxTips}}).}