surveillance (version 1.20.3)

algo.twins: Fit a Two-Component Epidemic Model using MCMC

Description

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

This is an experimental implementation that may be removed in future versions of the package.

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

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.

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.

Author

M. Hofmann and M. Höhle and D. Sabanés Bové

References

Held, L., Hofmann, M., Höhle, 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 (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)
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")])))
}