Today I
responded
to an R-sig-phylo
query by
posting some code to optimize Pagel's λ tree transformation for a discrete
character evolving by a continuous-time Markov chain. I did this by writing a very simple
wrapper around ape's ace function for ancestral character estimation, which
also fits this model:
## here's the wrapper function
fitLambda<-function(tree,x,model="ER"){
lik<-function(lambda,tree,x,model)
logLik(ace(x,rescale(tree,model="lambda",lambda),
type="discrete",model=model))
obj<-optimize(lik,c(0,1),tree=tree,x=x,model=model,maximum=TRUE)
fit<-ace(x,rescale(tree,model="lambda",lambda=obj$maximum),
type="discrete",model=model)
I<-fit$index.matrix
fitted.Q=matrix(fit$rates[I],dim(I)[1],dim(I)[2],
dimnames=list(dimnames(fit$lik.anc)[[2]],
dimnames(fit$lik.anc)[[2]]))
diag(fitted.Q)<--rowSums(fitted.Q,na.rm=TRUE)
list(Q=fitted.Q,lambda=obj$maximum,logLik=logLik(fit))
}
library(geiger)
library(phytools)
## simulate some data to test it
tree<-pbtree(n=200,scale=1)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
x<-sim.history(rescale(tree,model="lambda",lambda=0.7),Q)$states
## Done simulation(s).
fitLambda(tree,x)
## $Q
## a b
## a -0.7995 0.7995
## b 0.7995 -0.7995
##
## $lambda
## [1] 0.5154
##
## $logLik
## [1] -128.6
Note that the same model can also be fit using geiger's fitDiscrete
function, but the poster reported some issues with convergence which made me want to
post the second way. Here is fitDiscrete
fitDiscrete(tree,x,transform="lambda")
## DLSODA- Warning..Internal T (=R1) and H (=R2) are
## such that in the machine, T + H = T on the next step
## (H = step size). Solver will continue anyway.
## In above message, R1 = 0, R2 = 0
##
## DINTDY- T (=R1) illegal
## In above message, R1 = 2.8036e-222
##
## T not in interval TCUR - HU (= R1) to TCUR (=R2)
## In above message, R1 = 0, R2 = 0
##
## DINTDY- T (=R1) illegal
## In above message, R1 = 8.66551e-221
##
## T not in interval TCUR - HU (= R1) to TCUR (=R2)
## In above message, R1 = 0, R2 = 0
##
## DLSODA- Trouble in DINTDY. ITASK = I1, TOUT = R1
## In above message, I1 = 1
##
## In above message, R1 = 8.66551e-221
##
Here I've excluding a whole bunch of error messages....
##
## In above message, R1 = 8.66551e-221
##
## GEIGER-fitted comparative model of discrete data
## fitted Q matrix:
## a b
## a -0.7928 0.7928
## b 0.7928 -0.7928
##
## fitted 'lambda' model parameter:
## lambda = 0.508644
##
## model summary:
## log-likelihood = -128.569446
## AIC = 261.138891
## AICc = 261.199805
## free parameters = 2
##
## Convergence diagnostics:
## optimization iterations = 100
## failed iterations = 56
## frequency of best fit = NA
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimates
Also, fitting the λ model to discrete trait data has always struck me as a somewhat peculiar enterprise, so this is posted without any prejudice towards whether this is a good idea, a bad idea, or an average idea in the first place!
