Package 'tsiR' reference manual

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 <https://doi:10.1371/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.2
Built:	2024-11-09 04:22:31 UTC
Source:	https://github.com/adbecker/tsiR

corr

Description

Plot the correlation of the true data against the fitted resimulated data.

Usage

corr(sim)
corr(sim)

Arguments

sim

The dataframe or list produced by the 'runtsir' function.

derivative

Description

This function computes an 8 point derivative.

Usage

derivative(X, Y)
derivative(X, Y)

Arguments

`X`	The variable to differentiate with respect to.
`Y`	The function / vector to differentiate.

epitimes

Description

The times at which we declare a new outbreak has started based on the threshold parameter.

Usage

epitimes(data, threshold, epi.length = 3)
epitimes(data, threshold, epi.length = 3)

Arguments

`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.

estpars

Description

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.

Usage

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)

Arguments

`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.

Examples

## 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)

jagsfilter

Description

Used internally to filter jags results to give just the inference well use.

Usage

jagsfilter(mcmcresults)
jagsfilter(mcmcresults)

Arguments

mcmcresults

is the input from the jags model.

logcorr

Description

Plot the correlation of the true data against the fitted resimulated data.

Usage

logcorr(sim)
logcorr(sim)

Arguments

sim

The dataframe or list produced by the 'runtsir' function.

maxthreshold

Description

A function used to optimize the threshold parameter to give the best fit to the data. Optimizes the fit based on R squared.

Usage

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)

Arguments

`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.

Examples

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)

mcmcestpars

Description

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.

Usage

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)

Arguments

`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.

Usage

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)

Arguments

`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.

plotbeta

Description

Plots the inferred beta with confidence intervals (when they can be calculated)

Usage

plotbeta(dat)
plotbeta(dat)

Arguments

dat

the list produced from the runtsir, mcmctsir, and simulatetsir function.

plotbreaks

Description

Plots the cases data with a line whenever the forward simulation is seeded using the real data.

Usage

plotbreaks(data, threshold)
plotbreaks(data, threshold)

Arguments

`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.

plotcases

Description

Plots just the cases data.

Usage

plotcases(data)
plotcases(data)

Arguments

data

The data frame with cases.

plotcomp

Description

Plots just the comparison of the forward simulation fit to the data.

Usage

plotcomp(sim, errtype = "95", max.plot = 10)
plotcomp(sim, errtype = "95", max.plot = 10)

Arguments

`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.

plotdata

Description

Plots the cases data as well as birth and population dynamics.

Usage

plotdata(data)
plotdata(data)

Arguments

data

The dataframe with time, cases, births, and pop.

plotforward

Description

Plots the forward simulation from the TSIR model

Usage

plotforward(dat, inverse = F)
plotforward(dat, inverse = F)

Arguments

`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.

plotLLE

Description

Function to plot the Local Lyapunov Exponents. The output is of class ggplot2 so you can add standard ggplot2 options to it if desired.

Usage

plotLLE(LLE)
plotLLE(LLE)

Arguments

LLE

The output from TSIR_LLE

Examples



## 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)

plotregression

Description

Plots the cumulative cases - cumulative births data and regression fit

Usage

plotregression(dat)
plotregression(dat)

Arguments

dat

the list produced from the runtsir, mcmctsir, and simulatetsir function.

plotres

Description

Plots diagnostics and results of the runtsir function.

Usage

plotres(dat, max.plot = 10)
plotres(dat, max.plot = 10)

Arguments

`dat`	the list produced from the runtsir, mcmctsir, and simulatetsir function.
`max.plot`	the number of individual stochastic simulations to plot. Defaults to 10.

plotrho

Description

Plots the inferred reporting rate, rho

Usage

plotrho(dat)
plotrho(dat)

Arguments

dat

the list produced from the runtsir, mcmctsir, and simulatetsir function.

plotsbar

Description

Plots the profile log likelihood calculation for inferred sbar

Usage

plotsbar(dat)
plotsbar(dat)

Arguments

dat

the list produced from the runtsir, mcmctsir, and simulatetsir function.

predicttsir

Description

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.

Usage

predicttsir(times, births, beta, alpha, S0, I0, nsim, stochastic)
predicttsir(times, births, beta, alpha, S0, I0, nsim, stochastic)

Arguments

`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.

Examples

## 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)

residuals.births

Description

computes the residuals for when X is the cumulative births. Used internally.

Usage

residual.births(rho, Yhat, Y)
residual.births(rho, Yhat, Y)

Arguments

`rho`	The reporting rate, used to get units correct.
`Yhat`	The fitted regression line.
`Y`	The cumulative cases.

residuals.cases

Description

Computes the residuals for when X is the cumulative cases. Used internally.

Usage

residual.cases(Yhat, Y)
residual.cases(Yhat, Y)

Arguments

`Yhat`	The fitted regression line.
`Y`	The cumulative births.

runtsir

Description

This function runs the TSIR model.

Usage

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)

Arguments

`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.

Examples

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)

simulatetsir

Description

This function just simulates the forward prediction given the data and a parms list generated from estpars or mcmcestpars.

Usage

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)

Arguments

`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.

TSIR_LE

Description

A function to calculate the Lyapunov Exponennt (LE) from the TSIR model

Usage

TSIR_LE(time, S, I, alpha, beta, IP)
TSIR_LE(time, S, I, alpha, beta, IP)

Arguments

`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)

Examples


## 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)

TSIR_LLE

Description

A function to calculate the Local Lyapunov Exponennt (LLE) from the TSIR model

Usage

TSIR_LLE(LE, m = 1)
TSIR_LLE(LE, m = 1)

Arguments

`LE`	The output of TSIR_LE to pass the Jacobian elements
`m`	The window to sweep the time-varying Jacobian elements. Defaults to one.

Examples



## 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)

tsiRdata

Description

A function to take in time cases births and pop vectors (of any lengths) and interpolate them using the given infectious period.

Usage

tsiRdata(time, cases, births, pop, IP = 2)
tsiRdata(time, cases, births, pop, IP = 2)

Arguments

`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.

Measles incidence data from 20 cities

Description

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.

Usage

data("twentymeas")data("twentymeas")

Source

From Bryan Grenfell

Examples

names(twentymeas)
london <- twentymeas[["London"]]
plotdata(london)names(twentymeas)
london <- twentymeas[["London"]]
plotdata(london)

vplayout

Description

the function just breaks up the plot area into a grid. Called internally.

Usage

vplayout(x, y)
vplayout(x, y)

Arguments

`x`	is the x location of the plot
`y`	is the y lcoation of the ploy

Package 'tsiR'

Help Index

corr

Description

Usage

Arguments

derivative

Description

Usage

Arguments

epitimes

Description

Usage

Arguments

estpars

Description

Usage

Arguments

Examples

jagsfilter

Description

Usage

Arguments

logcorr

Description

Usage

Arguments

maxthreshold

Description

Usage

Arguments

Examples

mcmcestpars

Description

Usage

Arguments

mcmctsir

Description

Usage

Arguments

plotbeta

Description

Usage

Arguments

plotbreaks

Description

Usage

Arguments

plotcases

Description

Usage

Arguments

plotcomp

Description

Usage

Arguments

plotdata

Description

Usage

Arguments

plotforward

Description

Usage

Arguments

plotLLE

Description

Usage

Arguments

Examples

plotregression

Description

Usage

Arguments

plotres

Description

Usage

Arguments

plotrho

Description

Usage