Friday, December 26, 2014

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!

3 comments:

  1. Thanks for sharing, nice post!
    Máy đưa võng hay thiet bi dua vong tu dong đang ngày càng “được lòng” nhiều bậc phụ huynh bởi tính năng ưu việt của vong em be tu dong cùng sự tiện ích của vong tu dong. Những lợi ích mà may dua vong mang lại là vô cùng thiết thực. An Thái Sơn luôn cung cấp sản phẩm máy đưa võng hay vong ru tu dong cho be chất lượng nhất TP.HCM, tự hào là địa chỉ bán may dua vong tu dong gia re cho trẻ em.
    Chia sẽ các bạn những cách giúp trẻ nhanh biết đi hay bí quyết giúp trẻ nhanh biết đi, bí quyết làm thế nào để cho bé đi ngủ sớm hay những món ăn giảm cân cho trẻ béo phì, bí quyết giúp trẻ hết biếng ăn hiệu quả, trị chứng mất ngủ ở trẻ em hiệu quả, bí quyết cải thiện làn da cho bé hiệu quả, bí quyết giúp bé ngủ ngon giấc hay những cách chế biến đông trùng hạ thảo nguyên con, chia sẻ cách chữa trị bệnh rụng tóc ở trẻ em hiệu quả, những cách chống nắng cho bé hiệu quả hay thực phẩm giúp giải độc gan cho bạn, thực phẩm giúp nhanh liền sẹo hiệu quả, phòng trị bệnh viêm khớp ở trẻ em an toàn, những thực phẩm bổ não cho trẻ cải thiện trí nhớ, những thực phẩm không nên ăn khi thiếu máu não bạn nên lưu ý, chia sẻ thực phẩm cho người bị rối loạn tiền đình hoặc món ăn chữa bệnh mất ngủ hiệu quả
    Những thực phẩm giúp làm giảm tại http://thucphamlamgiam.blogspot.com/
    Những thực phẩm tốt cho tại http://thucphamtotcho.blogspot.com/
    Những thực phẩm tốt cho da tại http://thucphamtotchoda.blogspot.com/

    ReplyDelete
  2. Hi Liam

    I am getting the following error when I run fitLambda function


    Error in ace(x, rescale(tree, model = "lambda", lambda), type = "discrete", :
    length of phenotypic and of phylogenetic data do not match.

    Cany you plz help me with this error.


    Shiva

    ReplyDelete
  3. Hi Liam

    I am getting the following error when I run fitLambda function


    Error in ace(x, rescale(tree, model = "lambda", lambda), type = "discrete", :
    length of phenotypic and of phylogenetic data do not match.

    Cany you plz help me with this error.


    Shiva

    ReplyDelete