version 0.3

DESCRIPTION

@@ -1,17 +1,19 @@
 Package: coalescentMCMC
-Version: 0.2
-Date: 2013-01-02
+Version: 0.3
+Date: 2013-10-27
 Title: MCMC Algorithms for the Coalescent
-Author: Emmanuel Paradis
-Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
-Depends: ape, pegas, phangorn, adegenet, coda
+Authors@R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email = "Emmanuel.Paradis@ird.fr"))
+Depends: ape, coda
+Imports: phangorn, stats
 ZipData: no
 Description: coalescentMCMC provides a flexible framework for
   coalescent analyses in R. It includes a main function running the
-  MCMC algorithm, auxilliary functions for tree rearrangement, and some
+  MCMC algorithm, auxiliary functions for tree rearrangement, and some
   functions to compute population genetic parameters.
 License: GPL (>= 2)
-Packaged: 2013-07-19 10:01:33 UTC; paradis
+Packaged: 2013-10-27 04:34:21 UTC; paradis
+Author: Emmanuel Paradis [aut, cre, cph]
+Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
 NeedsCompilation: yes
 Repository: CRAN
-Date/Publication: 2013-07-19 16:23:13
+Date/Publication: 2013-10-27 07:08:27

MD5

@@ -1,26 +1,24 @@
 0632cc7589c46bbf4981f0a778f1ef2a *COPYING
-41db9f31b72a1a1d407dcdc7ef1f3442 *DESCRIPTION
-67aaf5cffbeea97b2eeaf843d51736a8 *NAMESPACE
-cbba597e045c358f3de88d65790d0ddf *NEWS
-76fc210f685a966a80af3d766bf2a5f9 *R/coalescentMCMC.R
+055c18b2e6e0a17502d01cd1f894166b *DESCRIPTION
+4c0e9e2412b04458ac4fe1f24368b785 *NAMESPACE
+d9e6cdd8fe023cb932c9202efbc3de92 *NEWS
+0f81220bc1e0b8d81bec4efdc50214a6 *R/coalescentMCMC.R
 a52942fca58804291fa4cca2a5e9c3dd *R/dcoal.R
-01edab916944e326887dcf69a832c144 *R/treeOperators.R
+5855c71436ebb147dfc8d5b038d500ef *R/treeOperators.R
 d3bc00b2363b4155824eff32034801ad *R/zzz.R
-e5199c70ce61c9060f806ed82cb289d4 *Thanks
+01319ab1a68cc304cbb0593461e10541 *build/vignette.rds
 dfa3ae7f0e0bc0918691ff5894094c03 *inst/doc/CoalescentModels.R
 38240cef46b41eeb490ec9514930a48c *inst/doc/CoalescentModels.Rnw
-0bcdbabfb7092b7640dfcfbef43d486a *inst/doc/CoalescentModels.pdf
-89de4cff6c4b5e549a113a2a9fbdb011 *inst/doc/Running_coalescentMCMC.R
-e23289a4d68e33f3c7c9a15ab6b41b7a *inst/doc/Running_coalescentMCMC.Rnw
-9a296a8f19c56174dc7e784bfaa19c79 *inst/doc/Running_coalescentMCMC.pdf
-f5d8bb0591f7b77f4f0741657243cce6 *man/coalesceMCMC.Rd
+8da39d8199433c96cf9ee1c1fa9ff5e1 *inst/doc/CoalescentModels.pdf
+0d5a83d0784207bb0911e0352170ab6b *inst/doc/Running_coalescentMCMC.R
+e2e9e5b527e95804201ce3d789cc2eed *inst/doc/Running_coalescentMCMC.Rnw
+90a4f0bf20a91e8348586a180be3c104 *inst/doc/Running_coalescentMCMC.pdf
+5776e0661b44e5213b0689cc176aa3d9 *man/coalesceMCMC.Rd
 9aa980aea36886b18ac6a09b823765da *man/coalescentMCMC-internal.Rd
-75cee3886d17bba64ac13c52d9e87beb *man/dcoal.Rd
-822affb72dc7c49b0b62f4d05bb47acf *man/treeOperators.Rd
+f768c542f46f407ec6cce19410f749cc *man/dcoal.Rd
+90b1945d199985a0b7162999b6f140d0 *man/treeOperators.Rd
 2a6f9e9e044a78154d3cfda5936d6f48 *src/Makevars
-412493cc6a8a250e740281af9e1c9086 *src/coalescentMCMC.c
+c3e22b2a119684f9a878926188ddd39e *src/coalescentMCMC.c
 38240cef46b41eeb490ec9514930a48c *vignettes/CoalescentModels.Rnw
-9d144704316ef46e6c8ec412bc1fde3e *vignettes/Running_coalescentMCMC-001.pdf
-8dc91bbd23b6a474511cec93c97a6534 *vignettes/Running_coalescentMCMC-002.pdf
-e23289a4d68e33f3c7c9a15ab6b41b7a *vignettes/Running_coalescentMCMC.Rnw
+e2e9e5b527e95804201ce3d789cc2eed *vignettes/Running_coalescentMCMC.Rnw
 8bc7717e10b96e671d31408560cf5981 *vignettes/coalescentMCMC.bib

NAMESPACE

@@ -1,8 +1,10 @@
-useDynLib(coalescentMCMC)
+useDynLib(coalescentMCMC, .registration = TRUE)
 
-exportPattern("^[^\\.]")
+export(.coalescentMCMCenv, coalescentMCMC, dcoal, dcoal.linear,
+       dcoal.step, dcoal.time, dcoal.time2, EdgeLengthJittering,
+       getMCMCtrees, NeighborhoodRearrangement, TipInterchange)
 
-import(ape)
-import(phangorn)
-import(pegas)
-import(coda)
+importFrom(ape, as.phylo, branching.times, dist.dna, Ntip, reorder.phylo)
+importFrom(coda, mcmc)
+importFrom(phangorn, phyDat, pml)
+importFrom(stats, hclust, rbinom, reorder, runif)

NEWS

@@ -1,6 +1,40 @@
+		CHANGES IN coalescentMCMC VERSION 0.3
+
+
+NEW FEATURES
+
+    o coalescentMCMC() has been extensively modified. It has a new
+      option (model) to select the coalescent model. Currently, only
+      two models can be fitted.
+
+    o The coalescent parameters are now output as part of the "coda"
+      object.
+
+    o Trees of the MCMCs are extracted with the new function
+      getMCMCtrees. The trees of successive chains are stored and can
+      be retrieved separately. If several lists of trees are stored,
+      getMCMCtrees() will call an interactive menu to select the list
+      to retrieve.
+
+    o The vignette "Running_coalescentMCMC" has been updated and now
+      presents a (almost) complete analysis.
+
+
+OTHER CHANGES
+
+    o The packages adegenet and pegas are no more required.
+
+    o Improved DESCRIPTION and NAMESPACE files.
+
+    o The arguments of NeighborhoodRearrangement() have been modified.
+
+
+
 		CHANGES IN coalescentMCMC VERSION 0.2
 
 
 NEW FEATURES
 
     o The main MCMC has been completely rewritten.
+
+    o The package coda is now used to analyse MCMC outputs.

R/coalescentMCMC.R

@@ -1,4 +1,4 @@
-## coalescentMCMC.R (2013-07-19)
+## coalescentMCMC.R (2013-10-18)
 
 ##   Run MCMC for Coalescent Trees
 
@@ -9,25 +9,50 @@
 
 coalescentMCMC <-
     function(x, ntrees = 3000, burnin = 1000, frequency = 1,
-             tree0 = NULL, quiet = FALSE)
+             tree0 = NULL, model = NULL, quiet = FALSE)
 {
     if (is.null(tree0)) {
         d <- dist.dna(x, "JC69")
-        X <- phangorn::phyDat(x)
         tree0 <- as.phylo(hclust(d, "average"))
     }
 
-    n <- Ntip(tree0)
+    X <- phyDat(x)
+    n <- length(tree0$tip.label)
     nodeMax <- 2*n - 1
-    nOut <- ntrees# + burnin
+    nOut <- ntrees
+    nOut2 <- ntrees * frequency + burnin
 
-    getlogLik <- function(phy, X) phangorn::pml(phy, X)$logLik #phangorn:::pml6(phy, X)
+    getlogLik <- function(phy, X) pml(phy, X)$logLik
 
     TREES <- vector("list", nOut)
-    LL <- numeric(nOut)
+    LL <- numeric(nOut2)
     TREES[[1L]] <- tree0
     lnL0 <- getlogLik(tree0, X)
-    LL[[1L]] <- lnL0
+    LL[1L] <- lnL0
+
+    if (is.null(model)) {
+        np <- 1L
+        para.nms <- "theta"
+        ## quantities to calculate THETA:
+        two2n <- 2:n
+        K4theta <- length(two2n)
+        tmp <- choose(two2n, 2)
+        getparams <- function(phy, bt) {
+            x4theta <- rev(diff(c(0, sort(bt))))
+            sum(x4theta * tmp)/K4theta
+        }
+        f.theta <- function(t, p) p
+    } else { # only "time" (ie, exponential model) for the moment
+        np <- 2L
+        para.nms <- c("theta0", "rho")
+        getparams <- function(phy, bt) { # 'bt' is not used but is needed to have the same arguments than above
+            out <- nlminb(c(0.02, 0),
+                          function(p) -dcoal.time(phy, p[1], p[2], log = TRUE))
+            out$par
+        }
+        f.theta <- function(t, p) p[1] * exp(p[2] * t)
+    }
+    params <- matrix(0, nOut2, np)
 
     i <- 2L
     j <- 0L # number of accepted trees
@@ -42,21 +67,34 @@ coalescentMCMC <-
         cat("Generation    Nb of accepted trees\n")
     }
 
+    ##
+    bt0 <- branching.times(tree0)
+    params[1L, ] <- para0 <- getparams(tree0, bt0)
+
+    nodesToSample <- (n + 2):nodeMax
+
     while (k < nOut) {
         if (!quiet)
             cat("\r  ", i, "                ", j, "           ")
 
-        tr.b <- NeighborhoodRearrangement(tree0, n, nodeMax)
+        ## select one internal node excluding the root:
+        target <- sample(nodesToSample, 1L) # target node for rearrangement
+        THETA <- f.theta(bt0[target - n], para0) # the value of THETA at this node
+
+        tr.b <- NeighborhoodRearrangement(tree0, n, nodeMax, target, THETA, bt0)
+        ## do TipInterchange() every 10 steps:
+        ## tr.b <-
+        ##     if (!i %% 10) TipInterchange(tree0, n)
+        ##     else NeighborhoodRearrangement(tree0, n, nodeMax, target, THETA, bt0)
         if (!(i %% frequency) && i > burnin) {
             k <- k + 1L
             TREES[[k]] <- tr.b
         }
-        ## do TipInterchange() every 10 steps:
-        ## tr.b <-
-        ##     if (!i %% 10) TipInterchange(TREES[[i]], n)
-        ##     else NeighborhoodRearrangement(TREES[[i]], n, nodeMax)
         lnL.b <- getlogLik(tr.b, X)
-        LL[[i]] <- lnL.b
+        LL[i] <- lnL.b
+        ## calculate theta for the proposed tree:
+        bt <- branching.times(tr.b)
+        params[i, ] <- para <- getparams(tr.b, bt)
         i <- i + 1L
         ACCEPT <- if (is.na(lnL.b)) FALSE else {
             if (lnL.b >= lnL0) TRUE
@@ -66,19 +104,40 @@ coalescentMCMC <-
             j <- j + 1L
             lnL0 <- lnL.b
             tree0 <- tr.b
+            para0 <- para
+            bt0 <- bt
         }
     }
 
-    dim(LL) <- c(i - 1, 1)
-    colnames(LL) <- "logLik"
+    #dim(LL) <- c(i - 1, 1)
+    LL <- cbind(LL, params)
+    colnames(LL) <- c("logLik", para.nms)
     LL <- mcmc(LL, start = 1, end = i - 1)
 
+    ## compress the list of trees:
+    attr(TREES, "TipLabel") <- TREES[[1L]]$tip.label
+    for (i in seq_len(nOut)) TREES[[i]]$tip.label <- NULL
     class(TREES) <- "multiPhylo"
-    assign(".TREES", TREES, envir = .coalescentMCMCenv)
-    ## TREES <- .compressTipLabel(TREES)
+
+    j <- 1
+    list.trees <- ls(envir = .coalescentMCMCenv)
+    if (l <- length(list.trees))
+        j <- 1 + as.numeric(sub("TREES_", "", list.trees[l]))
+    assign(paste("TREES", j, sep = "_"), TREES, envir = .coalescentMCMCenv)
     if (!quiet) cat("\nDone.\n")
-    #list(mcmc = LL, trees = TREES)
     LL
 }
 
-getMCMCtrees <- function() get(".TREES", envir = .coalescentMCMCenv)
+getMCMCtrees <- function() {
+    list.trees <- ls(envir = .coalescentMCMCenv)
+    l <- length(list.trees)
+    if (!l) return(NULL)
+    if (l == 1)
+        return(get(list.trees, envir = .coalescentMCMCenv))
+    ## l > 1:
+    cat("Several lists of MCMC trees are stored:\n\n")
+    for (i in 1:l) cat(i, ":", list.trees[i], "\n")
+    cat("\nReturn which number? ")
+    i <- as.numeric(readLines(n = 1))
+    get(paste("TREES", i, sep = "_"), envir = .coalescentMCMCenv)
+}

R/treeOperators.R

@@ -1,4 +1,4 @@
-## treeOperators.R (2013-01-10)
+## treeOperators.R (2013-10-18)
 
 ##   Trees Operators for Running MCMC
 
@@ -9,22 +9,18 @@
 
 getIndexEdge <- function(tip, edge)
     ## 'integer(1)' mustn't be substituted by '0L' except if 'DUP = TRUE':
-    .C("get_single_index_integer", as.integer(edge[, 2L]),
-       as.integer(tip), integer(1L), PACKAGE = "coalescentMCMC",
-       NAOK = TRUE, DUP = FALSE)[[3L]]
+    .C(get_single_index_integer, as.integer(edge[, 2L]),
+       as.integer(tip), integer(1L), NAOK = TRUE, DUP = FALSE)[[3L]]
 
 getIndexEdge2 <- function(node, edge)
-    .C("get_two_index_integer", as.integer(edge[, 1L]),
-       as.integer(node), integer(2L), PACKAGE = "coalescentMCMC",
-       NAOK = TRUE, DUP = FALSE)[[3L]]
+    .C(get_two_index_integer, as.integer(edge[, 1L]),
+       as.integer(node), integer(2L), NAOK = TRUE, DUP = FALSE)[[3L]]
 
-NeighborhoodRearrangement <- function(phy, n, nodeMax)
+NeighborhoodRearrangement <- function(phy, n, nodeMax, target, THETA, brtimes)
 {
-    THETA <- pegas::theta.tree(phy, 1)$theta
-    bt <- c(rep(0, n), branching.times(phy))
-    ## select one internal node excluding the root
-    target <- sample((n + 2):nodeMax, size = 1L)
-
+    ## pegas is no more needed:
+    ## THETA <- theta.tree(phy, 1)$theta
+    bt <- c(rep(0, n), brtimes)
     e <- phy$edge # local copy
 
 ### i1, i2, and i3 are edge indices
@@ -33,13 +29,11 @@ NeighborhoodRearrangement <- function(phy, n, nodeMax)
     ## i1 <- which(e2 == target)
     i1 <- getIndexEdge(target, e)
     anc <- e[i1, 1L] # the ancestor of 'target'
-    ## i2 <- which(e1 == target)
     i2 <- getIndexEdge2(target, e) # the 2 edges where 'target' is basal
-    ## i3 <- which(e1 == anc)
     i3 <- getIndexEdge2(anc, e) # this includes i1, so:
     i3 <- i3[i3 != i1]
     sister <- e[i3, 2L] # the sister-node of 'target'
-    sel <- sample(2L, size = 1L)
+    sel <- sample.int(2L, 1L)
     i2.move <- i2[sel]
     i2.stay <- i2[-sel]
     phy$edge[i3, 2L] <- child2move <- e[i2.move, 2L]
@@ -89,10 +83,10 @@ EdgeLengthJittering <- function(phy)
 ### all edge lengths are added to a random value on U[-MIN, MAX]
 ### (the ultrametric nature of the tree is kept)
 {
-   z <- range(phy$edge.length)
-   MIN <- z[1]
-   MAX <- z[2]
-   x <- runif(1, -MIN, MAX) # should be OK even if MIN=0
-   phy$edge.length <- phy$edge.length + x
-   phy
+    z <- range(phy$edge.length)
+    MIN <- z[1]
+    MAX <- z[2]
+    x <- runif(1, -MIN, MAX) # should be OK even if MIN=0
+    phy$edge.length <- phy$edge.length + x
+    phy
 }