[adegenet-commits] r455 - pkg/R

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

Author: jombart
Date: 2009-10-18 18:48:21 +0200 (Sun, 18 Oct 2009)
New Revision: 455

Added:
   pkg/R/PCtest.R
Log:
Added pctest


Added: pkg/R/PCtest.R
===================================================================
--- pkg/R/PCtest.R	                        (rev 0)
+++ pkg/R/PCtest.R	2009-10-18 16:48:21 UTC (rev 455)
@@ -0,0 +1,51 @@
+pctest <- function(x, nperm=99, center=TRUE, scale=TRUE, method=c("sigma", "binom"), quiet=FALSE, plot=TRUE){
+    ## check x
+    if(!is.genind(x) & !is.genpop(x)){
+        stop("x is not a genind or a genpop object")
+    }
+
+    ## a few general variables
+    N <- nrow(x at tab)
+    P <- ncol(x at tab)
+
+    ## make tables of allele frequencies
+    X <- scaleGen(x, center=center, scale=scale, method=method, missing="mean")
+    fac.loc <- factor(sub("[.][^.]*$","",colnames(X)))
+    lX <- lapply(levels(fac.loc), function(id) X[,fac.loc==id,drop=FALSE])
+
+    ## auxil function to compute the first eigenvalue
+    if(N > P){ # N > P
+        f1 <- function(A){
+            Z <- t(A) %*% A / N
+            return(eigen(Z, symmetric=TRUE, only.values=TRUE)$values[1])
+        }
+    } else { #p <= n
+        f1 <- function(A){
+            Z <- A %*% t(A) / N
+            return(eigen(Z, symmetric=TRUE, only.values=TRUE)$values[1])
+        }
+    }
+
+
+    ## Monte Carlo procedure
+    makeOnePerm <- function(listX){
+        return(as.matrix(data.frame(lapply(listX, function(e) e[sample(N),,drop=FALSE]))))
+    }
+
+    if(quiet){
+        sim <- sapply(1:nperm, function(i) f1(makeOnePerm(lX)))
+    } else {
+        cat("\n Computing", nperm, "simulated eigenvalues ")
+        sim <- sapply(1:nperm, function(i) {cat(ifelse(i%%10==0, i, "."));return(f1(makeOnePerm(lX)))} )
+        cat(" done.\n")
+    }
+    ini <- f1(X)
+
+    ## return res
+    myCall <- match.call()
+    res <- as.randtest(sim=sim, obs=ini, alter="greater", call=myCall)
+    if(plot) {
+        plot(res, nclass=NULL, main="1st eigenvalue vs simulated eigenvalues (histogram)")
+    }
+    return(res)
+}