algo.twins | R Documentation |
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.
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"))
disProgObj |
object of class |
control |
control object:
|
Returns an object of class atwins
with elements
control |
specified control object |
disProgObj |
specified |
logFile |
contains the returned samples of the parameters |
logFile2 |
contains the sample means of the variables |
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.
M. Hofmann and M. Höhle and D. Sabanés Bové
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.
# 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")])))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.