algo.twins {surveillance} | R Documentation |
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, logFile="twins.log"))
disProgObj |
object of class disProg |
control |
control object:
|
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.
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. |
M. Hofmann and M. Höhle
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.
## 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)