algo.twins {surveillance}R Documentation

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, alpha_xi=10, beta_xi=10, psiRWSigma=0.25,
   alpha_psi=1, beta_psi=0.1, 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
alpha_xi
Parameter α_xi of the hyperprior of the epidemic parameter λ
beta_xi
Parameter β_xi of the hyperprior of the epidemic parameter λ
psiRWSigma
Starting value for the tuning of the variance of the random walk proposal for the overdispersion parameter psi.
alpha_psi
Parameter α_psi of the prior of the overdispersion parameter psi
beta_psi
Parameter β_psi of the prior of the overdispersion parameter psi
logFile
Base file name for the output files. The function writes three output files in your current working directory (i.e. getwd()). If logfile = "twins.log" the results are stored in the three files "twins.log", "twins.log2" and "twins.log.acc". "twins.log" containes the returned samples of the parameters psi, gamma_0, gamma_1, gamma_2, K, xi_λ λ_1,...,λ_{n}, the predictive distribution of the number of cases at time n+1 and the deviance. "twins.log2" 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. "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.

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.

Currently, the example section is commented, because the underlying C++ code appears to crash under some UNIX distributions and it takes rather long time.

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_λ λ_1,...,λ_{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.

Author(s)

M. Hofmann and M. Höhle

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

## Not run: 
#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. Apparently, the underlying C++
#code can cause problems on some LINUX distributions (e.g. ubuntu):
# *** glibc detected *** /usr/lib/R/bin/exec/R: corrupted double-linked
# so until further notice this example code is put into dontrun mode
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")])))
## End(Not run)

[Package surveillance version 1.0-3 Index]