R/fit_fp_c.R

In FPRgroup/FPEMcountry: fpemlocal

#' Fit one country model and save output
#'
#' A wrapper for \code{\link{fit_fp_csub}} augmented to save the output
#'
#' @param runname \emph{\sQuote{Character}} The name given to the file generated by this wrapper.
#' @param ... Please see \code{\link{fit_fp_csub}} for additional arguments.
#'
#' @export
fit_fp_c_autosave <- function(
  runname = NULL,
  ...
) {
  runlist <- fit_fp_c(
    ...
  )
  if (!dir.exists("output")) dir.create("output")
  if (!dir.exists("output/runs")) dir.create("output/runs")
  pathout <- file.path("output/runs", paste0(runname, ".rds"))
  saveRDS(runlist, pathout)
  print(paste0("Your file was saved to ", pathout))
}



#' Fit one country model
#'
#' Fits the family planning estimation model and returns samples and bias adjusted observations for respective run. Please see \code{\link{fit_fp_csub}} for additional argument default settings
#'
#' @param surveydata_filepath \emph{\sQuote{Character}} Path to survey data. Survey data should be a .csv with the following format \code{\link{contraceptive_use}}.
#' @param population_data \emph{\sQuote{Data.frame}} Population count data such as \code{\link{population_counts}}.
#' @param division_numeric_code \emph{\sQuote{Numeric}} A number associated with the country. See the data from \code{\link{divisions}}.
#' @param is_in_union \emph{\sQuote{Character}} "Y" if women are in union.
#' @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 surveydata_filepath \emph{\sQuote{Character}} Path to survey data. Survey data should be a .csv with the following format \code{\link{contraceptive_use}}.
#' @param service_stats_filepath \emph{\sQuote{Character}}
#' @param subnational '\emph{\sQuote{Logical}} If TRUE runs the sub national model.
#' @param diagnostic '\emph{\sQuote{Logical}} If TRUE saves the full jags output.
#' @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{core_data}}.
#' }
#' @examples See the reposiotry url in references for detailed examples
#' @references \url{https://github.com/FPcounts/FPEM}
#' @export
fit_fp_c <- function(
  population_data = NULL,
  is_in_union,
                           ...) {
  if (is_in_union == "ALL") {
    if (is.null(population_data)){
      warning("No population data provided, UNPD population counts used.")
    }
    fits <-
      purrr::pmap(list(is_in_union = list(
        Y = "Y", N = "N"
      ),
      ...),
      fit_fp_csub)
    samples_all <-
      posterior_samples_all_women(
        in_union_posterior_samples = fits %>% purrr::chuck("Y", "posterior_samples"),
        not_in_union_posterior_samples = fits %>% purrr::chuck("N", "posterior_samples"),
        core_data = fits %>% purrr::chuck("Y", "core_data"),
        population_data = population_data
      )
    core_data <- fits %>% purrr::chuck("Y", "core_data")
    core_data$observations <- tibble::tibble()
    core_data$is_in_union <- is_in_union
    ALL <- list(posterior_samples = samples_all,
                   core_data = core_data)
    fitr <- list(
      Y =  fits %>% purrr::chuck("Y"),
      N =  fits %>% purrr::chuck("N"),
      ALL = ALL
    )
  } else {
    fitr <- list(fit_fp_csub(is_in_union = is_in_union, ...))
    names(fitr) <- is_in_union
  }
  return(fitr)

}




#' Fit one country model
#'
#' Fits the family planning estimation model and returns samples and bias adjusted observations for respective run.
#'
#' @param surveydata_filepath \emph{\sQuote{Character}} Path to survey data. Survey data should be a .csv with the following format \code{\link{contraceptive_use}}.
#' @param division_numeric_code \emph{\sQuote{Numeric}} A number associated with the country. See the data from \code{\link{divisions}}.
#' @param is_in_union \emph{\sQuote{Character}} "Y" if women are in union.
#' @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 surveydata_filepath \emph{\sQuote{Character}} Path to survey data. Survey data should be a .csv with the following format \code{\link{contraceptive_use}}.
#' @param service_stats_filepath \emph{\sQuote{Character}}
#' @param subnational '\emph{\sQuote{Logical}} If TRUE runs the sub national model.
#' @param diagnostic '\emph{\sQuote{Logical}} If TRUE saves the full jags output.
#' @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{core_data}}.
#' }
#' @examples See the reposiotry url in references for detailed examples
#' @references \url{https://github.com/FPcounts/FPEM}
#' @export
fit_fp_csub <- function(
  is_in_union = "Y",
  division_numeric_code,
  first_year = NULL,
  last_year,
  surveydata_filepath = NULL,
  service_stats_filepath = NULL,
  subnational = FALSE,
  diagnostic = FALSE,
  params_to_save = c("mod.ct", "unmet.ct", "trad.ct", "mu.in", "logitratio.yunmet.hat.i"),
  nchains = 10,
  niter = 2500,
  nburnin = 500,
  nthin = max(1, floor((niter - nburnin)*nchains / 2000))
) {
  # EMUs are not to be used for unmarried women
  if (is_in_union == "N") service_stats_filepath <- NULL
  # 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, core_data)
  list_bias <- list_bias_data(core_data)
  list_service_stats <-
    list_service_stats(
      service_stats_filepath = service_stats_filepath,
      model_seq_years = core_data$year_sequence_list$model_seq_y,
      division_numeric_code,
      first_year_observed = core_data$year_sequence_list$first_year_observed
    )
  # 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 = params_to_save,
    model.file = "model.txt",
    n.chains = nchains,
    n.iter = niter,
    n.burnin = nburnin,
    n.thin = nthin
  )
  # This is a hack to remove jags data inputs from the global environment
  # This due to poor functionality within R2jags, can be removed if you switch to another package for parallel computing
  temp <- list(list_auxiliary, list_global, list_bias, list_service_stats)
  names_to_remove_from_envr <- c(unlist(lapply(temp, names)), "Y")
  rm(list = names_to_remove_from_envr, envir = globalenv())
  diagnostic_check(diagnostic = diagnostic,
                   mod = mod,
                   core_data = core_data)
  # 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)
  dimnames(posterior_samples) <- list(division_numeric_code)
    return(list(
      posterior_samples = posterior_samples,
      core_data = core_data))
}