knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(echo = TRUE)
library(FPEMcountry)
devtools::load_all()
library(ggplot2)
library(grid)
library(gridExtra)

Table of Contents

  1. Introduction
  2. Run models
  3. Plot model results

Introduction

In this vignette we will run two models for each country available, one for women in union and another for women not in union. We will use the package default data fpemdata::contraceptive_use for this example.

Make directories for model files, import observed data and global model results

# vector of countries where we have data
global_divs <- FPEMcountry::index_m$index_area_df$division_numeric_code
divs <- dplyr::tibble(division_numeric_code = as.numeric(as.character(global_divs[global_divs %in% FPEMcountry::divisions$division_numeric_code])))
divs <- divs %>% dplyr::filter(division_numeric_code %in% FPEMcountry::population_counts$division_numeric_code)
 # dplyr::pull()
divnames <- divs %>% 
  dplyr::left_join(FPEMcountry::divisions) %>%
  dplyr::select(name_country) %>%
  dplyr::pull()
divs <- divs %>% dplyr::pull()


indicators <- c(
    "unmet_need_any",
    "contraceptive_use_modern",
    "contraceptive_use_traditional",
    "contraceptive_use_any"
    )

#stuck at div = 732

# divnamesordered <- divnames %>% sort()
# divsordered <- divs[match(divnamesordered, divnames)]
# divs <- divsordered
# divnames <- divnamesordered
# 
# global <- fpemdata::global_estimates_married %>%
#   dplyr::rename(division_numeric_code = Iso

run models

# devtools::load_all()
i <- 0
total <- length(divs)
for(div in divs) {
  fpem_one_country_autosave(
    runname = paste0(div),
    is_in_union = "ALL",
    service_stats = FALSE,
    division_numeric_code = div,
    first_year = 1970,
    last_year = 2030,
    nchains = 10,
    niter = 2500,
    nburnin = 500,
    nthin = 10
  )
  # fpem_results_autosave(runname = paste0(div))
  # fpem_plot_autosave(runname = paste0(div),
  #                         indicators = indicators,
  #                         compare_to_global = TRUE)
  i <- i + 1
  print(i/total)
}
files <- paste0("output/plots/", divs, ".pdf")
pdftools::pdf_combine(files, output = "joined.pdf")


FPRgroup/FPEMcountry documentation built on April 24, 2023, 4:32 p.m.