diff --git a/DESCRIPTION b/DESCRIPTION index 3bf90a6..8fedb0c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: diversitree -Version: 0.9-20 +Version: 0.9-23 Title: Comparative 'Phylogenetic' Analyses of Diversification Authors@R: c(person("Richard G.", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), diff --git a/man/make.bisse.Rd b/man/make.bisse.Rd index ab37604..b983ee0 100644 --- a/man/make.bisse.Rd +++ b/man/make.bisse.Rd @@ -197,105 +197,6 @@ plot(h, phy) lik <- make.bisse(phy, phy$tip.state) lik(pars) # -159.71 - -## Heuristic guess at a starting point, based on the constant-rate -## birth-death model (uses \link{make.bd}). -p <- starting.point.bisse(phy) - -\dontrun{ -## Start an ML search from this point. This takes some time (~7s) -fit <- find.mle(lik, p, method="subplex") -logLik(fit) # -158.6875 - -## The estimated parameters aren't too far away from the real ones, even -## with such a small tree -rbind(real=pars, - estimated=round(coef(fit), 2)) - -## Test a constrained model where the speciation rates are set equal -## (takes ~4s). -lik.l <- constrain(lik, lambda1 ~ lambda0) -fit.l <- find.mle(lik.l, p[-1], method="subplex") -logLik(fit.l) # -158.7357 - -## Despite the difference in the estimated parameters, there is no -## statistical support for this difference: -anova(fit, equal.lambda=fit.l) - -## Run an MCMC. Because we are fitting six parameters to a tree with -## only 50 species, priors will be needed. I will use an exponential -## prior with rate 1/(2r), where r is the character independent -## diversificiation rate: -prior <- make.prior.exponential(1 / (2 * (p[1] - p[3]))) - -## This takes quite a while to run, so is not run by default -tmp <- mcmc(lik, fit$par, nsteps=100, prior=prior, w=.1, print.every=0) - -w <- diff(sapply(tmp[2:7], range)) -samples <- mcmc(lik, fit$par, nsteps=1000, prior=prior, w=w, - print.every=100) - -## See \link{profiles.plot} for more information on plotting these -## profiles. -col <- c("blue", "red") -profiles.plot(samples[c("lambda0", "lambda1")], col.line=col, las=1, - xlab="Speciation rate", legend="topright") -} - -## BiSSE reduces to the birth-death model and Mk2 when diversification -## is state independent (i.e., lambda0 ~ lambda1 and mu0 ~ mu1). -lik.mk2 <- make.mk2(phy, phy$tip.state) -lik.bd <- make.bd(phy) - -## 1. BiSSE / Birth-Death -## Set the q01 and q10 parameters to arbitrary numbers (need not be -## symmetric), and constrain the lambdas and mus to be the same for each -## state. The likelihood function now has just two parameters and -## will be proprtional to Nee's birth-death based likelihood: -lik.bisse.bd <- constrain(lik, - lambda1 ~ lambda0, mu1 ~ mu0, - q01 ~ .01, q10 ~ .02) -pars <- c(.1, .03) -## These differ by -167.3861 for both parameter sets: -lik.bisse.bd(pars) - lik.bd(pars) -lik.bisse.bd(2*pars) - lik.bd(2*pars) - -## 2. BiSSE / Mk2 -## Same idea as above: set all diversification parameters to arbitrary -## values (but symmetric this time): -lik.bisse.mk2 <- constrain(lik, - lambda0 ~ .1, lambda1 ~ .1, - mu0 ~ .03, mu1 ~ .03) -## Differ by -150.4740 for both parameter sets. -lik.bisse.mk2(pars) - lik.mk2(pars) -lik.bisse.mk2(2*pars) - lik.mk2(2*pars) - -## 3. Sampled BiSSE / Birth-Death -## Pretend that the tree is only .6 sampled: -lik.bd2 <- make.bd(phy, sampling.f=.6) -lik.bisse2 <- make.bisse(phy, phy$tip.state, sampling.f=c(.6, .6)) -lik.bisse2.bd <- constrain(lik.bisse2, - lambda1 ~ lambda0, mu1 ~ mu0, - q01 ~ .01, q10 ~ .01) - -## Difference of -167.2876 -lik.bisse2.bd(pars) - lik.bd2(pars) -lik.bisse2.bd(2*pars) - lik.bd2(2*pars) - -## 4. Unresolved clade BiSSE / Birth-Death -unresolved <- data.frame(tip.label=I(c("sp25", "sp30", "sp40", "sp56", "sp20")), - Nc =c(10, 9, 6, 5, 2), - n0=0, n1=0) -unresolved.bd <- structure(unresolved$Nc, names=unresolved$tip.label) -lik.bisse3 <- make.bisse(phy, phy$tip.state, unresolved) -lik.bisse3.bd <- constrain(lik.bisse3, - lambda1 ~ lambda0, mu1 ~ mu0, - q01 ~ .01, q10 ~ .01) -lik.bd3 <- make.bd(phy, unresolved=unresolved.bd) - -## Difference of -167.1523 -lik.bisse3.bd(pars) - lik.bd3(pars) -lik.bisse3.bd(pars*2) - lik.bd3(pars*2) } \references{ diff --git a/man/make.mkn.Rd b/man/make.mkn.Rd index b91ae63..03a7a2f 100644 --- a/man/make.mkn.Rd +++ b/man/make.mkn.Rd @@ -82,9 +82,9 @@ make.mkn.meristic(tree, states, k, control=list()) and return both values. (Note that this will not give you a likelihood for use with ML or MCMC functions!). } - \item \code{root.p}{Vector of probabilities/weights to use when + \item \code{root.p}: Vector of probabilities/weights to use when \code{ROOT.GIVEN} is specified. Must be of length \code{k} (2 for - \code{make.mk2}).} + \code{make.mk2}). % TODO: Undocumented: \item \code{intermediates}: Add intermediates to the returned value as diff --git a/man/simulate.Rd b/man/simulate.Rd index 7562abf..69da8b6 100644 --- a/man/simulate.Rd +++ b/man/simulate.Rd @@ -56,7 +56,7 @@ prune(phy, to.drop=NULL) \arguments{ \item{pars}{Vector of parameters. The parameters must be in the same order as an unconstrained likelihood function returned by - \code{make.x}, for tree type {x}. The MuSSE simulator automatically + \code{make.x}, for tree type \code{x}. The MuSSE simulator automatically detects the appropriate number of states, given a parameter vector.} \item{type}{Type of tree to generate: May be "bisse" or "bd".} diff --git a/src/GslOdeBase.cpp b/src/GslOdeBase.cpp index 7b169c1..3d9b119 100644 --- a/src/GslOdeBase.cpp +++ b/src/GslOdeBase.cpp @@ -45,7 +45,7 @@ GslOdeBase::r_derivs(double t_, std::vector y_, SEXP pars_) { // equivalent) so that pars are reset if the R call fails. if (y_.size() != size()) Rf_error("Incorrect input length (expected %d, got %d)", - size(), y_.size()); + (int)size(), (int)y_.size()); set_pars(pars_); std::vector ret(size()); @@ -99,8 +99,8 @@ void GslOdeBase::set_state(double t0, const std::vector y_) { must_be_initialised(); if ( y_.size() != size() ) - Rf_error("Expected 'y' of size %d (recieved %d)", - size(), y_.size()); + Rf_error("Expected 'y' of size %d (recieved %d)", + (int)size(), (int)y_.size()); t = t0; y = y_; diff --git a/src/TimeMachine.cpp b/src/TimeMachine.cpp index 92852db..3f969da 100644 --- a/src/TimeMachine.cpp +++ b/src/TimeMachine.cpp @@ -36,7 +36,7 @@ TimeMachine::TimeMachine(std::vector names, void TimeMachine::set(std::vector pars) { if (pars.size() != np_in) - error("Expected %d parameters, recieved %d", np_in, pars.size()); + error("Expected %d parameters, recieved %d", (int)np_in, (int)pars.size()); // Only go through the extra effort below if the parameters differ. if ( pars == p_in ) return;