surveillance (version 1.2-1)

algo.twins: Model fit based on a two-component epidemic model

Description

Fits a negative binomial model (as described in Held et al. (2006) to an 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: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

Value

• Returns an object of class atwins with elements
• control
  specified control object
• disProgObj
  specified disProg-object
• logFile
  containes 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
  containes 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.

encoding

latin1

Details

Note that for the time being 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 your current working directory (i.e. getwd()): twins.log, twins.log.acc and twins.log2 Thus you need to have write permissions in the current getwd() directory.

References

Held, L., Hofmann, M., H�hle{Hoehle}, M. and Schmid V. (2006) A two-component model for counts of infectious diseases, Biostatistics, 7, pp. 422--437.

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.
otwins <- algo.twins(hepatitisA)

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

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