Nothing
### 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')
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.