Permalink
Browse files

changes brownieREML including object class

  • Loading branch information...
liamrevell committed Sep 21, 2015
1 parent efd9879 commit ee876147e9a41bdea4d6ccd566d078cd4f1d029c
Showing with 30 additions and 1 deletion.
  1. +1 −0 NAMESPACE
  2. +29 −1 R/brownieREML.R
View
@@ -67,6 +67,7 @@ S3method(print, anc.ML)
S3method(print, fitMk)
S3method(summary, fitMk)
S3method(logLik, fitMk)
S3method(print, brownieREML)
importFrom(animation, ani.options, ani.record, ani.replay, saveVideo)
importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, bind.tree, branching.times, collapse.singles)
View
@@ -17,6 +17,8 @@ brownieREML<-function(tree,x,maxit=2000){
}
sig1<-mean(pic(x,tree)^2)
logL1<--lik1(sig1,tree,x)
H<-optimHess(sig1,lik1,tree=tree,x=x)
v1<-1/H
# fit the multiple rate model
lik2<-function(sig2,tree,x){
tt<-scaleByMap(tree,sig2)
@@ -26,9 +28,16 @@ brownieREML<-function(tree,x,maxit=2000){
}
YY<-optim(setNames(rep(1,p)*runif(n=p),colnames(tree$mapped.edge)),lik2,tree=tree,x=x,method="L-BFGS-B",lower=rep(1e-8,p))
sig2<-YY$par
v2<-solve(optimHess(sig2,lik2,tree=tree,x=x))
logL2<--YY$value
if(YY$convergence==0) converged<-"Optimization has converged."
else converged<-"Optimization may not have converged. Consider increasing maxit."
convergence=(YY$convergence==0)
return(list(sig2.single=sig1,logL1=logL1,sig2.multiple=sig2,logL2=logL2,convergence=convergence))
obj<-list(sig2.single=sig1,var.single=v1,logL1=logL1,k1=1,sig2.multiple=sig2,vcv.multiple=v2,logL.multiple=logL2,
k2=length(sig2),convergence=converged)
class(obj)<-"brownieREML"
obj
}
# This function scales a mapped tree by sig2
@@ -39,3 +48,22 @@ scaleByMap<-function(mtree,sig2){
class(tree)<-"phylo"
return(tree)
}
## S3 print method for "brownieREML"
## S3 print method for object of class "brownie.lite"
## written by Liam J. Revell 2013
print.brownieREML<-function(x, ...){
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
x<-lapply(x,function(a,b) if(is.numeric(a)) round(a,b) else a,b=digits)
cat("REML single-rate model:\n")
cat("\ts^2\tse\tk\tlogL\n")
cat(paste("value",x$sig2.single,round(sqrt(x$var.single),digits),x$k1,x$logL1,"\n\n",sep="\t"))
cat("REML multi-rate model:\n")
cat(paste(c("",paste("s^2(",names(x$sig2.multiple),")","\tse(",names(x$sig2.multiple),")",sep=""),
"k","logL","\n"),collapse="\t"))
cat(paste(paste(c("value",paste(x$sig2.multiple,round(sqrt(diag(x$vcv.multiple)),digits),sep="\t"),x$k2,
x$logL.multiple),collapse="\t"),"\n\n",sep=""))
if(x$convergence[1]=="Optimization has converged.") cat("R thinks it has found the ML solution.\n\n")
else cat("Optimization may not have converged.\n\n")
}

0 comments on commit ee87614

Please sign in to comment.