Package 'twins'

Fit a Two-Component Epidemic Model using MCMC

Description

Fits a negative-binomial model as described in Held et al. (2006) to a univariate time series of counts.

Usage

algo.twins(disProgObj, control=list(burnin=1000, filter=10,
   sampleSize=2500, noOfHarmonics=1, alpha_xi=10, beta_xi=10,
   psiRWSigma=0.25,alpha_psi=1, beta_psi=0.1, nu_trend=FALSE,
   logFile="twins.log"))

Arguments

disProgObj

object of class disProg

control

control object: burnin Number of burn in samples. filter Thinning parameter. If filter = 10 every 10th sample is after the burn in is returned. sampleSize Number of returned samples. Total number of samples = burnin+filter*sampleSize noOfHarmonics Number of harmonics to use in the modelling, i.e. $L$ in (2.2) of Held et al (2006). alpha_xi Parameter $\alpha_{\xi}$ of the hyperprior of the epidemic parameter $\lambda$ beta_xi Parameter $\beta_{\xi}$ of the hyperprior of the epidemic parameter $\lambda$ psiRWSigma Starting value for the tuning of the variance of the random walk proposal for the overdispersion parameter $\psi$. alpha_psi Parameter $\alpha_{\psi}$ of the prior of the overdispersion parameter $\psi$ beta_psi Parameter $\beta_{\psi}$ of the prior of the overdispersion parameter $\psi$ nu_trend Adjust for a linear trend in the endemic part? (default: FALSE) logFile Base file name for the output files. The function writes three output files in the current working directory getwd(). If logfile = "twins.log" the results are stored in the three files ‘twins.log’, ‘twins.log2’ and ‘twins.log.acc’. ‘twins.log’ contains the returned samples of the parameters $\psi$, $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}$, K, $\xi_{\lambda}$ $\lambda_{1},...,\lambda{n}$, the predictive distribution of the number of cases at time $n+1$ and the deviance. ‘twins.log2’ contains the sample means of the variables $X_{t}, Y_{t}, \omega_{t}$ and the relative frequency of a changepoint at time t for t=1,...,n and the relative frequency of a predicted changepoint at time n+1. ‘twins.log.acc’ contains the acceptance rates of $\psi$, the changepoints and the endemic parameters $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}$ in the third column and the variance of the random walk proposal for the update of the parameter $\psi$ in the second column.

Value

Returns an object of class atwins with elements

control	specified control object
disProgObj	specified disProg-object
logFile	contains the returned samples of the parameters $\psi$, $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}$, K, $\xi_{\lambda}$ $\lambda_{1},...,\lambda{n}$, the predictive distribution and the deviance.
logFile2	contains the sample means of the variables $X_{t}, Y_{t}, \omega_{t}$ and the relative frequency of a changepoint at time t for t=1,...,n and the relative frequency of a predicted changepoint at time n+1.

Note

This function is not a surveillance algorithm, but only a modelling approach as described in the Held et. al (2006) paper.

Note also that the function writes three logfiles in the current working directory getwd(): ‘twins.log’, ‘twins.log.acc’ and ‘twins.log2’. Thus you need to have write permissions in the current working directory.

Author(s)

M. Hofmann and M. Hoehle with contributions by M. Paul, D. Sabanes Bove and S. Meyer

References

Held, L., Hofmann, M., Hoehle, M. and Schmid V. (2006) A two-component model for counts of infectious diseases. Biostatistics, 7, 422–437. doi:10.1093/biostatistics/kxj016.

Examples

# Load the data used in the Held et al. (2006) paper
data("hepatitisA")

# Fix seed - this is used for the MCMC samplers in twins
set.seed(123)

# Call algorithm and save result (use short chain without filtering for speed)
oldwd <- setwd(tempdir())  # where logfiles will be written
otwins <- algo.twins(hepatitisA,
                     control=list(burnin=500, filter=1, sampleSize=1000))
setwd(oldwd)

# This shows the entire output (use ask=TRUE for pause between plots)
if (requireNamespace("surveillance"))
  plot(otwins, ask=FALSE)

# Direct access to MCMC output
hist(otwins$logFile$psi,xlab=expression(psi),main="")
if (require("coda"))
  print(summary(mcmc(otwins$logFile[,c("psi","xipsi","K")])))

---

Hepatitis A in Germany

Description

Weekly number of reported hepatitis A infections in Germany 2001-2004.

Usage

data("hepatitisA")

Format

A "disProg" object containing $208\times 1$ observations starting from week 1 in 2001 to week 52 in 2004.

Source

Robert Koch-Institut: SurvStat: https://survstat.rki.de/; Queried on 11-01-2005.

Examples

data("hepatitisA")
str(hepatitisA)
if (requireNamespace("surveillance"))
  plot(hepatitisA)

---

Plots for Fitted algo.twins Models

Description

Plot results of fitting a twins model using MCMC output. Plots similar to those in the Held et al. (2006) paper are generated.

Usage

## S3 method for class 'atwins'
plot(x, which=c(1,4,6,7), ask=TRUE, ...)

Arguments

x	An object of class "atwins" as returned by algo.twins.
which	a vector containing the different plot types to show 1 A plot of the observed time series Z is shown together with posterior means for the number of endemic cases (X) and number of epidemic cases (Y). 2 This plot shows trace plots of the gamma parameters over all MCMC samples. 3 This shows a trace plot of psi, which controls the overdispersion in the model. 4 Autocorrelation functions for K and psi are shown in order to judge whether the MCMC sampler has converged. 5 Shows a plot of the posterior mean of the seasonal model nu[t] together with 95% credibility intervals based on the quantiles of the posterior. 6 Histograms illustrating the posterior density for K and psi. The first one corresponds to Fig. 4(f) in the paper. 7 Histograms illustrating the predictive posterior density for the next observed number of cases Z[n+1]. Compare with Fig.5 in the paper.
ask	Boolean indicating whether to ask for a newline before showing the next plot (only if multiple are shown).
...	Additional arguments for stsplot_time, used for plot type 1.

Details

For details see the plots in the paper. Basically MCMC output is visualized. This function is experimental, as is algo.twins.

Author(s)

M. Hofmann and M. Hoehle

References

Held, L., Hofmann, M., Hoehle, M. and Schmid V. (2006) A two-component model for counts of infectious diseases. Biostatistics, 7, 422–437. doi:10.1093/biostatistics/kxj016.

See Also

algo.twins (with an example)

---