Title: | An Implementation of the TSIR Model |
---|---|
Description: | An implementation of the time-series Susceptible-Infected-Recovered (TSIR) model using a number of different fitting options for infectious disease time series data. The manuscript based on this package can be found here <doi:10.1371/journal.pone.0185528>. The method implemented here is described by Finkenstadt and Grenfell (2000) <doi:10.1111/1467-9876.00187>. |
Authors: | Alexander D. Becker [aut, cre], Sinead E. Morris [ctb], Ottar N. Bjornstad [ctb] |
Maintainer: | Alexander D. Becker <[email protected]> |
License: | GPL-3 |
Version: | 0.4.3 |
Built: | 2024-10-31 22:10:32 UTC |
Source: | https://github.com/cran/tsiR |
Plot the correlation of the true data against the fitted resimulated data.
corr(sim)
corr(sim)
sim |
The dataframe or list produced by the 'runtsir' function. |
This function computes an 8 point derivative.
derivative(X, Y)
derivative(X, Y)
X |
The variable to differentiate with respect to. |
Y |
The function / vector to differentiate. |
The times at which we declare a new outbreak has started based on the threshold parameter.
epitimes(data, threshold, epi.length = 3)
epitimes(data, threshold, epi.length = 3)
data |
The inputed data frame with the cases vector. This is the same data you put into runtsir. |
threshold |
The required number of cases observed to declare it an outbreak. |
epi.length |
The required duration (in 52/IP weeks) to declare it an outbreak. |
This function computes the set up to run the TSIR model, i.e. reconstructs susecptibles and estimates beta and alpha. This can be plugged into simulatetsir.
estpars( data, xreg = "cumcases", IP = 2, seasonality = "standard", regtype = "gaussian", sigmamax = 3, family = "gaussian", link = "identity", userYhat = numeric(), alpha = NULL, sbar = NULL, printon = F )
estpars( data, xreg = "cumcases", IP = 2, seasonality = "standard", regtype = "gaussian", sigmamax = 3, family = "gaussian", link = "identity", userYhat = numeric(), alpha = NULL, sbar = NULL, printon = F )
data |
The data frame containing cases and interpolated births and populations. |
xreg |
The x-axis for the regression. Options are 'cumcases' and 'cumbirths'. Defaults to 'cumcases'. |
IP |
The infectious period in weeks. This should be the same as your timestep. Defaults to 2 weeks. |
seasonality |
The type of contact to use. Options are standard for 52/IP point contact or schoolterm for just a two point on off contact or none for a single contact parameter. Defaults to standard. |
regtype |
The type of regression used in susceptible reconstruction. Options are 'gaussian', 'lm' (linear model), 'spline' (smooth.spline with 2.5 degrees freedom), 'lowess' (with f = 2/3, iter = 1), 'loess' (degree 1), and 'user' which is just a user inputed vector. Defaults to 'gaussian' and if that fails then defaults to loess. |
sigmamax |
The inverse kernal width for the gaussian regression. Default is 3. Smaller, stochastic outbreaks tend to need a lower sigma. |
family |
The family in the GLM regression. One can use any of the GLM ones, but the options are essentially 'poisson' (with link='log'), 'gaussian' (with link='log' or 'identity'), or 'quasipoisson' (with link='log'). Default is 'gaussian'. |
link |
The link function used with the glm family. Options are link='log' or 'identity'. Default is 'identity'. to include some bayesian approaches. For 'bayesglm' we use a gaussian prior with mean 1e-4. |
userYhat |
The inputed regression vector if regtype='user'. Defaults to NULL. |
alpha |
The mixing parameter. Defaults to NULL, i.e. the function estimates alpha. |
sbar |
The mean number of susceptibles. Defaults to NULL, i.e. the function estimates sbar. |
printon |
Whether to show diagnostic prints or not, defaults to FALSE. |
## Not run: require(kernlab) London <- twentymeas[["London"]] parms <- estpars(London) names(parms) sim <- simulatetsir(London,parms=parms,inits.fit=FALSE) plotres(sim) ## End(Not run)
## Not run: require(kernlab) London <- twentymeas[["London"]] parms <- estpars(London) names(parms) sim <- simulatetsir(London,parms=parms,inits.fit=FALSE) plotres(sim) ## End(Not run)
Used internally to filter jags results to give just the inference well use.
jagsfilter(mcmcresults)
jagsfilter(mcmcresults)
mcmcresults |
is the input from the jags model. |
Plot the correlation of the true data against the fitted resimulated data.
logcorr(sim)
logcorr(sim)
sim |
The dataframe or list produced by the 'runtsir' function. |
A function used to optimize the threshold parameter to give the best fit to the data. Optimizes the fit based on R squared.
maxthreshold( data, nsim = 2, IP = 2, method = "deterministic", inits.fit = FALSE, parms, thresholdmin = 2, thresholdmax = 20, printon = FALSE )
maxthreshold( data, nsim = 2, IP = 2, method = "deterministic", inits.fit = FALSE, parms, thresholdmin = 2, thresholdmax = 20, printon = FALSE )
data |
The time, cases, births, pop data frame. |
nsim |
The number of simulations to do. |
IP |
The infectious period, which should the time step of the data. |
method |
The forward simulation method used, i.e. deterministic, negbin, pois. |
inits.fit |
Whether or not to fit initial conditions as well. Defaults to FALSE here. This parameter is more necessary in more chaotic locations. |
parms |
The estimated parameters from estpars or mcmcestpars. |
thresholdmin |
The minimum number of cases to be considered an outbreak. |
thresholdmax |
The max number of cases to be considered an outbreak. |
printon |
A T/F statement to print the progress. |
require(kernlab) Mold <- twentymeas[["Mold"]] plotdata(Mold) ## Not run: parms <- estpars(data=Mold,alpha=0.97) tau <- maxthreshold(data=Mold,parms=parms, thresholdmin=8,thresholdmax=12,inits.fit=FALSE) res <- simulatetsir(data=Mold,parms=parms, epidemics='break',threshold=tau,method='negbin',inits.fit=FALSE) plotres(res) ## End(Not run)
require(kernlab) Mold <- twentymeas[["Mold"]] plotdata(Mold) ## Not run: parms <- estpars(data=Mold,alpha=0.97) tau <- maxthreshold(data=Mold,parms=parms, thresholdmin=8,thresholdmax=12,inits.fit=FALSE) res <- simulatetsir(data=Mold,parms=parms, epidemics='break',threshold=tau,method='negbin',inits.fit=FALSE) plotres(res) ## End(Not run)
This function computes the set up to run the TSIR model, i.e. reconstructes susecptibles and estimates beta and alpha using MCMC computations. Used the same way as estpars.
mcmcestpars( data, xreg = "cumcases", IP = 2, regtype = "gaussian", sigmamax = 3, seasonality = "standard", userYhat = numeric(), update.iter = 10000, n.iter = 30000, n.chains = 3, n.adapt = 1000, burn.in = 100, sbar = NULL, alpha = NULL, printon = F )
mcmcestpars( data, xreg = "cumcases", IP = 2, regtype = "gaussian", sigmamax = 3, seasonality = "standard", userYhat = numeric(), update.iter = 10000, n.iter = 30000, n.chains = 3, n.adapt = 1000, burn.in = 100, sbar = NULL, alpha = NULL, printon = F )
data |
The data frame containing cases and interpolated births and populations. |
xreg |
The x-axis for the regression. Options are 'cumcases' and 'cumbirths'. Defaults to 'cumcases'. |
IP |
The infectious period in weeks. Defaults to 2 weeks. |
regtype |
The type of regression used in susceptible reconstruction. Options are 'gaussian', 'lm' (linear model), 'spline' (smooth.spline with 2.5 degrees freedom), 'lowess' (with f = 2/3, iter = 1), 'loess' (degree 1), and 'user' which is just a user inputed vector. Defaults to 'gaussian' and if that fails then defaults to loess. |
sigmamax |
The inverse kernal width for the gaussian regression. Default is 3. Smaller, stochastic outbreaks tend to need a lower sigma. |
seasonality |
The type of contact to use. Options are standard for 52/IP point contact or schoolterm for just a two point on off contact or none for a single contact parameter. Defaults to standard. |
userYhat |
The inputed regression vector if regtype='user'. Defaults to NULL. |
update.iter |
Number of MCMC iterations to use in the update aspect. Default is 10000. |
n.iter |
Number of MCMC iterations to use. Default is 30000. |
n.chains |
Number of MCMC chains to use. Default is 3. |
n.adapt |
Adaptive number for MCMC. Default is 1000. |
burn.in |
Burn in number. Default is 100. |
sbar |
The mean number of susceptibles. Defaults to NULL, i.e. the function estimates sbar. |
alpha |
The mixing parameter. Defaults to NULL, i.e. the function estimates alpha. |
printon |
Whether to show diagnostic prints or not, defaults to FALSE. |
This function runs the TSIR model using a MCMC estimation. The susceptibles are still reconstructed in the same way as the regular tsir model, however beta, alpha, and sbar (or whatever combination you enter) are estimated using rjargs.
mcmctsir( data, xreg = "cumcases", IP = 2, nsim = 100, regtype = "gaussian", sigmamax = 3, userYhat = numeric(), update.iter = 10000, n.iter = 30000, n.chains = 3, n.adapt = 1000, burn.in = 100, method = "deterministic", epidemics = "cont", pred = "forward", seasonality = "standard", inits.fit = FALSE, threshold = 1, sbar = NULL, alpha = NULL, add.noise.sd = 0, mul.noise.sd = 0, printon = F )
mcmctsir( data, xreg = "cumcases", IP = 2, nsim = 100, regtype = "gaussian", sigmamax = 3, userYhat = numeric(), update.iter = 10000, n.iter = 30000, n.chains = 3, n.adapt = 1000, burn.in = 100, method = "deterministic", epidemics = "cont", pred = "forward", seasonality = "standard", inits.fit = FALSE, threshold = 1, sbar = NULL, alpha = NULL, add.noise.sd = 0, mul.noise.sd = 0, printon = F )
data |
The data frame containing cases and interpolated births and populations. |
xreg |
The x-axis for the regression. Options are 'cumcases' and 'cumbirths'. Defaults to 'cumcases'. |
IP |
The infectious period in weeks. Defaults to 2 weeks. |
nsim |
The number of simulations to do. Defaults to 100. |
regtype |
The type of regression used in susceptible reconstruction. Options are 'gaussian', 'lm' (linear model), 'spline' (smooth.spline with 2.5 degrees freedom), 'lowess' (with f = 2/3, iter = 1), 'loess' (degree 1), and 'user' which is just a user inputed vector. Defaults to 'gaussian' and if that fails then defaults to loess. |
sigmamax |
The inverse kernal width for the gaussian regression. Default is 3. Smaller, stochastic outbreaks tend to need a lower sigma. |
userYhat |
The inputed regression vector if regtype='user'. Defaults to NULL. |
update.iter |
Number of MCMC iterations to use in the update aspect. Default is 10000. |
n.iter |
Number of MCMC iterations to use. Default is 30000. |
n.chains |
Number of MCMC chains to use. Default is 3. |
n.adapt |
Adaptive number for MCMC. Default is 1000. |
burn.in |
Burn in number. Default is 100. |
method |
The type of next step prediction used. Options are 'negbin' for negative binomial, 'pois' for poisson distribution, and 'deterministic'. Defaults to 'deterministic'. |
epidemics |
The type of data splitting. Options are 'cont' which doesn't split the data up at all, and 'break' which breaks the epidemics up if there are a lot of zeros. Defaults to 'cont'. |
pred |
The type of prediction used. Options are 'forward' and 'step-ahead'. Defaults to 'forward'. |
seasonality |
The type of contact to use. Options are standard for 52/IP point contact or schoolterm for just a two point on off contact or none for a single contact parameter. Defaults to standard. |
inits.fit |
Whether or not to fit initial conditions using simple least squares as well. Defaults to FALSE. This parameter is more necessary in more chaotic locations. |
threshold |
The cut off for a new epidemic if epidemics = 'break'. Defaults to 1. |
sbar |
The mean number of susceptibles. Defaults to NULL, i.e. the function estimates sbar. |
alpha |
The mixing parameter. Defaults to NULL, i.e. the function estimates alpha. |
add.noise.sd |
The sd for additive noise, defaults to zero. |
mul.noise.sd |
The sd for multiplicative noise, defaults to zero. |
printon |
Whether to show diagnostic prints or not, defaults to FALSE. |
Plots the inferred beta with confidence intervals (when they can be calculated)
plotbeta(dat)
plotbeta(dat)
dat |
the list produced from the runtsir, mcmctsir, and simulatetsir function. |
Plots the cases data with a line whenever the forward simulation is seeded using the real data.
plotbreaks(data, threshold)
plotbreaks(data, threshold)
data |
Data frame with the cases vector. |
threshold |
The epidemic threshold, i.e. the number of cases required to spark a new outbreak in the model. |
Plots just the cases data.
plotcases(data)
plotcases(data)
data |
The data frame with cases. |
Plots just the comparison of the forward simulation fit to the data.
plotcomp(sim, errtype = "95", max.plot = 10)
plotcomp(sim, errtype = "95", max.plot = 10)
sim |
is list produced by runtsir or mcmctsir |
errtype |
is the type of error bands to show. Defaults to '95' for 95 percent CI, the other option is 'sd' to standard deviation. |
max.plot |
the number of individual stochastic simulations to plot. Defaults to 10. |
Plots the cases data as well as birth and population dynamics.
plotdata(data)
plotdata(data)
data |
The dataframe with time, cases, births, and pop. |
Plots the forward simulation from the TSIR model
plotforward(dat, inverse = F)
plotforward(dat, inverse = F)
dat |
the list produced from the runtsir, mcmctsir, and simulatetsir function. |
inverse |
a TRUE or FALSE option to plot the forward simulate negative (TRUE) or positive (FALSE). Defaults to FALSE. |
Function to plot the Local Lyapunov Exponents. The output is of class ggplot2 so you can add standard ggplot2 options to it if desired.
plotLLE(LLE)
plotLLE(LLE)
LLE |
The output from TSIR_LLE |
## Not run: require(kernlab) require(ggplot2) require(kernlab) London <- twentymeas$London ## just analyze the biennial portion of the data London <- subset(London, time > 1950) ## define the interval to be 2 weeks IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log') ## look at beta and alpha estimate plotbeta(parms) ## simulate the fitted parameters sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2) ## now lets predict forward 200 years using the mean birth rate, ## starting from rough initial conditions times <- seq(1965,2165, by = 1/ (52/IP)) births <- rep(mean(London$births),length(times)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## take the last 10 years pred <- lapply(pred, function(x) tail(x, 52/IP * 20) ) ## now compute the Lyapunov Exponent for the simulate and predicted model simLE <- TSIR_LE( time=sim$res$time, S=sim$simS$mean, I=sim$res$mean, alpha=sim$alpha, beta=sim$contact$beta, IP=IP ) predLE <- TSIR_LE( time=pred$I$time, S=pred$S$X3, I=pred$I$X3, alpha=parms$alpha, beta=parms$contact$beta, IP=IP ) simLE$LE predLE$LE simLLE <- TSIR_LLE(simLE) predLLE <- TSIR_LLE(predLE) plotLLE(simLLE) plotLLE(predLLE) ## End(Not run)
## Not run: require(kernlab) require(ggplot2) require(kernlab) London <- twentymeas$London ## just analyze the biennial portion of the data London <- subset(London, time > 1950) ## define the interval to be 2 weeks IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log') ## look at beta and alpha estimate plotbeta(parms) ## simulate the fitted parameters sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2) ## now lets predict forward 200 years using the mean birth rate, ## starting from rough initial conditions times <- seq(1965,2165, by = 1/ (52/IP)) births <- rep(mean(London$births),length(times)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## take the last 10 years pred <- lapply(pred, function(x) tail(x, 52/IP * 20) ) ## now compute the Lyapunov Exponent for the simulate and predicted model simLE <- TSIR_LE( time=sim$res$time, S=sim$simS$mean, I=sim$res$mean, alpha=sim$alpha, beta=sim$contact$beta, IP=IP ) predLE <- TSIR_LE( time=pred$I$time, S=pred$S$X3, I=pred$I$X3, alpha=parms$alpha, beta=parms$contact$beta, IP=IP ) simLE$LE predLE$LE simLLE <- TSIR_LLE(simLE) predLLE <- TSIR_LLE(predLE) plotLLE(simLLE) plotLLE(predLLE) ## End(Not run)
Plots the cumulative cases - cumulative births data and regression fit
plotregression(dat)
plotregression(dat)
dat |
the list produced from the runtsir, mcmctsir, and simulatetsir function. |
Plots diagnostics and results of the runtsir function.
plotres(dat, max.plot = 10)
plotres(dat, max.plot = 10)
dat |
the list produced from the runtsir, mcmctsir, and simulatetsir function. |
max.plot |
the number of individual stochastic simulations to plot. Defaults to 10. |
Plots the inferred reporting rate, rho
plotrho(dat)
plotrho(dat)
dat |
the list produced from the runtsir, mcmctsir, and simulatetsir function. |
Plots the profile log likelihood calculation for inferred sbar
plotsbar(dat)
plotsbar(dat)
dat |
the list produced from the runtsir, mcmctsir, and simulatetsir function. |
function to predict incidence and susceptibles using the tsir model. This is different than simulatetsir as you are inputting parameters as vectors. The output is a data frame I and S with mean and confidence intervals of predictions.
predicttsir(times, births, beta, alpha, S0, I0, nsim, stochastic)
predicttsir(times, births, beta, alpha, S0, I0, nsim, stochastic)
times |
The time vector to predict the model from. This assumes that the time step is equal to IP |
births |
The birth vector (of length length(times) or a single element) where each element is the births in that given (52/IP) time step |
beta |
The length(52/IP) beta vector of contact. |
alpha |
A single numeric which acts as the homogeniety parameter. |
S0 |
The starting initial condition for S. This should be greater than one, i.e. not a fraction. |
I0 |
The starting initial condition for I. This should be greater than one, i.e. not a fraction. |
nsim |
The number of simulations to perform. |
stochastic |
A TRUE / FALSE argument where FALSE is the deterministic model, and TRUE is a negative binomial distribution. |
## Not run: require(kernlab) require(ggplot2) require(kernlab) require(tsiR) London <- twentymeas$London London <- subset(London, time > 1950) IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian') plotbeta(parms) ## now lets predict forward 20 years using the mean birth rate, ## starting from rough initial conditions births <- min(London$births) times <- seq(1965,1985, by = 1/ (52/IP)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## plot this prediction ggplot(pred$I,aes(time,mean))+geom_line()+geom_ribbon(aes(ymin=low,ymax=high),alpha=0.3) ## End(Not run)
## Not run: require(kernlab) require(ggplot2) require(kernlab) require(tsiR) London <- twentymeas$London London <- subset(London, time > 1950) IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian') plotbeta(parms) ## now lets predict forward 20 years using the mean birth rate, ## starting from rough initial conditions births <- min(London$births) times <- seq(1965,1985, by = 1/ (52/IP)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## plot this prediction ggplot(pred$I,aes(time,mean))+geom_line()+geom_ribbon(aes(ymin=low,ymax=high),alpha=0.3) ## End(Not run)
computes the residuals for when X is the cumulative births. Used internally.
residual.births(rho, Yhat, Y)
residual.births(rho, Yhat, Y)
rho |
The reporting rate, used to get units correct. |
Yhat |
The fitted regression line. |
Y |
The cumulative cases. |
Computes the residuals for when X is the cumulative cases. Used internally.
residual.cases(Yhat, Y)
residual.cases(Yhat, Y)
Yhat |
The fitted regression line. |
Y |
The cumulative births. |
This function runs the TSIR model.
runtsir( data, xreg = "cumcases", IP = 2, nsim = 10, regtype = "gaussian", sigmamax = 3, userYhat = numeric(), alpha = NULL, sbar = NULL, family = "gaussian", link = "identity", method = "deterministic", inits.fit = FALSE, epidemics = "cont", pred = "forward", threshold = 1, seasonality = "standard", add.noise.sd = 0, mul.noise.sd = 0, printon = F, fit = NULL, fittype = NULL )
runtsir( data, xreg = "cumcases", IP = 2, nsim = 10, regtype = "gaussian", sigmamax = 3, userYhat = numeric(), alpha = NULL, sbar = NULL, family = "gaussian", link = "identity", method = "deterministic", inits.fit = FALSE, epidemics = "cont", pred = "forward", threshold = 1, seasonality = "standard", add.noise.sd = 0, mul.noise.sd = 0, printon = F, fit = NULL, fittype = NULL )
data |
The data frame containing cases and interpolated births and populations. |
xreg |
The x-axis for the regression. Options are 'cumcases' and 'cumbirths'. Defaults to 'cumcases'. |
IP |
The infectious period in weeks. Defaults to 2 weeks. |
nsim |
The number of simulations to do. Defaults to 100. |
regtype |
The type of regression used in susceptible reconstruction. Options are 'gaussian', 'lm' (linear model), 'spline' (smooth.spline with 2.5 degrees freedom), 'lowess' (with f = 2/3, iter = 1), 'loess' (degree 1), and 'user' which is just a user inputed vector. Defaults to 'gaussian' and if that fails then defaults to loess. |
sigmamax |
The inverse kernal width for the gaussian regression. Default is 3. Smaller, stochastic outbreaks tend to need a lower sigma. |
userYhat |
The inputed regression vector if regtype='user'. Defaults to NULL. |
alpha |
The mixing parameter. Defaults to NULL, i.e. the function estimates alpha. |
sbar |
The mean number of susceptibles. Defaults to NULL, i.e. the function estimates sbar. |
family |
The family in the GLM regression. One can use any of the GLM ones, but the options are essentially 'poisson' (with link='log'), 'gaussian' (with link='log' or 'identity'), or 'quasipoisson' (with link='log'). Default is 'gaussian'. |
link |
The link function used with the glm family. Options are link='log' or 'identity'. Default is 'identity'. |
method |
The type of next step prediction used. Options are 'negbin' for negative binomial, 'pois' for poisson distribution, and 'deterministic'. Defaults to 'deterministic'. |
inits.fit |
Whether or not to fit initial conditions using simple least squares as well. Defaults to FALSE. This parameter is more necessary in more chaotic locations. |
epidemics |
The type of data splitting. Options are 'cont' which doesn't split the data up at all, and 'break' which breaks the epidemics up if there are a lot of zeros. Defaults to 'cont'. |
pred |
The type of prediction used. Options are 'forward' and 'step-ahead'. Defaults to 'forward'. |
threshold |
The cut off for a new epidemic if epidemics = 'break'. Defaults to 1. |
seasonality |
The type of contact to use. Options are standard for 52/IP point contact or schoolterm for just a two point on off contact, or none for a single contact parameter. Defaults to standard. |
add.noise.sd |
The sd for additive noise, defaults to zero. |
mul.noise.sd |
The sd for multiplicative noise, defaults to zero. |
printon |
Whether to show diagnostic prints or not, defaults to FALSE. |
fit |
Now removed but gives a warning. |
fittype |
Now removed but gives a warning. |
require(kernlab) London <- twentymeas[["London"]] ## Not run: plotdata(London) res <- runtsir(data=London,method='pois',nsim=10, IP=2,inits.fit=FALSE) plotres(res) ## End(Not run)
require(kernlab) London <- twentymeas[["London"]] ## Not run: plotdata(London) res <- runtsir(data=London,method='pois',nsim=10, IP=2,inits.fit=FALSE) plotres(res) ## End(Not run)
This function just simulates the forward prediction given the data and a parms list generated from estpars or mcmcestpars.
simulatetsir( data, nsim = 100, IP = 2, parms, method = "deterministic", epidemics = "cont", pred = "forward", threshold = 1, inits.fit = FALSE, add.noise.sd = 0, mul.noise.sd = 0 )
simulatetsir( data, nsim = 100, IP = 2, parms, method = "deterministic", epidemics = "cont", pred = "forward", threshold = 1, inits.fit = FALSE, add.noise.sd = 0, mul.noise.sd = 0 )
data |
The data frame containing cases and interpolated births and populations. |
nsim |
The number of simulations to do. Defaults to 100. |
IP |
The infectious period. Defaults to 2. |
parms |
Either the parameters estimated by estpars or mcmcestpars, or a list containing beta, rho, Z, sbar, alpha, X, Y, Yhat, contact, alphalow, alphahigh, loglik, pop vectors. |
method |
The type of next step prediction used. Options are 'negbin' for negative binomial, 'pois' for poisson distribution, and 'deterministic'. Defaults to 'deterministic'. |
epidemics |
The type of data splitting. Options are 'cont' which doesn't split the data up at all, and 'break' which breaks the epidemics up if there are a lot of zeros. Defaults to 'cont'. |
pred |
The type of prediction used. Options are 'forward' and 'step-ahead'. Defaults to 'forward'. |
threshold |
The cut off for a new epidemic if epidemics = 'break'. Defaults to 1. |
inits.fit |
Whether or not to fit initial conditions using simple least squares as well. Defaults to FALSE. This parameter is more necessary in more chaotic locations. |
add.noise.sd |
The sd for additive noise, defaults to zero. |
mul.noise.sd |
The sd for multiplicative noise, defaults to zero. |
A function to calculate the Lyapunov Exponennt (LE) from the TSIR model
TSIR_LE(time, S, I, alpha, beta, IP)
TSIR_LE(time, S, I, alpha, beta, IP)
time |
The time vector from the data or simulated data |
S |
The S output from the simulated or predicted TSIR model |
I |
The I output from the simulated or predicted TSIR model |
alpha |
The homogeneity parameter from the simulated or predicted TSIR model |
beta |
The inferred contact rate, use beta = contact$beta where contact is an output from runtsir or simulatetsir |
IP |
The generation interval of the pathogen (in weeks) |
## Not run: require(kernlab) require(ggplot2) require(kernlab) London <- twentymeas$London ## just analyze the biennial portion of the data London <- subset(London, time > 1950) ## define the interval to be 2 weeks IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log') ## look at beta and alpha estimate plotbeta(parms) ## simulate the fitted parameters sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2) ## now lets predict forward 200 years using the mean birth rate, ## starting from rough initial conditions times <- seq(1965,2165, by = 1/ (52/IP)) births <- rep(mean(London$births),length(times)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## take the last 10 years pred <- lapply(pred, function(x) tail(x, 52/IP * 20) ) ## now compute the Lyapunov Exponent for the simulate and predicted model simLE <- TSIR_LE( time=sim$res$time, S=sim$simS$mean, I=sim$res$mean, alpha=sim$alpha, beta=sim$contact$beta, IP=IP ) predLE <- TSIR_LE( time=pred$I$time, S=pred$S$X3, I=pred$I$X3, alpha=parms$alpha, beta=parms$contact$beta, IP=IP ) simLE$LE predLE$LE ## End(Not run)
## Not run: require(kernlab) require(ggplot2) require(kernlab) London <- twentymeas$London ## just analyze the biennial portion of the data London <- subset(London, time > 1950) ## define the interval to be 2 weeks IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log') ## look at beta and alpha estimate plotbeta(parms) ## simulate the fitted parameters sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2) ## now lets predict forward 200 years using the mean birth rate, ## starting from rough initial conditions times <- seq(1965,2165, by = 1/ (52/IP)) births <- rep(mean(London$births),length(times)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## take the last 10 years pred <- lapply(pred, function(x) tail(x, 52/IP * 20) ) ## now compute the Lyapunov Exponent for the simulate and predicted model simLE <- TSIR_LE( time=sim$res$time, S=sim$simS$mean, I=sim$res$mean, alpha=sim$alpha, beta=sim$contact$beta, IP=IP ) predLE <- TSIR_LE( time=pred$I$time, S=pred$S$X3, I=pred$I$X3, alpha=parms$alpha, beta=parms$contact$beta, IP=IP ) simLE$LE predLE$LE ## End(Not run)
A function to calculate the Local Lyapunov Exponennt (LLE) from the TSIR model
TSIR_LLE(LE, m = 1)
TSIR_LLE(LE, m = 1)
LE |
The output of TSIR_LE to pass the Jacobian elements |
m |
The window to sweep the time-varying Jacobian elements. Defaults to one. |
## Not run: require(kernlab) require(ggplot2) require(kernlab) London <- twentymeas$London ## just analyze the biennial portion of the data London <- subset(London, time > 1950) ## define the interval to be 2 weeks IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log') ## look at beta and alpha estimate plotbeta(parms) ## simulate the fitted parameters sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2) ## now lets predict forward 200 years using the mean birth rate, ## starting from rough initial conditions times <- seq(1965,2165, by = 1/ (52/IP)) births <- rep(mean(London$births),length(times)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## take the last 10 years pred <- lapply(pred, function(x) tail(x, 52/IP * 20) ) ## now compute the Lyapunov Exponent for the simulate and predicted model simLE <- TSIR_LE( time=sim$res$time, S=sim$simS$mean, I=sim$res$mean, alpha=sim$alpha, beta=sim$contact$beta, IP=IP ) predLE <- TSIR_LE( time=pred$I$time, S=pred$S$X3, I=pred$I$X3, alpha=parms$alpha, beta=parms$contact$beta, IP=IP ) simLE$LE predLE$LE simLLE <- TSIR_LLE(simLE) predLLE <- TSIR_LLE(predLE) plotLLE(simLLE) plotLLE(predLLE) ## End(Not run)
## Not run: require(kernlab) require(ggplot2) require(kernlab) London <- twentymeas$London ## just analyze the biennial portion of the data London <- subset(London, time > 1950) ## define the interval to be 2 weeks IP <- 2 ## first estimate paramters from the London data parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log') ## look at beta and alpha estimate plotbeta(parms) ## simulate the fitted parameters sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2) ## now lets predict forward 200 years using the mean birth rate, ## starting from rough initial conditions times <- seq(1965,2165, by = 1/ (52/IP)) births <- rep(mean(London$births),length(times)) S0 <- parms$sbar I0 <- 1e-5*mean(London$pop) pred <- predicttsir(times=times,births=births, beta=parms$contact$beta,alpha=parms$alpha, S0=S0,I0=I0, nsim=50,stochastic=T) ## take the last 10 years pred <- lapply(pred, function(x) tail(x, 52/IP * 20) ) ## now compute the Lyapunov Exponent for the simulate and predicted model simLE <- TSIR_LE( time=sim$res$time, S=sim$simS$mean, I=sim$res$mean, alpha=sim$alpha, beta=sim$contact$beta, IP=IP ) predLE <- TSIR_LE( time=pred$I$time, S=pred$S$X3, I=pred$I$X3, alpha=parms$alpha, beta=parms$contact$beta, IP=IP ) simLE$LE predLE$LE simLLE <- TSIR_LLE(simLE) predLLE <- TSIR_LLE(predLE) plotLLE(simLLE) plotLLE(predLLE) ## End(Not run)
A function to take in time cases births and pop vectors (of any lengths) and interpolate them using the given infectious period.
tsiRdata(time, cases, births, pop, IP = 2)
tsiRdata(time, cases, births, pop, IP = 2)
time |
The time vector. |
cases |
The cases vector. |
births |
The births vector. |
pop |
The population vector. |
IP |
The infectious period (in weeks) to discretize to. Defaults to 2. |
twentymeas is a list containing 20 dataframes with cases, births, populations. Each dataframe is a 22 year time series at biweekly (i.e. IP=2) intervals.
data("twentymeas")
data("twentymeas")
From Bryan Grenfell
names(twentymeas) london <- twentymeas[["London"]] plotdata(london)
names(twentymeas) london <- twentymeas[["London"]] plotdata(london)
the function just breaks up the plot area into a grid. Called internally.
vplayout(x, y)
vplayout(x, y)
x |
is the x location of the plot |
y |
is the y lcoation of the ploy |