emvolz

library( deSolve )
library( ggplot2 )
library( cowplot )

dy <- function( t, y, parms, ... )
{
	S = y[1] 
	I0 = y[2] 
	I1 = y[3] 
	R = y[4] 

emvolz

#~ '
#~ using sim_replicate3 -- includes breakthrough & updated models
#~ Based on situation 17 sep 2021
#~ Including hier Bayes 
#~ smaller rho0
#~ s = +25%
#~ '

selcoef <- .25 
VAR_IMPORTS <- 650

emvolz

'
v3:
- using june 19 r1 data
- deltrans clusters n > 40 
- fixing clock rate on rtt consensus

'
path = '2020-06-19.r1/deltrans_trees'
mdfn = '2020-06-19.r1/metadata/master.csv'

emvolz

# Note this depends on some convenience functions in the 'sarscov2' package available at 
# https://github.com/emvolz-phylodynamics/sarscov2Rutils
library( sarscov2 )
library( ape ) 
library( lubridate )

# copy over the ref sequence 
file.copy( system.file( package='sarscov2', 'extdata/ref.fasta' ) , '.' ) 

# align 

emvolz

library( sarscov2 )
library( ape )

res = 10 
tau0 = 30 / 365 / 36.5^2
ncpu = 6

# define some parameters 
## path to the GISAID alignment: 
fn = 'may18_filter1.fas' 

emvolz

library( sarscov2 )
library( ape )

# define some parameters 
## path to the GISAID alignment, distance file, metadata: 
fn = 'may10_filter1.fas' 
distfn = 'dist.txt'
mdfn = 'may10_dataC.csv'

## A simple name for the region:

emvolz

library( sarscov2 )
library( ape )

# define some parameters 
## path to the GISAID alignment: 
fn = 'large_alignment_may5_Erik_Brazil.fasta' 
## A simple name for the region:
region='brazil'
## A path to the file containing small genetic distance pairs
distfn = 'large_alignment_may5_Erik_Brazil.dst'

emvolz

<!--
transmission
-->
<matrixeq id="MatrixEquation.2" spec="phydyn.model.MatrixEquation" destination="E" origin="Il" type="birth">beta*Il*S / ( S + E + Il + Ih +R )</matrixeq>
<matrixeq id="MatrixEquation.3" spec="phydyn.model.MatrixEquation" destination="E" origin="Ih" type="birth">beta*Ih*tau * S / (S + E + Il + Ih +R)</matrixeq>
<!--
stage progression
-->
<matrixeq id="MatrixEquation.0" spec="phydyn.model.MatrixEquation" destination="Il" origin="E" type="migration">if ( t &gt; SEIR_START ) then gamma0*E*(1-p_h) else 0.0</matrixeq>
<matrixeq id="MatrixEquation.1" spec="phydyn.model.MatrixEquation" destination="Ih" origin="E" type="migration">if ( t &gt; SEIR_START ) then gamma0*E*p_h else 0.0</matrixeq>

emvolz

invisible('
Extinction probability as a function of 
- Reproduction number R
- Dispersion k
- Number imported infections

Assumes negative binomial offspring distribution
') 

# negative binomial probability generating function 

emvolz

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)