Chapter 3 Analysis Wealth

3.1 Processing - Wealth

rm(list = ls())
library(IC2)
library(codebook)
library(tidyverse)

dat <- read.csv("./_dat/dat_mal_ineq_v3.csv")
source("./mal_ineq_fun.R")

# epiDisplay::tab1(dat$rapid_test)

#codebook_browser(dat)


# ==== NOT RUN =======
# ao <- dat %>%
#   filter(country == "SN") %>%
#   dplyr::select(rapid_test, wealth, w) %>%
#   filter(complete.cases(.))
# 
# ao_ci <- IC2::calcSConc(x = ao$rapid_test, 
#                         y = ao$wealth, 
#                         w = ao$w)
# ==== NOT RUN =======

dat.n <- dat %>%
  group_by(country, REGNAME, cluster_number) %>%
  nest() %>%
  mutate(sample.size.N = map(.x = data, .f = ~dplyr::select(.x, rapid_test, wealth, w) %>%
                               count(name = "N") %>% unlist()),
         sample.size.m = map(.x = data, .f = ~dplyr::select(.x, rapid_test, wealth, w) %>%
                               filter(!is.na(rapid_test)) %>%
                               count(name = "N_m") %>% unlist()),
         dat_s = map(.x = data, .f = ~dplyr::select(.x, rapid_test, wealth, w) %>%
                       filter(complete.cases(.))),
         sample.size = map(.x = dat_s, .f = ~count(.x) %>% unlist()),
         prev = map(.x = dat_s, .f = ~mean(.x$rapid_test, na.rm=T)),
         wealth_cats = map(.x = dat_s, .f = ~distinct(.x, wealth) %>% nrow())) %>%
  # FILTERS =======================
  filter(prev > 0 & prev < 1) %>%
  filter(sample.size>=10) %>%
  filter(wealth_cats>1)

dat_ci.n <- dat.n  %>%
  mutate(ci = map(.x = dat_s, .f = ~IC2::calcSConc(x = .x$rapid_test,
                        y = .x$wealth,
                        w = .x$w)),
         ci_val = map(.x = ci, .f = ~.x$ineq$index),
         h_calc = map(.x = dat_s, .f = ~h_ineq(dat = .x, var_soc = wealth, var_outcome = rapid_test))
         )
  
dat_ci <- dat_ci.n %>%
  dplyr::select(country, cluster_number, sample.size.N, sample.size.m, sample.size, prev, wealth_cats, ci_val, h_calc) %>%
  unnest() %>%
  ungroup() 

# Checks ========
dat_ci %>%
  ggplot(aes(x=c_index, y = ci_val)) +
  geom_point(alpha = .4) +
  labs(x = "hand calculation", y = "package calculation")

# ==== NOT RUN =======
# unique(dat.n$country)
# unique(dat_ci$country)
# 
# range(dat$w)
# range(dat_ci$N_m)
# range(dat_ci$n)
# range(dat_ci$ci_val, na.rm = T)
# 
# hist(dat_ci$N_m)
# hist(dat_ci$n)
# hist(dat_ci$ci_val)
# ==== NOT RUN =======

#saveRDS(dat_ci, "./_dat/dat_ci.rds")

3.2 Spatial data - Wealth

## Reading layer `DHS_adm' from data source `/Users/gabrielcarrasco/Dropbox/Work/Tarik LAB/Malaria Ineq/mal_ineq/_dat/SHP/union/DHS_adm.shp' using driver `ESRI Shapefile'
## Simple feature collection with 147 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13.30198 ymin: -26.86819 xmax: 50.49459 ymax: 15.7047
## geographic CRS: WGS 84

3.3 Distributions - Wealth

library(cowplot)

a <- hist_ineq(dat_ci, sample.size, "Sample Size")
b <- hist_ineq(dat_ci, prev, "Prevalence")
c <- bi_hist_ineq(dat_ci, var_x = prev, var_y = sample.size, lab_x = "Prevalence", lab_y = "Sample Size")

(dist <- plot_grid(a,b,c, labels = c("A)", "B)", "C)"), nrow = 1))

3.4 CI - Wealth

3.4.1 Plots CI - Wealth

library(cowplot)

ci_box_10 <- dat_ci %>% filter(sample.size>=10) %>% ci_box(var = ci_val)
ci_box_20 <- dat_ci %>% filter(sample.size>=20) %>% ci_box(var = ci_val)
ci_box_30 <- dat_ci %>% filter(sample.size>=30) %>% ci_box(var = ci_val)

(ci_box1 <- plot_grid(ci_box_10, ci_box_20, ci_box_30, labels = c("A)", "B)", "C)"), ncol = 1))

ci_prev_20 <- dat_ci %>% filter(sample.size>=20) %>% ci_prev(var = ci_val)
ci_prev_10 <- dat_ci %>% filter(sample.size>=10) %>% ci_prev(var = ci_val)
ci_prev_30 <- dat_ci %>% filter(sample.size>=30) %>% ci_prev(var = ci_val)

(ci_prev1 <- plot_grid(ci_prev_10, ci_prev_20, ci_prev_30, labels = c("A)", "B)", "C)"), ncol = 1))

3.4.2 Maps CI - Wealth

library(mapview)

m1 <- dat_ci_map %>%
  filter(!is.na(ci_val)) %>%
  mapview(zcol = "ci_val", legend = TRUE, layer.name = "CI (psu)")

m2 <- dat_ci_adm %>%
  mapview(zcol = "ci_val", legend = TRUE, layer.name = "CI (Adm)")

m1 + m2