modelling.Rmd

---
title: 'Exposure and Fixed-Effects modelling '
output: html_notebook
---

# Introduction

This aim of this research is to build a series of models which include the following response variables: 

1. Male age-standardised mortality rates (all-cause)
2. Female age-standardised mortality rates (all-cause)
3. Relative difference in all-cause rates (1 / 2)
4. Absolute difference in all-cause rates (1 - 2) 
5. Male age-standardised mortality rates (NCDs)
6. Female age-standardised mortality rates (NCDs)
7. Relative difference in NCD rates (5 / 6)
8. Absolute differnece in NCD rates (5 - 6)

This will be done for each country (level 3 in location hierarchy). 

And a series of exposures from the World Bank will be used: 

1. GDP per cap by purchasing power parity 
2. Male literacy rates 
3. Female literacy rates 
4. etc etc 

# Loading GBD data 


- Source is [here](http://s3.healthdata.org/gbd-api-2016-production/3f1a784fe33bf6d9447c7b90e8828c3f_files/IHME-GBD_2016_DATA-3f1a784f-1.zip)


```{r}
require(tidyverse)
source("scripts/download_completed_request.R")
```

Will `read_csv_directly` work?

```{r, eval = F}
asr_dta <- read_csv_directly(
  url_loc = "http://s3.healthdata.org/gbd-api-2016-production/3f1a784fe33bf6d9447c7b90e8828c3f_files/IHME-GBD_2016_DATA-3f1a784f-1.zip",
  csv_name = "IHME-GBD_2016_DATA-3f1a784f-1.csv")

```

*Note:* Though `read_csv_directly()` works as a function, the links are not permanent, so downloading them is preferable. The aim should be to reproduce the above and check with the online webpage to confirm. 

```{r, eval = F}
download_completed_request("http://s3.healthdata.org/gbd-api-2016-production/67c84b9e938574bd56e66ab1c58653db_files",
"IHME-GBD_2016_DATA-67c84b9e-",
"data/age_standardised_allcause_ncd/"
                           )

```

```{r}
asr_dta <- read_csv("data/age_standardised_allcause_ncd/1.csv")
```

Now to link this to the hierarchy 


```{r}

gbd_locations <- readxl::read_excel("raw_data/IHME_GBD_2016_CODEBOOK/IHME_GBD_2016_GBD_LOCATION_HIERARCHY_Y2018M05D23.XLSX")

gbd_country_ids <- gbd_locations %>% 
  filter(location_id %in% 1:218) %>% 
  filter(level == 3) %>% 
  pull(location_id)

asr_dta %>% 
  filter(location_id %in% gbd_country_ids) %>% 
  filter(measure_name == "Deaths") %>% 
  select(location = location_name, year, sex = sex_name, cause = cause_name, std_mort_rate = val) -> mort_rate_data

asr_dta %>% 
  filter(location_id %in% gbd_country_ids) %>% 
  filter(measure_name == "DALYs (Disability-Adjusted Life Years)") %>% 
  select(location = location_name, year, sex = sex_name, cause = cause_name, std_mort_rate = val) -> daly_rate_data

```

And now to join this to the world bank measures 

```{r, cache = T}
wb_data <- read_rds("data/world_bank/extracted_indicators.RData")
nrow(wb_data)

gdp_percap_data <- wbstats::wb(indicator= "NY.GDP.PCAP.PP.CD") %>% as_tibble() 
nrow(gdp_percap_data)

wb_data <- bind_rows(wb_data, gdp_percap_data)
rm(gdp_percap_data)
nrow(wb_data)

```

join to mort rate

```{r} 
mort_rate_data_joined <- mort_rate_data %>%
  mutate(country = location) %>% 
  inner_join(
    wb_data %>% 
      mutate(year = as.integer(date)) %>% 
      select(year, country, indicator_id = indicatorID, indicator, indicator_value = value) 
  ) 

```

```{r} 
daly_rate_data_joined <- daly_rate_data %>%
  mutate(country = location) %>% 
  inner_join(
    wb_data %>% 
      mutate(year = as.integer(date)) %>% 
      select(year, country, indicator_id = indicatorID, indicator, indicator_value = value) 
  ) 

```



Models by gender and indicator



```{r}

pull_indicator_coeff <- function(x){x %>% filter(term == "indicator_value") %>% pull(estimate) }
pull_indicator_se <- function(x){x %>% filter(term == "indicator_value") %>% pull(std.error) }

simple_models <- mort_rate_data_joined %>% 
  filter(cause == "All causes") %>% 
  select(country, year, sex, std_mort_rate, indicator_id, indicator, indicator_value) %>% 
  filter(!is.na(indicator_value)) %>% 
  group_by(indicator_id) %>% 
  nest() %>% 
  mutate(male_model = map(
    data,
    function(x) {
      x %>%
      filter(sex == "Male") %>%
      lm(std_mort_rate ~ year + indicator_value, data = .) 
    }
      )
    ) %>%
  mutate(female_model = map(
    data, 
    function(x) {
      x %>% 
        filter(sex == "Female") %>% 
        lm(std_mort_rate ~ year + indicator_value, data = .) 
        }
      )
    ) %>% 
  mutate(
    male_model_tidy = map(male_model, broom::tidy),
    female_model_tidy = map(female_model, broom::tidy)     
  ) %>% 
  mutate(
    male_indicator_coeff        = map_dbl(male_model_tidy,   pull_indicator_coeff),
    female_indicator_coeff      = map_dbl(female_model_tidy, pull_indicator_coeff),
    male_indicator_se           = map_dbl(male_model_tidy,   pull_indicator_se),
    female_indicator_se         = map_dbl(female_model_tidy, pull_indicator_se)
    )



```

Note: for some of these indicators the log should probably be used 


```{r, results = "asis"}
simple_models %>% 
  select(indicator_id, male_indicator_coeff:female_indicator_se) %>% 
  left_join(wb_data %>% select(indicator_id = indicatorID, indicator) %>% distinct) %>%
  mutate(
    male_t = male_indicator_coeff / male_indicator_se, 
    female_t = female_indicator_coeff / female_indicator_se
    ) %>% 
  select(indicator_id, indicator, male_t, female_t) %>% 
  arrange(desc(abs(male_t))) %>%
  mutate(
    male_rank = 1 + length(male_t) - rank(abs(male_t)), 
    female_rank = 1 + length(female_t) - rank(abs(female_t))) %>% 
  filter(abs(male_t) > 2 | abs(female_t) > 2) %>%
  mutate(male_t = round(male_t, 2), female_t = round(female_t, 2)) %>% 
  mutate(diff = male_rank - female_rank) %>% 
  knitr::kable(caption = "Signal strength of indicators on male and female age-standardised all-cause mortality trends")

```

The above shows the 'signal strength' of the different indicators once year is controlled for. For both genders GDP per capital has the strongest signal, followed by gender parity in primary school enrollment. 
For males, age-standardised mortality rates are then most strongly predicted by the female urban population (% of total), followed by male urban population as % of total. 
For females, the third strongest signal is gender parity in primary and secondary school enrollment, followed by gender parity in secondary school (only) enrollment rates.   

# As above, but CVD only

```{r}
produce_table <- function(DTA, this_cause, this_caption){

  DTA %>% 
  filter(cause ==  this_cause) %>% 
  select(country, year, sex, std_mort_rate, indicator_id, indicator, indicator_value) %>%
  filter(!is.na(indicator_value)) %>%
  group_by(indicator_id) %>%
  nest() %>%
  mutate(male_model = map(
    data,
    function(x) {
      x %>%
      filter(sex == "Male") %>%
      lm(std_mort_rate ~ year + indicator_value, data = .)
    }
      )
    ) %>%
  mutate(female_model = map(
    data,
    function(x) {
      x %>%
        filter(sex == "Female") %>%
        lm(std_mort_rate ~ year + indicator_value, data = .)
        }
      )
    ) %>%
  mutate(
    male_model_tidy = map(male_model, broom::tidy),
    female_model_tidy = map(female_model, broom::tidy)
  ) %>%
  mutate(
    male_indicator_coeff        = map_dbl(male_model_tidy,   pull_indicator_coeff),
    female_indicator_coeff      = map_dbl(female_model_tidy, pull_indicator_coeff),
    male_indicator_se           = map_dbl(male_model_tidy,   pull_indicator_se),
    female_indicator_se         = map_dbl(female_model_tidy, pull_indicator_se)
    ) %>%
  select(indicator_id, male_indicator_coeff:female_indicator_se) %>%
  left_join(wb_data %>% select(indicator_id = indicatorID, indicator) %>% distinct) %>%
  mutate(
    male_t = male_indicator_coeff / male_indicator_se,
    female_t = female_indicator_coeff / female_indicator_se
    ) %>%
  select(indicator_id, indicator, male_t, female_t) %>%
  arrange(desc(abs(male_t))) %>%
  mutate(
    male_rank = 1 + length(male_t) - rank(abs(male_t)),
    female_rank = 1 + length(female_t) - rank(abs(female_t))) %>%
  filter(abs(male_t) > 2 | abs(female_t) > 2) %>%
  mutate(male_t = round(male_t, 2), female_t = round(female_t, 2)) %>%
  mutate(diff = male_rank - female_rank) %>%
  knitr::kable(caption = this_caption)

  
}

```


```{r, results = "asis"}

mort_rate_data_joined %>% 
  produce_table(
    this_cause = "Non-communicable diseases", 
    this_caption = "Signal strength of indicators on male and female age-standardised NCD mortality trends"
    )


```

Gender differences in the signal strength of these indicators appear much larger for NCD mortality than all cause mortality. Note also that the direction of effects goes in the opposite direction for some indicators, though only for the Mo Ibrahim index are the effects in opposite directions and the t values > 2. 

The rank importance of access to anti-retroviral drugs is increased, to 2nd and 3rd place after GDP per capita. 


# Signal strength, all cause, DALY

Now to look at the equivalent for DALYs 


```{r, results = "asis"}
daly_rate_data_joined %>% 
  produce_table(
    this_cause = "All causes", 
    this_caption = "Signal strength of indicators on male and female age-standardised all-cause DALY rate trends")



```

For DALY rate trends, the signal strength of school enrollment gender parity is stronger for females than that of GDP per capita. 

Now, finally, to look at DALY rates related to NCDs


```{r, results = "asis"}

daly_rate_data_joined %>% 
  produce_table(
    this_cause = "Non-communicable diseases", 
    this_caption = "Signal strength of indicators on male and female age-standardised NCD DALY rate trends")


```

Let's now try to regress this on abs or rel difference


```{r}
diffs_mort_joined <- mort_rate_data_joined %>% 
  select(location, year, sex, cause, std_mort_rate) %>% 
  distinct() %>% 
  spread(sex, std_mort_rate) %>% 
  mutate(rel = Male / Female, abs = Male - Female, logrel = log(rel)) %>% 
  select(location, year, cause, rel, abs, logrel) %>% 
  mutate(country = location) %>% 
  inner_join(
    wb_data %>% 
      mutate(year = as.integer(date)) %>% 
      select(year, country, indicator_id = indicatorID, indicator, indicator_value = value) 
  ) 

```


```{r}
diffs_daly_joined <- daly_rate_data_joined %>% 
  select(location, year, sex, cause, std_mort_rate) %>% 
  distinct() %>% 
  spread(sex, std_mort_rate) %>% 
  mutate(rel = Male / Female, abs = Male - Female, logrel = log(rel)) %>% 
  select(location, year, cause, rel, abs, logrel) %>% 
  mutate(country = location) %>% 
  inner_join(
    wb_data %>% 
      mutate(year = as.integer(date)) %>% 
      select(year, country, indicator_id = indicatorID, indicator, indicator_value = value) 
  ) 

```

Now produce similar tables as above 

```{r}
produce_diff_table <- function(DTA, this_cause, this_caption){

  DTA %>% 
    gather(rel:logrel, key = "outcome", value = "outcome_value") %>% 
    filter(cause ==  this_cause) %>% 
  select(country, year, outcome, outcome_value, indicator_id, indicator, indicator_value) %>%
  filter(!is.na(indicator_value)) %>%
  group_by(indicator_id) %>%
  nest() %>%
  mutate(rel_model = map(
    data,
    function(x) {
      x %>%
      filter(outcome == "rel") %>%
      lm(outcome_value ~ year + indicator_value, data = .)
    }
      )
    ) %>%
  mutate(abs_model = map(
    data,
    function(x) {
      x %>%
        filter(outcome == "abs") %>%
        lm(outcome_value ~ year + indicator_value, data = .)
        }
      )
    ) %>%
  mutate(
    rel_model_tidy = map(rel_model, broom::tidy),
    abs_model_tidy = map(abs_model, broom::tidy)
  ) %>%
  mutate(
    rel_indicator_coeff        = map_dbl(rel_model_tidy,   pull_indicator_coeff),
    abs_indicator_coeff      = map_dbl(abs_model_tidy, pull_indicator_coeff),
    rel_indicator_se           = map_dbl(rel_model_tidy,   pull_indicator_se),
    abs_indicator_se         = map_dbl(abs_model_tidy, pull_indicator_se)
    ) %>%
  select(indicator_id, rel_indicator_coeff:abs_indicator_se) %>%
  left_join(wb_data %>% select(indicator_id = indicatorID, indicator) %>% distinct) %>%
  mutate(
    rel_t = rel_indicator_coeff / rel_indicator_se,
    abs_t = abs_indicator_coeff / abs_indicator_se
    ) %>%
  select(indicator_id, indicator, rel_t, abs_t) %>%
  arrange(desc(abs(rel_t))) %>%
  mutate(
    rel_rank = 1 + length(rel_t) - rank(abs(rel_t)),
    abs_rank = 1 + length(abs_t) - rank(abs(abs_t))) %>%
  filter(abs(rel_t) > 2 | abs(abs_t) > 2) %>%
  mutate(rel_t = round(rel_t, 2), abs_t = round(abs_t, 2)) %>%
  knitr::kable(caption = this_caption)

  
}


```

```{r, results = 'asis'}
diffs_mort_joined %>% 
  produce_diff_table(
    this_cause = "All causes", 
    this_caption = "Signal strength of indicators: gender inequalities in health; mort due to all causes")


```



```{r, results = 'asis'}
diffs_mort_joined %>% 
  produce_diff_table(
    this_cause = "Non-communicable diseases", 
    this_caption = "Signal strength of indicators: gender inequalities in health; mort due to NCDs")


```


```{r, results = 'asis'}
diffs_daly_joined %>% 
  produce_diff_table(
    this_cause = "All causes", 
    this_caption = "Signal strength of indicators: gender inequalities in health; DALYs due to all causes")


```

```{r, results = 'asis'}
diffs_daly_joined %>% 
  produce_diff_table(
    this_cause = "Non-communicable diseases", 
    this_caption = "Signal strength of indicators: gender inequalities in health; DALYs due to NCDs")


```


Now a CART to group countries into quintile of gender inequality in outcomes

```{r}
tmp_dta <- diffs_mort_joined %>% 
  filter(cause == "All causes") %>% 
  mutate(quint_rel = cut(rel, breaks = quantile(rel, seq(0, 1, by = 0.2)), labels = F) %>% as.factor() ) %>% 
  select(location,  year, cause, quint_rel, indicator_id, indicator_value) %>% 
  spread(indicator_id, indicator_value) 

mod_simple <- tmp_dta %>% 
  rpart::rpart(quint_rel ~ year + NY.GDP.PCAP.PP.CD, data = ., method = "class")

# plot(model_01)
# text(model_01, cex = 0.8, use.n = T, xpd = T)

rpart.plot::prp(mod_simple, faclen = 0, cex = 0.8, varlen = 0, extra = 2, under = T)

mod_full <- tmp_dta %>% 
  select(-location) %>% 
  rpart::rpart(quint_rel ~ ., data = ., method = "class")

rpart.plot::prp(mod_full, faclen = 0, cex = 0.8, varlen = 0, extra = 2, under = T)


```

|indicator_id         |indicator                                                                             |
|:--------------------|:-------------------------------------------------------------------------------------|
| NY.GDP.PCAP.PP.CD | GDP per capita, PPP (current international $)  |             
| SE.ENR.SECO.FM.ZS | School enrollment, secondary (gross), gender parity index (GPI)   |    
| MO.INDEX.XQ |   Overall Mo Ibrahim index  |          
| SL.UEM.1524.FM.ZS | Ratio of female to male youth unemployment rate (% ages 15-24) (modeled ILO estimate) | 
| SP.URB.TOTL.MA.ZS | Urban population, male (% of total) | 
| SE.ENR.TERT.FM.ZS | School enrollment, tertiary (gross), gender parity index (GPI) | 

So, the above model seems to be suggesting:

* If GDP per cap < $3k, gender parity in 2ndry school enrolment is low, and Mo Ibrahim index high -> lowest/negative male gender disadvantage
* If GDP per cap < $3k, gender parity in 2ndry school enrolment is now, and Mo Ibrahim index low -> 2nd quintile 
* If GDP per cap < $3k, gender parity in 2ndry school enrolment is high -> 3rd quintile
* If GDP per cap < $13k -> 3rd quintile
* If GDP per cap > $13k, high female:male youth unemployment rate -> 3rd quintile
* If GDP per cap > $13k, < 38% of males in urban areas (?) -> 5th quintile
* If GDP per cap > $13k, >= 38% of males in urban areas (?) but low gender parity index in tertiary school enrolment -> 5th quintile
* If GDP per cap > $13k, >= 38% of males in urban areas, high gender parity in tertiary school enrolment, but female:male youth unemployment >= 110 -> 5th quintile
* Otherwise 4th quintile

Now the same for NCDs


```{r}
tmp_dta <- diffs_mort_joined %>% 
  filter(cause == "Non-communicable diseases") %>% 
  mutate(quint_rel = cut(rel, breaks = quantile(rel, seq(0, 1, by = 0.2)), labels = F) %>% as.factor() ) %>% 
  select(location,  year, cause, quint_rel, indicator_id, indicator_value) %>% 
  spread(indicator_id, indicator_value) 

mod_simple <- tmp_dta %>% 
  rpart::rpart(quint_rel ~ year + NY.GDP.PCAP.PP.CD, data = ., method = "class")

# plot(model_01)
# text(model_01, cex = 0.8, use.n = T, xpd = T)

rpart.plot::prp(mod_simple, faclen = 0, cex = 0.8, varlen = 0, extra = 2, under = T)

mod_full <- tmp_dta %>% 
  select(-location) %>% 
  rpart::rpart(quint_rel ~ ., data = ., method = "class")

rpart.plot::prp(mod_full, faclen = 0, cex = 0.8, varlen = 0, extra = 2, under = T)


```


But what are these thresholds?

```{r}
diffs_mort_joined %>% 
  filter(cause == "All causes") %>% 
  pull(rel) %>% 
  quantile(seq(0, 1, by = 0.2)) -> tmp1

diffs_mort_joined %>% 
  filter(cause == "Non-communicable diseases") %>% 
  pull(rel) %>% 
  quantile(seq(0, 1, by = 0.2)) -> tmp2

tab <- rbind(tmp1, tmp2); rm(tmp1, tmp2)
row.names(tab) <- c("All causes", "NCDs")

tab


```