emvolz

#' This script simulates genealogies for a "Binary state speciation and extinction model" (BiSSE)
#' Trees are simulated using a structured coalescent framework implemented in the phydynR package
#' In this model, a wild type mutates at a constant rate to the mutant type which has lower transmission fitness. 
#'
#' @author Erik M Volz <erik.volz@gmail.com>
#' @date 14 January 2020 

library( phydynR )
library( ape )
library( ggtree )

emvolz

# |A|
# | |
# | |--B-|
# | |    |
# | |C|  |
# | | |  |D|
# | | |  | |
# a b c  d e

emvolz

#' This version separetely simulates the model producing a trajectory object `tfgy`
#' The likelihood is then computed using this trajectory and the tree 
#' An additional likelihood term compares the simulated trajectory to surveillance data giving the proportion of infeced men who are MSM

#' @param tfgy simuation trajectory produced by demographic.process.model
#' @value proportion of infected men in 2010 that are MSM 
tfgy2prevalenceStat <- function(tfgy)
{
  times <- tfgy[[1]]

emvolz

compute.paf <- function(THETA){
  tfgy <- dm( THETA, x0 = X0 , t0 = 1978, t1 = 2014, integrationMethod='lsoda' )
  fs <- tfgy[[2]]
  times <- tfgy[[1]]
  pafs <- sapply( fs, function(f){
    paf <- rowSums( f)[1:3] ; paf <- paf /sum(paf )
  })
  X <- cbind( times, t(pafs) )
  X[ order(times), ]
}

emvolz

require(ape)
require(treedater)
load( 'Untitled.RData' ) 

s <- 2153922 # sequence length 
bounds  <- data.frame( lower = dates, upper = dates + 1 )
(td0 <- dater( unroot(t), dates, s, minblen = 7/365, estimateSampleTimes = bounds, ncpu = 6))
pb0 <- parboot.treedater( td0 , overrideTempConstraint=FALSE, ncpu = 6)

#~ > td0