README.Rmd

---
output:
  github_document:
    pandoc_args: --webtex
always_allow_html: true
bibliography: Accessibility-Vaccination-Sites-Hamilton/references.bib
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# An examination of the accessibility implications of a pilot COVID-19 vaccination program in Hamilton, Ontario

<!-- badges: start -->
<!-- badges: end -->

Antonio Páez (McMaster University)  
Christopher D. Higgins (University of Toronto Scarborough)  

Findings (2020) https://doi.org/10.32866/001c.24082

## Abstract

The Government of Ontario in Canada announced the pilot for a new vaccination program, with designated pharmacies across the province now able to offer COVID-19 vaccines. The accessibility of this program raises questions about travel times to vaccination sites and the distribution of these times among the population. In our examination of the City of Hamilton we find that selected sites do not serve rural and urban residents well; particularly, the associated cost of travel (in terms of travel time) is expected to be disproportionally borne by lower income urban populations and rural residents. Modest additions to the list of pilot sites in the city can substantially alleviate this inequity.

## Keywords

- COVID-19  
- Vaccination sites
- Accessibility  
- Hamilton  
- Ontario
- Reproducible research  



```{r install-data-package, include = FALSE}
# Run only once if needed to install data package `vaccHamilton`
if (!require("vaccHamilton", character.only = TRUE)) {
      devtools::install_github("https://github.com/paezha/Accessibility-Pharmacies-Hamilton-Vaccines", 
                               subdir = "vaccHamilton")
  }
```

```{r load-packages, include=FALSE, cache=FALSE}
#Load the packages needed to read and work with the data
library(grid)
library(gridExtra)
library(kableExtra)
library(patchwork)
library(png)
library(sf)
library(tidyverse)
library(units)
library(vaccHamilton) # Data package
```

```{r invoke-data, include=FALSE}
data("data_da_2016")
data("dwelling_network_points_2016")
data("hamilton_cma")
data("modes_55to69")
data("pharmacy_locations")
data("ttm_car")
data("ttm_transit")
data("ttm_walk")
data("urban_types")
```

```{r population-age-55-plus, include=FALSE}
# Data preparation
## Calculate number of people age 55+ per DA:
data_da_2016 <- data_da_2016 %>%
  filter(`Region Name` == "Hamilton") %>%
  mutate(Population_55to69 = (`v_CA16_208: 55 to 59 years` + 
                  `v_CA16_226: 60 to 64 years` + 
                  `v_CA16_247: 65 to 69 years`),
         Population_70plus = (`v_CA16_265: 70 to 74 years` +
                  `v_CA16_283: 75 to 79 years` +
                  `v_CA16_301: 80 to 84 years` + 
                  `v_CA16_319: 85 years and over`))
```

```{r select-pharmacies, include=FALSE}
# Data preparation
## Select pharmacy locations:
pharmacy_locations <- pharmacy_locations %>%
  slice(1:23)

## Also filter selected pharmacies from travel time tables:
ttm_car <- ttm_car %>%
  dplyr::filter(OBJECTID <= 23)

ttm_transit <- ttm_transit %>%
  dplyr::filter(OBJECTID <= 23)

ttm_walk <- ttm_walk %>%
  dplyr::filter(OBJECTID <= 23)
```

```{r clip-da-geometry, include=FALSE}
# Data preparation
## This is purely cosmetic: clip `data_da_2016`:
data_da_2016 <- data_da_2016 %>%
  st_intersection(hamilton_cma)
```

```{r number-of-residential-units-per-DA, include=FALSE}
# Data preparation
## Find number of residential units in each DA. This is needed to assign the population to parcels
residential_units <- dwelling_network_points_2016 %>%
  st_drop_geometry() %>%
  group_by(DAUID) %>%
  summarize(residential_units = sum(Dwellings))
```

```{r join-population-to-parcels, include=FALSE}
# Data preparation
## Join number of people age 55+ per DA to dwellings and total number of residential units:
dwelling_network_points_2016 <- dwelling_network_points_2016 %>%
  left_join(data_da_2016 %>%
              st_drop_geometry() %>%
              dplyr::transmute(DAUID = GeoUID, 
                               Population_55to69,
                               Population_70plus),
            by = "DAUID") %>%
  left_join(residential_units,
            by = "DAUID")
```

```{r calculate-population-per-parcel, include=FALSE}
# Data preparation
##Calculate population per point on the network. This is the population in the DA divided by the total number of residential units. This gives the average number of people per residential unit. Then, multiply by the number of dwellings (residential units) at the point:
dwelling_network_points_2016 <- dwelling_network_points_2016 %>%
  mutate(Population_55to69 = (Population_55to69/residential_units) * Dwellings,
         Population_70plus = (Population_70plus/residential_units) * Dwellings)
```

```{r join-region-to-parcels, include=FALSE}
dwelling_network_points_2016 <- dwelling_network_points_2016 %>%
  st_join(urban_types)
```

```{r calculate-population-statistics, include=FALSE}
population_statistics <- dwelling_network_points_2016 %>% 
  st_drop_geometry() %>% 
  dplyr::filter(!is.na(Type)) %>% 
  group_by(Type) %>% 
  summarize(Population_55to69 = sum(Population_55to69, 
                                    na.rm = TRUE), 
            Population_70plus = sum(Population_70plus, 
                                    na.rm = TRUE))
```

```{r join-population-data-to-travel-time-tables, include=FALSE}
# Data preparation
## Join the population data to the travel time matrices:
ttm_car <- ttm_car %>%
  left_join(dwelling_network_points_2016 %>%
              st_drop_geometry() %>%
              dplyr::transmute(ID, 
                            Population_55to69,
                            Population_70plus,
                            Urban_type = Type,
                            DAUID, 
                            TAZUID),
            by = c("UID" = "ID"))

ttm_transit <- ttm_transit %>%
  left_join(dwelling_network_points_2016 %>%
              st_drop_geometry() %>%
              dplyr::transmute(ID, 
                            Population_55to69,
                            Population_70plus,
                            Urban_type = Type,
                            DAUID, 
                            TAZUID),
            by = c("UID" = "ID"))

ttm_walk <- ttm_walk %>%
  left_join(dwelling_network_points_2016 %>%
              st_drop_geometry() %>%
              dplyr::transmute(ID, 
                            Population_55to69,
                            Population_70plus,
                            Urban_type = Type,
                            DAUID, 
                            TAZUID),
            by = c("UID" = "ID"))
```

```{r join-pharmacy-data-to-travel-time-tables, include=FALSE}
# Data preparation
## Join the pharmacy data to the travel time matrices:
ttm_car <- ttm_car %>%
  left_join(pharmacy_locations %>%
              st_drop_geometry() %>%
              dplyr::select(ID, 
                            Type),
            by = c("OBJECTID" = "ID"))

ttm_transit <- ttm_transit %>%
  left_join(pharmacy_locations %>%
              st_drop_geometry() %>%
              dplyr::select(ID, 
                            Type),
            by = c("OBJECTID" = "ID"))

ttm_walk <- ttm_walk %>%
  left_join(pharmacy_locations %>%
              st_drop_geometry() %>%
              dplyr::select(ID, 
                            Type),
            by = c("OBJECTID" = "ID"))
```

```{r drop-missing-values, include=FALSE}
# Data preparation
## Drop parcels with missing values:
ttm_car <- ttm_car %>%
  drop_na(TAZUID,
          Population_55to69)

ttm_transit <- ttm_transit %>%
  drop_na(TAZUID,
          Population_55to69)

ttm_walk <- ttm_walk %>%
  drop_na(TAZUID,
          Population_55to69)
```

```{r proportion-of-trips-by-mode, include=FALSE}
# Data preparation
## Calculate proportion of trips by mode per TAZ:
modes <- modes_55to69 %>%
  mutate(p_car = Driver/(Driver + Transit + Walk),
         p_transit = Transit/(Driver + Transit + Walk),
         p_walk = Walk/(Driver + Transit + Walk))
```

```{r complete-proportion-of-trips-by-mode, include=FALSE}
# Data preparation
## Assume that TAZ with NAs have a proportion of car trips of one:
modes <- modes %>%
  dplyr::mutate(p_car = replace_na(p_car, 1),
                p_transit = replace_na(p_transit, 0),
                p_walk = replace_na(p_walk, 0))
```

```{r join-proportion-of-trips-to-travel-time-tables, include=FALSE}
# Data preparation
## Join proportion of trips by mode to travel time tables:
ttm_car <- ttm_car %>%
  left_join(modes %>%
              st_drop_geometry() %>%
              transmute(TAZUID,
                        p_car,
                        p_transit,
                        p_walk),
            by = "TAZUID")

ttm_transit <- ttm_transit %>%
  left_join(modes %>%
              st_drop_geometry() %>%
              transmute(TAZUID,
                        p_car,
                        p_transit,
                        p_walk),
            by = "TAZUID")

ttm_walk <- ttm_walk %>%
  left_join(modes %>%
              st_drop_geometry() %>%
              transmute(TAZUID,
                        p_car,
                        p_transit,
                        p_walk),
            by = "TAZUID")
```

```{r income-quintiles, include=FALSE}
# Data preparation
## Calculate quintiles of median household income in Hamilton:
mhi_q <- data_da_2016 %>%
  dplyr::filter(`Region Name` == "Hamilton") %>%
  pull(`v_CA16_2397: Median total income of households in 2015 ($)`) %>%
  quantile(c(0, 0.2, 0.4, 0.6, 0.8, 1), 
           na.rm = TRUE)
```

```{r label-income-quintiles, include=FALSE}
# Data preparation
## Label the DAs according to income quintiles
data_da_2016 <- data_da_2016 %>%
  rename(Median_household_income = `v_CA16_2397: Median total income of households in 2015 ($)`) %>%
  mutate(income_quintile = case_when(Median_household_income < mhi_q[2] ~ "Bottom 20%",
                                     Median_household_income >= mhi_q[2] & Median_household_income < mhi_q[3] ~ "Fourth 20%",
                                     Median_household_income >= mhi_q[3] & Median_household_income < mhi_q[4] ~ "Third 20%",
                                     Median_household_income >= mhi_q[4] & Median_household_income < mhi_q[5] ~ "Second 20%",
                                     Median_household_income >= mhi_q[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile,
                                  levels = c("Top 20%",
                                             "Second 20%",
                                             "Third 20%",
                                             "Fourth 20%",
                                             "Bottom 20%"),
                                  ordered = TRUE))
```

```{r minimum-travel-time-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Retrieve the minimum travel time from a point to a pharmacy by each mode for people aged 55 to 69:
min_ttm_car <- ttm_car %>%
  dplyr::filter(Type == "Pilot") %>%
  group_by(UID) %>%
  summarize(travel_time = min(travel_time),
            Population = first(Population_55to69),
            Urban_type = first(Urban_type),
            DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            p_car = first(p_car),
            p_transit = first(p_transit),
            p_walk = first(p_walk),
            .groups = "drop")

min_ttm_transit <- ttm_transit %>%
  dplyr::filter(Type == "Pilot") %>%
  group_by(UID) %>%
  summarize(travel_time = min(travel_time),
            Population = first(Population_55to69),
            Urban_type = first(Urban_type),
            DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            p_car = first(p_car),
            p_transit = first(p_transit),
            p_walk = first(p_walk),
            .groups = "drop")

min_ttm_walk <- ttm_walk %>%
  dplyr::filter(Type == "Pilot") %>%
  group_by(UID) %>%
  summarize(travel_time = min(travel_time),
            Population = first(Population_55to69),
            Urban_type = first(Urban_type),
            DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            p_car = first(p_car),
            p_transit = first(p_transit),
            p_walk = first(p_walk),
            .groups = "drop")
```

```{r weighted-travel-time-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Calculate weighted travel time:
min_ttm_car <- min_ttm_car %>%
  mutate(weighted_travel_time = travel_time * p_car)

min_ttm_transit <- min_ttm_transit %>%
  mutate(weighted_travel_time = travel_time * p_transit)

min_ttm_walk <- min_ttm_walk %>%
  mutate(weighted_travel_time = travel_time * p_walk)
```

```{r bind-travel-time-tables-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Bind the three tables:
min_ttm <- rbind(min_ttm_car,
                 min_ttm_transit,
                 min_ttm_walk)
```

```{r weighted-travel-time-by-mode-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Bind the three tables and label them by mode:
min_ttm_mode_baseline <- rbind(data.frame(min_ttm_car, Mode = "Car"),
                               data.frame(min_ttm_transit, Mode = "Transit"),
                               data.frame(min_ttm_walk, Mode = "Walking"))
```

```{r aggregate-weighted-travel-times-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Aggregate weighted travel times by mode:
min_ttm <- min_ttm %>%
  group_by(UID) %>%
  summarize(DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            Population = first(Population),
            Urban_type = first(Urban_type),
            weighted_travel_time = sum(weighted_travel_time))
```

```{r calculate-person-hours-travel-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Calculate person-hours of travel:
min_ttm <- min_ttm %>%
  mutate(weighted_person_hours = Population * weighted_travel_time/60)
```

```{r summarize-results-by-taz-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Summarize by TAZ and join geometry:
min_ttm_taz <- min_ttm %>%
  group_by(TAZUID) %>%
  summarize(weighted_travel_time = mean(weighted_travel_time),
            weighted_person_hours = sum(weighted_person_hours)) %>%
  left_join(modes %>%
              dplyr::select(TAZUID),
            by = "TAZUID") %>%
  st_as_sf()
```

```{r weighted-travel-time-map-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Map weighted travel time:
ttm_map_baseline <- ggplot() +
  geom_sf(data = min_ttm_taz,
          aes(fill = weighted_travel_time),
          color = NA) +
  scale_fill_distiller(name = "Expected travel time (min)",
                       palette = "OrRd", 
                       direction = 1) +
  geom_sf(data = urban_types,
          aes(color = Type),
          fill = NA) +
  scale_color_manual(name = "",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen"))  +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))
```

```{r person-hours-traveled-map-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Map weighted person-hours of travel:
pht_map_baseline <- ggplot() +
  geom_sf(data = min_ttm_taz,
          aes(fill = weighted_person_hours),
          color = NA) +
  scale_fill_distiller(name = "Person-hours of travel",
                       palette = "OrRd", 
                       direction = 1) +
  geom_sf(data = urban_types,
          aes(color = Type),
          fill = NA) +
  scale_color_manual(name = "",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen"))  +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))
```

```{r join-income-information-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Join median household information to table:
min_ttm <- min_ttm %>%
  left_join(data_da_2016 %>%
              st_drop_geometry() %>%
              dplyr::transmute(DAUID = GeoUID,
                               income_quintile),
            by = "DAUID")
```

```{r plot-distribution-travel-time-by-income-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Distribution of weighted travel time by income:
ttm_income_baseline <- ggplot(data = min_ttm, 
       aes(x = weighted_travel_time, 
           color= income_quintile)) +
  geom_density(size = 0.5, adjust = 2.5) +
  labs(x = "Travel time (in minutes; weighted by mode)") +
  scale_color_manual(values = c("Top 20%" = "blue",
                                "Second 20%" = "steelblue1",
                                "Third 20%" = "gray40", 
                                "Fourth 20%" = "tomato",
                                "Bottom 20%" = "red")) +
  theme_minimal() +
  theme(text = element_text(size = 7),
        legend.text = element_text(size = 7))
```

```{r plot-distribution-travel-time-by-region-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Distribution of weighted travel time by region:
ttm_region_baseline <- ggplot(data = min_ttm %>%
                                drop_na(Urban_type), 
       aes(x = weighted_travel_time, 
           color= Urban_type)) +
  geom_density(size = 0.5, adjust = 3) +
  labs(x = "Travel time (in minutes; weighted by mode)") +
  scale_color_manual(name = "",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen"))  +
  theme_minimal() +
  theme(text = element_text(size = 7),
        legend.text = element_text(size = 7))
```

```{r total-weighted-person-hours-by-group-baseline, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Total weighted person-hours by income:
person_hours_income_baseline <- min_ttm %>%
  group_by(income_quintile) %>%
  summarize(total_weighted_person_hours = sum(weighted_person_hours),
            total_population = sum(Population),
            .groups = "drop") %>%
  mutate(ratio = total_weighted_person_hours/total_population) %>%
  rename(group = income_quintile)

## Total weighted person-hours by region:
person_hours_region_baseline <- min_ttm %>%
  group_by(Urban_type) %>%
  summarize(total_weighted_person_hours = sum(weighted_person_hours),
            total_population = sum(Population),
            .groups = "drop") %>%
  drop_na() %>%
  mutate(ratio = total_weighted_person_hours/total_population) %>%
  rename(group = Urban_type)

person_hours_baseline <- rbind(person_hours_income_baseline,
                               person_hours_region_baseline)
```

```{r minimum-travel-time, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Retrieve the minimum travel time from a point to a pharmacy by each mode for people aged 55 to 69:
min_ttm_car <- ttm_car %>%
  #dplyr::filter(Type == "Pilot") %>%
  group_by(UID) %>%
  summarize(travel_time = min(travel_time),
            Population = first(Population_55to69),
            Urban_type = first(Urban_type),
            DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            p_car = first(p_car),
            p_transit = first(p_transit),
            p_walk = first(p_walk),
            .groups = "drop")

min_ttm_transit <- ttm_transit %>%
  #dplyr::filter(Type == "Pilot") %>%
  group_by(UID) %>%
  summarize(travel_time = min(travel_time),
            Population = first(Population_55to69),
            Urban_type = first(Urban_type),
            DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            p_car = first(p_car),
            p_transit = first(p_transit),
            p_walk = first(p_walk),
            .groups = "drop")

min_ttm_walk <- ttm_walk %>%
  #dplyr::filter(Type == "Pilot") %>%
  group_by(UID) %>%
  summarize(travel_time = min(travel_time),
            Population = first(Population_55to69),
            Urban_type = first(Urban_type),
            DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            p_car = first(p_car),
            p_transit = first(p_transit),
            p_walk = first(p_walk),
            .groups = "drop")
```

```{r weighted-travel-time, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Calculate weighted travel time:
min_ttm_car <- min_ttm_car %>%
  mutate(weighted_travel_time = travel_time * p_car)

min_ttm_transit <- min_ttm_transit %>%
  mutate(weighted_travel_time = travel_time * p_transit)

min_ttm_walk <- min_ttm_walk %>%
  mutate(weighted_travel_time = travel_time * p_walk)
```

```{r bind-travel-time-tables, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Bind the three tables:
min_ttm <- rbind(min_ttm_car,
                 min_ttm_transit,
                 min_ttm_walk)
```

```{r weighted-travel-time-by-mode, include=FALSE}
# Baseline analysis (PILOT LOCATIONS ONLY)
## Bind the three tables and label them by mode:
min_ttm_mode <- rbind(data.frame(min_ttm_car, Mode = "Car"),
                      data.frame(min_ttm_transit, Mode = "Transit"),
                      data.frame(min_ttm_walk, Mode = "Walking"))
```

```{r aggregate-weighted-travel-times, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Aggregate weighted travel times by mode:
min_ttm <- min_ttm %>%
  group_by(UID) %>%
  summarize(DAUID = first(DAUID),
            TAZUID = first(TAZUID),
            Population = first(Population),
            Urban_type = first(Urban_type),
            weighted_travel_time = sum(weighted_travel_time))
```

```{r calculate-person-hours-travel, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Calculate person-hours of travel:
min_ttm <- min_ttm %>%
  mutate(weighted_person_hours = Population * weighted_travel_time/60)
```

```{r summarize-results-by-taz, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Summarize by TAZ and join geometry:
min_ttm_taz <- min_ttm %>%
  group_by(TAZUID) %>%
  summarize(weighted_travel_time = mean(weighted_travel_time),
            weighted_person_hours = sum(weighted_person_hours)) %>%
  left_join(modes %>%
              dplyr::select(TAZUID),
            by = "TAZUID") %>%
  st_as_sf()
```

```{r weighted-travel-time-map, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Map weighted travel time. This is the average of the weighted travel time per TAZ:
ttm_map <- ggplot() +
  geom_sf(data = min_ttm_taz,
          aes(fill = weighted_travel_time),
          color = NA) +
  geom_sf(data = hamilton_cma,
          fill = NA) +
  scale_fill_distiller(palette = "OrRd", 
                       direction = 1)
```

```{r total-person-hours-map, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Map weighted person-hours of travel. This is the total person-hours traveled per TAZ:
pht_map <- ggplot() +
  geom_sf(data = min_ttm_taz,
          aes(fill = weighted_person_hours),
          color = NA) +
  geom_sf(data = hamilton_cma,
          fill = NA) +
  scale_fill_distiller(palette = "OrRd", 
                       direction = 1)
```

```{r join-income-information, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Join median household information to table:
min_ttm <- min_ttm %>%
  left_join(data_da_2016 %>%
              st_drop_geometry() %>%
              dplyr::transmute(DAUID = GeoUID,
                               income_quintile),
            by = "DAUID")
```

```{r plot-distribution-travel-time-by-income, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Distribution of weighted travel time by income:
ttm_income <- ggplot(data = min_ttm, 
       aes(x = weighted_travel_time, 
           color= income_quintile)) +
  geom_density(size = 0.5, adjust = 2.5) +
  labs(x = "Travel time (in minutes; weighted by mode)") +
  scale_color_manual(name = "Income quintile",
                     values = c("Top 20%" = "blue",
                                "Second 20%" = "steelblue1",
                                "Third 20%" = "gray40", 
                                "Fourth 20%" = "tomato",
                                "Bottom 20%" = "red")) +
  theme_minimal() +
  theme(text = element_text(size = 7),
        legend.text = element_text(size = 7))
```

```{r plot-distribution-travel-time-by-region, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Distribution of weighted travel time by region:
ttm_region <- ggplot(data = min_ttm %>%
                                drop_na(Urban_type), 
       aes(x = weighted_travel_time, 
           color= Urban_type)) +
  geom_density(size = 0.5, adjust = 3) +
  labs(x = "Travel time (in minutes; weighted by mode)") +
  scale_color_manual(name = "Region",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen"))  +
  theme_minimal() +
  theme(text = element_text(size = 7),
        legend.text = element_text(size = 7))
```

```{r total-weighted-person-hours-by-group, include=FALSE}
# Scenario analysis (PILOT AND CANDIDATE LOCATIONS)
## Total weighted person-hours by income:
person_hours_income <- min_ttm %>%
  group_by(income_quintile) %>%
  summarize(total_weighted_person_hours = sum(weighted_person_hours),
            total_population = sum(Population),
            .groups = "drop") %>%
  mutate(ratio = total_weighted_person_hours/total_population) %>%
  rename(group = income_quintile)

## Total weighted person-hours by region:
person_hours_region <- min_ttm %>%
  group_by(Urban_type) %>%
  summarize(total_weighted_person_hours = sum(weighted_person_hours),
            total_population = sum(Population),
            .groups = "drop") %>%
  drop_na() %>%
  mutate(ratio = total_weighted_person_hours/total_population) %>%
  rename(group = Urban_type)

person_hours <- rbind(person_hours_income,
                      person_hours_region)
```

```{r table-comparison, echo=FALSE}
table_comparison <- cbind(person_hours_baseline %>%
        dplyr::select(group, total_population, total_weighted_person_hours, ratio), 
      person_hours %>%
        dplyr::select(total_weighted_person_hours, ratio))

colnames(table_comparison) <- c("group", "total_population", "total_weighted_person_hours_baseline", "ratio_baseline",
                                "total_weighted_person_hours", "ratio")

table_comparison <- table_comparison %>%
  mutate(total_weighted_person_hours_baseline = round(total_weighted_person_hours_baseline, 2),
         ratio_baseline = set_units(ratio_baseline, h) %>% 
           set_units(min) %>% 
           drop_units(),
         ratio_baseline = round(ratio_baseline, 2),
         total_weighted_person_hours = round(total_weighted_person_hours, 2),
         ratio = set_units(ratio, h) %>% 
           set_units(min) %>% 
           drop_units(),
         ratio = round(ratio, 2))
```

```{r comparison-travel-time-by-mode, echo=FALSE}
# Comparison of weighted travel times by mode between pilot and scenario
comparison_ttm_mode <- cbind(min_ttm_mode_baseline %>%
  group_by(Mode) %>%
  summarize(pht_mode_baseline = sum(Population * weighted_travel_time/60),
            .groups = "drop"),
  min_ttm_mode %>%
  group_by(Mode) %>%
  summarize(pht_mode = sum(Population * weighted_travel_time/60),
            .groups = "drop") %>%
    select(pht_mode)) %>% 
  mutate(diff = paste0(round((pht_mode - pht_mode_baseline)/pht_mode_baseline * 100, 2), "%", sep = ""))
```

Questions
==========================

Along with the provision of health care facilities to treat severe cases of COVID-19 [e.g., @Ghorbanzadeh2021spatial; @Pereira2021geographic], another front in the fight against the pandemic is the rolling out of vaccination programs. The Government of Ontario, in Canada, implemented a pilot program to offer vaccines in pharmacies for people 60 years and older in early 2021. In the first incarnation of this program, the province approved 325 pharmacies in Toronto, Windsor-Essex, and Kingston (the latter including the region of Frontenac-Lennox and Addington). The City of Hamilton, the third most populous urban center in the province, operated dedicated vaccination centers for people aged 70+ and mobile pop-up clinics for people aged 75+; however, it did not have any approved pharmacies under the pilot for other populations until April 1st 2021, with the announcement of 20 pharmacies in the city \footnote{\url{https://www.cbc.ca/news/canada/hamilton/astrazeneca-vaccine-hamilton-1.5972704}}. At the same time that the expansion of the pilot was announced, the province extended eligibility for this program to people 55 years and older.

Critics were swift to point out that the list of pharmacies approved for Hamilton by the province were mostly located in lower density parts of the city that are not well serviced by transit and are difficult to reach by foot\footnote{See inter alia: \url{https://twitter.com/RyanMcGreal/status/1378027149790224386?s=20} and \url{https://twitter.com/NrinderWard3/status/1378679195514060801?s=20}.}. Indeed, as seen in Figure \ref{fig:pharmacies-and-regions}, a vast majority of the pharmacies are in suburban Hamilton. Whether this is an issue is less clear-cut when we consider that Hamilton's older population skews suburban (see Figure \ref{fig:population-map}). Given the target demographic for the program, it is possible that suburban sites could be convenient for mature adults and the young old: the population aged 55 to 69 in Hamilton is approximately `r prettyNum(round((table_comparison[7,2]/1000), 2) * 1000, big.mark = ",")` suburban, `r prettyNum(round((table_comparison[8,2]/1000), 2) * 1000, big.mark = ",")` urban, and only `r prettyNum(round((table_comparison[6,2]/1000), 2) * 1000, big.mark = ",")` rural. Nevertheless, the selection of sites by the province raised some important questions\footnote{The decision-making process to select these sites appears to have been opaque, and the Mayor of the city was caught flat footed by the announcement; see: \url{https://twitter.com/FredEisenberger/status/1378350123114242053?s=20}}. As @Yu2021sustained note, good geographical coverage is a key element for a successful vaccination campaign; at the same time, siting vaccinations sites in car-oriented locations may lead to inequities in access by requiring people who do not own cars or who typically travel by other modes (mostly lower income urban residents) to spend more of their time reaching an approved pharmacy.

In this research, we investigate the accessibility implications of the sites selected for the pilot vaccination program. Concretely, we ask:

- What is the estimated travel time needed to reach the nearest approved pharmacy, assuming that every person requires a vaccine?
- What is the distribution of travel times across the population of the city?
- How does the cost of time needed for travel and its distribution change with the addition of candidate vaccination sites in urban Hamilton?

We concentrate on the 55 to 69 year old population segment because the older 70+ group have access to other dedicated facilities besides those in the provincial pharmacy pilot (including pop-up clinics for people 75+).

```{r pharmacies-map, echo=FALSE, fig.height=4, fig.align="center", fig.cap="\\label{fig:pharmacies-and-regions}Regions with the City of Hamilton; the location of pharmacies in pilot is shown (black triangles) along with urban locations for scenario analysis (white circles)."}
#map_pharmacies <-
ggplot() +
  geom_sf(data = urban_types,
          aes(fill = Type),
          color = "lightgray") +
  geom_sf(data = pharmacy_locations,
          aes(shape = Type,
              color = Type),
          size = 2.5,
          show.legend = FALSE) +
  scale_color_manual(values = c("Candidate" = "white",
                                "Pilot" = "black")) +
  scale_fill_manual(values = c("Urban" = "firebrick2",
                               "Suburban" = "cornflowerblue",
                               "Rural" = "forestgreen")) +
  guides(fill = guide_legend(title = "")) + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))
```

```{r save-pharmacies-map-1, include=FALSE, eval=FALSE}
ggplot() +
  geom_sf(data = urban_types,
          aes(fill = Type),
          color = "lightgray") +
  geom_sf(data = pharmacy_locations %>%
            filter(Type == "Pilot"),
          aes(shape = Type,
              color = Type),
          size = 2.5,
          show.legend = FALSE) +
  scale_shape_manual(values = c("Candidate" = 16,
                                "Pilot" = 17)) +
  scale_color_manual(values = c("Candidate" = "white",
                                "Pilot" = "black")) +
  scale_fill_manual(values = c("Urban" = "firebrick2",
                               "Suburban" = "cornflowerblue",
                               "Rural" = "forestgreen")) +
  guides(fill = guide_legend(title = "")) + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))

ggsave("images/pharmacies-map-1.jpeg")
```

```{r save-pharmacies-map-2, include=FALSE, eval=FALSE}
ggplot() +
  geom_sf(data = urban_types,
          aes(fill = Type),
          color = "lightgray") +
  geom_sf(data = pharmacy_locations,
          aes(shape = Type,
              color = Type),
          size = 2.5,
          show.legend = FALSE) +
  scale_shape_manual(values = c("Candidate" = 16,
                                "Pilot" = 17)) +
  scale_color_manual(values = c("Candidate" = "white",
                                "Pilot" = "black")) +
  scale_fill_manual(values = c("Urban" = "firebrick2",
                               "Suburban" = "cornflowerblue",
                               "Rural" = "forestgreen")) +
  guides(fill = guide_legend(title = "")) + 
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))

ggsave("images/pharmacies-map-2.jpeg")
```

```{r population-map, echo=FALSE, warning=FALSE, fig.height=8, fig.align="center", fig.cap="\\label{fig:population-map}Distribution of population aged 55+ in the City of Hamilton."}
ggplot() +
  geom_sf(data = data_da_2016 %>%
            filter(`Region Name` == "Hamilton") %>%
            dplyr::select(Population_55to69, Population_70plus) %>%
            pivot_longer(cols = -geometry, names_to = "Age", values_to = "Population") %>%
            mutate(Age = factor(Age,
                                levels = c("Population_55to69",
                                           "Population_70plus"),
                                labels = c("Age: 55 to 69",
                                           "Age: 70+"),
                                ordered = TRUE)) %>%
            st_as_sf(),
          aes(fill = Population),
          color = NA) +
  geom_sf(data = urban_types,
          aes(color = Type),
          fill = NA) +
  scale_fill_distiller(name = "Population",
                       breaks = c(20, 50, 150, 400),
                       palette = "OrRd", 
                       trans = "log",
                       direction = 1) +
  scale_color_manual(name = "",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen"))  +
  facet_wrap(~ Age, nrow = 2) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))
```

```{r save-population-map, include=FALSE, fig.height=8, dpi=300, eval=FALSE}
ggplot() +
  geom_sf(data = data_da_2016 %>%
            filter(`Region Name` == "Hamilton") %>%
            dplyr::select(Population_55to69, Population_70plus) %>%
            pivot_longer(cols = -geometry, names_to = "Age", values_to = "Population") %>%
            mutate(Age = factor(Age,
                                levels = c("Population_55to69",
                                           "Population_70plus"),
                                labels = c("Age: 55 to 69",
                                           "Age: 70+"),
                                ordered = TRUE)) %>%
            st_as_sf(),
          aes(fill = Population),
          color = NA) +
  geom_sf(data = urban_types,
          aes(color = Type),
          fill = NA) +
  scale_fill_distiller(name = "Population",
                       breaks = c(20, 50, 150, 400),
                       palette = "OrRd", 
                       trans = "log",
                       direction = 1) +
  scale_color_manual(name = "",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen"))  +
  facet_wrap(~ Age, nrow = 2) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))

ggsave("images/population-map.jpeg")

ggplot() +
  geom_sf(data = data_da_2016 %>%
            filter(`Region Name` == "Hamilton"),
          aes(fill = Population_55to69),
          color = NA) +
  geom_sf(data = urban_types,
          aes(color = Type),
          fill = NA) +
  scale_fill_distiller(name = "Population aged 55 to 69",
                       breaks = c(20, 50, 150, 400),
                       palette = "OrRd", 
                       trans = "log",
                       direction = 1) +
  scale_color_manual(name = "",
                     values = c("Urban" = "firebrick2",
                                "Suburban" = "cornflowerblue",
                                "Rural" = "forestgreen")) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text = element_text(size = 7),
        legend.position = "bottom",
        legend.text = element_text(size = 8))

ggsave("images/population-55-to-69-map.jpeg")
```

Methods
============

We use data from the following sources.

- Urban, suburban, and rural boundary definitions from the City of Hamilton\footnote{\url{https://open.hamilton.ca/}}
- Population and median total household income for 2016 by Dissemination Area (DA) boundary using the `cancensus` package [@vonBergmann2021cancensus]
- Modal split by age for traffic analysis zones (TAZ) from the 2016 Transportation Tomorrow Survey (TTS)\footnote{\url{http://dmg.utoronto.ca/}}
- Locations of pilot pharmacies from public records
- Locations of three additional major chain pharmacies for scenario testing
- Residential units by land parcel

Using the population aged 55 to 69 y.o. we first calculate the average number of people per dwelling and assign them proportionally to the dwellings by parcel. Second, median household incomes and modal splits across walking, transit, and car are joined to the parcels. Our data are from 2016 (the most recent available), but we assume they are largely representative of spatial, mobility, and demographic trends in the city. In the absence of travel behavior information during the pandemic we suggest this is an inevitable limitation of the analysis. Third, we use the `r5r` [@Pereira2021r5r] package to calculate the travel time from each parcel to all pharmacies by three modes using a cutoff value of 180 min and a maximum walking distance of 10,000 m. These thresholds are chosen to maximize the coverage of the travel time tables, but the analysis, as explained next, is based on travel to the nearest approved pharmacy. Once we obtained travel time tables with population, proportion of trips by mode, and income information, we calculated the expected travel time $ett$ from each parcel $i$ to a pharmacy $j$ as follows:
$$
ett_i = p^c_i \min(tt^c_{ij}) + p^t_i \min(tt^t_{ij}) + p^w_i \min(tt^w_{ij})
$$
\noindent where $p^k_i$ is the proportion of trips by mode $k$ in the TAZ of parcel $i$, and $tt^k_{ij}$ is the vector of travel times from parcel $i$ to the pharmacies. The expected travel time is thus the weighted sum of travel times to the _nearest_ pharmacy, with the weights given by the expected modal split in the TAZ.

The expected travel time $i$ was multiplied by the assigned population in parcel $i$ to obtain a measure of person-hours of travel ($PHT$) that gives a cumulative cost (in travel time) for a population group. This is as follows:
$$
PHT_i = P_i\cdot ett_i
$$

Please note that this paper is a reproducible research document [see @Brunsdon2020opening] conducted using open source tools for transportation analysis [@Lovelace2021open]. The code and data necessary to reproduce the analysis are available in a public repository\footnote{\url{https://github.com/paezha/Accessibility-Pharmacies-Hamilton-Vaccines}}.

Findings
===================

The top panel of Figure \ref{fig:maps-baseline} shows the average multimodal travel times (weighted by modal shares) by TAZ in Hamilton. It is apparent that travel times tend to be lower in much of suburban Hamilton and higher in the urban core and some rural parts of the city, particularly to the west. This is unsurprising, given the higher probability of travel by car and the predominantly suburban character of the vaccination sites. However, even accounting for the distribution of population, this leads to large disparities in the number of person-hours of travel across the city, with a concentration of the burden of travel in the urban core and the rural west (see bottom panel of Figure \ref{fig:maps-baseline}). 

The disparities are not trivial. 

As seen in Table \ref{tab:distribution-results}, under the pilot program approximately `r round((table_comparison[5,2] + table_comparison[4,2]) / (sum(table_comparison[1:5,2])) * 100, 2)`\% of people live in DAs in the bottom 40% of the median household income scale, but they account for `r round((table_comparison[5,3] + table_comparison[4,3]) / (sum(table_comparison[1:5,3])) * 100, 2)`\% of the total person-hours of travel. In contrast, `r round((table_comparison[1,2] + table_comparison[2,2]) / (sum(table_comparison[1:5,2])) * 100, 2)`\% of people aged 55 to 69 in DAs in the top 40% of the median household income scale accrue only `r round((table_comparison[1,3] + table_comparison[2,3]) / (sum(table_comparison[1:5,3])) * 100, 2)`\% of the total person-hours of travel. Where the mean travel time of residents of DAs with high median household income is `r table_comparison[1, 4] %>% set_units(h) %>% set_units(min) %>% round() %>% drop_units()` minutes, residents of lower income DAs average `r table_comparison[5, 4] %>% set_units(h) %>% set_units(min) %>% round() %>% drop_units()` minutes in travel time. In addition to greater average travel times, the top-left panel in Figure \ref{fig:results} suggests residents in lower income DAs also see substantially larger variability in travel times and some may face considerably longer trip durations.

There are also important disparities by region. As shown in Table \ref{tab:distribution-results}, the urban and rural populations in Hamilton are approximately `r round((table_comparison[6,2] + table_comparison[8,2])/ (sum(table_comparison[6:8,2])) * 100, 2)`\% of the population but they bear `r round((table_comparison[6,3] + table_comparison[8,3])/ (sum(table_comparison[6:8,3])) * 100, 2)`\% of the total person-hours of travel, with also much greater variability in expected travel times (Figure \ref{fig:results}, bottom-left panel). 

For comparison purposes we consider a scenario with some modest additions to the list of pharmacies in the provincial pilot. For the scenario we choose three plausible sites (i.e., pharmacies that belong to chains already used by the province for this program). We do not conduct a systematic search of candidate locations. We note that some papers propound the design of optimal accessibility landscapes [e.g., @Horner2008optimal; @Paez2013exploring], but such an endeavor is beyond the scope of a short findings paper. Instead, we repeat the analysis after including the three urban sites shown in white circles in Figure \ref{fig:pharmacies-and-regions}.

The results of this scenario appear in the last two columns of Table \ref{tab:distribution-results} and the two right panels of Figure \ref{fig:results}. While all income groups benefit from the addition of these three sites with shorter mean trip durations, the most remarkable difference is the large reduction in the disparities between residents in DAs with lower levels of income. The top-right panel of Figure \ref{fig:results} shows that the distribution of expected travel time is now more in line for all income groups, even if the bottom two income quintiles still have somewhat wider spreads. Unsurprisingly, the addition of three urban vaccination sites does not have a large impact for rural residents. The scenario also suggests that important reductions in the need for travel can be achieved. While this is true for all three modes of transportation (see Table \ref{tab:distribution-results-by-mode}), the reductions are more substantial for travel by transit and walking. At a time that non-pharmaceutical interventions call for reductions in mobility, it is clearly in the public interest to reduce the need to travel even for this essential purpose. 

The results indicate that the locations chosen by the province for the pilot vaccination program do not serve urban or rural residents of the city well, and there are some important questions regarding equity of access to the program. As shown, a disproportionate burden in the cost of travel (in terms of travel time) falls on lower income urban populations and rural populations. The selection of three sensible urban locations does much to alleviate disparities in the burden of transportation - and the fact that this is so despite the lack of a systematic search for optimal locations does not speak well of the effort put forth by the province in the choice of vaccination sites. The situation is different for rural populations, mainly because there are not many candidate locations in rural parts of the city; in this case, it is likely that increasing access for residents there will necessitate an expansion of the existing mobile vaccination pop-up clinic program\footnote{\url{https://www.hamilton.ca/government-information/news-centre/news-releases/hamiltons-covid-19-vaccination-program-expansion-1}}.

```{r figure-maps-baseline, echo=FALSE, warning=FALSE, fig.height=7, fig.align="center", fig.cap="\\label{fig:maps-baseline}Average expected travel time by TAZ (in minutes) and total person-hours of travel by TAZ."}
grid.arrange(ttm_map_baseline +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Pilot"),
                       shape = 17,
                       size = 1,
                       color = "black",
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 16,
                       color = "white",
                       size = 1,
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 1,
                       color = "black",
                       size = 1,
                       show.legend = FALSE) +
               theme(legend.title = element_text(size = 8),
                     legend.text = element_text(size = 7), 
                     axis.text = element_text(size = 7)),
             pht_map_baseline  +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Pilot"),
                       shape = 17,
                       size = 1,
                       color = "black",
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 16,
                       color = "white",
                       size = 1,
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 1,
                       color = "black",
                       size = 1,
                       show.legend = FALSE) +
               theme(legend.title = element_text(size = 8),
                     legend.text = element_text(size = 7), 
                     axis.text = element_text(size = 7)),
             nrow = 2)
```

```{r save-figure-maps-baseline, include=FALSE, fig.height=8, dpi=600, eval=FALSE}
ttm_plot <- ttm_map_baseline +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Pilot"),
                       shape = 17,
                       size = 1,
                       color = "black",
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 16,
                       color = "white",
                       size = 1,
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 1,
                       color = "black",
                       size = 1,
                       show.legend = FALSE) +
               theme(legend.title = element_text(size = 8),
                     legend.text = element_text(size = 7), 
                     axis.text = element_text(size = 7))

pht_plot <- pht_map_baseline  +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Pilot"),
                       shape = 17,
                       size = 1,
                       color = "black",
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 16,
                       color = "white",
                       size = 1,
                       show.legend = FALSE) +
               geom_sf(data = pharmacy_locations %>%
                         dplyr::filter(Type == "Candidate"),
                       shape = 1,
                       color = "black",
                       size = 1,
                       show.legend = FALSE) +
               theme(legend.title = element_text(size = 8),
                     legend.text = element_text(size = 7), 
                     axis.text = element_text(size = 7))

ttm_plot/pht_plot

ggsave("images/baseline-maps.jpeg")
```

```{r table-results, echo=FALSE}
table_comparison  %>%
  kable("html",
        booktabs = TRUE,
        col.names = c("Group", 
                      "Population", 
                      "Total PHT", 
                      "Minutes per person", 
                      "Total PHT",
                      "Minutes per person"),
        digits = 3,
        align = c("lccccc"),
        caption = "\\label{tab:distribution-results}Distribution of person-hours of travel (PHT) by median total household income and region: pilot locations only (baseline), and scenario with three urban locations added") %>% 
  add_header_above(c(" " = 2, "Pilot Program" = 2, "Scenario" = 2)) %>%
  kable_styling(bootstrap_options = c("striped")) %>%
  pack_rows("Income Quintile", 1, 5) %>%
  pack_rows("Region", 6, 8) %>%
  footnote(general = "The population totals differ due to small differences in the classification of the regions")
```

```{r table-results-by-mode, echo=FALSE}
comparison_ttm_mode %>%
  kable("html",
        booktabs = TRUE,
        col.names = c("Mode", 
                      "Baseline PHT", 
                      "Scenario PHT", 
                      "Reduction in PHT"),
        digits = 3,
        align = c("lccc"),
        caption = "\\label{tab:distribution-results-by-mode}Distribution of person-hours of travel (PHT) by mode of transportation: pilot locations only (baseline) and scenario with three urban locations added")  %>%
  kable_styling(bootstrap_options = c("striped"))
```

```{r figure-results, echo=FALSE, warning=FALSE, fig.height=4, fig.align="center", fig.cap="\\label{fig:results}Distribution of expected travel time for different population groups."}
grid.arrange(ttm_income_baseline + 
               ggtitle("By income: Pilot Program") + 
               labs(x = "") + 
               xlim(c(0,40)) +
               theme(title = element_text(size = 8), 
                     axis.text = element_text(size = 7), 
                     legend.position = "none"), 
             ttm_income + 
               ggtitle("By income: Scenario") + 
               labs(x = "", y = "") + 
               xlim(c(0,40)) +
               theme(title = element_text(size = 8), 
                     axis.text = element_text(size = 7)), 
             ttm_region_baseline + 
               ggtitle("By region: Pilot Program") + 
               xlim(c(0,40)) +
               theme(title = element_text(size = 8),
                     axis.title = element_text(size = 8), 
                     axis.text = element_text(size = 7), 
                     legend.position = "none"), 
             ttm_region + 
               ggtitle("By region: Scenario") + 
               labs(y = "") +
               xlim(c(0,40)) +
               theme(title = element_text(size = 8),
                     axis.title = element_text(size = 8), 
                     axis.text = element_text(size = 7)), 
             nrow = 2)
```

```{r save-figure-results, include=FALSE, fig.height=8, dpi=900, eval=FALSE}
ttm_income_baseline <- ttm_income_baseline + 
               ggtitle("By income: Pilot Program") + 
               labs(x = "") + 
               xlim(c(0,40)) +
               theme(title = element_text(size = 8), 
                     axis.text = element_text(size = 7), 
                     legend.position = "none")

ttm_income <- ttm_income + 
               ggtitle("By income: Scenario") + 
               labs(x = "", y = "") + 
               xlim(c(0,40)) +
               theme(title = element_text(size = 8), 
                     axis.text = element_text(size = 7)) 

ttm_region_baseline <- ttm_region_baseline + 
               ggtitle("By region: Pilot Program") + 
               xlim(c(0,40)) +
               theme(title = element_text(size = 8),
                     axis.title = element_text(size = 8), 
                     axis.text = element_text(size = 7), 
                     legend.position = "none")

ttm_region <- ttm_region + 
               ggtitle("By region: Scenario") + 
               labs(y = "") +
               xlim(c(0,40)) +
               theme(title = element_text(size = 8),
                     axis.title = element_text(size = 8), 
                     axis.text = element_text(size = 7))

(ttm_income_baseline + ttm_income) / (ttm_region_baseline + ttm_region)

ggsave("images/figure-results.jpeg")
```

References {#references .unnumbered}
==========