expoTree (version 1.0.1)

expoTree-package: Calculate density dependent likelihood of a phylogenetic tree

Description

Calculates the density dependent likelihood of a phylogenetic tree. It takes branching and sampling times as an argument and integrates the likelihoood function over the whole tree.

Details

Package:

expoTree

Type:

Package

Version:

1.0.1

Date:

2013-09-03

License:

BSD 3 clause

References

Leventhal, Guenthard, Bonhoeffer & Stadler, (2013) "Using an epidemiological model for phylogenetic inference reveals density-dependence in HIV transmission".

Al-Mohy & Higham (2011) "Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators". SIAM Journal on Scientific Computing, 33.

Cheng & Highami (2001) "Implementation for LAPACK of a Block Algorithm for Matrix 1-Norm Estimation". Numerical Analysis Report No. 393, Manchester Centre for Computational Mathematics, Manchester, England, August 2001, and LAPACK Working Note 152.

Blackford et al. (1997) "ScaLAPACK Users' Guide". Society for Industrial and Applied Mathematics, Philadelphia, PA.

Examples

  # simulate trees
  N <- 15
  beta <- 1
  mu <- 0.1
  psi <- 0.1
  rho <- 0
  nsamp <- 20
  epis <- lapply(1:2,function(x) sim.epi(N,beta,mu,psi,nsamp))
  forest <- lapply(epis,function(x) cbind(x$times,x$ttypes))

  # plot lineages-through-time
  plotLTT(forest)

  extant <- sapply(forest,function(t) sum(2*t[,2]-1))
  lineages <- lapply(forest,function(t) sum(2*t[,2]-1)+cumsum(1-2*t[,2]))
  max.lineages <- sapply(lineages,max)

  # calculate likelihood for the forest
  lik <- sapply(forest,function(tree) {
    runExpoTree(matrix(c(N,beta,mu,psi,rho),nrow=1),tree[,1],tree[,2])
  })
  cat("Likelihood = ",sum(lik),"\n")

  if (! is.nan(sum(lik))) {
    # Find optimal parameters
    (opt <- expoTree.optim(forest,pars=c(N,beta,psi),
                   lo=c(max(max.lineages),0,0),hi=c(50,3,1),
                   fix=c(FALSE,FALSE,TRUE,FALSE,TRUE),
                   fix.val=c(0,0,-psi/(psi+mu),0,0)))
  }