diff --git a/DESCRIPTION b/DESCRIPTION index 8a0c5c7a..b4492e4d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,12 @@ Package: phytools -Version: 0.6-26 -Date: 2017-09-18 +Version: 0.6-27 +Date: 2017-09-19 Title: Phylogenetic Tools for Comparative Biology (and Other Things) Author: Liam J. Revell Maintainer: Liam J. Revell Depends: R (>= 3.2.0), ape (>= 4.0), maps Imports: animation, clusterGeneration, coda, combinat, graphics, grDevices, - methods, mnormt, msm, nlme, numDeriv, phangorn (>= 2.1.1), plotrix, + MASS, methods, mnormt, msm, nlme, numDeriv, phangorn (>= 2.1.1), plotrix, scatterplot3d, stats, utils Suggests: geiger, rgl ZipData: no @@ -56,6 +56,6 @@ Description: Package contains various functions for phylogenetic analysis. research. License: GPL (>= 2) URL: http://github.com/liamrevell/phytools -Packaged: 2017-09-18 12:00:00 EST +Packaged: 2017-09-19 12:00:00 EST Repository: -Date/Publication: 2017-09-18 12:00:00 EST +Date/Publication: 2017-09-19 12:00:00 EST diff --git a/NAMESPACE b/NAMESPACE index 17ac1635..a13f62f7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -132,3 +132,4 @@ pdf, dev.off, col2rgb) importFrom(combinat, permn) importFrom(coda, HPDinterval) importFrom(nlme, gls, varFixed) +importFrom(MASS, ginv) diff --git a/R/ratebytree.R b/R/ratebytree.R index 8f0dd70d..bb255935 100644 --- a/R/ratebytree.R +++ b/R/ratebytree.R @@ -233,7 +233,7 @@ rbt.cont<-function(trees,x,...){ ## compute covariance matrix H.multi<-hessian(lik.multi,fit.multi$par,trees=trees,y=x, se=se,model=model,trace=FALSE) - Cov.multi<-solve(H.multi) + Cov.multi<-if(qr(H.multi)$rank!=ncol(H.multi)) ginv(H.multi) else solve(H.multi) ## now fit single-rate model lik.onerate<-function(theta,trees,y,se,model,trace=FALSE){ n<-sapply(trees,Ntip) @@ -289,7 +289,7 @@ rbt.cont<-function(trees,x,...){ ## compute covariance matrix H.onerate<-hessian(lik.onerate,fit.onerate$par,trees=trees,y=x, se=se,model=model,trace=FALSE) - Cov.onerate<-solve(H.onerate) + Cov.onerate<-if(qr(H.onerate)$rank!=ncol(H.onerate)) ginv(H.onerate) else solve(H.onerate) ## compare models: LR<-2*(-fit.multi$value+fit.onerate$value) km<-2*N+if(model=="BM") 0 else if(model%in%c("OU","EB")) N diff --git a/man/ratebytree.Rd b/man/ratebytree.Rd index c72cd747..94f013e5 100644 --- a/man/ratebytree.Rd +++ b/man/ratebytree.Rd @@ -27,11 +27,13 @@ posthoc.ratebytree(x, ...) } \details{ At present it is not possible to specify different models to fit for the different trees - although if (for instance) character evolution on tree 1 proceeded by a strong \emph{OU} process while character evolution on tree 2 was by \emph{BM}, we would probably reject a constant-process model and tree 2 should show a very low value of \code{alpha}. + + To compute the standard errors for each fitted paramater value, the function computes the negative inverse of the Hessian matrix at the MLEs; however, if this matrix is computationally singular generalized inverse (\code{\link{ginv}}) will be used instead without warning. The function also conducts a likelihood-ratio test to compare the two models. } \value{ - An object of class \code{"ratebytree"}. + An object of class \code{"ratebytree"} or an object of class \code{"posthoc.ratebytree"} in the case of the method \code{posthoc}. } \references{ Adams, D. C. (2012) Comparing evolutionary rates for different phenotypic traits on a phylogeny using likelihood. \emph{Syst. Biol.}, \bold{62}, 181-192.