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

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.

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

Finally, inspection of the C++ code using valgrind shows some memory leaks when running the old underlying C++ program. As we are unable to fix this impurity at the present time, we have instead put the example code in a 'dontrun' environment. The example code, however, works fine -- the measure is thus more aimed at reducing the number of CRAN problems with the package.

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

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

## Author

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

## Examples

if (FALSE) {
# 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)
otwins <- algo.twins(hepatitisA,
control=list(burnin=500, filter=1, sampleSize=1000))

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

hist(otwins$logFile$psi,xlab=expression(psi),main="")
if (require("coda")) {
print(summary(mcmc(otwins\$logFile[,c("psi","xipsi","K")])))
}
}