surveil-mcmc.R
In surveil: Time Series Models for Disease Surveillance

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  results = "hold", 
  collapse = TRUE, 
  eval = TRUE,
  fig.pos = 'h',
  fig.width = 5,
  fig.height = 3.5,
  fig.align = 'center'
)

## ----message = FALSE----------------------------------------------------------
library(surveil)

## -----------------------------------------------------------------------------
data(cancer)
cancer2 <- subset(cancer, Age == "10-14")
fit <- stan_rw(cancer2, time = Year, iter = 1500, chains = 4, refresh = 0)

## -----------------------------------------------------------------------------
samples <- fit$samples
class(samples)

## -----------------------------------------------------------------------------
phi <- as.matrix(samples, pars = "rate")
dim(phi)

## -----------------------------------------------------------------------------
phi <- 10e3 * phi
phi_1999 <- phi[,1] 
head(phi_1999)

## -----------------------------------------------------------------------------
hist(phi_1999, xlab = "Rate per 10,000, 1999", main = NA)

## -----------------------------------------------------------------------------
mean(phi_1999)


## -----------------------------------------------------------------------------
quantile(phi_1999, probs = c(0.025, 0.975))

## -----------------------------------------------------------------------------
plot(fit, style = "lines", scale = 10e3)

## -----------------------------------------------------------------------------
phi_2017 <- phi[,19]
diff <- (phi_2017 - phi_1999)
hist(diff, xlab = "Rate difference", main = NA)

## -----------------------------------------------------------------------------
mean(diff)

## -----------------------------------------------------------------------------
100 * mean( diff / phi_1999 )

## -----------------------------------------------------------------------------
print(samples)

## ----fig.width = 4, fig.height = 3.5, eval = FALSE----------------------------
#  # plot not shown
#  rstan::stan_rhat(samples)