## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4,
  fig.align = "center"
)

## ----setup--------------------------------------------------------------------
library(MAPCtools)
data("toy_data")

## ----plot-missing-data--------------------------------------------------------
plot_missing_data(toy_data, x = period, y = age)

## ----plot-missing-data-stratify-----------------------------------------------
plot_missing_data(
  data = toy_data,
  x = period,
  y = age,
  stratify_by = education
)

## ----plot-missing-data-for-each-----------------------------------------------
plot_missing_data(
  data = toy_data,
  x = period,
  y = age,
  stratify_by = education,
  for_each = sex
)

## ----plot-coints-1d-----------------------------------------------------------
plot_counts_1D(
  toy_data,
  x = age,
  stratify_by = education,
  for_each = sex
)

## ----plot-counts-2d-----------------------------------------------------------
plot_counts_2D(
  toy_data,
  x = age,
  y = period,
  stratify_by = education,
  for_each = sex
)

## ----plot-binnd-counts--------------------------------------------------------
plot_binned_counts(
  toy_data,
  x = period,
  bin_by = age,
  n_bins = 4,
  stratify_by = education,
  for_each = sex
)

## ----plot-mean-1d-------------------------------------------------------------
plot_mean_response_1D(
  toy_data,
  response = count,
  x = age,
  stratify_by = education,
  for_each = sex
)

## ----plot-mean-2d-------------------------------------------------------------
plot_mean_response_2D(
  toy_data,
  response = count,
  x = period,
  y = age,
  stratify_by = education
)

## ----plot-known-rate----------------------------------------------------------
require(ggplot2)

# Over age
ggplot(toy_data, aes(x = age, y = known_rate, color = education)) +
  stat_summary(fun=mean, geom="line") +
  facet_wrap(~ sex, ncol = 2) +
  labs(
    title = "Poisson rates by age and education level",
    x = "Age",
    y = "Rate",
    color = "Education"
  ) +
  scale_color_viridis_d() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")

# Over period
ggplot(toy_data, aes(x = period, y = known_rate, color = education)) +
  stat_summary(fun=mean, geom="line") +
  facet_wrap(~ sex, ncol = 2) +
  labs(
    title = "Poisson rates by period and education level",
    x = "Period",
    y = "Rate",
    color = "Education"
  ) +
  scale_color_viridis_d() + 
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")

# Over cohort
ggplot(toy_data, aes(x = cohort, y = known_rate, color = education)) +
  stat_summary(fun=mean, geom="line") +
  facet_wrap(~ sex, ncol = 2) +
  labs(
    title = "Poisson rates by cohort and education level",
    x = "Cohort",
    y = "Rate",
    color = "Education"
  ) +
  scale_color_viridis_d() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")

## ----filter-sex---------------------------------------------------------------
require(dplyr)
toy_data.f <- toy_data %>% filter(sex == "female") %>% subset(cohort > 1931)
toy_data.m <- toy_data %>% filter(sex == "male") %>% subset(cohort > 1931 & cohort < 1999)

## ----fit-apC, eval=FALSE, echo=TRUE-------------------------------------------
# apC_fit.f <- fit_MAPC(
#   data = toy_data.f,
#   response = count,
#   family = "poisson",
#   apc_format = "apC",
#   stratify_by = education,
#   reference_strata = 1,
#   age = age,
#   period = period
# )
# 
# apC_fit.m <- fit_MAPC(
#   data = toy_data.m,
#   response = count,
#   family = "poisson",
#   apc_format = "apC",
#   stratify_by = education,
#   reference_strata = 1,
#   age = age,
#   period = period
# )

## ----load-apC-fit, results = 'hide'-------------------------------------------
apC_fit.f <- readRDS(system.file("extdata", "quickstart-apC_fit_f.rds", package = "MAPCtools"))
apC_fit.m <- readRDS(system.file("extdata", "quickstart-apC_fit_m.rds", package = "MAPCtools"))

## ----print-apC-fit------------------------------------------------------------
print(apC_fit.f)      # Concise summary the model that was fit
# print(apC_fit.f)

## ----plot-apC-fit-f-----------------------------------------------------------
plot(apC_fit.f)       # Plots estimated cross-strata contrast trends

## ----plot-apC-fit-m-----------------------------------------------------------
plot(apC_fit.m)       # Plots estimated cross-strata contrast trends

## ----summary-apC-fit----------------------------------------------------------
# This doesn't print nice in a rmd/qmd file
# summary(apC_fit)       # Detailed posterior summaries

## ----fit-all-mapc, eval=FALSE, echo=TRUE--------------------------------------
# all_fits.f <- fit_all_MAPC(
#   data = toy_data.f,
#   response = count,
#   family = "poisson",
#   stratify_by = education,
#   reference_strata = 1,
#   age = age,
#   period = period,
#   include.random = TRUE
# )
# 
# all_fits.m <- fit_all_MAPC(
#   data = toy_data.m,
#   response = count,
#   family = "poisson",
#   stratify_by = education,
#   reference_strata = 1,
#   age = age,
#   period = period,
#   include.random = TRUE
# )

## ----load-all-fits, results = 'hide'------------------------------------------
all_fits.f <- readRDS(system.file("extdata", "quickstart-all_fits_f.rds", package = "MAPCtools"))
all_fits.m <- readRDS(system.file("extdata", "quickstart-all_fits_m.rds", package = "MAPCtools"))

## ----print-all_fits-----------------------------------------------------------
print(all_fits.f)    # concise summary of each model
# print(all_fits.m)

## ----plot-all-fits-f----------------------------------------------------------
plot(all_fits.f)     # model comparison plots (DIC/WAIC/log-score)

## ----plot-all-fits-m----------------------------------------------------------
plot(all_fits.m)     # model comparison plots (DIC/WAIC/log-score)

## ----summary-all-fits---------------------------------------------------------
# summary(all_fits.f)  # detailed posterior summaries for each fit
# summary(all_fits.m)