#########################################################################
##       R PROGRAM: Prior_Sensitivity_Analysis_Plots.R
##
##         PROJECT: Evaluation of multiple-encounter population size
##                  estimators
##
##     DESCRIPTION: Plot estimates from the sensitivity-analysis
##                  estimates
##
##      WRITTEN BY: Steve Gutreuter, CDC/CGH/DGHT Statistics, Estimation
##                                   and Modeling Team
##                  sgutreuter@cdc.gov
##
##            DATE: 2021-10-20
##
##      DISCLAIMER: Although this program has been used by the Centers
##                  for Disease Control & Prevention (CDC), no warranty,
##                  expressed or implied, is made by the CDC or the U.S.
##                  Government as to the accuracy and functioning of the
##                  program and related program material nor shall the
##                  fact of distribution constitute any such warranty,
##                  and no responsibility is assumed by the CDC in
##                  connection therewith.
##
#########################################################################
library(tidyverse)
library(ggplot2)
basepath <- file.path(Sys.getenv("PROJ"), "PSE/PSEsim")
workpath <- file.path(basepath, "R")
datapath <- file.path(basepath, "data")
outpath <- file.path(basepath, "output")
setwd(workpath)
#########################################################################
## Get the estimates
#########################################################################
ests <- readRDS(file.path(datapath,
                          "Priors_sensitivity_analysis_estimates.rds"))
#########################################################################
## Plot effect of Nmax on BMA estimates
#########################################################################
e1 <- ests %>%
    filter(Model == "BMA") %>%
    select(Events, Rep, Nest, HPDCIwidth, bmaNmax) %>%
    mutate(Lists = ordered(Events),
           bmaNmax = bmaNmax * 10)
p1 <- ggplot(e1, aes(x = bmaNmax, y = Nest, group = Rep, colour = Lists)) +
    theme_classic() +
    facet_grid(rows = vars(Lists)) +
    geom_line() + geom_point() +
    scale_color_manual(values = c("red", "blue")) +
    labs(x = "Nmax (thousands)",
         y = expression(hat(N))) +
    xlim(0, 500)
ggsave(file.path(outpath, "graphs/BMApriors_N_vs_Nmax.png"),
       device = "png")
#########################################################################
## Plot effect of Nmax on BMA HPD interval width
#########################################################################
p2 <- ggplot(e1, aes(x = bmaNmax, y = HPDCIwidth, group = Rep,
                     colour = Lists)) +
    theme_classic() +
    facet_grid(rows = vars(Events)) +
    geom_line() + geom_point() +
    scale_color_manual(values = c("red", "blue")) +
    labs(x = "Nmax (thousands)",
         y = expression(hat(N))) +
    xlim(0, 500)
ggsave(file.path(outpath, "graphs/BMApriors_HPDCIwidth_vs_Nmax.png"),
       device = "png")
#########################################################################
## Plot effect max number of latent classes on LCM estimates
#########################################################################
e2 <- ests %>%
    filter(Model == "LCMCR" & !is.na(lcmK)) %>%
    select(Events, Rep, Nest, HPDCIwidth, lcmK) %>%
    mutate(Lists = ordered(Events))
p2 <- ggplot(e2, aes(x = lcmK, y = Nest, group = Lists, colour = Lists)) +
    theme_classic() +
    geom_line() + geom_point() +
    scale_color_manual(values = c("red", "blue")) +
    labs(x = "Prior maximum number of latent classes",
         y = expression(hat(N)))
ggsave(file.path(outpath, "graphs/LCMpriors_N_vs_Kprior.png"),
       device = "png")
#########################################################################
## Plot effect alpha hyperprior parameters on LCM estimates
#########################################################################
e3 <- ests %>%
    filter(Model == "LCMCR" & !is.na(lcmAlpha)) %>%
    select(Events, Rep, Nest, HPDCIwidth, lcmAlpha) %>%
    mutate(Lists = ordered(Events),
           lcmAlpha = ordered(lcmAlpha))
p3 <- ggplot(e3, aes(x = lcmAlpha, y = Nest, group = Rep, colour = Lists)) +
    theme_classic() +
    facet_grid(rows = vars(Lists)) +
    geom_line() + geom_point() +
    scale_color_manual(values = c("red", "blue")) +
    labs(x = "Gamma hyperprior shape and rate parameters",
         y = expression(hat(N)))
ggsave(file.path(outpath, "graphs/LCMpriors_N_vs_Gamma_parms.png"),
       device = "png")
################################   END of FILE   ########################
