Permalink
Browse files

update to ratebytree to use ginv if Hessian nearly singular

  • Loading branch information...
liamrevell committed Sep 19, 2017
1 parent 4407c72 commit a510bb2420d28a893f1df62f8d287735b38b65ab
Showing with 11 additions and 8 deletions.
  1. +5 −5 DESCRIPTION
  2. +1 −0 NAMESPACE
  3. +2 −2 R/ratebytree.R
  4. +3 −1 man/ratebytree.Rd
View
@@ -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 <liam.revell@umb.edu>
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
View
@@ -132,3 +132,4 @@ pdf, dev.off, col2rgb)
importFrom(combinat, permn)
importFrom(coda, HPDinterval)
importFrom(nlme, gls, varFixed)
importFrom(MASS, ginv)
View
@@ -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
View
@@ -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.

0 comments on commit a510bb2

Please sign in to comment.