algo.twins: Model fit based on a two-component epidemic model

In jimhester/surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

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, 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 α_ξ of the hyperprior of the epidemic parameter λ beta_xi Parameter β_ξ 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 ψ. alpha_psi Parameter α_ψ of the prior of the overdispersion parameter ψ beta_psi Parameter β_ψ of the prior of the overdispersion parameter ψ 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 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 ψ, γ_0, γ_1, γ_2, K, ξ_λ λ_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, ω_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 ψ, the changepoints and the endemic parameters γ_0, γ_1, γ_2 in the third column and the variance of the random walk proposal for the update of the parameter ψ 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.

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	containes the returned samples of the parameters ψ, γ_0, γ_1, γ_2, K, ξ_λ λ_1,...,λ_{n}, the predictive distribution and the deviance.
logFile2	containes the sample means of the variables X_t, Y_t, ω_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 and D. Sabanés Bové

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 (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")])))
}

## End(Not run)

jimhester/surveillance documentation built on May 19, 2019, 10:33 a.m.