injury_tables.Rmd

---
title: "Injrury Outcomes - WHW matrices"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    toc: true
    toc_depth: 5
    toc_float: true
---

```{r setup, include=FALSE}
library(knitr)
library(pracma)
library(summarytools)
library(dplyr)
library(purrr)
library(flextable)
library(ggplot2)
library(tidyverse)
knitr::opts_chunk$set(comment=NA, prompt=FALSE, cache=FALSE, echo=F, message = F, warning = F, results='asis')

smodes <- read_csv('data/global/modes/standardized_modes.csv')

```

```{js echo=FALSE}
window.location.href='#Who_Hit_Whom_matrix_(normalized_to_annual_figures)';
```


```{r loadLibraries, echo = F, message = F}

st_options(bootstrap.css     = FALSE,       # Already part of the theme so no need for it
           plain.ascii       = FALSE,       # One of the essential settings
           style             = "rmarkdown", # Idem.
           dfSummary.silent  = TRUE,        # Suppresses messages about temporary files
           footnote          = NA,          # Keeping the results minimalistic
           subtitle.emphasis = FALSE)       # For the vignette theme, this gives
# much better results. Your mileage may vary

```

```{r load_objects = "asis"}
io <- readRDS("results/multi_city/io.rds")
# Assumes that multi_city_script.R has been run till 
# Get names of cities from the io object
cities <- names(io)[!names(io) %in% 'scen_prop']

round_to <- 1

# Set plot_mmets to F
plot_mmets <- F

sum_and_round_and_print <- function(data,text=''){
  data <- lapply(data, function(x)round(x,round_to))
  data <- lapply(data,function(x)rbind(x,Total=colSums(x)))
  for(city in cities) {
    print(kable(data[[city]], caption = paste(text, city)))
    cat("\n")
  }
}
round_and_print <- function(data,text=''){
  data <- lapply(data, function(x)round(x,round_to))
  for(city in cities) {
    print(kable(data[[city]], caption = paste(text, city)))
    cat("\n")
  }
}

```


```{r preprocessing}


collapse_ages <- function(data,ages=c('15-49','50-69'),age_lab='age_cat'){
  target_min_age <- as.numeric(sapply(ages,function(x)strsplit(x,'-')[[1]][1]))
  target_max_age <- as.numeric(sapply(ages,function(x)strsplit(x,'-')[[1]][2]))
  min_ages <- as.numeric(sapply(data[[age_lab]],function(x)strsplit(x,'-')[[1]][1]))
  max_ages <- as.numeric(sapply(data[[age_lab]],function(x)strsplit(x,'-')[[1]][2]))
  genders <- unique(data$sex)
  if(ncol(data)>3) {
    reformatted <- do.call(rbind,lapply(1:length(ages),
                                        function(x) t(sapply(genders,
                                                             function(y) 
                                                               if(ncol(data)>3) colSums(data[min_ages<=target_max_age[x]&max_ages>=target_min_age[x]&data$sex==y,-c(1,2),drop=F])
                                                             else sum(data[min_ages<=target_max_age[x]&max_ages>=target_min_age[x]&data$sex==y,-c(1,2),drop=F])
                                        ))
    ))
    data.frame(age_cat=rep(ages,each=2),sex=rep(genders,2),reformatted,stringsAsFactors=F)
  }else{
    reformatted <- do.call(c,lapply(1:length(ages),
                                    function(x) sapply(genders,
                                                       function(y) 
                                                         sum(data[min_ages<=target_max_age[x]&max_ages>=target_min_age[x]&data$sex==y,-c(1,2),drop=F])
                                    )
    ))
    data.frame(age=rep(ages,each=2),sex=rep(genders,2),population=as.numeric(reformatted),stringsAsFactors=F)
  }
}

scen_prop <- io$scen_prop
io <- io[-2]
for(i in 1:length(io)) {
  io[[i]]$demographic$age[io[[i]]$demographic$age=='5-9'] <- '05-9'
  io[[i]]$outcomes$hb$deaths$age_cat[io[[i]]$outcomes$hb$deaths$age_cat=='5-9'] <- '05-9'
  io[[i]]$outcomes$hb$ylls$age_cat[io[[i]]$outcomes$hb$ylls$age_cat=='5-9'] <- '05-9'
  io[[i]]$outcomes$pathway_hb$deaths$age_cat[io[[i]]$outcomes$pathway_hb$deaths$age_cat=='5-9'] <- '05-9'
  io[[i]]$outcomes$pathway_hb$ylls$age_cat[io[[i]]$outcomes$pathway_hb$ylls$age_cat=='5-9'] <- '05-9'
  io[[i]]$outcomes$hb$deaths <- io[[i]]$outcomes$hb$deaths[,!sapply(names(io[[i]]$outcomes$hb$deaths),function(x)grepl('ac|neo',as.character(x)))]
  io[[i]]$outcomes$hb$ylls <- io[[i]]$outcomes$hb$ylls[,!sapply(names(io[[i]]$outcomes$hb$ylls),function(x)grepl('ac|neo',as.character(x)))]
  io[[i]]$outcomes$pathway_hb$deaths <- io[[i]]$outcomes$pathway_hb$deaths[,!sapply(names(io[[i]]$outcomes$pathway_hb$deaths),function(x)grepl('ac|neo',as.character(x)))]
  io[[i]]$outcomes$pathway_hb$ylls <- io[[i]]$outcomes$pathway_hb$ylls[,!sapply(names(io[[i]]$outcomes$pathway_hb$ylls),function(x)grepl('ac|neo',as.character(x)))]
}
for(city in cities)
  for(type in c('hb','pathway_hb'))
    for(out in c('deaths','ylls'))
      io[[city]]$outcomes[[type]][[out]] <- collapse_ages(io[[city]]$outcomes[[type]][[out]])
for(city in cities) io[[city]]$demographic <- collapse_ages(io[[city]]$demographic,age_lab='age')
pop_by_age <- lapply(io,function(x)sapply(unique(x$demographic$age),function(y)sum(subset(x$demographic,age==y)$population)))
pop_by_gender <- lapply(io,function(x)sapply(unique(x$demographic$sex),function(y)sum(subset(x$demographic,sex==y)$population)))
injury_col <- which(colnames(io[[1]]$outcomes$hb$deaths)=='scen1_deaths_inj')
ap_cols <- which(sapply(colnames(io[[1]]$outcomes$pathway_hb$deaths),function(x)grepl('ap',as.character(x))))
pa_cols <- which(sapply(colnames(io[[1]]$outcomes$pathway_hb$deaths),function(x)grepl('pa',as.character(x))))
scen_names <- rownames(scen_prop)
# scen_names <- c('Walking', 'Bicycling', 'Driving', 'Public Transport', 'Motorcycling')


```

# Who Hit Whom matrix (normalized to annual figures)

Who hit whom matrix for all scenarios for each case study

## Who hit whom matrix of each case study (by baseline + scenarios) {#whw}

```{r whw, results = 'asis', fig.width=7, fig.height=4, echo=FALSE}

scen_names_w_baseline <- c('bl', 'w_sc', 'bi_sc', 'car_sc', 'mc_sc', 'bus_sc')
scen_lt <- data.frame("scen_name" = names(io[[1]]$outcomes$whw), "val" = scen_names_w_baseline, stringsAsFactors = FALSE)

all_whw <- list()

for(city in cities){
  {
    cat('\n')
    
    cat('###', stringr::str_to_title(stringr::str_replace(city, '_', ' ')), '\n')

    for (cs in names(io[[1]]$outcomes$whw)){
      
      if (is.null(io[[city]]$outcomes$whw[[cs]]$whw)){
        if (!is.null(io[[city]]$outcomes$whw[[cs]]$nov)){
          td2 <- t(io[[city]]$outcomes$whw[[cs]]$nov) %>% as.data.frame()
          td2$mode <- 'NOV'
          
          td2 <- td2 %>% dplyr::select(mode, names(.))
          
          td3 <- td2
        }
        
      }else{
        td1 <- (io[[city]]$outcomes$whw[[cs]]$whw) %>% as.data.frame() %>% tibble::rownames_to_column("mode")
        td3 <- td1
        if (!is.null(io[[city]]$outcomes$whw[[cs]]$nov)){
          td2 <- t(io[[city]]$outcomes$whw[[cs]]$nov) %>% as.data.frame()
          td2$mode <- 'NOV'
          
          td3 <- plyr::rbind.fill(td1, td2)
        }
      }
      td3 <- td3 %>% mutate(rowSum = rowSums(.[2:ncol(td3)], na.rm = T))
      td3 <- td3 %>% janitor::adorn_totals("row")
      td3[, 2:ncol(td3)] <- round(td3[, 2:ncol(td3)], 2)
      
      var <- scen_lt %>% filter(scen_name == cs) %>% dplyr::select(val) %>% as.character()
      
      # scen_names <- c('Walking', 'Bicycling', 'Driving', 'Public Transport', 'Motorcycling')
      
      
      if (var == 'bl'){
        qualified_scen_name <- 'Baseline'
        scen <- ''
      }else if(var == "w_sc"){
        qualified_scen_name <- 'Walking'
      }else if(var == "bi_sc"){
        qualified_scen_name <- 'Bicycling'
      }else if(var == "car_sc"){
        qualified_scen_name <- 'Driving'
      }else if(var == "mc_sc"){
        qualified_scen_name <- 'Motorcycling'
      }else if(var == "bus_sc"){
        qualified_scen_name <- 'Public Transport'
      }
      
      cat('\n')
      
      cat('####', qualified_scen_name, ' WHW', '\n')
      
      cat('\n')
      
      td3 <- td3[,c(1, na.omit(match(smodes$exhaustive_list, colnames(td3))), ncol(td3))]
      
      td3 <- td3[c(na.omit(match(smodes$exhaustive_list, td3$mode)), nrow(td3)),]
      
      temp <- td3 %>% mutate_if(is.numeric, .funs = funs(case_when( . < 1 ~ round(., 2) ,  . >= 1 ~  as.numeric(round(.))))) %>%  mutate_if(is.numeric, ~as.character(.))
      
      ft <- flextable(temp)
      ft <- add_header_row(ft, values = c('Strike Mode', 'Casualty Mode'), colwidths = c(1, ncol(temp) - 1))
      
      cat(knit_print(ft))
      
      cat("\n")
      
      td <- td3
      names(td)[2:ncol(td)] <- paste(names(td)[2:ncol(td)], city, sep = "_")
      
      all_whw[[var]][[city]] <- td
    }
    
  }
}

```


```{r whw_for_all_cs, results = 'asis', fig.width=7, fig.height=4, echo=FALSE}

st <- list()

file_list <- list()

td <- NULL

whw_list <- list()

for (cs in names(all_whw)){
  # cs <- names(all_whw)[1]
  
  td <- all_whw[[cs]] %>% purrr::reduce(full_join, by = "mode") %>% as.data.frame() %>% dplyr::select(mode, sort(names(.)))
  
  td[is.na(td)] <- 0
  
  td <- (td[!duplicated(td), ])
  
  td <- rbind(td %>% dplyr::filter(!mode %in% c("Total", "NOV")) %>% arrange(mode), td %>% filter(mode == 'NOV'), td %>% filter(mode == 'Total'))
  
  td[is.na(td)] <- 0
  
  td <- td %>% mutate_if(is.numeric, round, 2)
  
  backup_td <- td
  
  colnames(td) = gsub("auto_rickshaw", "ar", colnames(td))
  
  whw_lng <- reshape2::melt(td)
  
  col_split <- stringr::str_split(whw_lng$variable, "_", simplify = TRUE, n = 2)
  
  whw_lng <- cbind(whw_lng, col_split)
  names(whw_lng)[4] <- 'strike_mode'
  names(whw_lng)[5] <- 'city'
  whw_lng <- whw_lng %>% dplyr::select(-variable)
  whw_lng <- whw_lng %>% dplyr::filter(strike_mode != "rowSum")
  whw_lng <- whw_lng %>% rename(str_mode = mode, cas_mode = strike_mode)
  
  if (nrow(whw_lng[whw_lng$cas_mode == 'ar',]) > 0){
    whw_lng[whw_lng$cas_mode == 'ar',]$cas_mode <- "auto_rickshaw"
  }
  
  if (nrow(whw_lng[whw_lng$str_mode == 'ar',]) > 0){
    whw_lng[whw_lng$str_mode == 'ar',]$str_mode <- "auto_rickshaw"
  }
  
  whw_lng$str_mode <- factor(whw_lng$str_mode, levels = unique(whw_lng$str_mode))
  
  qualified_scen_name <- 'Public Transport'
  scen <- 'Scenario'
  
  if (cs == 'bl'){
    qualified_scen_name <- 'Baseline'
    scen <- ''
  }else if(cs == "w_sc"){
    qualified_scen_name <- 'Walking'
  }else if(cs == "bi_sc"){
    qualified_scen_name <- 'Bicycling'
  }else if(cs == "car_sc"){
    qualified_scen_name <- 'Driving'
  }else if(cs == "mc_sc"){
    qualified_scen_name <- 'Motorcycling'
  }
  
  whw_lng$scenario <- qualified_scen_name
  
  if (length(whw_list) == 0)
    whw_list <- whw_lng
  else
    whw_list <- rbind(whw_list, whw_lng)
  
  f <- ggplot(data = whw_lng) +
    aes(x = str_mode, fill = cas_mode, weight = value) +
    geom_bar(position = "stack") +
    scale_fill_brewer(palette = "Dark2") +
    labs(x = "Strike Mode", y = "# of incidents", title = paste("WHW - ", qualified_scen_name, scen)) +
    coord_flip() +
    facet_wrap(vars(city), scales = "free_y")
  
  file_list[[cs]] <- f
  
  
  #ggsave(f, file = paste0("results/multi_city/whw_matrices/plots/whw_", cs,".png"), limitsize = F, dpi = 300, scale = 1.5)
  
  d <- paste0("results/multi_city/whw_matrices/interactive_plots/whw_", cs,".html")
  
  # htmlwidgets::saveWidget(plotly::ggplotly(f), file.path(normalizePath(dirname(d)), basename(d)))
  #ggsave(f, file=paste0("results/multi_city/whw_matrices/whw_", cs,".png"), width = 14, height = 10, units = "cm")
  
  qual_name <- paste(qualified_scen_name, scen)
  st[[cs]] <- format(td, scientific = F)
  
}

# Remove total from whw long form, and write as csv
readr::write_csv(whw_list %>% filter(str_mode != 'Total'), paste0('results/multi_city/whw_matrices/whw_lng.csv'))

```

## Who hit whom matrix for all case cities {#whw_all}


```{r whw_print, results = 'asis', fig.width=7, fig.height=4, echo=FALSE}

list_temp <- list()
index <- 1

scen_names <- unique(whw_list$scenario)
for (i in 1:length(scen_names)){
  sn <- scen_names[i]
  cat('###', sn, ' scenario', '\n')
  
  temp <- whw_list %>% filter(scenario == sn) %>% dplyr::select(cas_mode)
  
  cas_mode <- unique(temp$cas_mode)
  
  for (j in 1:length(cas_mode)){

    cm <- cas_mode[j]
    temp1 <- whw_list %>% filter(scenario == sn & cas_mode == cm) %>% spread(value = value, key = str_mode)
    
    temp1 <- temp1[,c(1, 2, 3, na.omit(match(smodes$exhaustive_list, colnames(temp1))), ncol(temp1))]
    
    temp1 <- temp1 %>% mutate_if(is.numeric, .funs = funs(case_when( . < 1 ~ round(., 2) ,  . >= 1 ~  as.numeric(round(.))))) %>%  mutate_if(is.numeric, ~as.character(.))
    
    list_temp[[index]] <- temp1
    
    index <- index + 1
    cat('\n')
    cat('#### Casualty mode: ', stringr::str_to_title(cm), '\n')
    #print(kable(format(temp1, scientific = F), caption = paste('Who Hit Whom (WHW) for ', scen_names[i], 'scenario and ', cas_mode[j], ' mode')))
    
    ft <- flextable(temp1)
    ft <- add_header_row(ft, values = c(' ', 'strike mode'), colwidths = c(3, ncol(temp1) - 3))
    ft <- merge_v(ft, j = 1:3, part = 'body')
    
    cat('\n')
    
    cat(knit_print(ft))
    
    cat('\n')
    
  }
  
  cat('\n')
  
}

```

## Who hit whom matrix for all case cities as CSV {#whw_all_csv}

[Download from](https://raw.githubusercontent.com/ITHIM/ITHIM-R/master/results/multi_city/whw_matrices/whw_lng.csv)



### Change in deaths due to injury per billion km

```{r inj_100k = "asis", echo = F, message = F}
require(tibble)
overall_el <- list() # To save all cities results
overall_dd <- list() # To save all cities results
for (city in cities){
  
  #city <- 'delhi'
  
  # print(city)
  
  el <- list() # To save city specific results
  
  #if (!city %in% c('buenos_aires', 'mexico_city', 'bogota', 'vizag')){
  
  for (cs in names(io[[city]]$outcomes$whw)){
    if (length(names(io[[city]]$outcomes$whw$Baseline)) == 2){
      # td1: Number of fatalities in nov
      td1 <- round(io[[city]]$outcomes$whw[[cs]]$nov) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      # td2: Number of fatalities in whw
      td2 <- colSums(round(io[[city]]$outcomes$whw[[cs]]$whw)) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      # td3: Number of fatalities: sum of nov and whw
      td3 <- full_join(td2, td1, by = 'mode') %>% mutate(count = rowSums(.[2:3], na.rm = T)) %>% dplyr::select(-c('count.x', 'count.y'))
      #td3[td3$mode == 'pedestrian', ]$mode <- 'walking'
      
    }else if(length(names(io[[city]]$outcomes$whw$Baseline)) == 1 && names(io[[city]]$outcomes$whw$Baseline) == 'whw'){
      td3 <- colSums(round(io[[city]]$outcomes$whw[[cs]]$whw)) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      #td3[td3$mode == 'pedestrian', ]$mode <- 'walking'
    }else if(length(names(io[[city]]$outcomes$whw$Baseline)) == 1 && names(io[[city]]$outcomes$whw$Baseline) == 'nov'){
      
      td3 <- round(io[[city]]$outcomes$whw[[cs]]$nov) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      #td3[td3$mode == 'pedestrian', ]$mode <- 'walking'
      
    }
    # td4: distance travelled per mode
    td4 <- io[[city]]$true_dist %>% filter(stage_mode %in% td3$mode) %>% dplyr::select(stage_mode, cs) %>% as.data.frame()
    
    if (length(el) == 0){
      el <- td4 %>% dplyr::select(stage_mode)
      dd <- el
    }
    # td4: merge total fatalities to distance
    td4 <- full_join(td4, td3 %>% dplyr::select(mode, count) %>% rename(stage_mode = mode), by = 'stage_mode') 
    
    var <- scen_lt %>% filter(scen_name == cs) %>% dplyr::select(val) %>% as.character()
    
    names(td4)[2] <- var
    
    td5 <- td4
    
    td4[, 2] <- as.numeric(td4[,2])
    # Compute risk per bn km travelled
    td4[, 2] <- round((td4[,3] / ( td4[,2] * 365)) * 
                        1000000000, 4)
    
    names(td4)[3] <- paste(names(td4)[2], names(td4)[3], sep = "_")
    
    el <- full_join(el, td4, by = 'stage_mode')
    
    # City's population
    pop <- sum(io[[city]]$demographic$population)
    # Compute distance travelled per capita
    td6 <- td5
    td6[[paste0(var,"_risk")]] <- unlist(td4[2])
    td6[[paste0(var,"_percapita")]] <- td6[[var]] * 365 / pop
    names(td6)[3] <- paste(names(td6)[2], names(td6)[3], sep = "_")
    dd <- full_join(dd, td6, by = 'stage_mode')
  }
  print(kable(el, caption = city))
  
  td <- el %>% dplyr::select(-contains('count'))
  names(td)[2:ncol(td)] <- paste(names(td)[2:ncol(td)], city, sep = "_")
  
  overall_el[[city]] <- td
  overall_dd[[city]] <- dd
  cat("\n")
}
```


### Change in deaths due to injury per Billion km across cities

```{r inj_100k_all_cities = "asis"}
require(tibble)

td <- overall_el %>% purrr::reduce(full_join, by = "stage_mode") %>% as.data.frame()
td[is.na(td)] <- 0
td <- td %>% dplyr::select(stage_mode, sort(names(.)))
#readr::write_csv(td, 'results/multi_city/whw_matrices/injury_risks_per_100k_kms.csv')
colnames(td) = gsub("bi_sc", "bisc", colnames(td))
colnames(td) = gsub("bus_sc", "bssc", colnames(td))
colnames(td) = gsub("car_sc", "csc", colnames(td))
colnames(td) = gsub("mc_sc", "mcsc", colnames(td))
colnames(td) = gsub("w_sc", "wsc", colnames(td))

injury_risks_b <- td
# # Convert it into Billion by multiping the values by 10000
# injury_risks_b[, 2:ncol(injury_risks_b)] <- injury_risks_b[, 2:ncol(injury_risks_b)] * 10000

#readr::write_csv(injury_risks_b, 'results/multi_city/whw_matrices/injury_risks_per_billion_kms.csv')
colnames(injury_risks_b) = gsub("bi_sc", "bisc", colnames(injury_risks_b))
colnames(injury_risks_b) = gsub("bus_sc", "bssc", colnames(injury_risks_b))
colnames(injury_risks_b) = gsub("car_sc", "csc", colnames(injury_risks_b))
colnames(injury_risks_b) = gsub("mc_sc", "mcsc", colnames(injury_risks_b))
colnames(injury_risks_b) = gsub("w_sc", "wsc", colnames(injury_risks_b))
injury_risks_lng <- reshape2::melt(injury_risks_b)
col_split <- stringr::str_split(injury_risks_lng$variable, "_", simplify = TRUE, n = 2)
injury_risks_lng <- cbind(injury_risks_lng, col_split)
names(injury_risks_lng)[4] <- 'scenario'
names(injury_risks_lng)[5] <- 'city'
injury_risks_lng <- injury_risks_lng %>% dplyr::select(-variable)
cols <- c("#1b9e77", "#d95f02", "#7570b3", "#e7298a", "#66a61e", "#e6ab02")
injury_risks_lng$scenario <- as.character(injury_risks_lng$scenario)
injury_risks_lng[injury_risks_lng$scenario == 'bl',]$scenario <- 'Baseline'
injury_risks_lng[injury_risks_lng$scenario == 'bisc',]$scenario <- 'Bicycling'
injury_risks_lng[injury_risks_lng$scenario == 'bssc',]$scenario <- 'Public Transport'
injury_risks_lng[injury_risks_lng$scenario == 'csc',]$scenario <- 'Driving'
injury_risks_lng[injury_risks_lng$scenario == 'mcsc',]$scenario <- 'Motorcycling'
injury_risks_lng[injury_risks_lng$scenario == 'wsc',]$scenario <- 'Walking'

rd <- rename(injury_risks_lng, mode = stage_mode)
rd <- rd %>% filter(mode != 'Total')

readr::write_csv(rd, 'results/multi_city/whw_matrices/injury_risks_per_billion_kms_lng.csv')

temp <- injury_risks_lng %>% group_by(stage_mode, scenario) %>% summarise(mv = mean(value), sd = sd(value))
# injury_risks_lng <- left_join(injury_risks_lng, td, by = c('stage_mode', 'scenario'))

cols_cbrewer <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#c51b7d','#ff7f00','#cab2d6','#6a3d9a','#543005','#b15928')

d <- ggplot(injury_risks_lng) +
  aes(x = stage_mode, fill = scenario, weight = value) +
  geom_bar(position = "dodge", colour = 'black') +
  scale_color_manual(values = cols) +
  scale_fill_discrete(guide = guide_legend(reverse=TRUE)) +
  coord_flip() +
  theme_minimal() +
  facet_wrap(vars(city)) +
  labs(x = "Strike Mode", y = "# of incidents", title = paste("Injury risk per billion km")) 
df <- injury_risks_lng %>% spread(value = value, key = city) 
df[df == 0] <- NA
df <- df %>% mutate(mv = rowMeans(dplyr::select(., -c(stage_mode, scenario)))) %>% reshape2::melt()
df1 <- df %>% dplyr::filter(variable == 'mv')
df <- df %>% dplyr::filter(variable != 'mv')
d1 <- ggplot(df, aes(x = stage_mode, fill = variable, weight = value)) +
   geom_bar(position = "dodge", colour = 'black') +
   scale_fill_manual(values = cols_cbrewer) +
   theme_minimal() + 
   geom_point(data = df1, aes(x = stage_mode, fill = variable, y = value), colour = "black") + 
   facet_wrap(vars(scenario)) + 
   #scale_y_continuous(breaks = seq(0, max(df$value, na.rm = T), 25)) +
   labs(x = "Strike Mode", y = "# of incidents", title = paste("Injury risk per billion km")) +
   theme(axis.text.x = element_text(face="bold", color="#993333", size=8, angle=90))
d2 <- ggplot(df %>% filter(!stage_mode %in% c('auto_rickshaw', 'bus', 'truck')), aes(x = stage_mode, fill = variable, weight = value)) +
    geom_bar(position = "dodge", colour = 'black') +
    scale_fill_manual(values = cols_cbrewer) +
    theme_minimal() + 
    geom_point(data = df1 %>% filter(!stage_mode %in% c('auto_rickshaw', 'bus', 'truck')) , aes(x = stage_mode, fill = variable, y = value), colour = "black") + 
    facet_grid(vars(scenario), scales='free_x') + 
    # scale_y_continuous(breaks = seq(0, max(df$value, na.rm = T), 50)) +
    labs(x = "Strike Mode", y = "# of incidents", title = paste("Injury risk per billion km")) +
    theme(axis.text.x = element_text(face="bold", color="#993333", size=14))
# ggsave(d, file = paste0("results/multi_city/whw_matrices/plots/injury_risk_billion_km.png"), limitsize = F, dpi = 300, scale = 1.5)
fp <- paste0("results/multi_city/whw_matrices/interactive_plots/injury_risk_billion_km.html")
# htmlwidgets::saveWidget(plotly::ggplotly(d), file.path(normalizePath(dirname(fp)), basename(fp)))
# ggsave(d1, file = paste0("results/multi_city/whw_matrices/plots/injury_risk_billion_km_by_mode.png"),  width = 8, height = 5, dpi = 300, units = "in", scale = 2, device='png')
 fp <- paste0("results/multi_city/whw_matrices/interactive_plots/injury_risk_billion_km_by_mode.html")
 # htmlwidgets::saveWidget(plotly::ggplotly(d1), file.path(normalizePath(dirname(fp)), basename(fp)))
# print(kable(format(injury_risks_b, scientific = F), caption = paste('Injury per billion km across all cities')))
# ggsave(d2, file = paste0("results/multi_city/whw_matrices/plots/injury_risk_billion_km_by_specific_mode.png"), limitsize = F, dpi = 300, scale = 1.5)
fp <- paste0("results/multi_city/whw_matrices/interactive_plots/injury_risk_billion_km_by_specific_mode.html")
# htmlwidgets::saveWidget(plotly::ggplotly(d2), file.path(normalizePath(dirname(fp)), basename(fp)))
```

### # of deaths due to injury per 100,000 people

```{r inj_100k = "asis"}
require(tibble)
overall_el_normalized <- list()
for (city in cities){
  
  print(city)
  
  el <- list()
  
  for (cs in names(io[[city]]$outcomes$whw)){
    if (length(names(io[[city]]$outcomes$whw$Baseline)) == 2){
      td1 <- round(io[[city]]$outcomes$whw[[cs]]$nov) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      td2 <- colSums(round(io[[city]]$outcomes$whw[[cs]]$whw)) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      td3 <- full_join(td2, td1, by = 'mode') %>% mutate(count = rowSums(.[2:3], na.rm = T)) %>% dplyr::select(-c('count.x', 'count.y'))
      
    }else if(length(names(io[[city]]$outcomes$whw$Baseline)) == 1 && names(io[[city]]$outcomes$whw$Baseline) == 'whw'){
      td3 <- colSums(round(io[[city]]$outcomes$whw[[cs]]$whw)) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
    }else if(length(names(io[[city]]$outcomes$whw$Baseline)) == 1 && names(io[[city]]$outcomes$whw$Baseline) == 'nov'){
      
      td3 <- round(io[[city]]$outcomes$whw[[cs]]$nov) %>% as.data.frame() %>% rownames_to_column() %>% rename(mode = rowname) %>% rename_at(2, ~"count")
      
    }
    
    td4 <- io[[city]]$true_dist %>% filter(stage_mode %in% td3$mode) %>% dplyr::select(stage_mode)
    
    if (length(el) == 0){
      el <- td4 %>% dplyr::select(stage_mode)
    }
    
    td4 <- full_join(td4, td3 %>% dplyr::select(mode, count) %>% rename(stage_mode = mode), by = 'stage_mode')
    
    var <- scen_lt %>% filter(scen_name == cs) %>% dplyr::select(val) %>% as.character()
    
    names(td4)[2] <- var
    
    td5 <- td4
    
    td4[, 2] <- round(td4[, 2] / sum(io[[city]]$demographic$population) * 100000, 2)
    
    el <- inner_join(el, td4, by = 'stage_mode')
  }
  
  el <- el %>% ungroup() %>% janitor::adorn_totals(c('row', 'col'))
  
  print(kable(el, caption = city))
  td <- el
  names(td)[2:ncol(td)] <- paste(names(td)[2:ncol(td)], city, sep = "_")
  
  overall_el_normalized[[city]] <- td
  
  cat("\n")
  
  
  
}

td <- overall_el_normalized %>% purrr::reduce(full_join, by = "stage_mode") %>% as.data.frame()
td[is.na(td)] <- 0
td <- td %>% dplyr::select(stage_mode, sort(names(.)))

td <- rename(td, mode = stage_mode)
colnames(td) = gsub("bi_sc", "bisc", colnames(td))
colnames(td) = gsub("bus_sc", "bssc", colnames(td))
colnames(td) = gsub("car_sc", "csc", colnames(td))
colnames(td) = gsub("mc_sc", "mcsc", colnames(td))
colnames(td) = gsub("w_sc", "wsc", colnames(td))
injury_risks_per_100k <- reshape2::melt(td)
col_split <- stringr::str_split(injury_risks_per_100k$variable, "_", simplify = TRUE, n = 2)
injury_risks_per_100k <- cbind(injury_risks_per_100k, col_split)
names(injury_risks_per_100k)[4] <- 'scenario'
names(injury_risks_per_100k)[5] <- 'city'
injury_risks_per_100k <- injury_risks_per_100k %>% dplyr::select(-variable)
injury_risks_per_100k$scenario <- as.character(injury_risks_per_100k$scenario)
injury_risks_per_100k[injury_risks_per_100k$scenario == 'bl',]$scenario <- 'Baseline'
injury_risks_per_100k[injury_risks_per_100k$scenario == 'bisc',]$scenario <- 'Bicycling'
injury_risks_per_100k[injury_risks_per_100k$scenario == 'bssc',]$scenario <- 'Public Transport'
injury_risks_per_100k[injury_risks_per_100k$scenario == 'csc',]$scenario <- 'Driving'
injury_risks_per_100k[injury_risks_per_100k$scenario == 'mcsc',]$scenario <- 'Motorcycling'
injury_risks_per_100k[injury_risks_per_100k$scenario == 'wsc',]$scenario <- 'Walking'

readr::write_csv(injury_risks_per_100k %>% filter(mode != 'Total' & scenario != 'Total'), 'results/multi_city/whw_matrices/injury_risks_per_100k_pop.csv')


```