R/do_1country_run.R

In FPRgroup/fpemmodeling: fpemmodeling

#' Run one country model
#'
#' Runs the family planning estimation model and returns samples and bias adjusted observations for respective run.
#'
#' @param is_in_union \emph{\sQuote{Character}} "Y" if women are in union.
#' @param surveydata_filepath \emph{\sQuote{Character}} If NULL package data is used
#' @param service_stats \emph{\sQuote{Logical}} If FALSE service statistics are not used. Service stats are not public. Requires FPcounts/fpemservicestat package
#' @param service_stats_filepath \emph{\sQuote{Character}}If NULL private package data is used. Requires FPcounts/fpemservicestat package
#' @param division_numeric_code \emph{\sQuote{Numeric}} A number associated with the country. See the data from \code{\link[fpemdata:divisions]{fpemdata::divisions}}
#' @param first_year \emph{\sQuote{Numeric}} The first year of model estimates in output. The model will be fit to all data, including dates before this date if available
#' @param last_year \emph{\sQuote{Numeric}} The last year of model estimates in output. The model will be fit to all data, including dates after this date if available
#' @param subnational '\emph{\sQuote{Logical}} If FALSE runs the national model
#'
#'
#' @return \emph{\sQuote{List}}
#' \enumerate{
#'   \item \strong{posterior_samples}  \emph{\sQuote{Numeric array}} An array of samples of dimension chains x samples x years x proportions
#'   \item \strong{core_data}          \emph{\sQuote{Data.frame}} The processed data associated with the model run from \code{\link[fpemmodeling:core_data]{fpemmodeling::core_data}}
#' }
#'
#' @examples See the reposiotry url in references for detailed examples
#' @references \url{https://github.com/FPcounts/FPEM}
#' @export
do_1country_run <- function(
                            is_in_union = "Y",
                            surveydata_filepath = NULL,
                            service_stats = FALSE,
                            service_stats_filepath = NULL,
                            division_numeric_code,
                            first_year = NULL,
                            last_year,
                            subnational = FALSE
) {
  # check inputs to this wrapper
  check_inputs(
    surveydata_filepath = surveydata_filepath,
    subnational = subnational,
    division_numeric_code = division_numeric_code
  )
  # core data consists of imported data which is filtered and settings for the run
  core_data <- core_data(
    is_in_union = is_in_union,
    surveydata_filepath = surveydata_filepath,
    division_numeric_code = division_numeric_code,
    first_year = first_year,
    last_year = last_year,
    subnational = subnational
  )
  # processing for jags model including index
  list_auxiliary <- list_auxiliary_data(core_data)
  list_global <- list_global_data(is_in_union = is_in_union,
                                  core_data = core_data)
  list_bias <- list_bias_data(core_data$observations)
  list_service_stats <-
    list_service_stats(
      service_stats = service_stats,
      service_stats_filepath = service_stats_filepath,
      seq_years = core_data$year_sequence_list$model_seq_years,
      division_numeric_code,
      first_year
    )
  # write and run model
  write_jags_model(
    old_dm = TRUE,
    #is_in_union = is_in_union,
    include_ss_data = !is.null(list_service_stats),
    nulldata = nrow(core_data$observations) == 0,
    is_in_union = is_in_union
  )
  mod <- jags.parallel(
    data = c(list_auxiliary, list_global, list_bias, list_service_stats, Y = 1),
    parameters.to.save = c("mod.ct", "unmet.ct", "trad.ct", "mu.jn", "logitratio.yunmet.hat.j"),
    model.file = "model.txt",
    n.chains = 3,
    n.iter = 2000,
    n.burnin = 500
    )
  # bias adjusted data added to core data based on model estimates of bias
  if (nrow(core_data$observations) > 0) {
    core_data$observations <- bias_adj(
      core_data = core_data,
      list_auxiliary = list_auxiliary,
      list_global = list_global,
      mod = mod)
  }
  # reformat samples
  posterior_samples <- posterior_samples_array_format(mod, core_data)
  return(list(
    posterior_samples = posterior_samples,
    core_data = core_data
    )
  )
}