summary_tables_PA_AP.Rmd

---
title: "AP/PA Tables ITHIM Global"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document: 
    toc: true
    toc_depth: 5
    toc_float: true
  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(tidyverse)
library(readxl)
#library(plotly)
#library(htmltools)

```


# Boxplots of Individual-level PM2.5 Concentrations

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

scen_names <- c("Baseline", "Bicycling", "Driving", "Public Transport", "Motorcycling")

input_parameter_file <<- "InputParameters_v18.0.xlsx"
# scenario definition
scenario_name <- "GLOBAL"
reference_scenario <- 'Baseline'

compute_mode <- 'constant' # constant parameters from the given parameters


all_inputs <- read_excel(input_parameter_file, sheet = "all_city_parameter_inputs")
all_inputs[is.na(all_inputs)] <- ""
all_inputs <- as.data.frame(all_inputs)

# get input parameters into correct format
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()

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,ifelse(is.numeric(val), as.numeric(val), as.character(val)))
    })
    names(parameter_list[[parameter_names[i]]]) <- cities
  }else if(all_inputs[parameter_index,2]=='constant'){
    if (compute_mode != 'sample'){
      indices <- 0
      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))
      })
    }
    if(compute_mode=='sample'){ # if sampling from distribution, check that distribution parameters exist
      parameter_list[[parameter_names[i]]] <- lapply(cities,function(x) {
        indices <- 1:2
        city_index <- which(colnames(all_inputs)==x)  
        val <- all_inputs[parameter_index+indices,city_index] 
        if (val[1] == '' & val[2]==''){  # if no distribution parameters given in input file, read in constant value instead
          indices <-0
          city_index <- which(colnames(all_inputs)==x) 
          val <- all_inputs[parameter_index+indices,city_index]} 
        val <- 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
  }
}

list2env(parameter_list, environment()) 

# read in global parameters

all_global_inputs <- read_excel(input_parameter_file, sheet = "all_global_parameter_inputs")
all_global_inputs[is.na(all_global_inputs)] <- ""
all_global_inputs <- as.data.frame(all_global_inputs)

# get input parameters into correct format
global_parameter_names <- all_global_inputs$parameter
global_parameter_starts <- which(global_parameter_names!='')
global_parameter_stops <- c(global_parameter_starts[-1] - 1, nrow(all_global_inputs)) 
global_parameter_names <- global_parameter_names[global_parameter_names!='']
global_parameter_list <- list()

for(i in 1:length(global_parameter_names)){
  global_parameter_list[[global_parameter_names[i]]] <- list()
  global_parameter_index <- which(all_global_inputs$parameter==global_parameter_names[i]) 
  if(all_global_inputs[global_parameter_index,2]=='')  { 
    
    global_parameter_list[[global_parameter_names[i]]] <- all_global_inputs[global_parameter_index,'global']
    
  }else if(all_global_inputs[global_parameter_index,2]=='constant'){
    if (compute_mode != 'sample'){
      global_parameter_list[[global_parameter_names[i]]] <- ifelse(all_global_inputs[global_parameter_index,'global']=='',
                                                                   0,as.numeric(all_global_inputs[global_parameter_index,'global']))
    }
    else if(compute_mode=='sample'){ # if sampling from distribution, check that distribution parameters exist
      indices <- 1:2
      val <- all_global_inputs[global_parameter_index+indices,'global'] 
      if (val[1] == '' & val[2]==''){  # if no distribution parameters given in input file, read in constant value instead
        val <- all_global_inputs[global_parameter_index,'global']} 
      val <- as.numeric(val)
      global_parameter_list[[global_parameter_names[i]]] <- val
    }
  }
}

list2env(global_parameter_list, environment()) 

dist_cat <- unlist(strsplit(gsub(" ", "", dist_cat, fixed = TRUE), "\\,"))

outcome_age_min <- as.numeric(unlist(strsplit(gsub(" ", "", outcome_age_min, fixed = TRUE), "\\,")))
outcome_age_max <- as.numeric(unlist(strsplit(gsub(" ", "", outcome_age_max, fixed = TRUE), "\\,")))
outcome_age_groups <- unlist(strsplit(gsub(" ", "", outcome_age_groups, fixed = TRUE), "\\,"))

min_age <- as.numeric(min_age)
max_age <- as.numeric(max_age)

day_to_week_scalar <- as.numeric(day_to_week_scalar)

```

```{r}
# high<- c(70, 25, 150,60)
# low<- c(40,15,100,40)
# l <- htmltools::tagList()
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")
                                                         
  n_col <- ncol(io[[cities[x]]]$outcomes$pm_conc_pp)
  names(io[[cities[x]]]$outcomes$pm_conc_pp)[(n_col - 4):n_col] <- scen_names
  
  data_long <- gather(io[[cities[x]]]$outcomes$pm_conc_pp, scenario, pm_conc,
                      Baseline:`Motorcycling`, factor_key = TRUE)
  y <- ggplot(data_long, aes(x = scenario, y = pm_conc, fill = scenario)) +
    geom_boxplot(outlier.shape = 8) + ggtitle(cities[x]) +
    labs(y = "Daily PM2.5 Concentration", x = "Scenarios") #+
  # +ylim(low[x],high[x])
  # yy <- ggplot() +
  #   geom_density(data_long, aes(x = scenario, y = pm_conc, fill = scenario)) +
  #   ggtitle(cities[x]) +
  #   labs(y = "Daily PM2.5 Concentration", x = "Scenarios")
  print(y)
  #print(yy)
  # l[[i]] <- ggplotly(ggplot(data_long, aes(x = scenario, y = pm_conc, fill = scenario)) +
  #   geom_boxplot(outlier.shape = 8) + ggtitle(cities[x]) +
  #   labs(y = "Daily PM2.5 Concentration", x = "Scenarios"))
  # print(ggplotly(yy))
}
```




# Descriptive tables of Individual-level PM2.5 Concentrations

```{r, message=FALSE, warning=FALSE, echo=FALSE}
data_long <- NA
summary <- NA
for (x in 1:length(cities)) {
  #names(io[[cities[x]]]$outcomes$pm_conc_pp)[6:11]<-c("Baseline","Walking", "Bicycling", "Driving", "Motorcycling", "Public Transport")
  #
  n_col <- ncol(io[[cities[x]]]$outcomes$pm_conc_pp)
  names(io[[cities[x]]]$outcomes$pm_conc_pp)[(n_col - 4):n_col] <- scen_names
  data_long <- gather(io[[cities[x]]]$outcomes$pm_conc_pp, scenario, pm_conc,
                      Baseline:`Motorcycling`, 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$scenario_PM <- round(io[[cities[x]]]$outcomes$scenario_pm, digits=2)
  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))
{
  print(cities[x])
  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}

co2_conc <- list()
td <- NA
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(co2_conc) == 0)
    co2_conc <- td
  else
    co2_conc <- left_join(co2_conc, td)
  
}

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

print(kable(co2_conc))

```



# 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")
  #
  n_col <- ncol(io[[cities[x]]]$outcomes$mmets)
  names(io[[cities[x]]]$outcomes$mmets)[(n_col - 4):n_col] <- scen_names
  
  data_long <- gather(io[[cities[x]]]$outcomes$mmets, scenario, mmet, Baseline:`Motorcycling`, 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")
  n_col <- ncol(io[[cities[x]]]$outcomes$mmets)
  names(io[[cities[x]]]$outcomes$mmets)[(n_col - 4):n_col] <- scen_names
  
  data_long <- gather(io[[cities[x]]]$outcomes$mmets, scenario, mmet, Baseline:Motorcycling, 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 table for Physical Activity (MMETs)
```{r, message=FALSE, warning=FALSE, echo=FALSE}

cl <- list()
for (city in cities)
{
  #td <- round(colMeans(io[[city]]$outcomes$mmets[,5:10], na.rm = T), 1) %>%
  #  as.data.frame() %>% tibble::rownames_to_column()
  #  
  n_col <- ncol(io[[cities[x]]]$outcomes$mmets)

  td <- round(colMeans(io[[city]]$outcomes$mmets[,(n_col - 4):n_col], 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 <- scen_names

print(kable(cl))

```

<!-- # 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])) -->

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



# **2nd option**

# PM 2.5 summary and CO2
```{r message=FALSE, warning=FALSE, echo=FALSE}

l <- list()
# I have to keep the first chunk of code to run this part
# high<- c(70, 25, 150,60)
# low<- c(40,15,100,40)
# l <- htmltools::tagList()
for (x in 1:length(cities)) {
  # print(cities[x])
  cat("\n \n")
  cat( paste0('## **', cities[x], '** \n' )) 
  # Emission inventory
  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
  
   
  # PM 2.5 table and plot
  n_col <- ncol(io[[cities[x]]]$outcomes$mmets)
  
  names(io[[cities[x]]]$outcomes$pm_conc_pp)[(n_col - 4):n_col] <- scen_names
  
  data_long <- gather(io[[cities[x]]]$outcomes$pm_conc_pp, scenario, pm_conc,
                      Baseline:`Motorcycling`, factor_key = TRUE)
  
  summary2 <- 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)))
  summary2 <- cbind(summary2$scenario, round(summary2[,2:5], digits = 1))
  summary2$scenario_PM <- round(io[[cities[x]]]$outcomes$scenario_pm, 2)
  summary2$change_PM <- round(io[[cities[x]]]$outcomes$scenario_pm - io[[cities[x]]]$outcomes$scenario_pm[1], digits = 2)
  names(summary2)[1] <- "Scenario"
  
  # Plot
  y <- ggplot(data_long, aes(x = scenario, y = pm_conc, fill = scenario)) +
    geom_boxplot(outlier.shape = 8) + ggtitle(cities[x]) +
    labs(y = "Daily PM2.5 Concentration", x = "Scenarios")
  
  # CO2
  td <- round(colSums(io[[cities[x]]]$outcomes$co2_conc, na.rm = T), 1) %>%
    as.data.frame() %>% tibble::rownames_to_column()
  names(td) <- c('Scenario', 'CO2')
  td$Scenario <- scen_names
  
  cat("\n \n")
  cat( paste('#### Boxplot of Individual-level PM2.5 concentrations \n' ))
  print(y)
  
  cat("\n \n")
  cat( paste('#### Emission inventory \n' )) 
  print(kable(summary))
  
  cat("\n \n")
  cat( paste('#### Descriptive table of Individual-level PM2.5 concentrations and CO2 emissions\n' ))
  l[[cities[x]]] <- left_join(summary2, td, by = "Scenario")
  l[[cities[x]]]$city <- cities[x]
  print(kable(cbind(summary2, td), caption = cities[x]))
  
}

output_file <- "results/multi_city/AP/desc_stats.xlsx"
writexl::write_xlsx(list(PM_2.5_emission_inventory = sl,
                         CO2_conc = co2_conc,
                         summary_stats_PM2.5_CO2 = bind_rows(l)), path = output_file)


write_csv(bind_rows(l), "results/multi_city/AP/desc_stats.csv")


```

# 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")
  n_col <- ncol(io[[cities[x]]]$outcomes$mmets)
  names(io[[cities[x]]]$outcomes$mmets)[(n_col - 4):n_col] <- scen_names
  data_long <- gather(io[[cities[x]]]$outcomes$mmets, scenario, mmet, Baseline:`Motorcycling`, 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)
  
}
```