inst/doc/sir_vignette.R

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

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

###################################################
### code chunk number 1: sir_vignette.Rnw:25-27
###################################################
library(rcolgem)
tree <- read.tree(system.file( 'extdata/sirModel0.nwk', package='rcolgem'))


###################################################
### code chunk number 2: sir_vignette.Rnw:30-32
###################################################
library(rjson)
epidata <- fromJSON(file=system.file('extdata/sirModel0.json', package='rcolgem')) 


###################################################
### code chunk number 3: sir_vignette.Rnw:71-72
###################################################
parms_truth <- list( beta = .00020002, gamma = 1, S0 = 9999, t0 = 0 )


###################################################
### code chunk number 4: sir_vignette.Rnw:77-81
###################################################
sampleTimes <- rep(12, 75)
names(sampleTimes) <- tree$tip.label
bdt <- binaryDatedTree( tree, sampleTimes=sampleTimes) 
bdt


###################################################
### code chunk number 5: sir_vignette.Rnw:89-92
###################################################
births <- c( I = 'parms$beta * S * I' )
deaths <- c( I = 'parms$gamma * I' )
nonDemeDynamics <- c(S = '-parms$beta * S * I')


###################################################
### code chunk number 6: sir_vignette.Rnw:102-104
###################################################
x0 <- c(I=1, S= unname(parms_truth$S0) )
t0 <- bdt$maxSampleTime - max(bdt$heights) -1 


###################################################
### code chunk number 7: sir_vignette.Rnw:109-120
###################################################
print( 
system.time( 
print(
  coalescent.log.likelihood(
    bdt
    , births, deaths, nonDemeDynamics
    , t0, x0
    , parms=parms_truth
    , fgyResolution = 1000
    , integrationMethod = 'rk4')
)))


###################################################
### code chunk number 8: sir_vignette.Rnw:125-126
###################################################
	library(bbmle)


###################################################
### code chunk number 9: sir_vignette.Rnw:129-146
###################################################
obj_fun <- function(lnbeta, lnI0)
{
	beta <- exp(lnbeta)
	I0 <- exp(lnI0)
	parms <- parms_truth
	parms$beta <- beta
	x0 <- c(I=unname(I0), S = unname(parms$S0) )
	mll <- -coalescent.log.likelihood(
		bdt
		, births, deaths, nonDemeDynamics
		,  t0, x0
		, parms=parms
		, fgyResolution = 1000
		, integrationMethod = 'rk4')
	print(paste(mll, beta, I0))
	mll
}


###################################################
### code chunk number 10: sir_vignette.Rnw:151-158 (eval = FALSE)
###################################################
## fit <- mle2(
##   obj_fun
##   , start=list(lnbeta=log(parms_truth$beta*.75), lnI0=log(1))
##   , method='Nelder-Mead'
##   , optimizer='optim' 
##   , control=list(trace=6, reltol=1e-8)
## )


###################################################
### code chunk number 11: sir_vignette.Rnw:163-170
###################################################
	load( system.file('extdata/sirModel0-fit.RData', package='rcolgem') )
	AIC(fit)
	logLik(fit)
	coef(fit)
	exp(coef(fit))
	# how biased is the estimate? 
	exp(coef(fit)['lnbeta']) - parms_truth$beta


###################################################
### code chunk number 12: results0plot
###################################################
beta <- exp(coef(fit)['lnbeta'])
I0 <- exp(coef(fit)['lnI0'])
parms <- parms_truth
parms$beta <- beta
x0 <- c(I=unname(I0), S = unname(parms$S0) )

o <- solve.model.unstructured(t0,bdt$maxSampleTime, x0
   , births
   ,  deaths
   , nonDemeDynamics, parms)
otruth <- solve.model.unstructured(t0, bdt$maxSampleTime, x0
  , births
  ,  deaths
  , nonDemeDynamics, parms_truth)
plot(epidata$t, epidata$I, type='line'
  , ylim=c(0, 100+max(max(o[,2]),max(epidata$I))))
lines(o[,1], o[,2], col='red' )
lines(otruth[,1], otruth[,2], col='blue' )


###################################################
### code chunk number 13: results0
###################################################
beta <- exp(coef(fit)['lnbeta'])
I0 <- exp(coef(fit)['lnI0'])
parms <- parms_truth
parms$beta <- beta
x0 <- c(I=unname(I0), S = unname(parms$S0) )

o <- solve.model.unstructured(t0,bdt$maxSampleTime, x0
   , births
   ,  deaths
   , nonDemeDynamics, parms)
otruth <- solve.model.unstructured(t0, bdt$maxSampleTime, x0
  , births
  ,  deaths
  , nonDemeDynamics, parms_truth)
plot(epidata$t, epidata$I, type='line'
  , ylim=c(0, 100+max(max(o[,2]),max(epidata$I))))
lines(o[,1], o[,2], col='red' )
lines(otruth[,1], otruth[,2], col='blue' )


###################################################
### code chunk number 14: sir_vignette.Rnw:205-207 (eval = FALSE)
###################################################
## profbeta <- profile(fit, which=1, alpha=.05
##   , std.err=1, trace=TRUE, tol.newmin=1 )


###################################################
### code chunk number 15: sir_vignette.Rnw:210-211
###################################################
load( system.file('extdata/sirModel0-profbeta.RData',  package='rcolgem') )


###################################################
### code chunk number 16: sir_vignette.Rnw:214-215
###################################################
c( exp( confint( profbeta ) ), TrueVal=parms_truth$beta )


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


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