README.Rmd

---
output:
  github_document:
    pandoc_args: --webtex
always_allow_html: true
bibliography: [Bike-share-spatial-equity/packages.bib, Bike-share-spatial-equity/references.bib]
---

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

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

# Examining spatial equity of bike share: A balanced floating catchment area approach

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

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

Desjardins, E., Higgins, C.D., Paez, A. (2022) Examining spatial equity of bike share: A balanced floating catchment area approach, Transportation Research Part D (102) 103091 (doi.org/10.1016/j.trd.2021.103091)

## Abstract

Public bicycle share programs (PBSPs) can play a role in advancing transportation equity if they make bicycling more accessible to disadvantaged populations. In Ontario, Hamilton Bike Share expanded their program in 2018 by adding twelve "equity" stations with the explicit objective of increasing access for under-serviced neighbourhoods. In this case study, we investigate differentials in accessibility to stations using a balanced floating catchment area approach and compare accessibility with and without the equity stations. We analyze population interpolated to small cells to better reflect walking to a station and conduct a sensitivity analysis at several walking time thresholds. We then reaggregate the estimated accessibility by income groups for further analysis. Our findings indicate that equity stations increased accessibility for the serviced population at every threshold examined, but the increase was relatively modest especially for population in the bottom 20% of median total household income.

## Keywords

- Public Bike Sharing Programs  
- Transport equity  
- Horizontal equity  
- Vertical equity  
- Accessibility  
- Balanced floating catchment areas  
- Hamilton  
- Ontario
- Reproducible research  

```{r load-packages-1, include=FALSE}
library(glue)
library(here)
```

```{r install-data-package, include = FALSE, renv.ignore=TRUE}
# Run only once if needed to install data package `sobiEquity`
if (!require("sobiEquity", character.only = TRUE)) {
      install.packages(glue(here::here(),'/sobiEquity_0.1.0.tar.gz'), repos=NULL, type='source')
  }
```

```{r load-packages-2, include=FALSE}
library(biscale)
library(cowplot)
library(data.table)
#library(disk.frame)
library(dplyr)
library(gdistance)
library(ggmap)
library(ggplot2)
library(ggspatial)
library(glue)
library(gridExtra)
library(here)
library(kableExtra)
library(pycno)
#library(r5r) # the r5r package requires Java Development Kit version 11, which can be downloaded from https://www.oracle.com/java/technologies/javase-jdk11-downloads.htmllibrary(readr)
library(raster)
#library(rgdal)
library(sf)
library(sobiEquity) # Data package
library(tidyr)
#library(tinytex)
library(units)
```

```{r generate-package-bibliography, include=FALSE}
knitr::write_bib(c(.packages(), "knitr", "rticles"), "packages.bib")
```

```{r invoke-data, include=FALSE}
# The data is part of data package `sobiEquity` 
data("community_downtown")
data("hamilton_cma")
data("hamilton_da_2016")
data("population_50x50")
data("sobi_service")
data("sobi_hubs")
data("ttm_walk")
```

```{r crop-das, include=FALSE}
# Crop DAs using bounding box of service area
sobi_service_bbox <- sobi_service %>% 
  st_bbox() %>%
  st_as_sfc()

cropped_da <- hamilton_da_2016 %>%
  st_crop(sobi_service_bbox) %>%
  dplyr::select(DA, geometry)
```

\newpage

1. Introduction
===================

The potential of public bicycle share programs (PBSPs) to increase bicycling levels is but one of many reasons for implementing such programs in urban areas [see @fishmanBikeshareImpactActive2015; @hosfordEvaluationImpactPublic2018; @hosfordEvaluatingImpactImplementing2019]. As a healthy, inexpensive, and convenient form of public transportation, shared bicycles can encourage individuals to take up bicycling for short local trips or first and last mile trips to other public transportation instead of using personal vehicles. These programs can also play a role in advancing transportation equity if they make bicycling more accessible to disadvantaged populations. Although PBSPs are available to the general public in over 800 cities worldwide [@fishmanBikeshareReviewRecent2016], and ought to be accessible to any individual who wishes to use them, research on PBSPs indicates that inequities persist with respect to who can use and access them.

Many PBSPs in North America now offer specific programs to address equity that primarily focus on removing cost barriers and increasing access for groups that are under-represented among existing users [@trec2019]. The locations of docking stations is a common consideration for increasing access to improve equity [@howlandCurrentEffortsMake2017]. Hamilton Bike Share (HBS), located in Hamilton, Ontario, was the only Canadian PBSP included in a North American scan of bike share equity programs [see @trec2019]. HBS was launched in 2015 and currently has over 900 operational bicycles and 130 docking stations. An equity program, _Everyone Rides Initiative_, was implemented in 2018 which expanded the program by introducing twelve "equity" stations to more disadvantaged neighbourhoods in the core service area. 

This paper examines accessibility to Hamilton Bike Share using a balanced floating catchment area (BFCA) approach. We also conduct a comparative analysis to the conventional two-step floating catchment area (2SFCA) approach to highlight the benefits of this method which, to our knowledge, has not yet been used in the cycling literature. The paper also assesses the contribution of the equity stations to reducing inequities in accessibility for different groups according to median total household income, and provides policy recommendations to further improve equity.

This paper is an example of open and reproducible research that uses only open software for transportation and statistical analysis [@bivand2020progress; @lovelace2021open]. All data were obtained from publicly available sources and organized in the form of a data package. Following best practices in spatial data science [@brunsdon2020opening], an open data product [@arribas2021open] along with the code needed to reproduce, modify or extend the analysis are available for download.^[https://github.com/paezha/Accessibility-Sobi-Hamilton] 

2. Literature Review
===================

## 2.1. Public Bicycle Share Programs

Public bicycle share programs have been implemented in over 800 cities worldwide and a great deal has been learned about their typical users to date [@fishmanBikeshareReviewRecent2016]. In many cities, males use bike share more than females [@breyWantRideMy2017; @nickkarSpatialtemporalGenderLand2019; @ogilvieInequalitiesUsagePublic2012; @reillyGenderDisparitiesNew2020; @wangGenderGapGenerators2019; @wintersWhoAreSuperusers2019] as do younger age cohorts [@breyWantRideMy2017; @buckAreBikeshareUsers2013; @fullerUseNewPublic2011]. However, one study found that bike share users in Washington, DC were more likely to be female [@buckAreBikeshareUsers2013], which suggests that the gender gap among bicyclists who use PBSPs is less disparate than the gap for personal bicycle use [@fishmanBikeshareReviewRecent2016]. There is some evidence that bike share users are less likely to own a car [@buckAreBikeshareUsers2013; @reillyNoncyclistsFrequentCyclists2020]. However, the relationship between income or education and bike share use is less clear-cut. Stations in disadvantaged communities in Chicago have been found to generate most of the average annual trips [@qianBikesharingEquityDisadvantaged2020] and individuals from minority or lower socioeconomic status neighborhoods in Minneapolis-St. Paul also used the city's PBSP more [@wangNeighborhoodSociodemographicCharacteristics2019]. Similar findings were reported in London [@ogilvieInequalitiesUsagePublic2012]. Being university educated was a significant correlate of bike share use in Montreal, Canada [@fullerUseNewPublic2011]. Not coincidentally, financial savings have been found to motivate those on a low income to use bike share [@fishmanBikeshareReviewRecent2016]. 

## 2.2. Equity of PBSPs

The introduction of PBSPs has been accompanied by a flurry of research focusing on who benefits from them. Many studies have examined differences in demographics or socioeconomic status between those who use or have access to PBSPs and those who don't, while other research has explored spatial inequities in where stations are located [e.g., @hosfordWhoArePublic2018; @mooneyFreedomStationSpatial2019; @hullgrassoBikeShareEquity2020; @qianBikesharingEquityDisadvantaged2020; @qianBikeshareDestinationChoices2021; @smith2015exploring]. In Chicago, Qian and Jaller [-@qianBikesharingEquityDisadvantaged2020] estimated ridership in the city's PBSP and found that a minority of docking stations were located in disadvantaged communities, while annual members from such areas had a lower share of trips compared to other areas in the city. Similar results were found in Philadelphia. Despite efforts to increase equity within the city's PBSP, census block groups with lower median income generated fewer trips [@caspiBikesharingPhiladelphiaLowerincome2019]. Trips from stations in such areas were utilitarian (e.g., commuting to work), which points to the importance of ensuring equitable access [@caspiBikesharingPhiladelphiaLowerincome2019]. In the case of Seattle, all neighbourhoods were found to have some level of access to dockless bicycles but those with higher incomes and more residents of higher education had more bicycles [@mooneyFreedomStationSpatial2019]. Babagoli et al. [-@babagoliExploringHealthSpatial2019] also found that neighbourhoods in New York City with higher affluence had the greatest proportion of Citi Bike stations.

Chen et al. [-@chenExploringEquityPerformance2019] distinguish two types of equity related to bike share, horizontal and vertical, based on the work of Delbosc and Currie [@delboscUsingLorenzCurves2011]. Horizontal equity leads to balanced or equal distribution of accessibility and costs for all similar groups, while vertical equity would involve greater and targeted access for under-represented or disadvantaged populations [@chenExploringEquityPerformance2019]. Both are of interest to researchers and transportation planners since they are often linked in that advantage, or conversely disadvantage, has spatial patterns. Bike share equity programs [see @trec2019; @howlandCurrentEffortsMake2017] can be considered efforts to improve vertical equity since they favour groups that have benefited less from PBSPs through the placement of stations or equitable fee structures.

On the whole, existing studies highlight the need for PBSPs to be more accessible for equity populations in order to increase use beyond the "typical" users. This has been the focus of recent research [see, among others, @auchinclossDesignBaselineDescription2020; @hullgrassoBikeShareEquity2020; @macarthurAdaptiveBikeShare2020]. Offering more people the option of using sustainable and active transportation, particularly those who have lower socioeconomic status and might benefit the most, is a worthy policy goal for cities with PBSPs. However, exploring transportation equity by investigating where docking stations are located, often using neighborhoods or census tracts as the geographical unit of analysis, can ignore or miss the benefits that may be derived from adjacent zones. Meaning that, stations may be lacking in certain neighborhoods but there may be stations accessible within a reasonable walking time. This is where geographical accessibility becomes an important consideration.

## 2.3. Accessibility Approaches

Accessibility has been applied in both a positive and normative way to inform transportation planning [@paezMeasuringAccessibilityPositive2012], but its utility to this field has evolved over the past century and has increasingly become linked with recent planning interests in prioritizing modes that are suitable for local trips like walking and cycling [@levineCenturyEvolutionAccessibility2020]. Beyond the utility derived from using shared bicycles to access destinations of value, docking stations themselves are amenities because they offer a transportation service. Therefore, the ease of reaching these stations, which are spread spatially in a given area, can affect use of the PBSP. 

The location and size of docking stations are two factors that are relevant to accessibility. Since the time or distance needed to reach a docking station decreases the potential of accessing the program, the location matters. Kabra et al. [-@kabraBikeShareSystemsAccessibility2020] found that the majority of bike share usage in Paris, France comes from areas within 300 m of stations, which amounts to 2-4 minutes walking by an adult who does not have a disability. Other studies have adopted  [-@hosfordEvaluationImpactPublic2018; -@hosfordEvaluatingImpactImplementing2019] Similar to other public amenities affected by crowding, the utility of docking stations is also limited by the _maximum_ number of bicycles that they can hold (e.g., their size). Accessibility analyses for PBSPs constitute a positive and evaluation-based approach that also has the potential to inform equity efforts. For instance, Wang and Lindsey [-@wangNewBikeShare2019] investigated whether new or relocated bike share stations increased accessibility and use, which offered important insights to improve the performance of the program.

Several approaches have been commonly used for measuring place-based accessibility, including cumulative opportunities, gravity, and utility-based measures [@geursAccessibilityEvaluationLanduse2004; @handyMeasuringAccessibilityExploration1997]. The gravity-based approach involves weighting destination opportunities, such as the quantity of bike share stations, by the time required to reach them from an origin using an impedance function [@handyMeasuringAccessibilityExploration1997; @kwanSpaceTimeIntegral1998]. While such measures are suitable for capturing the potential for reaching destinations from a given location, they do not take demand or congestion effects into account which is an important consideration when calculating accessibility for amenities such as bike share stations. 

In contrast, floating catchment area (FCA) methods have been widely employed in health care accessibility research. The benefit of this approach is that the method incorporates information on capacity and demand in calculating accessibility. Geurs and van Wee [-@geursAccessibilityEvaluationLanduse2004] note "competition for activities with restricted capacity" should be taken into account for land-use components of accessibility. The FCA approach is more appropriate and informative than calculating provider-to-population ratios (PPR) that simply divide the level of supply of a service (e.g., the number of bicycle racks at a station) by the population who have access to the service [@paezDemandLevelService2019]. In particular, the Two-Step Floating Catchment Area (2SFCA) method [@luoMeasuresSpatialAccessibility2003; @radkeSpatialDecompositionsModeling2000] produces flexible catchment areas instead of using rigid boundaries like PPR. In the first step, a ratio of supply to demand at service locations is calculated, such as the number of beds at a hospital divided by the number of people within the catchment area of the hospital, weighted by the distance involved in reaching the facility. Next, these service level ratios are allocated back to the population centers and summarized as a measure of congested accessibility. Thus, this model does a good job of considering potential crowding or competition for services. However, overlapping catchment areas from the conventional FCA approach lead to inflation of population totals and deflation of service levels across a study area which generates inaccurate or misleading accessibility estimates [@paezDemandLevelService2019]. While there have been many methodological innovations in FCA methods [e.g., @delamaterSpatialAccessibilitySuboptimally2013; @luoMeasuresSpatialAccessibility2003; @luoEnhancedTwostepFloating2009; @radkeSpatialDecompositionsModeling2000; @wanThreestepFloatingCatchment2012], a recent improvement to this approach was achieved through a simple and intuitive balancing of the impedance that addressed the effects of demand and service inflation found in earlier FCA approaches [see @paezDemandLevelService2019]. This provides more useful information because it does not assume that people are limited to service within pre-defined boundaries [@paezDemandLevelService2019].

When measuring accessibility, researchers have also taken different approaches when it comes to the aggregation of data, either by using the individual or household as the smallest unit of analysis or larger spatial zones. Previous research on bike share equity has typically used a meso- or macro-level approach with aggregated data from entire neighborhoods or census tracts [@babagoliExploringHealthSpatial2019; @mooneyFreedomStationSpatial2019; @qianBikesharingEquityDisadvantaged2020; @wangNeighborhoodSociodemographicCharacteristics2019], although there are recent exceptions [@chenExploringEquityPerformance2019; @chenUnobservedHeterogeneityTransportation2021]. This is also true for studies examining correlates of bike share demand [@wangNewBikeShare2019]. Handy and Niemeier [-@handyMeasuringAccessibilityExploration1997] note that using disaggregated data in accessibility analyses provides a more accurate estimate for individuals. Chen et al. [-@chenExploringEquityPerformance2019, pp. 530] are in favour of using disaggregated data, which they did in their recent analysis of Tampa's PBSP, because they note "the use of aggregated data might hinder our understanding of the equity impacts since individual disparities are absorbed after aggregation".

## 2.4 Previous Research

Our analysis builds upon a previous study [@hosfordWhoArePublic2018], which found that areas in Hamilton with less advantage are better served by the city's PBSP compared to other Canadian cities (i.e., Toronto, Vancouver, Montreal, and Ottawa-Gatineau) where areas that are less deprived have greater access. Hosford and Winters [-@hosfordWhoArePublic2018, pp. 47] acknowledge that "Hamilton stands out in that the lower income neighborhoods are located near the city center and wealthier neighborhoods are in the surrounding suburban areas". Therefore, the core service area for the PBSP in Hamilton by default covers more of the disadvantaged areas in the city. However, there is also a great deal of variation in income in the core service area because of the local university and increasing gentrification. Hosford and Winters [-@hosfordWhoArePublic2018] took a macro-level approach in their analysis by using dissemination areas _across_ the city as the unit of analysis. They did not focus specifically on the core service area and did not differentiate between conventional and equity stations.

Therefore, this paper also contributes to the cycling literature by expanding the analysis of Hosford and Winters [-@hosfordWhoArePublic2018] to assess equity in accessibility to HBS and explore the contribution of the equity stations added to Hamilton's PBSP. 

3. Case Study {#sec:study}
======================

## 3.1. Original System

The case study of this paper is the city of Hamilton, Ontario, Canada. Before June 2020, the program was known as Social Bicycles or SoBi Hamilton, but is now called Hamilton Bike Share (HBS). The core service area spans 40 sq.km of the city, although it was planned to be 20 sq.km [@hamiltonsobi2015], and roughly 138,000 people are within 30 minutes walking of a station [see Figure \ref{fig:hamilton-and-sobi-service-area}]. This represents roughly one fifth of the total population of the Hamilton Census Metropolitan Area according to the 2016 Canadian Census. The City of Hamilton undertook a large public engagement campaign to validate the locations of stations that had been selected and to crowdsource potential locations for additional locations [@hamiltonsobi2014]. Most of the potential locations that were suggested were in the east end of the core service area that lack stations or in neighborhoods not serviced by the PBSP. The program was enthusiastically welcomed in the city in 2015 - within three weeks of launching, 10,000 trips had been made [@hamiltonHamiltonBikeShare2015], however inadequate coverage given the size of the program's service area was identified as a problem early on, and transportation planners noted that the small size (i.e., supply of bicycle racks) and low quantity of stations would lead to challenges in balancing the system [@hamiltonsobi2015].

## 3.2. Equity Initiative

Hamilton Bike Share Inc., the non-profit operator, launched an equity program, Everyone Rides Initiative (ERI), in 2018 with the objective of reducing barriers that may prevent individuals from accessing the program. Additional bicycles and twelve "equity" stations were added to the core service area in more disadvantaged and under-serviced neighbourhoods. The equity program also offers subsidized memberships to individuals who identify as low income, and complements this service with cycle skills education. A comparable program can be found in Philadelphia [see @caspiBikesharingPhiladelphiaLowerincome2019].

## 3.3. Current System

As of June 2020, HBS has 900 bikes, 130 stations [see Figure \ref{fig:sobi-stations-in-hamilton}], and over 26,000 active memberships [@hamiltonHamiltonBikeShare2015]. The core service area remains 40 sq.km. The program has twelve "equity" stations and 118 "conventional" stations (i.e., that were not explicitly added to address inequities in the program). The City of Hamilton has positioned stations between 300 and 600 m apart [@scottRouteChoiceBike2021], but anticipates that they will service those living within a 250 m buffer [@hamiltonsobi2015]. The latter constitutes a normative statement: people ought to be able to access a station in less than 600 m if they live in the core service area, with most usage coming from 250 m around. However, it is not known how far people are actually willing to travel to reach a station. It would be reasonable to assume that people are willing to walk beyond this threshold to access other stations if the ones nearest them have no supply of bicycles.

```{r hamilton-and-sobi-service-area, echo=FALSE, message=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "The core service area of Hamilton Bike Share is outlined in blue. Hamilton Census Metropolitan Area is shown in grey."}
geo_context <- ggplot() + 
  geom_sf(data = hamilton_da_2016,
          color = NA,
          fill = "gray") + 
  geom_sf(data = sobi_service,
          fill = NA,
          colour = "blue") + 
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  #annotate("text", x = -80, y = 43.5, label = "Lake Ontario") +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(location = "tl", which_north = "true", 
        pad_x = unit(0.25, "in"), pad_y = unit(0.25, "in"),
        style = north_arrow_fancy_orienteering)

# Inset map
geo_context_inset <- get_stamenmap(bbox = c(left = -79,
                                    bottom = 43,
                                    right = -81,
                                    top = 44),
          maptype = "toner-lite", 
          crop = FALSE,
          zoom = 8)
# plot map
geo_context_inset <- ggmap(geo_context_inset)  +
    coord_sf(crs = st_crs(4326))  +
    #xlab("Longitude") + ylab("Latitude") +
  annotation_scale(location = "bl", width_hint = 0.4) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.border = element_rect(color = "black",
                                    fill = NA),
        axis.title = element_blank())

#geo_context + geo_context_inset
ggdraw() + 
  draw_plot(geo_context, 0.0, 0.0, 0.8, 1) + 
  draw_plot(geo_context_inset, 0.5, 0.4, 0.6, 0.6) +
  draw_plot_label(x = c(.25,.5), y = c(0.45,.6), 
                  label =  c("Hamilton Bikeshare \nCore Service Area",  "Lake Ontario"),
                  hjust = 0, 
                  size = 8) +
  draw_line(x = c(0.34, 0.4),
    y = c(0.43, 0.5),
    color = "blue", 
    size = 1,
    lineend = "round",
    linejoin = "mitre",
    arrow = arrow(length = unit(0.1, "inches")))
```

```{r sobi-stations-in-hamilton, echo=FALSE, fig.align = 'center', out.width = "100%", fig.cap = "The spatial distribution of bike share docking stations in Hamilton, Ontario. Conventional stations are in purple and equity (ERI) stations are in orange. The service area of Hamilton Bike Share is outlined in blue and the city's downtown core is outlined in dark green. Hamilton Census Metropolitan Area is shown in grey."}
ggplot() + 
  geom_sf(data = cropped_da,
          color = NA,
          fill = "gray") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          size = 1,
          color = "darkgreen") +
  geom_sf(data = sobi_hubs %>% 
            dplyr::filter(hub_status == "Active"), # Remove deactivated hubs
          aes(size = RACKS_AMOU, fill = hub_type),
          shape = 21) +
  scale_size(name = "Number of racks", 
             range = c(0.5, 3.5)) +
  scale_fill_manual(name = "Hub Type",
                    values = c("Conventional" = "purple", "ERI" = "orange")) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        legend.position = "bottom",
        axis.ticks = element_blank()) +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in"),) +
  annotation_north_arrow(location = "tl", which_north = "true", 
        pad_x = unit(0.25, "in"), pad_y = unit(0.20, "in"),
        style = north_arrow_fancy_orienteering)
```

4. Methods and Data {#sec:methods}
======================

## 4.1. Balanced Floating Catchment Area

The BFCA method was developed to address issues with demand and supply inflation that result from the overlapping catchment areas produced by earlier FCA methods [@paezDemandLevelService2019]. Páez et al. [-@paezDemandLevelService2019] adjusted the impedance weights so that both supply and demand are proportionally allocated. The result is a FCA method that balances the population and level of service by eliminating the over-counting of population and level of service that leads to distortions in demand and supply.

The first step in the FCA method is to allocate the population to be serviced by each docking station:

$$
P_j = {\sum_{i = 1}^{n} P_i{w_{ij}}}
$$

As seen in the equation above, the population allocated to station $j$ is the weighted sum of the population in the region; a spatial weight $w_{ij}$ represents the friction that the population at $i$ faces when reaching station $j$, and is usually given by a distance-decay function, so that each station is assumed to service only a segment of the population within a limited geographical range. The level of service of station $j$ in bicycle racks per person is the supply at each station (i.e., the maximum number of bicycle racks) divided by the population within the established catchment area:
$$
L_j = \frac {S_j}{P_j} = \frac {S_j}{{\sum_{i = 1}^{n} P_i{w_{ij}}}}
$$

In the second step, the accessibility of population unit $i$ is calculated as the weighted sum of the level of service of all stations that can be reached from there according to the spatial weights:
$$
A_i = {\sum_{j = 1}^{J} L_j{w_{ij}}} = {\sum_{j = 1}^{J} \frac {S_j{w_{ij}}}{\sum_{i = 1}^{n} P_i{w_{ij}}}}
$$

The balanced approach of Páez et al. [-@paezDemandLevelService2019] replaces the spatial weights with normalized versions as follows. In the first step, the population is weighted with:
$$
{w_{ij}^{i} = \frac {w_{ij}}{\sum_{j = 1}^{J} {w_{ij}}}}
$$

\noindent and in the second step, the level of service is weighted using:
$$
{w_{ij}^{j} = \frac {w_{ij}}{\sum_{i = 1}^{n} {w_{ij}}}}
$$
These weights satisfy the following properties:
$$
\sum_{j = 1}^{J} {w^i_{ij}} = 1
$$

\noindent and:
$$
\sum_{i = 1}^{n} {w^j_{ij}} = 1
$$

With these weights, accessibility can be calculated without risk of demand or supply inflation:
$$
A_i = {\sum_{j = 1}^{J} \frac {S_j{w^j_{ij}}}{\sum_{i = 1}^{n} P_i{w^i_{ij}}}}
$$

By allocating the population and level of service proportionally, this method preserves the values of the population and level of service since:
$$
\sum_{i=1}^n A_i = \sum_{j=1}^J L_j
$$

In fact, since the proportional allocation procedure means that any proportion of the population allocated to a station is never allocated to other stations, and conversely any level of service allocated to a population is never re-allocated elsewhere, this property is replicated for any level of aggregation. 

## 4.2. Pycnophylactic Interpolation

Since walking trips to a docking station likely happen at a level lower than even the smallest census geography, this requires a more disaggregated approach than the use of census geographies. For this reason, we implemented our analysis parting from small population cells to better reflect the friction of walking to a docking station, which is an important component of a bike share trip [@chenExploringEquityPerformance2019].

To obtain population at sub-census geography levels, we used pycnophylactic interpolation [@toblerSmoothPycnophylacticInterpolation1979]. We obtained population data from the 2016 Census of Canada for dissemination areas (DA), which is the smallest publicly available census geography in Canada. These zonal values of the population were interpolated to smaller polygons that are 50-by-50 m in size. Pycnophylactic interpolation involves smoothing out the population from each DA while preserving total volume. When interpolating the population at this high level of resolution, it is important to ensure that population numbers are not allocated to areas where people do not live in Hamilton (for example, to parks or large institutional buildings, etc.). To do so, we retrieved shapefiles for various geographic features from Open Hamilton. Next, we removed these features from the PBSP core service area and used pycnophylactic interpolation to disaggregate and reallocate population within the remaining area [see Figure \ref{fig:interpolated-population}].

```{r interpolated-population, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Interpolated population in Hamilton Bike Share's core service area (outlined in blue) and within 30 minutes of walking to the core service area. Each population cell is 50-by-50 m in size. The downtown area is outlined in dark green."}
ggplot() + 
  geom_sf(data = cropped_da,
          fill = "gray",
          color = NA) + 
  geom_tile(data = cbind(st_drop_geometry(population_50x50), 
                         st_coordinates(population_50x50)),
            aes(x = X, 
                y = Y, 
                fill = population)) +
  scale_fill_distiller(name = "Population",
                       palette = "YlOrRd",
                       direction = 1) +
  geom_sf(data = sobi_service,
          fill = NA,
          colour = "blue") +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank()) +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in"),) +
  annotation_north_arrow(location = "tl", which_north = "true", 
        pad_x = unit(0.25, "in"), pad_y = unit(0.20, "in"),
        style = north_arrow_fancy_orienteering)
```

## 4.3. Travel Time Matrix

To calculate walking times from the centroid of the population cells to docking stations, we extracted OpenStreetMap data for HBS's service area from [BBBike](https://download.bbbike.org/osm/bbbike/), an online cycle route planner that interfaces with OpenStreetMap. OpenStreetMap data provides the networks for calculating walking times from each population cell to nearby docking stations, using a maximum walking distance of 10 km and walking time of 30 minutes as thresholds. A travel time matrix was created with the origins as the coordinates of the population cells and the destinations as the coordinates of the docking stations within the maximum threshold. This process provides a more realistic measure of the friction of reaching stations by taking infrastructure into account network travel times, rather than using the Euclidean distance from population cell to station. Routing and travel time calculations were completed using the `R` package `r5r`, used for rapid realistic routing operations [@Pereira2021r5r].

## 4.4. Data

All data for this research were accessed from publicly available Census of Canada sources, from OpenStreetMaps, and from Open Hamilton^[https://open.hamilton.ca/], a public online repository of data curated by the City of Hamilton. Median total household income statistics were drawn from the 2016 Canadian Census.

5. Results
===================

## 5.1. Accessibility by Distance Thresholds

Consensus regarding the distance that individuals are willing to walk to access a docking station is lacking, but the literature on walking behaviour provided some guidelines to determine the thresholds for our sensitivity analysis. Previous studies have found that living within 250 m [@fullerUseNewPublic2011] and 300 m [@kabraBikeShareSystemsAccessibility2020] of a docking station is correlated with bike share use. Some studies have defined "close proximity" to a docking station as 500 metres or less [@fullerImpactEvaluationPublic2013; @hosfordEvaluationImpactPublic2018; @hosfordEvaluatingImpactImplementing2019], while other research has found that walking trips are less than 600 m and rarely more than 1200 m [@millwardActivetransportWalkingBehavior2013] or a median distance of 650 m [@larsenQuarterMileReexamining2010]. HBS will often depict a map at some docking stations to show the locations of the other nearest stations within a five minute walk, which suggests that this is an average distance that people are expected to walk. The National Association of City Transportation Officials (NACTO) has a similar normative guide [@nactowalkingstation2015]. 

In the present case, we found that accessibility calculated using the BFCA method increased with a threshold between two and four minutes, but was then maximized at five minutes. Accessibility decreased substantially after eight minutes, which is intuitive given that demand on a limited supply increases as more people can reach each station.

We experimented with various thresholds by conducting a sensitivity analysis to calculate accessibility at different walking times from population cell to docking stations: three minutes, five minutes, ten minutes, and fifteen minutes. We categorized these thresholds as minimum, average, maximum, and extreme, respectively. At each threshold, we compared accessibility between the current system and the original system to examine the contribution of the twelve equity stations. When considering the results reported below, it is important to remember that accessibility is technically a form of smoothing [@okelly2003aggregate, pp. 7-8]: smaller thresholds produce less smoothing (which can result in "spiky" accessibility landscapes), while larger thresholds produce more smoothing and fewer spikes.

In the analysis we use the following distance decay function: 
$$
{w_{ij} = \begin{cases} 
      1 & \text{when } t_{ij}\leq \gamma \\
      0 & \text{otherwise} 
   \end{cases}}
$$
\noindent where $\gamma$ is the relevant threshold. The weights are standardized as discussed above.

```{r join-population-to-travel-time-data, include=FALSE}
ttm_walk <- left_join(ttm_walk, 
                      population_50x50 %>%
                        st_drop_geometry(),
                      by = "UID")
```

```{r join-sobi-station-to-travel-time-data, include=FALSE}
ttm_walk <- ttm_walk %>% 
  left_join(sobi_hubs %>%
              st_drop_geometry() %>%
              dplyr::select(OBJECTID, 
                            RACKS_AMOU),
            by = c("OBJECTID" = "OBJECTID"))
```

```{r, include=FALSE}
names(ttm_walk) <- c("UID", 
                     "hub", 
                     "travel_time",
                     "hub_type",
                     "hub_status",
                     "population",
                     "racks")
```

```{r accessibility-sensibility, eval=FALSE, include=FALSE}
# Testing different thresholds to check how accessibility changes
accessibility_test <- data.frame(threshold = c(3:15), 
                                 accessibility = numeric(13))
for(i in 1:13){
  results_test <- ttm_walk %>%
    dplyr::filter(hub_type == "Conventional" &
                    hub_status == "Active") %>%
    b2sfca(threshold = i +2 )
  accessibility_test$accessibility[i] <- sum(results_test$accessibility$accessibility)
}

ggplot(accessibility_test, 
       aes(x = threshold, y = accessibility)) +
  geom_line()
```

<!-- Calculate accessibility for different thresholds -->

```{r 3-min-conventional, include=FALSE}
results_3min_conventional <- ttm_walk %>%
  dplyr::filter(hub_type == "Conventional" &
                  hub_status == "Active") %>%
  b2sfca(threshold = 3)
```

```{r 3-min-all, include=FALSE}
results_3min_all <- ttm_walk %>%
  dplyr::filter(hub_status == "Active") %>%
  b2sfca(threshold = 3)
```

```{r 3-min-sum-conventional, include=FALSE}
sum(results_3min_conventional$los$los)
sum(results_3min_conventional$accessibility$accessibility)
```

```{r 3-min-sum-all, include=FALSE}
sum(results_3min_all$los$los)
sum(results_3min_all$accessibility$accessibility)
```

```{r 5-min-conventional, include=FALSE}
results_5min_conventional <- ttm_walk %>%
  dplyr::filter(hub_type == "Conventional" &
                  hub_status == "Active") %>%
  b2sfca(threshold = 5)
```

```{r 5-min-all, include=FALSE}
results_5min_all <- ttm_walk %>%
  dplyr::filter(hub_status == "Active") %>%
  b2sfca(threshold = 5)
```

```{r, 5-min-sum-conventional, include=FALSE}
sum(results_5min_conventional$los$los)
sum(results_5min_conventional$accessibility$accessibility)
```

```{r 5-min-sum-all, include=FALSE}
sum(results_5min_all$los$los)
sum(results_5min_all$accessibility$accessibility)
```

```{r 10-min-conventional, include=FALSE}
results_10min_conventional <- ttm_walk %>%
  dplyr::filter(hub_type == "Conventional" &
                  hub_status == "Active") %>%
  b2sfca(threshold = 10)
```

```{r 10-min-all, include=FALSE}
results_10min_all <- ttm_walk %>%
  dplyr::filter(hub_status == "Active") %>%
  b2sfca(threshold = 10)
```

```{r 10-min-sum-conventional, include=FALSE}
sum(results_10min_conventional$los$los)
sum(results_10min_conventional$accessibility$accessibility)
```

```{r, 10-min-sum-all, include=FALSE}
sum(results_10min_all$los$los)
sum(results_10min_all$accessibility$accessibility)
```

```{r 15-min-conventional, include=FALSE}
results_15min_conventional <- ttm_walk %>%
  dplyr::filter(hub_type == "Conventional" &
                  hub_status == "Active") %>%
  b2sfca(threshold = 15)
```

```{r 15-min-all, include=FALSE}
results_15min_all <- ttm_walk %>%
  dplyr::filter(hub_status == "Active") %>%
  b2sfca(threshold = 15)
```

```{r 15-min-sum-conventional, include=FALSE}
sum(results_15min_conventional$los$los)
sum(results_15min_conventional$accessibility$accessibility)
```

```{r 15-min-sum-all, include=FALSE}
sum(results_15min_all$los$los)
sum(results_15min_all$accessibility$accessibility)
```

```{r, bind-all-results, include=FALSE}
los <- rbind(data.frame(results_3min_all$los, system = "All stations", threshold = "3 min"),
             data.frame(results_3min_conventional$los, system = "Conventional", threshold = "3 min"),
             data.frame(results_5min_all$los, system = "All stations", threshold = "5 min"),
             data.frame(results_5min_conventional$los, system = "Conventional", threshold = "5 min"),
             data.frame(results_10min_all$los, system = "All stations", threshold = "10 min"),
             data.frame(results_10min_conventional$los, system = "Conventional", threshold = "10 min"),
             data.frame(results_15min_all$los, system = "All stations", threshold = "15 min"),
             data.frame(results_15min_conventional$los, system = "Conventional", threshold = "15 min"))

accessibility <- rbind(data.frame(results_3min_all$accessibility, system = "All stations", threshold = "3 min"),
                       data.frame(results_3min_conventional$accessibility, system = "Conventional", threshold = "3 min"),
                       data.frame(results_5min_all$accessibility, system = "All stations", threshold = "5 min"),
                       data.frame(results_5min_conventional$accessibility, system = "Conventional", threshold = "5 min"),
                       data.frame(results_10min_all$accessibility, system = "All stations", threshold = "10 min"),
                       data.frame(results_10min_conventional$accessibility, system = "Conventional", threshold = "10 min"),
                       data.frame(results_15min_all$accessibility, system = "All stations", threshold = "15 min"),
                       data.frame(results_15min_conventional$accessibility, system = "Conventional", threshold = "15 min"))
```

```{r, add-geometry-to-results, include=FALSE}
los <- los %>%
  left_join(sobi_hubs %>%
              dplyr::select(OBJECTID),
            by = c("hub" = "OBJECTID")) %>%
  st_as_sf()

accessibility <- accessibility %>%
  left_join(population_50x50 %>%
              dplyr::select(UID),
            by = c("UID" = "UID")) %>%
  st_as_sf()
```

### 5.1.1. Minimum Threshold

With a walking distance of three minutes, we found that the total level of service is 25.2 bicycle racks per person system-wide in the original system configuration (i.e., without equity stations). The addition of equity stations increases this slightly to 25.4 bicycles per person. This total level of service is allocated to the population to obtain the levels of accessibility, which turn out to be relatively uniform overall, with the exception of two small areas where accessibility is slightly higher. This high level of system-wide accessibility occurs because the population that can reach a docking station when travel time is three minutes or less is very limited, and accessibility is strongly shaped by a few locations that concentrate population and stations. For this reason, accessibility is highly concentrated in small areas around those stations. The map is not shown, since it is less informative, but can be recreated using the source code and data.

```{r figure-6, eval = FALSE, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Accessibility at 3 minutes walk (minimum threshold) compared between current system with equity stations and the original system without equity stations."}
#accessibility plot for current system with equity stations at 3 minutes walk
accessibility_3min <- accessibility %>% 
  dplyr::filter(threshold == "3 min")

ggplot() +
  geom_sf(data = population_50x50,
          color = "gray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_3min), st_coordinates(accessibility_3min)),
            aes(x = X, y = Y, fill = accessibility)) +
  scale_fill_distiller(palette = "PuRd", direction = 1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ threshold)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in"),) +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering)
  
```

### 5.1.2. Medium Threshold

With a walking distance of five minutes, we found that there are 68.6 bicycle racks per person system-wide in the original system configuration (i.e., without equity stations). With the addition of equity stations, there are now 68.8 bicycle racks per person. At this threshold, there are more bicycle racks per person than at the minimum threshold. System-wide accessibility has in fact increased: the population that can reach the stations has grown, but not to the the point that congestion effects begin to take place. Figure \ref{fig:figure-7} presents a comparison of accessibility between the systems. Again, accessibility is fairly uniform, with the exception of one very small area. The equity stations noticeably increase accessibility in the east end of the core service area by filling gaps in PBSP coverage. 

```{r figure-7-a, include=FALSE}
#accessibility plot for current system with equity stations at 5 minutes walk
accessibility_5min <- accessibility %>% 
  dplyr::filter(threshold == "5 min",
                system == "All stations") %>%
  mutate(accessibility = cut(accessibility,
                             breaks = c(-Inf, 0.001, 0.02, 0.4, 8, Inf), 
                             labels=c("<0.001", "0.001-0.02", "0.02-0.4", "0.4-8", ">8")),
         accessibility = factor(accessibility,
                                levels = c(">8", "0.4-8", "0.02-0.4", "0.001-0.02", "<0.001"),
                                ordered = TRUE))


acc_5_all <- ggplot() +
  geom_sf(data = population_50x50,
          color = "gray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_5min), 
                         st_coordinates(accessibility_5min)),
            aes(x = X, y = Y, fill = accessibility)) +
  scale_fill_brewer(name = "Accessibility\n(racks/person)",
                    palette = "PuRd", 
                    direction = -1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ .)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering)

# Extract the legend
acc_5_legend <- cowplot::get_legend(acc_5_all)
```

```{r figure-7-b, include=FALSE}
#accessibility plot for current system with equity stations at 5 minutes walk
accessibility_5min <- accessibility %>% 
  dplyr::filter(threshold == "5 min",
                system == "Conventional") %>%
  mutate(accessibility = cut(accessibility,
                             breaks = c(-Inf, 0.001, 0.02, 0.4, 8, Inf), 
                             labels = c("<0.001", "0.001-0.02", "0.02-0.4", "0.4-8", ">8")),
         accessibility = factor(accessibility,
                                levels = c(">8", "0.4-8", "0.02-0.4", "0.001-0.02", "<0.001"),
                                ordered = TRUE))

acc_5_conv <- ggplot() +
  geom_sf(data = population_50x50,
          color = "gray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_5min),
                         st_coordinates(accessibility_5min)),
            aes(x = X, 
                y = Y, 
                fill = accessibility)) +
  scale_fill_brewer(name = "Accessibility\n(racks/person)",
                    palette = "PuRd", 
                    direction = -1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ .)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in"))
```

```{r figure-7, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Accessibility at 5 minutes walk (average threshold) compared between current system with equity stations and the original system without equity stations."}
ggdraw() + 
  draw_plot(acc_5_all + 
              theme(legend.position = "none"), 
            0.00, 0.5, 1, 0.5) + 
  draw_plot(acc_5_conv + 
              theme(legend.position = "none"), 0.00, 0.0, 1, 0.5) + 
  draw_plot(acc_5_legend, 0.01, 0.4, 0.2, 0.2)
```

### 5.1.3. High Threshold

With a walking distance of ten minutes, we found that there are 3.61 bicycle racks per person system-wide in the original system configuration (i.e., without equity stations). With the addition of equity stations, there are now 3.74 bicycles per person. Figure \ref{fig:figure-8} presents a comparison of accessibility between the systems. Differences in accessibility across the service area are now apparent, with users near the university and its adjacent neighborhoods, as well as neighborhoods north of the downtown area (the latter is outlined in green), having slightly higher accessibility. While the differences are modest, they are more apparent at this threshold than at shorter walking distances, especially in the more disadvantaged neighbourhoods in the east end of the core service area where equity stations were added.

```{r figure-8-a, include=FALSE}
#accessibility plot for current system with equity stations at 10 minutes walk
accessibility_10min <- accessibility %>% 
  dplyr::filter(threshold == "10 min",
                system == "All stations")  %>%
  mutate(accessibility = cut(accessibility,
                             breaks = c(-Inf, 0.00004, 0.0003, 0.002, Inf), 
                             labels = c("<0.00004", "0.00004-0.0003", "0.0003-0.002", ">0.002")),
         accessibility = factor(accessibility,
                                levels = c( ">0.002", "0.0003-0.002", "0.00004-0.0003", "<0.00004"),
                                ordered = TRUE))

acc_10_all <- ggplot() +
  geom_sf(data = population_50x50,
          color = "lightgray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_10min), 
                         st_coordinates(accessibility_10min)),
            aes(x = X, 
                y = Y, 
                fill = accessibility)) +
  scale_fill_brewer(name = "Accessibility\n(racks/person)",
                    palette = "PuRd",
                    direction = -1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ .)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering)

# Extract the legend
acc_10_legend <- cowplot::get_legend(acc_10_all)
```

```{r figure-8-b, include=FALSE}
#accessibility plot for current system without equity stations at 10 minutes walk
accessibility_10min <- accessibility %>% 
  dplyr::filter(threshold == "10 min",
                system == "Conventional")  %>%
  mutate(accessibility = cut(accessibility,
                             breaks = c(-Inf, 0.00004, 0.0003, 0.002, Inf), 
                             labels = c("<0.00004", "0.00004-0.0003", "0.0003-0.002", ">0.002")),
         accessibility = factor(accessibility,
                                levels = c( ">0.002", "0.0003-0.002", "0.00004-0.0003", "<0.00004"),
                                ordered = TRUE))

acc_10_conv <- ggplot() +
  geom_sf(data = population_50x50,
          color = "lightgray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_10min), 
                         st_coordinates(accessibility_10min)),
            aes(x = X, 
                y = Y, 
                fill = accessibility)) +
  scale_fill_brewer(name = "Accessibility\n(racks/person)",
                    palette = "PuRd",
                    direction = -1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ .)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in"))
```

```{r figure-8, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Accessibility at 10 minutes walk (maximum threshold) compared between current system with equity stations and the original system without equity stations."}
ggdraw() + 
  draw_plot(acc_10_all + 
              theme(legend.position = "none"), 
            0.00, 0.5, 1, 0.5) + 
  draw_plot(acc_10_conv + 
              theme(legend.position = "none"), 0.00, 0.0, 1, 0.5) + 
  draw_plot(acc_10_legend, 0.01, 0.4, 0.2, 0.2)
```

### 5.1.4. Extreme Threshold

With a walking distance of fifteen minutes, we found that there are 2.44 bicycle racks per person system-wide in the original system configuration (i.e., without equity stations). With the addition of equity stations, there are now 2.55 bicycles per person. Figure \ref{fig:figure-9} presents a comparison of accessibility between the systems. Users near the university and the neighborhoods north of the downtown area (the latter is outlined in green) have the highest accessibility, followed by those who live in the city's downtown area. Accessibility in the east end, where equity stations were implemented, of the core service area remains lower than other areas.

```{r figure-9-a, include=FALSE}
#accessibility plot for current system with equity stations at 15 minutes walk
accessibility_15min <- accessibility %>% 
  dplyr::filter(threshold == "15 min",
                system == "All stations")  %>%
  mutate(accessibility = cut(accessibility,
                             breaks = c(-Inf, 0.000006, 0.00005, 0.0003, Inf), 
                             labels = c("<0.000006", "0.000006-0.00005", "0.00005-0.0003", ">0.0003")),
         accessibility = factor(accessibility,
                                levels = c( ">0.0003", "0.00005-0.0003", "0.000006-0.00005", "<0.000006"),
                                ordered = TRUE))

acc_15_all <- ggplot() +
  geom_sf(data = population_50x50,
          color = "lightgray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_15min), st_coordinates(accessibility_15min)),
            aes(x = X, y = Y, fill = accessibility)) +
  scale_fill_brewer(name = "Accessibility\n(racks/person)",
                    palette = "PuRd", 
                    direction = -1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ .)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering)

# Extract the legend
acc_15_legend <- cowplot::get_legend(acc_15_all)
```

```{r figure-9-b, include=FALSE}
#accessibility plot for current system with equity stations at 15 minutes walk
accessibility_15min <- accessibility %>% 
  dplyr::filter(threshold == "15 min",
                system == "Conventional")  %>%
  mutate(accessibility = cut(accessibility,
                             breaks = c(-Inf, 0.000006, 0.00005, 0.0003, Inf), 
                             labels = c("<0.000006", "0.000006-0.00005", "0.00005-0.0003", ">0.0003")),
         accessibility = factor(accessibility,
                                levels = c( ">0.0003", "0.00005-0.0003", "0.000006-0.00005", "<0.000006"),
                                ordered = TRUE))

acc_15_conv <- ggplot() +
  geom_sf(data = population_50x50,
          color = "lightgray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_15min), st_coordinates(accessibility_15min)),
            aes(x = X, y = Y, fill = accessibility)) +
  scale_fill_brewer(name = "Accessibility\n(racks/person)",
                    palette = "PuRd", 
                    direction = -1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  facet_grid(system ~ .)  +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in"))
```

```{r figure-9, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Accessibility at 15 minutes walk (extreme threshold) compared between current system with equity stations and the original system without equity stations."}
ggdraw() + 
  draw_plot(acc_15_all + 
              theme(legend.position = "none"), 
            0.00, 0.5, 1, 0.5) + 
  draw_plot(acc_15_conv + 
              theme(legend.position = "none"), 0.00, 0.0, 1, 0.5) + 
  draw_plot(acc_15_legend, 0.01, 0.4, 0.2, 0.2)
```

## 5.2. Comparative Analysis

To highlight the benefit of using the BFCA approach, we conducted a comparative analysis with the 2SFCA method. The latter method does not allocate population and level of service proportionally, which means that there are often substantial amounts of multiple counting each. Figure \ref{fig:figure-10} shows the level of inflation and deflation at five minutes that results when the balanced and the conventional FCA methods are compared. The inflation and deflation is defined as the ratio of the level of service and accessibility calculated using the conventional approach to the balanced approach. This effect has been documented in other studies [@chenMeasuringAccessibilityHealth2021; @paezDemandLevelService2019]. Figure \ref{fig:figure-10} shows that in some locations the level of service is deflated by up to 80%. This is evident for the equity stations in the east end of the core service area. Noticeably, the level of service is least deflated in the downtown area, likely due to the density of docking stations. Figure \ref{fig:figure-11} shows that accessibility is inflated by as much as 50% in many areas across the city.

```{r conventional-5-min-non-equity, include=FALSE}
conventional_5min_non_equity <- ttm_walk %>%
  dplyr::filter(hub_type == "Conventional" &
                  hub_status == "Active") %>%
  c2sfca(threshold = 5)
```

```{r conventional-5-min-all, include=FALSE}
conventional_5min_all <- ttm_walk %>%
  dplyr::filter(hub_status == "Active") %>%
  c2sfca(threshold = 5)
```

```{r conventional-5-min-sum-non-equity, include=FALSE}
sum(conventional_5min_non_equity$los$los)
sum(conventional_5min_non_equity$accessibility$accessibility)
```

```{r conventional-5-min-sum-all, include=FALSE}
sum(conventional_5min_non_equity$los$los)
sum(conventional_5min_non_equity$accessibility$accessibility)
```

```{r balanced-5-min-non-equity, include=FALSE}
balanced_5min_non_equity <- ttm_walk %>%
  dplyr::filter(hub_type == "Conventional" &
                  hub_status == "Active") %>%
  b2sfca(threshold = 5)
```

```{r balanced-5-min-all, include=FALSE}
balanced_5min_all <- ttm_walk %>%
  dplyr::filter(hub_status == "Active") %>%
  b2sfca(threshold = 5)
```

```{r balanced-5-min-sum-non-equity, include=FALSE}
sum(balanced_5min_non_equity$los$los)
sum(balanced_5min_non_equity$accessibility$accessibility)
```

```{r balanced-5-min-sum-all, include=FALSE}
sum(conventional_5min_all$los$los)
sum(conventional_5min_all$accessibility$accessibility)
```

```{r los-inflation-5-min, include=FALSE}
los_inflation <- data.frame(hub = conventional_5min_non_equity$los$hub, 
                       los_inflation = conventional_5min_non_equity$los$los/balanced_5min_non_equity$los$los) %>% # Ratio of conventional los to balanced los
  left_join(sobi_hubs %>%
              dplyr::select(OBJECTID),
            by = c("hub" = "OBJECTID")) %>%
  st_as_sf()

accessibility_inflation <- data.frame(UID = conventional_5min_non_equity$accessibility$UID,
                            accessibility_inflation = conventional_5min_non_equity$accessibility$accessibility/balanced_5min_non_equity$accessibility$accessibility) %>%
  left_join(population_50x50 %>%
              dplyr::select(UID),
            by = c("UID" = "UID")) %>%
  st_as_sf()
```

```{r figure-10, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Deflation of level of service at 5 minutes walk (average threshold) in the current system with equity stations."}
ggplot() +
  geom_sf(data = population_50x50,
          color = "gray") +
  geom_sf(data = los_inflation,
            aes(size = los_inflation, color = los_inflation)) +
  scale_color_distiller(name = "Deflation of Level of Service",
                        palette = "PuRd", 
                        direction = 1) +
  guides(size = guide_legend(title = "")) +
  geom_sf(data = sobi_service,
          fill = NA) +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_scale(location = "bl", 
                   width_hint = 0.4,
                   pad_x = unit(0.4, "in"), 
                   pad_y = unit(0.2, "in"),) +
  annotation_north_arrow(location = "tr", 
                         which_north = "true", 
        pad_x = unit(0.0, "in"), 
        pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering)
```

```{r figure-11, echo=FALSE, fig.align = 'center', out.width = "90%", fig.cap = "Inflation of accessibility at 5 minutes walk (average threshold) in the current system with equity stations."}
ggplot() +
  geom_sf(data = population_50x50,
          color = "gray") +
  geom_tile(data = cbind(st_drop_geometry(accessibility_inflation), st_coordinates(accessibility_inflation)),
            aes(x = X, y = Y, fill = accessibility_inflation)) +
  scale_fill_distiller(name = "Inflation of accessibility",
                       palette = "PuRd", 
                       direction = 1) +
  geom_sf(data = sobi_service,
          fill = NA) +
  geom_sf(data = community_downtown %>%
            dplyr::filter(NAME == "Hamilton"),
          fill = NA,
          color = "dark green") +
  theme(panel.background = element_rect(fill = "white"),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())  +
  annotation_scale(location = "bl", 
                   width_hint = 0.4,
                   pad_x = unit(0.4, "in"), 
                   pad_y = unit(0.2, "in"),) +
  annotation_north_arrow(location = "tr", 
                         which_north = "true", 
        pad_x = unit(0.0, "in"), 
        pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering)
```

## 5.3. Accessibility by Median Total Household Income

To examine whether accessibility to HBS was increased for groups with lower income, we estimated what level of accessibility accrues to how many people in different income strata. To this end, we took our small population cells and aggregated the accessibility and population by the median total household income as imputed from the dissemination areas. A unique property of the BFCA method, which does not hold for the 2SFCA approach because of the inflation/deflation issues discussed above, is that accessibility and population data can be re-aggregated using income as an aggregator criterion, while preserving the total population and the level of service. This avoids demand and supply inflation, and also enables us to present findings in a way that is more intuitive to interpret. Figures \ref{fig:figure-bi-map-threshold-3}, \ref{fig:figure-bi-map-threshold-5}, \ref{fig:figure-bi-map-threshold-10}, and \ref{fig:figure-bi-map-threshold-15} depict bivariate choropleth maps that combine the spatial distribution of accessibility and median total household income, using tertiles for the coloring scheme.  

```{r, include=FALSE}
# drop geometry
test <- accessibility

# join accessibility to population cells
test <- test %>% 
  st_join(population_50x50 %>% 
            dplyr::select(-UID)) %>%
  rename(population50x50 = population)
```

```{r, include=FALSE}
test <- test %>%
  st_join(hamilton_da_2016)
```

```{r, include=FALSE}
agg_accessibility_da <- test %>%
  group_by(DA, 
           system,
           threshold) %>%
  summarize(accessibility = sum(accessibility),
            median_total_income = first(median_total_income), 
            .groups = "drop")
```

```{r figure-12, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "3 min", 
                system == "All stations") %>%
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_3_all <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("All stations") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  #annotation_scale(location = "bl", width_hint = 0.4,
  #                 pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))

legend <- bi_legend(pal = "GrPink",
                    dim = 3,
                    xlab = "Higher accessibility ",
                    ylab = "Higher income",
                    size = 6) + 
  theme(plot.background = element_rect(color = "lightgrey"))
```

```{r figure-13, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "3 min", 
                system == "Conventional") %>%
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_3_conv <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("Conventional") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  #annotation_north_arrow(location = "tr", which_north = "true", 
  #      pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
  #      style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-bi-map-threshold-3, echo = FALSE, out.width="1.2\\linewidth", fig.cap="\\label{fig-bivariate-map-threshold-3}Bivariate map of accessibility and income at the minimum threshold of three minutes with equity stations (top panel) and without equity stations (bottom panel)."}
ggdraw() + 
  draw_plot(bi_map_3_all, 0.00, 0.5, 1, 0.5) + 
  draw_plot(bi_map_3_conv, 0.00, 0.0, 1, 0.5) + 
  draw_plot(legend, 0.06, 0.4, 0.2, 0.2)
```

```{r figure-14, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "5 min", 
                system == "All stations") %>%
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_5_all <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("All stations") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  #annotation_scale(location = "bl", width_hint = 0.4,
  #                 pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-15, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "5 min", 
                system == "Conventional") %>%
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_5_conv <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("Conventional") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  #annotation_north_arrow(location = "tr", which_north = "true", 
  #      pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
  #      style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-bi-map-threshold-5, echo = FALSE, out.width="1.2\\linewidth", fig.cap="\\label{fig-bivariate-map-threshold-5}Bivariate map of accessibility and income at the average threshold of five minutes with equity stations (top panel) and without equity stations (bottom panel)."}
ggdraw() + draw_plot(bi_map_5_all, 0.00, 0.5, 1, 0.5) + 
  draw_plot(bi_map_5_conv, 0.00, 0.0, 1, 0.5) + 
  draw_plot(legend, 0.06, 0.4, 0.2, 0.2)
```

```{r figure-16, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "10 min", 
                system == "All stations") %>%
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_10_all <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("All stations") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  #annotation_scale(location = "bl", width_hint = 0.4,
  #                 pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-17, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "10 min", 
                system == "Conventional") %>%
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_10_conv <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("Conventional") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  #annotation_north_arrow(location = "tr", which_north = "true", 
  #      pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
  #      style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-bi-map-threshold-10, echo = FALSE, out.width="1.2\\linewidth", fig.cap="\\label{fig-bivariate-map-threshold-10}Bivariate map of accessibility and income at the maximum threshold of ten minutes with equity stations (top panel) and without equity stations (bottom panel)."}
ggdraw() + draw_plot(bi_map_10_all, 0.00, 0.5, 1, 0.5) + 
  draw_plot(bi_map_10_conv, 0.00, 0.0, 1, 0.5) + 
  draw_plot(legend, 0.06, 0.4, 0.2, 0.2)
```

```{r figure-18, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "15 min", 
                system == "All stations") %>%
  drop_na() %>% 
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_15_all <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("All stations") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  #annotation_scale(location = "bl", width_hint = 0.4,
  #                 pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  annotation_north_arrow(location = "tr", which_north = "true", 
        pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
        style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-19, echo=FALSE, fig.align = 'center', out.width = "90%"}
data_map <- agg_accessibility_da %>%
  dplyr::filter(threshold == "15 min", 
                system == "Conventional") %>%
  drop_na() %>% 
  bi_class(x = accessibility, 
           y = median_total_income, 
           style = "quantile", 
           dim = 3) %>%
  st_as_sf() 

bi_map_15_conv <- ggplot() +
  geom_sf(data = data_map %>%
            st_buffer(dist = 25, endCapStyle = "SQUARE"),
          mapping = aes(fill = bi_class), 
          color = NA,
          show.legend = FALSE) +
  geom_sf(data = cropped_da,
          fill = NA,
          color = "lightgrey") +
  geom_sf(data = sobi_service, 
          fill = NA,
          color = "blue") +
  ggtitle("Conventional") +
  bi_scale_fill(pal = "GrPink", 
                dim = 3) +
  annotation_scale(location = "bl", width_hint = 0.4,
                   pad_x = unit(0.4, "in"), pad_y = unit(0.2, "in")) +
  #annotation_north_arrow(location = "tr", which_north = "true", 
  #      pad_x = unit(0.0, "in"), pad_y = unit(0.0, "in"),
  #      style = north_arrow_fancy_orienteering) +
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        plot.title = element_text(hjust=0.5))
```

```{r figure-bi-map-threshold-15, echo = FALSE, out.width="1.2\\linewidth", fig.cap="\\label{fig-bivariate-map-threshold-15}Bivariate map of accessibility and income at the extreme threshold of fifteen with equity stations (top panel) and without equity stations (bottom panel)."}
ggdraw() + draw_plot(bi_map_15_all, 0.00, 0.5, 1, 0.5) + 
  draw_plot(bi_map_15_conv, 0.00, 0.0, 1, 0.5) + 
  draw_plot(legend, 0.06, 0.4, 0.2, 0.2)
```

```{r, include=FALSE}
# Find the income quantiles of the population served by SoBi
iq <- hamilton_da_2016 %>% 
  dplyr::filter(sobi_service == "Yes") %>%
  pull(median_total_income) %>%
  quantile(c(0, 0.2, 0.4, 0.6, 0.8, 1), na.rm = TRUE)
```

```{r population-core-service-by-income, include=FALSE}
# Find the population by income quantile
population_by_income <- population_50x50 %>%
  dplyr::select(-population) %>%
  st_join(hamilton_da_2016) %>%
  st_drop_geometry() %>%
  distinct(DA,
           .keep_all = TRUE) %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE)) %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population)) %>%
  drop_na()
```

```{r income-all-3, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_all_3 <- test %>% 
  dplyr::filter(system == "All stations", threshold == "3 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_all_3 <- agg_accessibility_all_3 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "All stations", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-con-3, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_con_3 <- test %>% 
  dplyr::filter(system == "Conventional", threshold == "3 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_con_3 <- agg_accessibility_con_3 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "Conventional", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-all-5, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_all_5 <- test %>% 
  dplyr::filter(system == "All stations", threshold == "5 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_all_5 <- agg_accessibility_all_5 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "All stations", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-con-5, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_con_5 <- test %>% 
  dplyr::filter(system == "Conventional", threshold == "5 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_con_5 <- agg_accessibility_con_5 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "Conventional", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-all-10, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_all_10 <- test %>% 
  dplyr::filter(system == "All stations", threshold == "10 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_all_10 <- agg_accessibility_all_10 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "All stations", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-con-10, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_con_10 <- test %>% 
  dplyr::filter(system == "Conventional", threshold == "10 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_con_10 <- agg_accessibility_con_10 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "Conventional", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-all-15, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_all_15 <- test %>% 
  dplyr::filter(system == "All stations", threshold == "15 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_all_15 <- agg_accessibility_all_15 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "All stations", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-con-15, include=FALSE}
# Add a variable to identify the income 
agg_accessibility_con_15 <- test %>% 
  dplyr::filter(system == "Conventional", threshold == "15 min") %>%
  mutate(income_quintile = case_when(median_total_income <= iq[2] ~ "Bottom 20%",
                                     median_total_income > iq[2] & median_total_income <= iq[3] ~ "Second 20%",
                                     median_total_income > iq[3] & median_total_income <= iq[4] ~ "Third 20%",
                                     median_total_income > iq[4] & median_total_income <= iq[5] ~ "Fourth 20%",
                                     median_total_income > iq[5] ~ "Top 20%"),
         income_quintile = factor(income_quintile, 
                                  levels = c("Bottom 20%", "Second 20%", "Third 20%", "Fourth 20%", "Top 20%"),
                                  ordered = TRUE))

inc_stats_con_15 <- agg_accessibility_con_15 %>%
  st_drop_geometry() %>%
  #dplyr::filter(system == "Conventional", threshold == "10 min") %>%
  group_by(income_quintile) %>%
  summarize(population = sum(population50x50),
            accessibility = sum(accessibility))
```

```{r income-all-hubs, include=FALSE}
all_hubs <- rbind(inc_stats_all_3,
                  inc_stats_all_5,
                  inc_stats_all_10,
                  inc_stats_all_15) %>%
  drop_na() # there is a very small number of population cells without income data

all_hubs
```

```{r income-conventional-hubs, include=FALSE}
conventional_hubs <- rbind(inc_stats_con_3,
                           inc_stats_con_5,
                           inc_stats_con_10,
                           inc_stats_con_15) %>%
  drop_na() # there is a very small number of population cells without income data

conventional_hubs
```

Table \ref{tab:accessibility-income} presents accessibility levels by income strata. As expected, the extreme threshold of fifteen minutes is associated with the largest number of people who are within the assumed service area of docking stations. We found that stations added to HBS expanded the spatial coverage of HBS, leading to more balance across the core service area. This is particularly evident at the minimum and average thresholds of three and five minutes, respectively, where the equity stations fill a number of gaps in program coverage. 

In line with Hosford and Winters' [-@hosfordWhoArePublic2018] findings, we found that accessibility to HBS was high for individuals in the bottom and second quintiles of median total household income who live in the centre of the core service area. This may well be an artifact of the spatial socioeconomic and demographic profile of Hamilton, where the most dense parts of the city (where a PBSP is most easily launched) are also those with relatively lower incomes [@hosfordWhoArePublic2018]. On the other hand, we found the levels of accessibility to Hamilton's PBSP are generally lower for populations in the bottom 20% of median total household income, compared to populations in the top 20%.

Our analysis revealed that the addition of equity stations increased accessibility to HBS by growing the population serviced irrespective of the walking threshold. The largest gains were made for dissemination areas in the second 20% median total household income, where an additional 3,073 and 5,395 people could reach a docking station within three and five minutes walk, respectively, after the addition of equity stations. However, we found that there were only small increases in the population in the bottom 20% of median total household income who are serviced by the equity stations, and the accessibility gains are also quite modest and smaller than for populations in the second and third quintiles of median total household income. This suggests that inequities in accessibility to HBS persist according to income, albeit this depends on the walking time thresholds. With and without equity stations, people in the top 20% of income have the highest level of access at a threshold of ten and fifteen minutes. Although dissemination areas in the second 20% have the highest level of access by a significant amount at lower distance thresholds, the bottom 20%, who may benefit more from increased access to HBS, have the lowest accessibility at three minutes threshold and the second lowest access at all other thresholds.

```{r accessibility-income, echo=FALSE}
accessibility_income <- cbind(
  conventional_hubs,
  all_hubs[, 2:3])

colnames(accessibility_income) <- c("iq", "population_conv", "accessibility_conv", "population_all", "accessibility_all")

accessibility_income <- accessibility_income %>%
  transmute(iq,
            population = rep(population_by_income$population, 4),
            population_conv = round(population_conv),
            accessibility_conv = round(accessibility_conv, 3), 
            population_all = round(population_all),
            accessibility_all = round(accessibility_all, 3),
            population_diff = population_all - population_conv, 
            accessibility_diff = accessibility_all - accessibility_conv)

kable(accessibility_income,
      "html",
      booktabs = TRUE,
      col.names = c("Income Quintile",
                    "Total Population",
                    "Population",
                    "Accessibility",
                    "Population",
                    "Accessibility",
                    "Population",
                    "Accessibility"),
      align = c("l", "c", "c", "c", "c", "c", "c"),
      caption = "\\label{tab:accessibility-by-income}Accessibility and population serviced by income quintile and between systems (with and without equity stations). Total population is the population by income quintile in the DAs that have any PBSP service at all.") %>%
  kable_styling(latex_options = c("striped", "scale_down")) %>%
  pack_rows("Threshold - 3 minutes", 1, 5, 
            label_row_css = "background-color: #8f9498; color: #fff;") %>%
  pack_rows("Threshold - 5 minutes", 6, 10, 
            label_row_css = "background-color: #8f9498; color: #fff;") %>%
  pack_rows("Threshold - 10 minutes", 11, 15, 
            label_row_css = "background-color: #8f9498; color: #fff;") %>%
  pack_rows("Threshold - 15 minutes", 16, 20, 
            label_row_css = "background-color: #8f9498; color: #fff;") %>%
  add_header_above(c(" " = 2, "Without Equity Stations" = 2, "With Equity Stations" = 2, "Difference" = 2)) %>%
  #column_spec(2, color = "white",
  #background = spec_color(accessibility_income$accessibility[1:40], end = 1)) %>%
  footnote(general = " ",
           alphabet = c("With equity stations = Hamilton Bike Share current system (118 conventional stations, 12 equity stations)", "Without equity stations = Hamilton Bike Share original system (118 conventional stations, no equity stations)"))
```

6. Discussion
===================

Using disaggregated population data, we examined equity in accessibility to Hamilton Bike Share (HBS), with a focus on assessing the contribution of the program's equity stations. The BFCA approach, combined with pycnophylactic interpolation, enabled us to measure accessibility on a micro scale which better reflects the scale at which walking takes place and avoids the "absorption of disparities", as articulated by Chen et al. [-@chenExploringEquityPerformance2019] and elsewhere in the transportation literature [@nationalacademiesofsciencesEffectiveMethodsEnvironmental2004; @rowangouldIdentifyingEnvironmentalJustice2016]. Our method also considers potential demand and supply, as well as congestion effects, in the calculation of accessibility, which is important for services like PBSPs but not accounted for in other common accessibility models. Unlike the 2SFCA method, there was no inflation or deflation in our accessibility estimates because population and level of service were preserved (see Figure \ref{fig:figure-10} and Figure \ref{fig:figure-11}). This differentiates our analysis from similar papers exploring equity in PBSPs that use larger geographical units of analysis or that focus only on station location instead of level of service. In this way, our paper has made an important contribution by applying an intuitive and useful approach not previously used in the cycling literature to measure accessibility to a PBSP. This open and reproducible research could be expanded for additional analysis using other individual-level data in Hamilton (e.g., age, gender, household size, education, etc.) or applied to PBSPs in other North American cities.

The sensitivity analysis revealed that accessibility to docking stations is maximized at five minutes and decreases significantly by eight minutes [see Table \ref{tab:accessibility-income}]. This reflects the normative guide advertised on some docking stations in Hamilton showing other stations within a five minute walk, as well as the directive of NACTO [@nactowalkingstation2015]. We found that over 118,000 people can access a bike share station within a 15 minute walk, which represents roughly 85% of the total population in the core service area [see Table \ref{tab:accessibility-income}]. At a minimum threshold of three minutes, too few people can reach stations which leads to relatively high levels of service since there is little crowding for the stations. However, accessibility is at its lowest after eight minutes whereby congestion effects due to increased potential demand kick in. Chen et al. [-@chenExploringEquityPerformance2019] also found that accessibility decreases as larger geographical units are used in the estimation, since more people can access the service.

The City of Hamilton has recognized from the launch of HBS that substantially more stations and bicycles are needed to service the area [@hamiltonsobi2015]. With a service area of 40 sq.km, it is estimated that Hamilton should have between 380 and 440 stations instead of 130, and 1,500 bicycles instead of 900 [@hamiltonsobi2015]. The City acknowledged that reduced capacity within the system led to gaps in coverage in some areas of the city "with some areas not having the recommended station density of 300m between stations or 10 stations per square km" [@hamiltonsobi2015]. Figures \ref{fig:figure-bi-map-threshold-3}, \ref{fig:figure-bi-map-threshold-5}, \ref{fig:figure-bi-map-threshold-10}, and \ref{fig:figure-bi-map-threshold-15} demonstrate how gaps in accessibility were filled by the additional equity stations implemented by Hamilton Bike Share. The improvement in accessibility is most noticeable at the average and maximum thresholds in the east end of the core service area, which corresponds to the neighbourhoods where stations were added. However, the contribution of equity stations to reducing inequities in accessibility was modest.

The equity stations were added to HBS as a targeted initiative to increase access for equity users in specific under-serviced neighbourhoods. However, we found that inequities persisted as evidenced by differences in accessibility according to income quintile (Table \ref{tab:accessibility-income}). While the addition of equity stations seems to modestly increase accessibility for all income groups at all thresholds, they did not increase accessibility substantially for any single income group. Most importantly, individuals in the bottom 20% of median total household income have the second lowest level of access to HBS at most thresholds (average, maximum, and extreme). At the minimum threshold, the bottom 20% have the lowest level of access. This suggests that additional equity stations are needed to more sufficiently address vertical equity. While previous research found that neighborhoods with more disadvantage are better serviced by HBS, the authors used the Pampalon Deprivation Index to determine the level of disadvantage for dissemination areas _across_ the city not just within the core service area [@hosfordWhoArePublic2018]. Instead, we used median total household income for each dissemination area _within_ the core service area. We conclude that Hamilton's PBSP, while by default located in areas with more deprivation compared to other cities as noted by Hosford and Winters [-@hosfordWhoArePublic2018], continues to have disparities in accessibility between income groups. With equity stations, many areas with low median total household income in the east end of the service area have zero or low accessibility at the average threshold of five minutes. At the maximum and extreme thresholds, the top 20% have the highest level of access to HBS. These findings align broadly with other studies from Tampa [@chenExploringEquityPerformance2019], Philadelphia [@caspiBikesharingPhiladelphiaLowerincome2019] and Seattle [@mooneyFreedomStationSpatial2019], which have found disparities in station location, annual trips, or access to bicycles, respectively, between levels of income and education. With the addition of equity stations, there were large gains in accessibility for the second 20% at the average threshold, amounting to over 5,000 more people, but much smaller gains for the bottom 20% with only 326 more people who are able to access a HBS station.

This paper highlights the importance of both the location and size of docking stations for increasing equity in accessibility to PBSPs. Congestion effects for services like bike share that offer a potential maximum supply of bicycles (i.e., number of racks) at each docking station must be taken into account when measuring accessibility. Several studies have examined who can access PBSPs based on the location of stations [among others, see @babagoliExploringHealthSpatial2019; @qianBikesharingEquityDisadvantaged2020; @qianEnhancingEquitableService2020; @wangNewBikeShare2019], but the size of docking stations in disadvantaged communities is not considered in the literature. An exception is a study on bike share ridership in New York City that found "the number of docks at a bike share station is statistically significant and positively associated with trip attractions for both men and women" [@wangGenderGapGenerators2019, p.5]. Station size was also not identified by PBSP owners and operators as a factor that needed to be addressed in equity efforts [@howlandCurrentEffortsMake2017]. PBSPs may not necessarily reduce inequities in access if stations are easy to reach for equity users but offer only a small number of bicycles. Likewise, more people may not opt to use the program if the supply of bicycles available at nearby stations is insufficient to meet demand. Babgoli et al. [-@babagoliExploringHealthSpatial2019] found a slight but not statistically significant increase in the proportion of neighborhoods with the highest levels of poverty that had stations after the Citi Bike expansion in 2015. Although the Citi Bike expansion was not specifically driven by a desire to reduce inequity in access, 16% of neighborhoods with the highest levels of poverty had stations compared to 12% before. This leaves open the question as to whether the stations held enough bicycles to meet demand, and whether congestion effects in neighbourhoods with more poverty are reducing the potential of bike share trips. 

The following policy and planning implications can inform efforts for increasing equity in HBS, and potentially other PBSPs in North America. Our analysis has identified specific areas that have both low accessibility and low median total household income. Figures \ref{fig:figure-bi-map-threshold-3}, \ref{fig:figure-bi-map-threshold-5}, \ref{fig:figure-bi-map-threshold-10}, and \ref{fig:figure-bi-map-threshold-15} highlight potential locations for new equity stations to increase access for low-income groups. Additional stations are needed to further reduce inequities in accessibility. This information will be useful to transportation planners in Hamilton and HBS operators. Our findings also suggest that PBSPs should consider increasing the size of stations in targeted communities. Notably, figure \ref{fig:sobi-stations-in-hamilton} highlights how the equity stations in Hamilton have a lower number of racks compared to docking stations in the city's downtown core. Docking stations added to increase vertical equity should have an adequate supply of public bicycles, particularly more in areas that have a strong interest in adopting bike share, to increase levels of service and potentially attract new users. Finally, distance to stations is associated with bike share use [@fullerUseNewPublic2011; @wangNewBikeShare2019], especially when the trip starts at the residence [@fishmanBikeShareSynthesis2013]. Distance can also be a barrier to using PBSPs [@fishmanBarriersBikesharingAnalysis2014]. To address this issue, transportation planners and PBSP operators should consult with groups and communities who are under-represented among bike share users to determine how far they are willing to walk to reach a docking station and to get their input on the location of new docking stations.

7. Study Limitations
===================

This paper did not examine or compare ridership data between conventional and equity stations. Therefore, further research is needed to determine whether the addition of equity stations encouraged more bicycling for equity users living near them. Other studies have specifically looked at differences in trip type, frequency, or length among users from disadvantaged neighborhoods [@caspiBikesharingPhiladelphiaLowerincome2019; @qianBikesharingEquityDisadvantaged2020; @wangNeighborhoodSociodemographicCharacteristics2019], but our analysis is limited by the lack of publicly available route and individual user data to conduct similar analyses for HBS. 

An additional limitation is the lack of publicly available information about the average number of bicycles at each station on a daily basis. HBS works to balance the number of bicycles across the core service area, but it is reasonable to expect that the number of bicycles will not match exactly the number of racks at every station. Ideally, instead of number of bicycle racks as our measure of supply, we would have liked to use the average number of bicycles at stations, perhaps at different times during the day across many seasons. Should this data become available, it would be worthwhile to revisit the study to examine how well the operation of the system (including balancing of bicycles across stations) works to maintain the nominal levels of accessibility examined in this paper.

\textcolor{blue}{Finally, it is important to note that this study examined accessibility to docking stations, that is, the entry point to the bikeshare system. There is a literature that examines transportation modes from this perspective, including accessibility to transit [e.g.,zhaoforecasting2003; @olszewskiusing2005; zieltracomparative2011] and more recently accessibility to bikeshare systems [e.g., @hosfordWhoArePublic2018; @mooneyFreedomStationSpatial2019; @caggianiequity2020; @duranfair2020]. These focus of these studies is to assess who has access to a transportation service. However, as argued by @moniruzzamanaccessibility2012, usage of the system likely depends as well on the accessibility to destinations that the system generates. An extension of the present study would be to investigate accessibility by the bikeshare system, and perhaps compare that to usage patterns.}

8. Conclusion
===================

The addition of specific stations to improve equity within the PBSP in Hamilton had the net effect of increasing accessibility and reducing vertical inequities to some extent. In particular, accessibility improved the most for those in the second 20% median total household income at all thresholds, but the gains were only modest for all income groups. Dissemination areas in the bottom 20% had the lowest accessibility at three minutes, and second lowest levels of accessibility at five, ten, and fifteen minutes. Congestion effects were observed at higher thresholds, with accessibility decreasing significantly once the catchment area is increased to ten minutes walking.  

Wang and Lindsey [-@wangNewBikeShare2019] have noted that there is a lack of research that examines how bike share users' behaviour changes as a result of program changes to station locations or improvements in accessibility. As such, a logical next step to this research is to examine whether the equity stations increased ridership or resulted in new memberships for equity users. An examination of the types of trips undertaken by residents in these areas would also be informative, such as the study undertaken by Caspi and Norland [-@caspiBikesharingPhiladelphiaLowerincome2019] after docking stations were implemented in low-income Philadelphia neighborhoods. The bulk of bicycle infrastructure that has been built in Hamilton to date is located in the core service area particularly near the conventional stations. It would be worthwhile to explore the route choice of bike share trips departing or ending at the equity stations and to identify factors that specifically influence trips from these stations, which would extend existing studies conducted by Scott and colleagues [@luUnderstandingBikeShare2018; @scottWhatFactorsInfluence2019; @scottRouteChoiceBike2021]. Perceived accessibility using bike share according to low-income neighbourhoods could also be studied in the future. Lättman et al. [-@lattmanDevelopmentTestPerceived2016] have developed the _Perceived Accessibility Scale_ that would be useful for this purpose. This paper, combined with additional studies such as those conceptualized above, would serve as a valuable case study for Hamilton and other cities with PBSPs that wish to evaluate and address inequities in accessibility to active transportation options in urban areas.

Acknowledgments
===================

We wish to thank the anonymous reviewers who provided constructive feedback to improve the quality of our paper. This research was completed using open software, and we also wish to acknowledge the developers of the following `R` packages: `biscale` [@R-biscale], `cowplot` [@R-cowplot], `data.table` [@R-data.table], `disk.frame` [@R-disk.frame], `gdistance` [@R-gdistance], `gridExtra` [@R-gridExtra], `kableExtra` [@R-kableExtra], `knitr` [@R-knitr], `pycno` [@R-pycno], `r5r` [@R-r5r], `raster` [@R-raster], `rgdal` [@R-rgdal], `rticles` [@R-rticles], `sf` [@R-sf], `tidyverse` [@R-tidyverse], `tinytex` [@R-tinytex],
`units` [@R-units].

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