[Adephylo-commits] r94 - in pkg: R man

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

Author: jombart
Date: 2008-12-13 15:20:53 +0100 (Sat, 13 Dec 2008)
New Revision: 94

Modified:
   pkg/R/moran.R
   pkg/R/ppca.R
   pkg/man/moranIdx.Rd
Log:
moran.idx can return more info.


Modified: pkg/R/moran.R
===================================================================
--- pkg/R/moran.R	2008-12-13 13:49:45 UTC (rev 93)
+++ pkg/R/moran.R	2008-12-13 14:20:53 UTC (rev 94)
@@ -1,9 +1,8 @@
 ############
 # moran.idx
 ############
-moran.idx <- function(x, prox){
+moran.idx <- function(x, prox, addInfo=FALSE){
 
-
     ## handle arguments
     if(any(is.na(x))) stop("NA entries in x")
     if(!is.numeric(x)) stop("x is not numeric")
@@ -26,6 +25,15 @@
 
     res <- num/denom
 
+    if(addInfo){
+        I0 <- -1/(n-1)
+        matToDiag <- .5 * (t(W) + W)
+        rangeI <- range(eigen(matToDiag)$values)
+        attr(res, "I0") <- I0
+        attr(res, "Imin") <- rangeI[1]
+        attr(res, "Imax") <- rangeI[2]
+    }
+
     return(res)
 
 } # end moran.idx

Modified: pkg/R/ppca.R
===================================================================
--- pkg/R/ppca.R	2008-12-13 13:49:45 UTC (rev 93)
+++ pkg/R/ppca.R	2008-12-13 14:20:53 UTC (rev 94)
@@ -225,6 +225,108 @@
 
 
 
+###############
+# summary.ppca
+###############
+summary.ppca <- function (object, ..., printres=TRUE) {
+
+    ## some checks
+    if (!inherits(object, "ppca"))stop("to be used with 'ppca' object")
+    if(!require(ade4)) stop("The package ade4 is not installed.")
+    
+    
+    norm.w <- function(X, w) {
+        f2 <- function(v) sum(v * v * w)/sum(w)
+        norm <- apply(X, 2, f2)
+        return(norm)
+    }
+    
+    resfin <- list()
+
+    if(printres) {
+        cat("\nphylogenetic Principal Component Analysis\n")
+        cat("\nCall: ")
+        print(object$call)
+    }
+
+    appel <- as.list(object$call)
+    ## compute original pca
+    X <- obj$tab # transformed data
+    
+ 
+    nfposi <- object$nfposi
+    nfnega <- object$nfnega
+
+    dudi <- dudi.pca(X, center=TRUE, scale=FALSE, scannf=FALSE, nf=nfposi+nfnega)
+    ## end of pca
+
+    lw <- object$lw
+
+                                        # I0, Imin, Imax
+    n <- nrow(X)
+    I0 <- -1/(n-1)
+    L <- listw2mat(lw)
+    eigL <- eigen(0.5*(L+t(L)))$values
+    Imin <- min(eigL)
+    Imax <- max(eigL)
+    Ival <- data.frame(I0=I0,Imin=Imin,Imax=Imax)
+    row.names(Ival) <- ""
+    if(printres) {
+        cat("\nConnection network statistics:\n")
+        print(Ival)
+    }
+
+    Istat <- c(I0,Imin,Imax)
+    names(Istat) <- c("I0","Imin","Imax")
+    resfin$Istat <- Istat
+
+
+                                        # les scores de l'analyse de base
+    nf <- dudi$nf
+    eig <- dudi$eig[1:nf]
+    cum <- cumsum(dudi$eig)[1:nf]
+    ratio <- cum/sum(dudi$eig)
+    w <- apply(dudi$l1,2,lag.listw,x=lw)
+    moran <- apply(w*as.matrix(dudi$l1)*dudi$lw,2,sum)
+    res <- data.frame(var=eig,cum=cum,ratio=ratio, moran=moran)
+    row.names(res) <- paste("Axis",1:nf)
+    if(printres) {
+        cat("\nScores from the centred PCA\n")
+        print(res)
+    }
+
+    resfin$pca <- res
+
+
+                                        # les scores de l'analyse spatiale
+                                        # on recalcule l'objet en gardant tous les axes
+    eig <- object$eig
+    nfposimax <- sum(eig > 0)
+    nfnegamax <- sum(eig < 0)
+
+    ms <- multispati(dudi=dudi, listw=lw, scannf=FALSE,
+                     nfposi=nfposimax, nfnega=nfnegamax)
+
+    ndim <- dudi$rank
+    nf <- nfposi + nfnega
+    agarder <- c(1:nfposi,if (nfnega>0) (ndim-nfnega+1):ndim)
+    varspa <- norm.w(ms$li,dudi$lw)
+    moran <- apply(as.matrix(ms$li)*as.matrix(ms$ls)*dudi$lw,2,sum)
+    res <- data.frame(eig=eig,var=varspa,moran=moran/varspa)
+    row.names(res) <- paste("Axis",1:length(eig))
+
+    if(printres) {
+        cat("\nppca eigenvalues decomposition:\n")
+        print(res[agarder,])
+    }
+
+    resfin$ppca <- res
+
+    return(invisible(resfin))
+} # end summary.ppca
+
+
+
 ### testing
 ## obj <- phylo4d(read.tree(text=mjrochet$tre),mjrochet$tab)
 ## x at edge.length= rep(1,length(x at edge.label))

Modified: pkg/man/moranIdx.Rd
===================================================================
--- pkg/man/moranIdx.Rd	2008-12-13 13:49:45 UTC (rev 93)
+++ pkg/man/moranIdx.Rd	2008-12-13 14:20:53 UTC (rev 94)
@@ -6,18 +6,33 @@
   variable and a matrix of proximities among observations.
 }
 \usage{
-moran.idx(x, prox)
+moran.idx(x, prox, addInfo=FALSE)
 }
 \arguments{
   \item{x}{a numeric vector whose autocorrelation is computed.}
   \item{prox}{a matrix of proximities between observations, as computed
    by the \code{\link{proxTips}}. Off-diagonal terms must be positive or
    null, while diagonal terms must all equal zero.}
+ \item{addInfo}{a logical indicating whether supplementary info (null
+   value, minimum and maximum values) should be returned (TRUE) or not
+   (FALSE, default); if computed, these quantities are returned as
+   attributes.}
 }
 \value{
   The numeric value of Moran's index.
 }
 \author{ Thibaut Jombart \email{jombart at biomserv.univ-lyon1.fr} }
+\references{
+  Moran, P.A.P. (1948) The interpretation of statistical
+  maps. \emph{Journal of the Royal Statistical Society, B}
+  \bold{10}, 243--251.
+
+  Moran, P.A.P. (1950) Notes on continuous stochastic
+  phenomena. \emph{Biometrika}, \bold{37}, 17--23.
+
+  de Jong, P. and Sprenger, C. and van Veen, F. (1984) On extreme values
+  of Moran's I and Geary's c. \emph{Geographical Analysis}, \bold{16}, 17--24.
+}
 \seealso{\code{\link{proxTips}} which computes phylogenetic proximities
   between tips of a phylogeny.
 }