algo.twins.Rd
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"))
object of class disProg
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.
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.
Returns an object of class atwins
with elements
specified control object
specified disProg
-object
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.
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.
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.
M. Hofmann and M. Höhle and D. Sabanés Bové
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)
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")])))
}
}