Skip to content

Commit

Permalink
Merge pull request #61 from annecori/new-version
Browse files Browse the repository at this point in the history
New version
  • Loading branch information
jstockwin authored Jun 4, 2019
2 parents 7f8363b + 7d40276 commit 1537bff
Show file tree
Hide file tree
Showing 49 changed files with 5,376 additions and 2,640 deletions.
3 changes: 2 additions & 1 deletion .Rbuildignore
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
fitRotavirus.R
^.*\.Rproj$
^\.Rproj\.user$
^\.travis.yml
^\.travis\.yml$
^codecov\.yml$
35 changes: 34 additions & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
@@ -1,2 +1,35 @@
language: r
cache: packages
cache: packages

matrix:
include:
- os: linux
r: release
env:
- R_CODECOV=true
- os: linux
r: devel
- os: linux
r: oldrel
- os: osx
osx_image: xcode8.3

warnings_are_errors: true

notifications:
email:
on_success: change
on_failure: change

after_success:
- if [[ "${R_CODECOV}" ]]; then Rscript -e 'covr::codecov()'; fi

env:
global:
- NOT_CRAN: true
- R_BUILD_ARGS="--resave-data --compact-vignettes=gs+qpdf"
- R_CHECK_ARGS="--as-cran --timings"
- R_CHECK_TIME="TRUE"
- R_CHECK_TESTS="TRUE"
- _R_CHECK_TIMINGS_="0"
- _R_CHECK_SYSTEM_CLOCK_="0"
9 changes: 4 additions & 5 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
Package: EpiEstim
Version: 2.0-0
Date: 2012/01/17
Version: 2.1-0
Title: EpiEstim: a package to estimate time varying reproduction
numbers from epidemic curves
Author: Anne Cori <[email protected]>
Expand All @@ -24,12 +23,12 @@ Imports:
scales,
grDevices
Remotes:
reconhub/incidence,
nickreich/coarseDataTools@hackout3
nickreich/coarseDataTools
Suggests:
testthat,
compare,
utils
utils,
covr
License: GPL (>=2)
LazyLoad: yes
Packaged: 2012-05-30 10:35:03 UTC; acori
Expand Down
6 changes: 4 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,16 +1,18 @@
# Generated by roxygen2: do not edit by hand

S3method(plot,estimate_R)
export(DiscrSI)
export(EstimateR)
export(OverallInfectivity)
export(WT)
export(check_cdt_samples_convergence)
export(coarse2estim)
export(discr_si)
export(estimate_r)
export(estimate_R)
export(init_mcmc_params)
export(make_config)
export(make_mcmc_control)
export(overall_infectivity)
export(plots)
export(wallinga_teunis)
import(grid)
import(gridExtra)
Expand Down
86 changes: 45 additions & 41 deletions R/check_cdt_samples_convergence.R
Original file line number Diff line number Diff line change
@@ -1,66 +1,70 @@
#######################################################################################################################
# check_cdt_samples_convergence runs a Gelman Rubin test to check convergence of the MCMC chain in coarseDataTools #
#######################################################################################################################
################################################################################
# check_cdt_samples_convergence runs a Gelman Rubin test to check convergence of
# the MCMC chain in coarseDataTools #
################################################################################

#' Checking convergence of an MCMC chain by using the Gelman-Rubin algorithm
#'
#' \code{check_cdt_samples_convergence} Checking convergence of an MCMC chain by using the Gelman-Rubin algorithm
#'
#' @param cdt_samples the \code{@sample} slot of a \code{cd.fit.mcmc} S4 object (see package \code{coarseDataTools})
#' @return TRUE if the Gelman Rubin test for convergence was successful, FALSE otherwise
#'
#' \code{check_cdt_samples_convergence} Checking convergence of an MCMC chain by
#' using the Gelman-Rubin algorithm
#'
#' @param cdt_samples the \code{@sample} slot of a \code{cd.fit.mcmc} S4 object
#' (see package \code{coarseDataTools})
#' @return TRUE if the Gelman Rubin test for convergence was successful, FALSE
#' otherwise
#' @details{
#' This function splits an MCMC chain in two halfs and uses the Gelman-Rubin algorithm to assess convergence of the chain by comparing its two halves.
#' This function splits an MCMC chain in two halfs and uses the Gelman-Rubin
#' algorithm to assess convergence of the chain by comparing its two halves.
#' }
#' @seealso \code{\link{estimate_r}}
#' @seealso \code{\link{estimate_R}}
#' @author Anne Cori
#' @importFrom coda gelman.diag
#' @export
#' @examples
#' \dontrun{
#' ## Note the following examples use an MCMC routine
#' ## to estimate the serial interval distribution from data,
#' ## Note the following examples use an MCMC routine
#' ## to estimate the serial interval distribution from data,
#' ## so they may take a few minutes to run
#'
#'
#' ## load data on rotavirus
#' data("MockRotavirus")
#'
#' ## estimate the serial interval from data
#' SI_fit <- coarseDataTools::dic.fit.mcmc(dat = MockRotavirus$si_data,
#' dist="G",
#' init_pars=init_mcmc_params(MockRotavirus$si_data, "G"),
#' burnin = 1000,
#' n.samples = 5000)
#'
#'
#' ## estimate the serial interval from data
#' SI_fit <- coarseDataTools::dic.fit.mcmc(dat = MockRotavirus$si_data,
#' dist="G",
#' init_pars=init_mcmc_params(MockRotavirus$si_data, "G"),
#' burnin = 1000,
#' n.samples = 5000)
#'
#' ## use check_cdt_samples_convergence to check convergence
#' converg_diag <- check_cdt_samples_convergence(SI_fit@samples)
#' converg_diag
#'
#'
#' }
#'
check_cdt_samples_convergence <- function(cdt_samples)
{

#############################################################################################################################
# checking convergence of the MCMC by using the Gelman-Rubin algorithm between the first and second half of the MCMC sample
spl1 <- cdt_samples[1:floor(nrow(cdt_samples)/2),]
spl2 <- cdt_samples[(ceiling(nrow(cdt_samples)/2)+1):nrow(cdt_samples),]
#'
check_cdt_samples_convergence <- function(cdt_samples) {
## checking convergence of the MCMC by using the Gelman-Rubin algorithm
## between the first and second half of the MCMC sample
spl1 <- cdt_samples[seq_len(floor(nrow(cdt_samples) / 2)), ]
spl2 <- cdt_samples[seq(ceiling(nrow(cdt_samples) / 2) + 1, nrow(cdt_samples)), ]
GRD <- gelman.diag(as.mcmc.list(list(as.mcmc(spl1), as.mcmc(spl2))))
# Is any of the potential scale reduction factors >1.1 (looking at the upper CI)?
# If so this would suggest that the MCMC has not converged well.
if(any(GRD$psrf[,"Upper C.I."]>1.1))
{
warning("The Gelman-Rubin algorithm suggests the MCMC may not have converged within the number of iterations (MCMC.burnin + n1) specified.
# Is any of the potential scale reduction factors >1.1
# (looking at the upper CI)?
# If so this would suggest that the MCMC has not converged well.
if (any(GRD$psrf[, "Upper C.I."] > 1.1)) {
warning("The Gelman-Rubin algorithm suggests the MCMC may not have converged
within the number of iterations (MCMC.burnin + n1) specified.
You can visualise the full MCMC chain using: \n
> par(mfrow=c(2,1))
> plot(res$SI.Moments[,'Mean'], type='l', xlab='Iterations', ylab='Mean SI')
> plot(res$SI.Moments[,'Std'], type='l', xlab='Iterations', ylab='Std SI'),
where res is the output of estimate_r
> plot(res$SI.Moments[,'Mean'], type='l', xlab='Iterations',
ylab='Mean SI')
> plot(res$SI.Moments[,'Std'], type='l', xlab='Iterations',
ylab='Std SI'),
where res is the output of estimate_R
and decide whether to rerun for longer.")
return(FALSE)
}else
{
} else {
cat("\nGelman-Rubin MCMC convergence diagnostic was successful.")
return(TRUE)
}

}
Loading

0 comments on commit 1537bff

Please sign in to comment.