inst/doc/hiv_simtree_vignette.R

In rcolgem: statistical inference and modeling of genealogies generated by epidemic and ecological processes

### R code from vignette source 'hiv_simtree_vignette.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: hiv_simtree_vignette.Rnw:58-59
###################################################
library(rcolgem)


###################################################
### code chunk number 2: hiv_simtree_vignette.Rnw:75-76
###################################################
INFECTEDNAMES <- c('I0', 'I1', 'I2')


###################################################
### code chunk number 3: hiv_simtree_vignette.Rnw:82-88
###################################################
births <- rbind( 
  c('beta0 * S * I0 / (S + I0 + I1 + I2)', '0', '0'), 
  c('beta1 * S * I1 / (S + I0 + I1 + I2)', '0', '0'), 
  c('beta2 * S * I2 / (S + I0 + I1 + I2)', '0', '0')
)
rownames(births)=colnames(births)<- INFECTEDNAMES


###################################################
### code chunk number 4: hiv_simtree_vignette.Rnw:94-99
###################################################
## births <- rbind( 
##   c('parms$beta0 * S * I0 / (S + I0 + I1 + I2)', '0', '0'), 
##   c('parms$beta1 * S * I1 / (S + I0 + I1 + I2)', '0', '0'), 
##   c('parms$beta2 * S * I2 / (S + I0 + I1 + I2)', '0', '0')
## )


###################################################
### code chunk number 5: hiv_simtree_vignette.Rnw:105-111
###################################################
migrations <- rbind(
  c('0', 'gamma0 * I0', '0'),
  c('0', '0', 'gamma1 * I1'),
  c('0', '0', '0')
)
rownames(migrations)=colnames(migrations) <- INFECTEDNAMES


###################################################
### code chunk number 6: hiv_simtree_vignette.Rnw:116-122
###################################################
deaths <- c(
 'mu*I0'
 , 'mu*I1'
 , 'mu*I2 + gamma2 * I2'
)
names(deaths) <- INFECTEDNAMES


###################################################
### code chunk number 7: hiv_simtree_vignette.Rnw:127-130
###################################################
nonDemeDynamics <- c( S = '-mu*S + mu*(S + I0 + I1 + I2) -
   S * (beta0*I0+beta1*I1+beta2*I2) / (S + I0 + I1 + I2)'
)


###################################################
### code chunk number 8: hiv_simtree_vignette.Rnw:135-152
###################################################
demo.model <- build.demographic.process( 
  births
  , nonDemeDynamics
  , migrations=migrations
  , deaths=deaths
  , parameterNames = c(
     'beta0'
     , 'beta1'
     , 'beta2'
     , 'gamma0'
     , 'gamma1'
     , 'gamma2'
     , 'mu')
  , rcpp = TRUE
  , sde=TRUE
)
class(demo.model)


###################################################
### code chunk number 9: hiv_simtree_vignette.Rnw:163-164
###################################################
## demo.model(theta, x0, t0, t1, res = 1e3, integrationMethod='adams')


###################################################
### code chunk number 10: hiv_simtree_vignette.Rnw:184-195
###################################################
theta <- c(gamma0 = 1
 , gamma1 = 1/7
 , gamma2 = 1/2
 , mu = 1/30
 , beta0 = 12./10
 , beta1=3./100
 , beta2=9./100
)
t0 <- 0
t1 <- 50
x0 <- c(S = 999, I0 = 1, I1 =.1, I2 = .1)


###################################################
### code chunk number 11: hiv_simtree_vignette.Rnw:200-201
###################################################
show.demographic.process(demo.model, theta, x0, t0, t1)


###################################################
### code chunk number 12: hiv_simtree_vignette.Rnw:207-209
###################################################
n <- 100 
sampleTimes <- seq( 15, 25, length.out = n)


###################################################
### code chunk number 13: hiv_simtree_vignette.Rnw:213-215
###################################################
sampleStates <- t(rmultinom( n, size = 1, prob = c(.025, .9, .075) ))
head(sampleStates)


###################################################
### code chunk number 14: hiv_simtree_vignette.Rnw:219-221
###################################################
tree <- sim.co.tree (theta, demo.model, x0, t0, sampleTimes, sampleStates, res = 1e3)
tree


###################################################
### code chunk number 15: hiv_simtree_vignette.Rnw:225-226
###################################################
plot.phylo(tree)


###################################################
### code chunk number 16: hiv_simtree_vignette.Rnw:228-229
###################################################
ltt.plot(tree)