inst/doc/hiv_vignette.R

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

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

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


###################################################
### code chunk number 2: hiv_vignette.Rnw:63-74
###################################################
parms_truth <- list(gamma0 = 1
 , gamma1 = 1/7
 , gamma2 = 1/2
 , mu = 1/30
 , b=.036
 , beta0 = 12./10
 , beta1=3./100
 , beta2=9./100
 , S_0=3000
 , I0_0=1, I1_0=0.01, I2_0=0.01
 , m=3, mm=1) 


###################################################
### code chunk number 3: hiv_vignette.Rnw:88-89
###################################################
INFECTEDNAMES <- c('I0', 'I1', 'I2')


###################################################
### code chunk number 4: hiv_vignette.Rnw:95-101
###################################################
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')
)
rownames(births)=colnames(births)<- INFECTEDNAMES


###################################################
### code chunk number 5: hiv_vignette.Rnw:109-115
###################################################
migrations <- rbind(
  c('0', 'parms$gamma0 * I0', '0'),
  c('0', '0', 'parms$gamma1 * I1'),
  c('0', '0', '0')
)
rownames(migrations)=colnames(migrations) <- INFECTEDNAMES


###################################################
### code chunk number 6: hiv_vignette.Rnw:120-126
###################################################
deaths <- c(
 'parms$mu*I0'
 , 'parms$mu*I1'
 , 'parms$mu*I2 + parms$gamma2 * I2'
)
names(deaths) <- INFECTEDNAMES


###################################################
### code chunk number 7: hiv_vignette.Rnw:131-136
###################################################
nonDemeDynamics <- paste(sep='',
'-parms$mu*S + parms$mu*(S + I0 + I1 + I2)'
 ,'- S * (parms$beta0*I0+parms$beta1*I1+parms$beta2*I2) / (S + I0 + I1 + I2)'
)
names(nonDemeDynamics) <- 'S'


###################################################
### code chunk number 8: hiv_vignette.Rnw:143-155
###################################################
# read the tree
tree <- read.tree(
   system.file('extdata/hivSimulation.nwk', package='rcolgem'))
#~ the sample times are the same, because it is a homochronous sample at 50 years
sampleTimes <- rep(50, length(tree$tip.label))
names(sampleTimes) <- tree$tip.label
# create a tree with dated tips and internal nodes, 
# will infer the sample states from tip labels
bdt <- binaryDatedTree(tree
  , sampleTimes
  , sampleStatesAnnotations=INFECTEDNAMES)
bdt


###################################################
### code chunk number 9: hiv_vignette.Rnw:162-173
###################################################
print(system.time( print(
	coalescent.log.likelihood( bdt
	 , births,  deaths, nonDemeDynamics
	 , t0 = 0
	 , x0=c(I0=1, I1=.01, I2=.01, S = parms_truth$S_0)
	 , migrations = migrations
	 , parms=parms_truth
	 , fgyResolution=1000
	 , integrationMethod='euler'
	)
)))


###################################################
### code chunk number 10: hiv_vignette.Rnw:178-179
###################################################
library(bbmle)


###################################################
### code chunk number 11: hiv_vignette.Rnw:185-204
###################################################
obj_fun <- function(lnbeta0, lnbeta1, lnbeta2, t0)
{
	parms <- parms_truth
	parms$beta0 <- exp(lnbeta0)
	parms$beta1 <- exp(lnbeta1)
	parms$beta2 <- exp(lnbeta2)
	mll <- -coalescent.log.likelihood( bdt
		 , births, deaths, nonDemeDynamics
		 , t0 = t0
		 , x0=c(I0=1, I1=.01, I2=.01, S = parms$S_0)
		 , migrations = migrations
		 , parms=parms
		 , fgyResolution = 1000
		 , integrationMethod = 'rk4'
	)
	# track progress: 
	print(c(mll, exp(c(lnbeta0, lnbeta1, lnbeta2) ), t0) )
	mll
}


###################################################
### code chunk number 12: hiv_vignette.Rnw:209-213 (eval = FALSE)
###################################################
## fit <- mle2(obj_fun
##   , start=list(lnbeta0=log(.6), lnbeta1=log(.2), lnbeta2=log(.05), t0=0)
##   , method='Nelder-Mead', optimizer='optim' 
##   ,  control=list(trace=6, reltol=1e-8))


###################################################
### code chunk number 13: hiv_vignette.Rnw:218-223
###################################################
load( system.file('extdata/hivModel0-fit.RData', package='rcolgem') )
AIC(fit)
logLik(fit)
coef(fit)
exp(coef(fit))


###################################################
### code chunk number 14: results0
###################################################
parms <- parms_truth
parms[c('beta0', 'beta1', 'beta2', 't0')] <- unname(exp(coef(fit)))
parms$t0 <- log(parms$t0)
ox <- solve.model(t0=coef(fit)['t0'], t1=bdt$maxSampleTime, x0=c(I0=1, I1=.01, I2=.01, S = parms_truth$S_0), births,  deaths, nonDemeDynamics, parms, migrations=migrations, integrationMethod = 'rk4')
matplot( rtimes, rY, xlab='Time', ylab='Number infected by stage') 
lines(ox[,1] , ox[, 'I0'], col='black')
lines(ox[,1] , ox[, 'I1'], col='red')
lines(ox[,1] , ox[, 'I2'], col='green')


###################################################
### code chunk number 15: hiv_vignette.Rnw:248-250 (eval = FALSE)
###################################################
## profbeta <- profile(fit, which=1, alpha=.05
## 	, std.err=.5, trace=TRUE, tol.newmin=1 )


###################################################
### code chunk number 16: hiv_vignette.Rnw:253-254
###################################################
load( system.file('extdata/hivModel0-profbeta.RData',  package='rcolgem') )


###################################################
### code chunk number 17: hiv_vignette.Rnw:257-258
###################################################
c( exp( confint( profbeta ) ), TrueVal=parms_truth$beta0 )


###################################################
### code chunk number 18: profplot
###################################################
plot(profbeta)
abline( v = log( parms_truth$beta0) , col='red')


###################################################
### code chunk number 19: proffigure
###################################################
plot(profbeta)
abline( v = log( parms_truth$beta0) , col='red')