summary_tables_PA_AP.Rmd

---
title: "AP/PA Tables ITHIM Global"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document: default
  word_document: default
  pdf_document: default
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(comment=NA, prompt=FALSE, cache=FALSE, echo=F, results='asis')

```


# Introduction
These are the summary tables of the following items

1) Individual-level PM2.5 concentrations for baseline and scenarios
2) Overall CO2 emissions for baseline and scenarios
3) Individual-level physical activity for baseline and scenarios



```{r, echo=FALSE}
#setwd('C:/Users/rg574/Dropbox/spatiotemporal analysis fatalities inida/Rajasthan tourism road deaths') #to create pretty tables
```

```{r, message=FALSE, warning=FALSE}
#library(INLA)     #loading the INLA package
library(ggplot2)  #loading ggplot package for plotting graphs
library(knitr)     
library(tidyr)
library(dplyr)

```


# Boxplots of Individual-level PM2.5 Concentrations

```{r, message=FALSE, warning=FALSE, echo=FALSE}
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']


all_inputs <- read.csv('all_city_parameter_inputs.csv',stringsAsFactors = F)

#all_inputs$cape_town <- all_inputs$accra
#all_inputs$vizag <- all_inputs$sao_paulo

parameter_names <- all_inputs$parameter
parameter_starts <- which(parameter_names!='')
parameter_stops <- c(parameter_starts[-1] - 1, nrow(all_inputs))
parameter_names <- parameter_names[parameter_names!='']
parameter_list <- list()
compute_mode <- 'constant'
for(i in 1:length(parameter_names)){
  parameter_list[[parameter_names[i]]] <- list()
  parameter_index <- which(all_inputs$parameter==parameter_names[i])
  if(all_inputs[parameter_index,2]=='')  {
    parameter_list[[parameter_names[i]]] <- lapply(cities,function(x) {
      city_index <- which(colnames(all_inputs)==x)
      val <- all_inputs[parameter_index,city_index]
      ifelse(val%in%c('T','F'),val,as.numeric(val))
    })
    names(parameter_list[[parameter_names[i]]]) <- cities
  }else if(all_inputs[parameter_index,2]=='constant'){
    indices <- 0
    if(compute_mode=='sample') indices <- 1:2
    parameter_list[[parameter_names[i]]] <- lapply(cities,function(x) {
      city_index <- which(colnames(all_inputs)==x)
      val <- all_inputs[parameter_index+indices,city_index]
      ifelse(val=='',0,as.numeric(val))
    })
    names(parameter_list[[parameter_names[i]]]) <- cities
  }else{
    parameter_list[[parameter_names[i]]] <- lapply(cities,function(x) {
      city_index <- which(colnames(all_inputs)==x)
      if(any(all_inputs[parameter_starts[i]:parameter_stops[i],city_index]!='')){
        sublist_indices <- which(all_inputs[parameter_starts[i]:parameter_stops[i],city_index]!='')
        thing <- as.list(as.numeric(c(all_inputs[parameter_starts[i]:parameter_stops[i],city_index])[sublist_indices]))
        names(thing) <- c(all_inputs[parameter_starts[i]:parameter_stops[i],2])[sublist_indices]
        thing
      }
    }
    )
    names(parameter_list[[parameter_names[i]]]) <- cities
  }
}

for(i in 1:length(parameter_list)) assign(names(parameter_list)[i],parameter_list[[i]])

# high<- c(70, 25, 150,60)
# low<- c(40,15,100,40)
for (x in 1: length(cities))
{
  # print(cities[x])
  names(io[[cities[x]]]$outcomes$pm_conc_pp)[6:11]<-c("Baseline","Walking", "Bicycling", "Driving", "Motorcycling", "Public Transport")
  data_long <- gather(io[[cities[x]]]$outcomes$pm_conc_pp, scenario, pm_conc, Baseline:`Public Transport`, factor_key=TRUE)
  y<-ggplot(data_long, aes(x=scenario, y=pm_conc, fill=scenario)) + geom_boxplot(outlier.shape = NA)+ ggtitle(cities[x])+ labs(y="Daily PM2.5 Concentration", x = "Scenarios")
  #+ylim(low[x],high[x])
  print(y)
  
}
```




# Descriptive tables of Individual-level PM2.5 Concentrations

```{r, message=FALSE, warning=FALSE, echo=FALSE}

for (x in 1: length(cities))
{
  names(io[[cities[x]]]$outcomes$pm_conc_pp)[6:11]<-c("Baseline","Walking", "Bicycling", "Driving", "Motorcycling", "Public Transport")
  data_long <- gather(io[[cities[x]]]$outcomes$pm_conc_pp, scenario, pm_conc, Baseline:`Public Transport`, factor_key=TRUE)
  summary<- as.data.frame(data_long %>% group_by(scenario) %>% summarise('mean'=mean(pm_conc),'5th'=quantile(pm_conc, 0.05),'50th'=quantile(pm_conc, 0.5),'95th'=quantile(pm_conc, 0.9)))
  summary<- cbind(summary$scenario ,round(summary[,2:5], digits=1))
  summary$change_PM<- round(io[[cities[x]]]$outcomes$scenario_pm - io[[cities[x]]]$outcomes$scenario_pm[1], digits=2)
  names(summary)[1]<-"Scenario"
  print(kable(summary, caption= cities[x]))
  
}
```




<!-- # Descriptive tables of emission inventory -->

<!-- ```{r, message=FALSE, warning=FALSE, echo=FALSE} -->

<!-- trans_share<-c(22.5, 40.0, 22.5, 28.1)  -->
<!-- pm_conc<- c(50, 18, 122, 47 ) -->
<!-- city<-cities -->
<!-- shares<- as.data.frame(cbind(city,as.numeric(trans_share), pm_conc)) -->
<!-- names(shares)[2]<- "trans_share" -->
<!-- for (x in 1:length(cities)) -->
<!-- { -->
<!--   modes<-names(unlist(io[[cities[x]]]$emission_inventory)) -->
<!--   emissions<-as.data.frame(unlist(io[[cities[x]]]$emission_inventory)) -->
<!--   city_emissions<-cbind(as.data.frame(modes),  as.data.frame(emissions$`unlist(io[[cities[x]]]$emission_inventory)`)) -->
<!--   names(city_emissions)[2]<- "emissions" -->
<!--   select<- c("car", "motorcycle", "bus_driver", "truck", "big_truck") -->
<!--   city_emissions$modes<- as.character(city_emissions$modes) -->
<!--   city_emissions$modes[!(city_emissions$modes %in% select)]<- "other" -->
<!--   summary<-city_emissions %>% group_by(modes) %>% summarise(sum(emissions)) -->
<!--   names(summary)[2]<- "emissions" -->
<!--   summary$emissions <- round(summary$emissions*100/sum(summary$emissions), digits=1) -->
<!--   summary<-as.data.frame(summary) -->
<!--   summary[nrow(summary)+1,1]<-"Transport share" -->
<!--   summary[nrow(summary),2]<-as.character(shares$trans_share[x]) -->
<!--   summary[nrow(summary)+1,1]<-"PM2.5 Conc" -->
<!--   summary[nrow(summary),2]<-as.character(shares$pm_conc[x]) -->
<!--   print(kable(summary, caption= cities[x])) -->

<!-- } -->

<!-- ``` -->


# Descriptive tables of PM 2.5 emission inventory

```{r, message=FALSE, warning=FALSE, echo=FALSE}

#trans_share<-c(22.5, 40.0, 22.5, 28.1) 
#pm_conc<- c(50, 18, 122, 47 )
city<-cities
#shares<- as.data.frame(cbind(city,as.numeric(io), pm_conc))
#names(shares)[2]<- "trans_share"

sl <- list()
for (x in 1:length(cities))
{
  modes<-names(unlist(io[[cities[x]]]$PM_emission_inventory))
  emissions<-as.data.frame(unlist(io[[cities[x]]]$PM_emission_inventory))
  city_emissions<-cbind(as.data.frame(modes),  as.data.frame(emissions$`unlist(io[[cities[x]]]$PM_emission_inventory)`))
  names(city_emissions)[2]<- "emissions"
  select<- c("car", "motorcycle", "bus_driver", "truck", "big_truck")
  city_emissions$modes<- as.character(city_emissions$modes)
  city_emissions$modes[!(city_emissions$modes %in% select)]<- "other"
  summary<-city_emissions %>% group_by(modes) %>% summarise(sum(emissions))
  names(summary)[2]<- cities[x]
  summary[[cities[x]]] <- round(summary[[cities[x]]]*100/sum(summary[[cities[x]]]), digits=1)
  summary<-as.data.frame(summary)
  summary[nrow(summary)+1,1] <- "Transport share"
  summary[nrow(summary),2]<- as.character(round(pm_trans_share[[cities[x]]] * 100))
  summary[nrow(summary)+1,1] <- "PM2.5 Conc"
  summary[nrow(summary),2] <- as.character(pm_conc_base[[cities[x]]])
  io[[cities[x]]]$summary_emission<- summary
  
  if (length(sl) == 0){
    sl <- summary
  }else{
    sl <- left_join(sl , summary)
  }
  
}

# summary_all<- cbind(io[[cities[1]]]$summary_emission,io[[cities[2]]]$summary_emission,io[[cities[3]]]$summary_emission,io[[cities[4]]]$summary_emission)
# summary_all<- summary_all[-c(3,5,7)]
# summary_all$modes[2]<- "bus"
# 
# 
print(kable(sl))
```

# Descriptive table for CO2 Emissions
```{r, message=FALSE, warning=FALSE, echo=FALSE}

cl <- list()
for (city in cities)
{
  
  td <- round(colSums(io[[city]]$outcomes$co2_conc, na.rm = T), 1) %>% as.data.frame() %>% tibble::rownames_to_column()
  names(td) <- c('Scenario', city)
  
  if (length(cl) == 0)
    cl <- td
  else
    cl <- left_join(cl, td)
  
}

cl$Scenario <- c("Baseline","Walking", "Bicycling", "Driving", "Motorcycling", "Public Transport")

print(kable(cl))

```



# Boxplots of Individual-level Physical activity (MMETs)

```{r, message=FALSE, warning=FALSE, echo=FALSE}
limit=100
for (x in 1: length(cities))
{
  names(io[[cities[x]]]$outcomes$mmets)[5:10]<-c("Baseline","Walking", "Bicycling", "Driving", "Motorcycling", "Public Transport")
  data_long <- gather(io[[cities[x]]]$outcomes$mmets, scenario, mmet, Baseline:`Public Transport`, factor_key=TRUE)
  y<-ggplot(data_long, aes(x=scenario, y=mmet, fill=scenario)) + geom_boxplot(outlier.shape = NA)+ ggtitle(cities[x])+ labs(y="Weekly MMET", x = "Scenarios")+  ylim(0, limit)
  print(y)
  
}
```

<!-- # Descriptive tables of Individual-level Physical activity (MMETs) -->

<!-- ```{r, message=FALSE, warning=FALSE, echo=FALSE} -->
<!-- for (x in 1: length(cities)) -->
<!-- { -->
<!--   names(io[[cities[x]]]$outcomes$mmets)[5:10]<-c("baseline","walk_scen", "bike_scen", "car_scen", "MC_scen", "bus_scen") -->
<!--   data_long <- gather(io[[cities[x]]]$outcomes$mmets, scenario, mmet, baseline:bus_scen, factor_key=TRUE) -->
<!--   summary<- as.data.frame(data_long %>% group_by(scenario) %>% summarise('mean'=mean(mmet),'5th'=quantile(mmet, 0.05),'50th'=quantile(mmet, 0.5),'95th'=quantile(mmet, 0.9))) -->
<!--   summary<- cbind(summary$scenario ,round(summary[,2:5], digits=1)) -->
<!--   names(summary)[1]<-"Scenario" -->
<!--   print(kable(summary, caption= cities[x])) -->



<!-- } -->
<!-- ``` -->

<!-- # Descriptive tables of injury outcomes -->

<!-- ```{r, message=FALSE, warning=FALSE, echo=FALSE} -->
<!-- for (x in 1: length(cities)) -->
<!-- { -->

<!--   summary<- as.data.frame(io[[cities[x]]]$outcomes$injuries %>% group_by(scenario) %>% summarise('Ped'=sum(pedestrian,na.rm=T),'Bike'= sum(cycle), 'Car'=sum(car), "Motorcycle"=sum(motorcycle,na.rm=T), "All modes"=sum(Deaths,na.rm=T))) -->
<!--   summary$scenario<- c("baseline","walk_scen", "bike_scen", "car_scen", "MC_scen", "bus_scen") -->
<!--   summary<- cbind(summary$scenario ,round(summary[,2:6], digits=0)) -->
<!--   names(summary)[1]<-"Scenario" -->
<!--   print(kable(summary, caption= cities[x])) -->

<!-- } -->
<!-- ``` -->

<!-- # Striking vehicles in scenarios -->

<!-- ```{r, message=FALSE, warning=FALSE, echo=FALSE} -->
<!-- for (x in 1: length(cities)) -->
<!-- { -->

<!-- scenarios<- c('Baseline', 'Scenario 1', 'Scenario 2', 'Scenario 3', 'Scenario 4', 'Scenario 5') -->
<!-- for (i in 1: length(scenarios)) -->
<!-- { -->
<!-- whw<-as.data.frame(rowSums(io[[cities[x]]]$outcomes$whw[[scenarios[i]]]$whw)) -->
<!-- noov<-as.data.frame(sum(io[[cities[x]]]$outcomes$whw[[scenarios[i]]]$nov)) -->
<!-- names(whw)[1]<- as.character(scenarios[i]) -->
<!-- names(noov)[1]<-as.character(scenarios[i]) -->

<!-- if (i ==1 ) -->
<!-- { -->
<!--  summary<-round(rbind(whw,nov=noov), digits = 1) -->
<!-- }  -->
<!-- else -->
<!-- { -->
<!-- summary<- cbind(summary, round(rbind(whw, nov=noov), digits=1) )  -->
<!-- } -->

<!-- } -->

<!-- names(summary)<-c("base", "walk", "bike", "car", "MC", "bus") -->
<!-- print(kable(summary, caption= cities[x])) -->

<!-- } -->

<!-- ``` -->



<!-- # Distance share by scenarios in percentages (excludes city-specific minority modes-- auto rickshaws, subway, etc.) -->

<!-- ```{r, message=FALSE, warning=FALSE, echo=FALSE} -->
<!-- for (x in 1: length(cities)) -->
<!-- { -->

<!--   select<- c("car", "motorcycle", "walking", "cycle", "bus") -->

<!--   io[[cities[x]]]$trip_scen_sets$trip_mode[!io[[cities[x]]]$trip_scen_sets$trip_mode %in% select]<- "other" -->

<!--   dist_scen<- io[[cities[x]]]$trip_scen_sets %>% group_by(trip_mode, scenario) %>% summarise(sum(trip_distance)) -->
<!--   dist_scen<- as.data.frame (dist_scen) -->
<!--   names(dist_scen)[3]<-"distance" -->
<!--   dist_scen$distance<- round(dist_scen$distance, digits=0) -->
<!--     dist_scen<-spread(dist_scen, trip_mode, distance) -->
<!--   dist_scen$sum<-rowSums (dist_scen[2:(ncol(dist_scen))], na.rm = FALSE, dims = 1) -->
<!--   dist_scen_prop<- round(dist_scen[,2:(ncol(dist_scen))]*100/dist_scen$sum, digits=1) -->
<!--   dist_scen_prop<- cbind(dist_scen[,1], dist_scen_prop) -->

<!--   names(dist_scen_prop)[1]<-"Scenario" -->
<!--   dist_scen_prop$Scenario<- c("baseline","walk_scen", "bike_scen", "car_scen", "MC_scen", "bus_scen") -->
<!--   print(kable(dist_scen_prop, caption= cities[x])) -->
<!--   print(kable(dist_scen, caption= cities[x])) -->

<!-- } -->
<!-- ``` -->