################################################################################
##       R PROGRAM: get_prior_sensitivity_estimates.R
##
##         PROJECT: Evaluation of multiple-encounter population size
##                  estimators
##
##     DESCRIPTION: Extract and consolidate PSE estimates into a single
##                  data frame and examine
##
##           INPUT: ~data/Priors_sensitivity_analysis.Rdata
##
##          OUTPUT: ~data/Priors_sensitivity_analysis_estimates.rds
##
##      WRITTEN BY: Steve Gutreuter, CDC/CGH/DGHT Statistics, Estimation
##                                   and Modeling Team
##                  sgutreuter@cdc.gov
################################################################################
library(tidyverse)
library(data.table)
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)
source(file.path(workpath, "PSE_sim_functions.R"))
################################################################################
## Create vector of file names
################################################################################
est.objs <- "Priors_sensitivity_analysis.Rdata"
################################################################################
## Combine all estimates
################################################################################
combineRdata <- function(Rdatanames, datapath){
    combined <- data.frame(NULL)
    counts <- data.frame(NULL)
    for(i in seq_along(Rdatanames)){
        efname <- file.path(datapath, Rdatanames[i])
        attach(efname)
        dfnames <- cgwtools::lsdata(efname)
        for(j in seq_along(dfnames)){
            dfj <- get(dfnames[j])
            dfj$fname <- dfnames[j]
            combined <- rbind(combined, dfj)
            x <- data.frame(with(dfj, table(Events)))
            counts <- rbind(counts,
                            cbind(rep(dfnames[j], nrow(x)), x))
        }
        detach()
    }
    names(counts) <- c("Estimates frame", "Events", "Estimates")
    invisible(list(combined = combined, counts = counts))
}
comb.list <- combineRdata(est.objs, datapath = datapath)
comb <- comb.list$combined
################################################################################
## Create prior parameter variables and save
################################################################################
rcomb <- comb %>%
    mutate(bmaNmax = NA,
           bmaNmax = case_when(
               fname == "bma.Mbht.10K.5xprior" ~ 5,
               fname == "bma.Mbht.10K.10xprior" ~ 10,
               fname == "bma.Mbht.10K.50xprior" ~ 50
           ),
           lcmK = NA,
           lcmK = case_when(
               fname == "lcm.Mbht.10K.K1prior" ~ 1,
               fname == "lcm.Mbht.10K.K2prior" ~ 2,
               fname == "lcm.Mbht.10K.K3prior" ~ 3,
               fname == "lcm.Mbht.10K.K4prior" ~ 4,
               fname == "lcm.Mbht.10K.K5prior" ~ 5,
               fname == "lcm.Mbht.10K.K10prior" ~ 10,
               fname == "lcm.Mbht.10K.K15prior" ~ 15,
               fname == "lcm.Mbht.10K.K20prior" ~ 20,
               fname == "lcm.Mbht.10K.K30prior" ~ 30
           ),
           lcmAlpha = NA,
           lcmAlpha = case_when(
               fname == "lcm.Mbht.10K.hyp.025" ~ 0.025,
               fname == "lcm.Mbht.10K.hyp.25" ~ 0.25,
               fname == "lcm.Mbht.10K.hyp1" ~ 1.0
           )
           )
saveRDS(rcomb,
        file = file.path(datapath, "Priors_sensitivity_analysis_estimates.rds"))
################################   END of FILE   ###############################
