Comparing likelihood functions for estimating TMRCA assuming star phylogeny

starTreeTMRCA.R

invisible('
Comparing likelihood functions for estimating TMRCA assuming star phylogeny
') 

library( lubridate) 

# The data: 
d = data.frame( 
 date = ymd( c("2020/01/16", "2020/01/17", "2019/12/30", "2019/12/30", "2019/12/30", "2019/12/30", "2019/12/30", "2020/01/13", "2020/01/08", "2019/12/30", "2019/12/30", "2019/12/30", "2019/12/24", "2019/12/26", "2019/12/30", "2019/12/30", "2019/12/30", "2020/01/17", "2020/01/15", "2020/01/14") )
 , snps = c(2, 0, 1, 2, 2, 0, 2, 0, 0, 0, 2, 0, 3, 0, 0, 1, 0, 1, 1, 3)
 , x = c(2020.04098360656, 2020.04371584699, 2019.99452054795, 2019.99452054795, 2019.99452054795, 2019.99452054795, 2019.99452054795, 2020.03278688525, 2020.01912568306, 2019.99452054795, 2019.99452054795, 2019.99452054795, 2019.97808219178, 2019.98356164384, 2019.99452054795, 2019.99452054795, 2019.99452054795, 2020.04371584699, 2020.03825136612, 2020.03551912568)
)


#' @param tmrca numeric 
#' @param rate scalar or vector or rates in units of subst / genome / year 
loglik_treelength <- function( tmrca, rate  ){
	blen <- d$x - tmrca 
	if ( any (  blen < 0  ))
		return ( NA )
	trelen <- sum( blen )
	mean( dpois( sum( d$snps ), lambda =  trelen * rate, log = TRUE ) )
}

#' @param tmrca numeric 
#' @param rate scalar or vector or rates in units of subst / genome / year 
loglik_branch <- function( tmrca, rate  ){
	blen <- d$x - tmrca 
	if ( any (  blen < 0  ))
		return ( NA )
	
	mean( sapply( rate, function(r) sum( dpois(  d$snps, lambda = r*blen , log = TRUE ) ) ) )
}


#~ ---------
s <- 29903 # genome length 
rates <- s * seq( 8.3e-4, 0.00109, length = 1e3 )

daterange <- seq( as.Date( '2019-12-05'), as.Date( '2019-12-24') , by = .001) 
xrange <- decimal_date( daterange )

# likelihoods 
l_tl <- sapply( xrange, function(x) loglik_treelength( x, rates ))
l_b <-  sapply( xrange, function(x) loglik_branch ( x, rates ))

# convert to densiy 
d_tl <- exp( l_tl ) / sum( exp( l_tl ))
d_b <- exp( l_b ) / sum( exp( l_b ))

# compute ci 
ci <- function(  dens  ){
	i <- c( max( which( cumsum( dens ) < .025 ))
	  , min( which( cumsum( dens ) >.975 )) )
	daterange[ i ] 
}


plot( daterange, d_tl, type = 'l', col = 'black' , ylab = 'Density', main = 'Black: Tree length likelihood, Red: Branch length likelihood')
abline( v = ci(  d_tl) , col = 'black',  lty =3 )
lines( daterange, d_b, type = 'l' , col = 'red', lty = 1)
abline( v = ci(  d_b) , col = 'red',  lty =3 )

invisible( 
'
The confidence intervals: 
> ci( d_tl ) 
[1] "2019-12-14" "2019-12-23"
> ci( d_b ) 
[1] "2019-12-11" "2019-12-21"

')