#########################################################################
##       R PROGRAM: Priors_sensitivity_analysis.R
##
##         PROJECT: Evaluation of multiple-encounter population size
##                  estimators
##
##     DESCRIPTION: A brief sensitivity analysis for the BMA and LCM
##                  priors.
##
##      WRITTEN BY: Steve Gutreuter, CDC/CGH/DGHT Statistics, Estimation
##                                   and Modeling Team
##                  sgutreuter@cdc.gov
##
##            DATE: 2021-10-16
##
##      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(HDInterval)
library(tidyr)
library(dga)
library(LCMCR)

basepath <- file.path(Sys.getenv("PROJ"), "PSE/PSEsim")
workpath <- file.path(basepath, "R")
datapath <- file.path(basepath, "data")
setwd(workpath)

source(file.path(workpath, "PSE_sim_functions.R"))
dataname <- file.path(datapath, "PSEsimSamples.Rdata")
attach(dataname)

cgwtools::lsdata(file.path(datapath, "PSEsimSamples.Rdata"))

#########################################################################
## BMA loglinear model fitting to Mbht
#########################################################################
data(graphs3)
data(graphs4)
data(graphs5)
## priorNmax = 5 x N
bma.Mbht.10K.5xprior <- bma.estimate(data = Mbht.10K.100, events = c(3, 5),
                               priorNmax = 50000, nreps = 3,
                               progress = FALSE)
## priorNmax = 10 x N
bma.Mbht.10K.10xprior <- bma.estimate(data = Mbht.10K.100, events = c(3, 5),
                               priorNmax = 100000, nreps = 3,
                               progress = FALSE)
## priorNmax = 50 x N
bma.Mbht.10K.50xprior <- bma.estimate(data = Mbht.10K.100, events = c(3, 5),
                               priorNmax = 500000, nreps = 3,
                               progress = FALSE)
#########################################################################
## LCM model fitting to Mbht varying K
#########################################################################
seed_  <- 321
lcm.Mbht.10K.K1prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 1, seed = seed_)
lcm.Mbht.10K.K2prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 2, seed = seed_)
lcm.Mbht.10K.K3prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 3, seed = seed_)
lcm.Mbht.10K.K4prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 4, seed = seed_)
lcm.Mbht.10K.K5prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 5, seed = seed_)
lcm.Mbht.10K.K10prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 10, seed = seed_)
lcm.Mbht.10K.K15prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 15, seed = seed_)
lcm.Mbht.10K.K20prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 20, seed = seed_)
lcm.Mbht.10K.K30prior <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 1, a_alpha = 0.25,
                               b_alpha = 0.25, K = 30, seed = seed_)
#########################################################################
## LCM model fitting to Mbht varying Gamma hyperparameters alpha and beta
#########################################################################
lcm.Mbht.10K.hyp.025 <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 3, a_alpha = 0.25,
                               b_alpha = 0.25, K = 10, seed = seed_)
lcm.Mbht.10K.hyp.25 <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 3, a_alpha = 0.25,
                               b_alpha = 0.25, K = 10, seed = seed_)
lcm.Mbht.10K.hyp1 <- lcm.estimate(data = Mbht.10K.100, events = c(3,5),
                               buffer_size = 1000000, thinning = 100,
                               burnin = 500000, samples = 50000,
                               nreps = 3, a_alpha = 1,
                               b_alpha = 1, K = 10, seed = seed_)
#########################################################################
## Save the results
#########################################################################
save(bma.Mbht.10K.5xprior, bma.Mbht.10K.10xprior, bma.Mbht.10K.50xprior,
     bma.Mbht.10K.100xprior, lcm.Mbht.10K.K1prior, lcm.Mbht.10K.K2prior,
     lcm.Mbht.10K.K3prior, lcm.Mbht.10K.K4prior, lcm.Mbht.10K.K5prior,
     lcm.Mbht.10K.K10prior, lcm.Mbht.10K.K15prior, lcm.Mbht.10K.K20prior,
     lcm.Mbht.10K.K30prior, lcm.Mbht.10K.hyp.025, lcm.Mbht.10K.hyp.25,
     lcm.Mbht.10K.hyp.1,
     file = file.path(datapath, "Priors_sensitivity_analysis.Rdata"))
################################   END of FILE   ########################
