Wrapper function to optimize the λ tree transformation for a discrete trait

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!