[adegenet-commits] r725 - pkg/R

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

Author: jombart
Date: 2010-11-24 18:31:03 +0100 (Wed, 24 Nov 2010)
New Revision: 725

Modified:
   pkg/R/inbreeding.R
Log:
added the awesome, the marvelous, inbreeding.ml


Modified: pkg/R/inbreeding.R
===================================================================
--- pkg/R/inbreeding.R	2010-11-24 15:18:03 UTC (rev 724)
+++ pkg/R/inbreeding.R	2010-11-24 17:31:03 UTC (rev 725)
@@ -77,3 +77,115 @@
     return(res)
 
 } # end inbreeding
+
+
+
+
+
+
+
+
+
+###############
+## inbreeding.ml
+###############
+inbreeding.ml <- function(x, pop=NULL, truenames=TRUE, res.type=c("sample","function"), N=500, M=N*10, ...){
+    ## CHECKS ##
+    if(!is.genind(x)) stop("x is not a valid genind object")
+    checkType(x)
+    res.type <- match.arg(res.type)
+
+    if(x$ploidy != 2) stop("this inbreeding coefficient is designed for diploid genotypes only")
+
+    if(!is.null(pop)) pop(x) <- pop
+    if(is.null(x at pop) && is.null(pop)) {
+        pop(x) <- factor(rep(1, nrow(x at tab)))
+    }
+
+      ## COMPUTATIONS ##
+
+    ## get allele frequencies and \sum p_i^2 by pop and loc ##
+    tabfreq2 <- (makefreq(x = genind2genpop(x, quiet = TRUE), quiet=TRUE, truenames=truenames)$tab) ^2
+    sumpi2 <- t(apply(tabfreq2, 1, tapply, x$loc.fac, sum))
+
+    ## function to check a 1-locus genotype for homozigosity
+    ## returns 1 if homoz, 0 otherwise
+    ## !!! NOTE : reverse the values returned by f1 to obtain a strange thing !!!
+    f1 <- function(gen){
+        if(any(is.na(gen))) return(NA)
+        if(any(round(gen, 10)==1)) return(1)
+        return(0)
+    }
+
+    ## get the table of binary hetero/homo data
+    if(truenames) {
+        X <- truenames(x)$tab
+    } else
+    X <- x$tab
+
+    homotab <- t(apply(X, 1, tapply, x at loc.fac, f1))
+
+
+    ## get pi2 for the appropriate pop
+    if(truenames){
+    popx <- pop(x)
+    } else {
+        popx <- x$pop
+    }
+
+    popx <- as.character(popx)
+    tabpi2 <- sumpi2[popx, , drop=FALSE]
+
+
+
+    ## function returning a likelihood function - multi-locus ##
+    ## x is a vector of 1/0 for each locus, with 1=homoz
+    LIK <- function(x, sumpi2){
+        ## return( sum(log(   x*(F+(1-F)*sumpi2) + (1-x)*(1-(F+(1-F)*sumpi2)) ) ,na.rm=TRUE) ) # returns the log-likelihood
+        myEnv <- new.env()
+        assign("x", x, envir=myEnv)
+        assign("sumpi2", sumpi2, envir=myEnv)
+        res <- function(F) {
+            ## cat("\nx used:\n") # debugging
+            ## print(x)
+            ## cat("\nsumpi2 used:\n") # debugging
+            ## print(sumpi2)
+
+            return(exp(sum(log(   x*(F+(1-F)*sumpi2) + (1-x)*(1-(F+(1-F)*sumpi2)) ) ,na.rm=TRUE)))
+        }
+        environment(res) <- myEnv
+
+        return(res)
+    }
+
+
+    ## get likelihood functions for all individuals
+    res <- lapply(1:nrow(homotab), function(i) LIK(homotab[i,], tabpi2[i,]) )
+    res <- lapply(res, Vectorize)
+
+    ## name output
+    if(truenames) {
+        names(res) <- x$ind.names
+    } else {
+        names(res) <- names(x$tab)
+    }
+
+
+    ## IF WE RETURN FUNCTIONS ##
+    if(res.type=="function"){
+          return(res)
+    }
+
+
+    ## IF WE RETURN SAMPLES ##
+    ## function to get one sample
+    getSample <- function(f){ # f is a vectorized density function
+        x <- runif(M, 0, 1)
+        fx <- f(x)
+        fx <- fx/sum(fx)
+        return(sample(x, size=N, prob=fx))
+    }
+
+    res <- lapply(res, getSample)
+    return(res)
+} # end inbreeding.ml