Notation: R codes and a part of this brief are done with my teammate Anna Duan because we had a misunderstanding about how much we should collaborate.

Motivation

Transit Oriented Development (TOD) is a popular planning method which uses the relationship between transportation and land use to achieve city redevelopment. It’s commonly believed that TOD can promote business,increase the land value around transit and help form a more environmental-friendly lifestyle for people living in TOD area. However, it is also commonly related with gentrification because people with higher education level and income who prefer the city life and commuting by transit will be attracted to move back to TOD area and thus replace those city inhabitants who are economically-inferior.

Boston, Massachusetts has one of the country’s most developed public transportation systems developed in late 20th century. In the past decade, the Massachusetts Bay Transportation Authority (MBTA) has worked with private developers and the city to support more than 50 TOD projects near its stations[1]. These facts make Boston an interesting case for studying the effects of TOD.

knitr::opts_chunk$set(echo = TRUE)
options(scipen=999)

# LIBRARIES
#install.packages("tidyverse")          
library(tidyverse)
#install.packages("tidycensus")
library(tidycensus)
library(sf)
library(kableExtra)

options(scipen=999)
options(tigris_class = "sf")


#Plot Formatting
mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 16,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    strip.text.x = element_text(size = 10))
}
plotTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 10,colour = "black"),
    plot.subtitle = element_text(face="italic"),
    plot.caption = element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line("grey80", size = 0.1),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2),
    strip.background = element_rect(fill = "grey80", color = "white"),
    strip.text = element_text(size=10),
    axis.title = element_text(size=10),
    axis.text = element_text(size=10),
    plot.background = element_blank(),
    legend.background = element_blank(),
    legend.title = element_text(colour = "black", face = "italic"),
    legend.text = element_text(colour = "black", face = "italic"),
    strip.text.x = element_text(size = 10)
  )
}
# Quintile Map Breaks
qBr <- function(df, variable, rnd) {
  if (missing(rnd)) {
    as.character(quantile(round(df[[variable]],0),
                          c(.01,.2,.4,.6,.8), na.rm=T))
  } else if (rnd == FALSE | rnd == F) {
    as.character(formatC(quantile(df[[variable]]), digits = 3),
                 c(.01,.2,.4,.6,.8), na.rm=T)
  }
}
q5 <- function(variable) {as.factor(ntile(variable, 5))}
# Palette
palette5 <- c("#b3cde0","#6497b1","#005b96","#03396c","#011f4b")


#ACS DATA


# Input API key
census_api_key("0e3cc3910723434685f1e4c5df2ac519c0790ba5", install = TRUE, overwrite = TRUE)

# Variables: Median Rent, Median HH Income, Population, Bachelor's, No Vehicle (home owner, renter), Households (owner, renter-occupied)
# Projection: (NAD 1983 StatePlane Massachusetts Mainland FIPS 2001 Feet)
#2010
oritracts10 <-  
  get_acs(geography = "tract", variables = c("B25058_001E", "B19013_001E", "B01003_001E", "B23006_023E", 
                                             "B25044_003E", "B25044_010E", "B07013_002E", "B07013_003E"), 
          year=2010, state=25, county=025, geometry=T) %>% 
  st_transform('ESRI:102686')

oritracts18 <-  
  get_acs(geography = "tract", variables = c("B25058_001E", "B19013_001E", "B01003_001E", "B23006_023E", 
                                             "B25044_003E", "B25044_010E", "B07013_002E", "B07013_003E"), 
          year=2018, state=25, county=025, geometry=T) %>% 
  st_transform('ESRI:102686')
##############################Long form to wide form###############################
#2010
tracts10 <- 
  oritracts10 %>%
  dplyr::select( -NAME, -moe) %>%
  spread(variable, estimate) %>%
  dplyr::select(-geometry) %>%
  rename(Rent = B25058_001, 
         medHHInc = B19013_001,
         population = B01003_001, 
         bachelor = B23006_023,
         noVehicle_hmow = B25044_003, 
         noVehicle_hmre = B25044_010,
         Households_hmow = B07013_002,
         Households_hmre = B07013_003)

st_drop_geometry(tracts10)[1:3,]
#2018
tracts18 <- 
  oritracts18 %>%
  dplyr::select( -NAME, -moe) %>%
  spread(variable, estimate) %>%
  dplyr::select(-geometry) %>%
  rename(Rent = B25058_001, 
         medHHInc = B19013_001,
         population = B01003_001, 
         bachelor = B23006_023,
         noVehicle_hmow = B25044_003, 
         noVehicle_hmre  = B25044_010,
         Households_hmow = B07013_002,
         Households_hmre = B07013_003)

st_drop_geometry(tracts18)[1:3,]

###################################Mutate##########################################
#2010
tracts10 <- 
  tracts10 %>%
  mutate(pctBach = ifelse(population > 0, bachelor / population, 0),
         pctNoVehicle = ifelse(Households_hmow + Households_hmre > 0, 
                               (noVehicle_hmow + noVehicle_hmre) / 
                                  (Households_hmow + Households_hmre),0),
         year = "2010") %>%
  dplyr::select(-Households_hmow,-Households_hmre,-noVehicle_hmow,-noVehicle_hmre,-bachelor)
#2018
tracts18 <- 
  tracts18 %>%
  mutate(pctBach = ifelse(population > 0, bachelor / population, 0),
         pctNoVehicle = ifelse(Households_hmow + Households_hmre > 0, 
                               (noVehicle_hmow + noVehicle_hmre) / 
                                 (Households_hmow + Households_hmre),0),
         year = "2018") %>%
  dplyr::select(-Households_hmow,-Households_hmre,-noVehicle_hmow,-noVehicle_hmre,-bachelor)

###################################Bind 2010 and 2018#############################
allTracts <- rbind(tracts10,tracts18)

##################################Remove non-Boston tracts######################
bosTracts <- c("25025010405", "25025010404", "25025010801", "25025010702", "25025010204", 
"25025010802", "25025010104", "25025000703", "25025000504", "25025000704", "25025010103", 
"25025000803", "25025980300", "25025120201", "25025120104", "25025110607", "25025000302", 
"25025000301", "25025140400", "25025140300", "25025140201", "25025140202", "25025140102", 
"25025130402", "25025130300", "25025130200", "25025130100", "25025120700", "25025120600", 
"25025120500", "25025120400", "25025110601", "25025110502", "25025110501", "25025110401", 
"25025101102", "25025101101", "25025101002", "25025101001", "25025100900", "25025100800", 
"25025100601", "25025100500", "25025100400", "25025100300", "25025981300", "25025981201", 
"25025990101", "25025981501", "25025981700", "25025981800", "25025100200", "25025100100", 
"25025092400", "25025092300", "25025092200", "25025092000", "25025091900", "25025091800", 
"25025091700", "25025091600", "25025981100", "25025140105", "25025980700", "25025120105", 
"25025120301", "25025071201", "25025091001", "25025091500", "25025091400", "25025091300", 
"25025091200", "25025091100", "25025090700", "25025090600", "25025090400", "25025090300", 
"25025090200", "25025980101", "25025040801", "25025010203", "25025110403", "25025110201", 
"25025981000", "25025090100", "25025082100", "25025082000", "25025081900", "25025081800", 
"25025081700", "25025081500", "25025081400", "25025081300", "25025110103", "25025110301", 
"25025140106", "25025010701", "25025010408", "25025000503", "25025081200", "25025081100", 
"25025080900", "25025080801", "25025080601", "25025080500", "25025080401", "25025080300", 
"25025071101", "25025070900", "25025140107", "25025130404", "25025130406", "25025120103", 
"25025100700", "25025100603", "25025092101", "25025061101", "25025070800", "25025070700", 
"25025070600", "25025070500", "25025070300", "25025070200", "25025070101", "25025061200", 
"25025061000", "25025060800", "25025060700", "25025080100", "25025060301", "25025090901", 
"25025060101", "25025981502", "25025060600", "25025060400", "25025060200", "25025051200", 
"25025050700", "25025050600", "25025050500", "25025050400", "25025050300", "25025050200", 
"25025040300", "25025040200", "25025030500", "25025981600", "25025051101", "25025051000", 
"25025030400", "25025030300", "25025030200", "25025030100", "25025020200", "25025010600", 
"25025010500", "25025010300", "25025000802", "25025000701", "25025050101", "25025050901", 
"25025060501", "25025981202", "25025040100", "25025040600", "25025000602", "25025000601", 
"25025000502", "25025000402", "25025000401", "25025000202", "25025000201", "25025040401", 
"25025020303", "25025070402", "25025020302", "25025020301", "25025020101", "25025081001", 
"25025010403", "25025000100")

allTractsBos <- 
  subset(allTracts, GEOID %in% bosTracts)

#################################Transit Data####################################
mbtaNode <- st_read("E:/Upenn/CPLN508/TOD-Assignment/mbta_node.geojson") %>% st_transform(st_crs(allTractsBos)) 

#######################Exclude MBTA stops outside Boston#######################
#Filter out non-Boston stations
bosStations <-
  mbtaNode %>%
  filter(station != "Alewife" & station != "Assembly" & station != "Beachmont" & station !=  "Beaconsfield" & station !=  "Bellingham Square" & station !=  "Box District" & station !=  
         "Braintree" & station !=  "Brandon Hall" & station !=  "Brookline Village" & station !=  "Brookline Hills" & station !=  "Capen Street" & station !=  "Central" & station !=  
         "Central Avenue" & station !=  "Chelsea" & station !=  "Chestnut Hill" & station !=  "Coolidge Corner" & station !=  "Davis" & station !=  "Dean Road" & station !=  "Eastern Avenue" & station !=  
         "Eliot" & station !=  "Englewood Avenue" & station !=  "Fairbanks Street" & station !=  "Harvard" & station !=  "Hawes Street" & station !=  "Kendall/MIT" & station !=  "Kent Street" & station !=  "Longwood" & station !=  
         "Malden Center" & station !=  "Milton" & station !=  "Newton Centre" & station !=  "Newton Highlands" & station !=  "North Quincy" & station !=  "Oak Grove" & station !=  "Porter" & station !=  
         "Quincy Adams" & station !=  "Quincy Center" & station !=  "Revere Beach" & station !=  "Riverside" & station !=  "Saint Marys Street" & station !=  "Saint Paul Street" & station !=  
         "Summit Avenue" & station !=  "Tappan Street" & station !=  "Valley Road" & station !=  "Waban" & station !=  "Washington Square" & station !=  "Wellington" & station !=  "Wollaston" & station !=  
         "Wonderland" & station !=  "Woodland" & station !=  "Bellingham Sq"  & station !=  "Eastern Ave"  & station !=  "Brandon Hall" & station != "Summit Ave/Winchester St" & station != "Lechmere")

####################################Set TOD buffer##################################
bosBuffers <- 
  rbind(
    st_buffer(bosStations, 2640) %>% #in feet
      mutate(Legend = "Buffer") %>%
      dplyr::select(Legend),
    st_union(st_buffer(bosStations, 2640)) %>% #union buffer
      st_sf() %>%
      mutate(Legend = "Unioned Buffer"))

#############################Select by Centroids################################
buffer <- filter(bosBuffers, Legend=="Unioned Buffer")

selectCentroids <-
  st_centroid(allTractsBos)[buffer,] %>%
    st_drop_geometry() %>%
    left_join(dplyr::select(allTractsBos, GEOID)) %>%
    st_sf() %>%
    dplyr::select(year, population, Rent) %>%
    mutate(Selection_Type = "Select by Centroids")
#2010
selectCentroids10 <-
  st_centroid(tracts10)[buffer,] %>%
    st_drop_geometry() %>%
    left_join(dplyr::select(tracts10, GEOID)) %>%
    st_sf() %>%
    dplyr::select(year, population, Rent) %>%
    mutate(Selection_Type = "Select by Centroids")
#2018
selectCentroids18 <-
  st_centroid(tracts18)[buffer,] %>%
  st_drop_geometry() %>%
  left_join(dplyr::select(tracts18, GEOID)) %>%
  st_sf() %>%
  dplyr::select(year, population, Rent) %>%
  mutate(Selection_Type = "Select by Centroids")

allTracts.group <- 
  rbind(
    st_centroid(allTractsBos)[buffer,] %>%
      st_drop_geometry() %>%
      left_join(allTractsBos) %>%
      st_sf() %>%
      mutate(TOD = "TOD"),
    st_centroid(allTractsBos)[buffer, op = st_disjoint] %>%
      st_drop_geometry() %>%
      left_join(allTractsBos) %>%
      st_sf() %>%
      mutate(TOD = "Non-TOD")) %>%
  mutate(Rent.inf = ifelse(year == "2010", Rent * 1.1521, Rent))

1. Does TOD attract growth in population and rent?

If the hypothesis that transit-rich neighborhoods are desirable, it is then expected that neighborhoods close to a transit station should experience growth in population and rent as more people are drawn to them and willing to pay a premium to live near transit. The maps below document changes in rent and population in TOD tracts (within 0.5 miles of a transit station) between 2010 and 2018. In figure 1.1, we see some population growth in census tracts in the downtown area, near the intersections of the subway lines. Similarly, in figure 1.2, there are distinguish increases in TOD-area’s rent, with the most marked change occurring in the downtown area. Based on this, it appears that over time TOD neighborhoods in Boston are growing in population and experiencing increases in land value.

centectroi10_pop <- st_centroid(selectCentroids10[,1])
centectroi18_pop <- st_centroid(selectCentroids18[,1])
troi_pop <- rbind(selectCentroids10[,1:2], selectCentroids18[,1:2]) 
oitroi_pop.group <- st_centroid(troi_pop)
centectroi10_re <- st_centroid(selectCentroids10[,2])
centectroi18_re <- st_centroid(selectCentroids18[,2])
troi_rent <- rbind(selectCentroids10[,c(1,3)], selectCentroids18[,c(1,3)]) 
oitroi_rent.group <- st_centroid(troi_rent)
#population
ggplot() + 
  geom_sf(data=st_union(allTracts.group)) +
  geom_sf(data = buffer, fill = "white") + 
  geom_sf(data = oitroi_pop.group, aes(size = population), shape = 21, 
      fill = "lightblue", alpha = 0.8, show.legend = "point") + 
  scale_size_continuous(range = c(0.05, 4))+
  facet_wrap(~year)+
  labs(title = "TOD population change in 2010 and 2018", subtitle = "Boston, MA", caption="Figure 1.1") +
  mapTheme() + 
  theme(plot.title = element_text(size=22))

#rent
ggplot() + 
  geom_sf(data=st_union(allTractsBos)) +
  geom_sf(data = buffer, fill = "white") + 
  geom_sf(data = oitroi_rent.group, aes(size = Rent), shape = 21, color = "transparent",
          fill = "pink", alpha = 1, show.legend = "point") + 
  scale_size_continuous(range = c(0.05, 4))+
  facet_wrap(~year)+
  labs(title = "TOD land value change in 2010 and 2018", subtitle = "Boston, MA", caption="Figure 1.2") +
  mapTheme() + 
  theme(plot.title = element_text(size=22))

2. How does TOD affect other neighborhood characteristics?

In this section, other indicators are added to our analysis to measure the changes of demography, income, education and mode preference of people living in TOD and non-TOD neighborhoods from 2010 to 2018.

2.1 Change in Population

The majority of the tracts with highest population (the 5th quintile) are allocated within or near the downtown area, where all the subway lines intersect. This area also experienced a remarkable population growth from 2010 to 2018, while those TOD tracts further from downtown experienced less growth than the south non-TOD area. This indicates that TOD may not be as important as downtown development in driving Boston’s population growth.

ggplot(allTracts.group)+
  geom_sf(data = st_union(allTractsBos))+
  geom_sf(aes(fill = q5(population))) +
  geom_sf(data = buffer, fill = "transparent", color = "red")+
  scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.group, "population"),
                    name = "Population\n(Quintile Breaks, TOD in Red)") +
  labs(title = "Total Population in 2010 and 2018", subtitle = "Boston, MA", caption="Figure 2.1") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))

2.2 Change in Share of Adults with Bachelor’s Degrees or Higher

Tracts that have the highest share of adults with Bachelor’s Degrees or higher are concentrated in TOD area as well as its close vicinity. With regards to change, some non-TOD tracts saw decline and nearly all TOD tracts saw improvement between 2010 and 2018. This result suggests that TOD area may have a strong appeal for more civilized population compared to non-TOD area, which may result from easy access to transit, but may also because of the multiplier effect brought by downtown development.

ggplot(allTracts.group)+
  geom_sf(data = st_union(allTractsBos))+
  geom_sf(aes(fill = q5(pctBach))) +
  geom_sf(data = buffer, fill = "transparent", color = "red")+
  scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.group %>% mutate(pctBach_100 = pctBach *100), "pctBach_100"),
                    name = "% >25 Years with Bachelor's\n(Quintile Breaks, TOD in Red)") +
  labs(title = "Share of Bachelor's Degree or higher in 2010 and 2018", subtitle = "Boston, MA", caption="Figure 2.2") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))

2.3 Change in Median Household Income

From 2010 to 2018, the whole Boston area experienced increase in household income. Yet, TOD neighborhoods tend to have a higher growth than non-TOD ones, and some of the most significant growth happened in downtown area. This growth could mean that TOD areas attract wealthier households, although this growth could also be the result of people making more money due to the rich economic activities and opportunities in downtown as well as other TOD area.

ggplot(allTracts.group)+
  geom_sf(data = st_union(allTractsBos))+
  geom_sf(aes(fill = q5(medHHInc))) +
  geom_sf(data = buffer, fill = "transparent", color = "red")+
  scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.group, "medHHInc"),
                    name = "Median HH Income\n(Quintile Breaks, TOD in Red)") +
  labs(title = "Median Household Income in 2010 and 2018", subtitle = "Boston, MA (adjusted US Dollars)", caption="Figure 2.3") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))

2.4 Chanege in share of Households without Vehicles

Due to the easy accessibility to public transportation, it’s expected that TOD neighborhoods have a higher without-vehicle share. However, according to Figure 2.4, this effect is most apparent in downtown and its vicinity, some TOD neighbors even had a decrease for this indicator. This indicates that while TOD may help people shift from driving to more environmentally-friendly travel modes, the density of transit lines may be a critical factor in this shifting.

ggplot(allTracts.group)+
  geom_sf(data = st_union(allTractsBos))+
  geom_sf(aes(fill = q5(pctNoVehicle))) +
  geom_sf(data = buffer, fill = "transparent", color = "red")+
  scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.group %>% mutate(pctNoVehicle_100 = pctNoVehicle *100), "pctNoVehicle_100"),
                    name = "% Households\n(Quintile Breaks, TOD in Red)") +
  labs(title = "Share of Households without Vehicle 2010 - 2018", subtitle = "Boston, MA", caption="Figure 2.4") +
  facet_wrap(~year)+
  mapTheme() + 
  theme(plot.title = element_text(size=22))

## 2.5 Comprehensive Presentation Below, the previous trends are summarized in table 2.1. Although for most variables, TOD area tends to have a higher value than the non-TOD area, it’s noted that the growth effect brought by TOD in land value, population and population’s educational level is not distinct.

#Summarize allTracts.group
allTracts.Summary <- 
  st_drop_geometry(allTracts.group) %>%
  group_by(year, TOD) %>%
  summarize(Rent = mean(Rent, na.rm = T),
            Population = mean(population, na.rm = T),
            pctBach = mean(pctBach, na.rm = T),
            pctNoVehicle = mean(pctNoVehicle, na.rm = T),
            MedHHInc = mean(medHHInc, na.rm = T))

allTracts.Summary %>%
  unite(year.TOD, year, TOD, sep = ": ", remove = T) %>%
  gather(Variable, Value, -year.TOD) %>%
  mutate(Value = round(Value, 2)) %>%
  spread(year.TOD, Value) %>%
  kable() %>%
  kable_styling() %>%
  footnote(general_title = "\n",
           general = "Table 2.1")
Variable 2010: Non-TOD 2010: TOD 2018: Non-TOD 2018: TOD
MedHHInc 56146.43 52550.29 69593.97 70423.53
pctBach 0.19 0.31 0.23 0.35
pctNoVehicle 0.09 0.22 0.09 0.20
Population 3238.05 3398.96 3873.59 3668.81
Rent 964.49 1124.75 1293.12 1475.46

Table 2.1

By plotting the data in Table 2.1, it’s clear that: First, other than the percent of household without Vehicle, all variables have increased citywide from 2010 to 2018n for reasons in session 2.2. Second, population is lower in TOD neighborhoods in 2018, and these neighborhoods see less population growth on average than non-TOD areas, contradicting the hypothesis that transit-rich neighborhoods see the most population growth. Third, the comparison of the Percent of households without Vehicle between TOD and non-TOD neighborhoods is most notable, showing that people living in TOD are less likely to own a car.

###########################TOD Indicator Plots##################################
allTracts.Summary %>%
  gather(Variable, Value, -year, -TOD) %>%
  ggplot(aes(year, Value, fill = TOD)) +
    geom_bar(stat = "identity", position = "dodge") +
    facet_wrap(~Variable, scales = "free", ncol=5) +
    scale_fill_manual(values = c("#bae4bc", "#0868ac")) +
    labs(title = "Indicator differences across time and space", caption="Figure 2.5") +
    plotTheme() + theme(legend.position="bottom")

3. How does the results vary among subway lines and downtown?

Figure 3.1 shows the six sub-markets within the Boston subway system.

##########################Examining six submarkets ###########################
downtown <-
  st_intersection(
    st_buffer(filter(bosStations, line == "ORANGE/RED"), 2640) %>% st_union()) %>%
  st_sf() %>%
  mutate(Submarket = "downtown")

blue <-
    st_buffer(filter(bosStations, line %in% c("BLUE", "BLUE/ORANGE", "BLUE/GREEN")), 2640) %>% st_union() %>%
  st_sf() %>%
  st_difference(downtown) %>%
  mutate(Submarket = "blue")

orange <-
  st_buffer(filter(bosStations, line %in% c("ORANGE", "BLUE/ORANGE", "GREEN/ORANGE", "ORANGE/RED")), 2640) %>% st_union() %>%
  st_sf() %>%
  st_difference(downtown) %>%
  mutate(Submarket = "orange")

red <-
  st_buffer(filter(bosStations, line %in% c("RED", "ORANGE/RED", "GREEN/RED")), 2640) %>% st_union() %>%
  st_sf() %>%
  st_difference(downtown) %>%
  mutate(Submarket = "red")

green <-
  st_buffer(filter(bosStations, line %in% c("GREEN", "BLUE/GREEN", "GREEN/ORANGE", "GREEN/RED")), 2640) %>% st_union() %>%
  st_sf() %>%
  st_difference(downtown) %>%
  mutate(Submarket = "green")


sixMarkets <- rbind(red, green, orange, blue, downtown)

ggplot() + 
  geom_sf(data=st_union(allTractsBos)) +
  geom_sf(data=sixMarkets, aes(colour = Submarket), fill = "transparent", show.legend = "point", size= 1) +
  geom_sf(data=downtown, aes(colour = Submarket), fill = "transparent", show.legend = "point", size= 1) +
  scale_colour_manual(values = c("blue","black","green","orange","red")) +
  labs(title="Six submarkets", subtitle="Boston, MA", caption="Figure 3.1") +
  mapTheme()

Figure 3.2 and 3.3 show a clear spatial distinction of rent and population among submarkets. If we take a look at Table 3.1 and Figure 3.4, it’s clearly indicated that for all the variables, downtown area possesses a more distinct effect, compared to submarkets that do not or have fewer intersections. However, when downtown factor is controlled, it’s still noticeable that differences exist between TOD and non-TOD area for variables other than population and median household income.

# Bind buffers to tracts and map them or make small multiple plots

allTracts.sixMarkets <-
  st_join(st_centroid(allTractsBos), sixMarkets) %>%
  st_drop_geometry() %>%
  left_join(allTractsBos) %>%
  mutate(Submarket = replace_na(Submarket, "Non-TOD")) %>%
  st_sf() 

ggplot() +
  geom_sf(data=st_union(allTracts.sixMarkets)) +
  geom_sf(data = allTracts.sixMarkets, aes(fill = q5(Rent))) +
  facet_wrap(~Submarket + year) + 
  scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.sixMarkets, "Rent"),
                    name = "Rent\n(Quintile Breaks)") +
  labs(title="Rent variance among six submarkets", subtitle="Boston, MA", caption="Figure 3.2") +
  mapTheme() + theme(plot.title = element_text(size=18))

ggplot() +
  geom_sf(data=st_union(allTracts.sixMarkets)) +
  geom_sf(data = allTracts.sixMarkets, aes(fill = q5(population))) +
  facet_wrap(~Submarket + year) + 
  scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.sixMarkets, "population"),
                    name = "population\n(Quintile Breaks)") +
  labs(title="Population variance among six submarkets", subtitle="Boston, MA", caption="Figure 3.3") +
  mapTheme() + theme(plot.title = element_text(size=18))

###########################INDICATOR MAPS#######################################
#$100 in 2010 is equivalent in purchasing power to about $115.21 in 2018
#Select tracts within and outside of 0.5mile buffer
allTracts.group <- 
  rbind(
    st_centroid(allTractsBos)[buffer,] %>%
      st_drop_geometry() %>%
      left_join(allTractsBos) %>%
      st_sf() %>%
      mutate(TOD = "TOD"),
    st_centroid(allTractsBos)[buffer, op = st_disjoint] %>%
      st_drop_geometry() %>%
      left_join(allTractsBos) %>%
      st_sf() %>%
      mutate(TOD = "Non-TOD")) %>%
  mutate(Rent.inf = ifelse(year == "2010", Rent * 1.1521, Rent))
#spread goes long to wide, gather opposite
allTracts.sixMarkets.Summary <- 
  st_drop_geometry(allTracts.sixMarkets) %>%
  group_by(year, Submarket) %>%
  summarize(Rent = mean(Rent, na.rm = T),
            Population = mean(population, na.rm = T),
            pctBach = mean(pctBach, na.rm = T),
            pctNoVehicle = mean(pctNoVehicle, na.rm = T),
            MedHHInc = mean(medHHInc, na.rm = T))

#kable(allTracts.sixMarkets.Summary) %>%
  #kable_styling() %>%
  #footnote(general_title = "\n",
           #general = "Table 3.1")

########change to long form######
allTracts.sixMarkets.Summary %>%
  unite(year.Submarket, year, Submarket, sep = ": ", remove = T) %>%
  gather(Variable, Value, -year.Submarket) %>%
  mutate(Value = round(Value, 2)) %>%
  spread(year.Submarket, Value) %>%
  kable() %>%
  kable_styling() %>%
  footnote(general_title = "\n",
           general = "Table 3.1")
Variable 2010: blue 2010: downtown 2010: green 2010: Non-TOD 2010: orange 2010: red 2018: blue 2018: downtown 2018: green 2018: Non-TOD 2018: orange 2018: red
MedHHInc 50144.20 64431.60 53228.77 55361.51 57709.79 44483.70 72428.53 92884.00 71258.67 69130.29 79781.05 63752.60
pctBach 0.24 0.38 0.37 0.20 0.36 0.21 0.31 0.39 0.41 0.23 0.42 0.29
pctNoVehicle 0.24 0.28 0.29 0.10 0.23 0.16 0.19 0.24 0.26 0.09 0.23 0.16
Population 2911.00 3335.33 3478.79 3125.54 3300.18 3575.50 3124.53 3997.17 3617.31 3718.46 3585.31 4060.05
Rent 1116.93 1407.40 1257.43 967.27 1146.95 966.30 1538.67 1924.80 1643.75 1294.68 1509.31 1286.15

Table 3.1 
###########################TOD Indicator Plots##################################

allTracts.sixMarkets.Summary %>%
  gather(Variable, Value, -year, -Submarket) %>%
  ggplot(aes(year, Value, fill = Submarket)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~Variable, scales = "free", ncol=5) +
  scale_fill_manual(values = c("#335BFF", "#000000", "#34D539", "#F9FD0A", "#FF7000", "#FF0000", "#CDCDCD")) +
  labs(title = "Indicator variance across time and space", subtitle="Boston, MA", caption="Figure 3.4") +
  plotTheme() + theme(legend.position="bottom") +
  theme(plot.title = element_text(size=18))

4. How does Rent Vary with Distance from Transit?

To further analyze the relationship between TOD and neighborhood desirability, we will look at the relationship between rent and distance from subway stations. In both 2010 and 2018, the highest rent is at the shortest distance from transit stations, within 0.5miles, and it drops until 1 mile from stations. After this, however, both plots display two rent price peaks at roughly 1.5 and 3.2 miles. This aligns with the previous maps and group plots, suggesting that while transit rich neighborhoods have the highest rent prices due to their desirability, transit is not the dominant factor influencing rent prices, particularly outside of the 0.5 mile TOD buffer.

# 2010
rent_dis10 <- data.frame(
  distance = double(),
  Mean_Rent = double(),
  stringsAsFactors=FALSE
)

for (i in seq(0.4,4, by = 0.4)){
  
  A <- st_union(st_buffer(bosStations, i*5280)) %>% st_sf()
  a <- st_union(st_buffer(bosStations, (i-0.4)*5280)) %>% st_sf()
  std <- st_difference(A, a)
  b <-
    st_centroid(tracts10)[std,] %>%
    st_drop_geometry() %>%
    left_join(dplyr::select(tracts10, GEOID)) %>%
    st_sf() %>%
    dplyr::select(Rent)
  b$Rent <- as.numeric(b$Rent)
  c <- mean(b[["Rent"]],na.rm=TRUE)
  d <- data.frame(distance=i, Mean_Rent=c, stringsAsFactors=FALSE)
  rent_dis10 <- rbind(rent_dis10, d)
  i <- i + 0.4
}
rent_dis10$year=2010

# for 2018
rent_dis18 <- data.frame(
  distance = double(),
  Mean_Rent = double(),
  stringsAsFactors=FALSE
)

for (i in seq(0.4,4, by = 0.4)){
  
  A <- st_union(st_buffer(bosStations, i*5280)) %>% st_sf()
  a <- st_union(st_buffer(bosStations, (i-0.4)*5280)) %>% st_sf()
  std <- st_difference(A, a)
  b <-
    st_centroid(tracts18)[std,] %>%
    st_drop_geometry() %>%
    left_join(dplyr::select(tracts18, GEOID)) %>%
    st_sf() %>%
    dplyr::select(Rent)
  b$Rent <- as.numeric(b$Rent)
  c <- mean(b[["Rent"]],na.rm=TRUE)
  d <- data.frame(distance=i, Mean_Rent=c, year=2018, stringsAsFactors=FALSE)
  rent_dis18 <- rbind(rent_dis18, d)
  i <- i + 0.4
}
rent_dis18$year=2018

###combine all together and map
rent_dis <- rbind (rent_dis10, rent_dis18)

ggplot(data=rent_dis, aes(x=distance, y=Mean_Rent, group=as.factor(year), shape=as.factor(year), colour=as.factor(year))) +
  geom_line(aes(linetype=as.factor(year)), size=1)+
  geom_point(size=5)+
  scale_colour_hue(name="Year", l=40)+
  scale_shape_manual(name="Year", values=c(21,22))+
  scale_linetype_discrete(name="Year")+
  xlab("Distance") + ylab("Average rent")+
  labs(title = "Rent as a function of distance to subway stations", caption="Figure 4.1", subtitle="Boston, MA (in US Dollars and Miles)") +
  theme_bw()

5. Are TOD Neighborhoods Safer?

Another proposed benefit of TOD is neighborhood safety. With denser, walkable, and transit-rich neighborhoods, there are more people on the street. This is based in the notion of “eyes on the street,” an idea termed by the great urbanist Jane Jacobs, which suggests that dense, lively neighborhoods where residents are out on the streets are safer because of the many people who witness the activity in the area. To test this idea, we plot the incidences of a common street crime category, robbery.

At a glance, it appears that TOD tracts experience a higher rate of robberies, disproving the “eyes on the street” idea. However, it is important to consider another factor that may affect crime: rent. Research has suggested that lower rent is associated with higher crime rates. On the map it appears that TOD and non-TOD tracts that have rents in the lower categories have similar amounts of crime. However, it appears that many TOD tracts in the highest rent category do not have incidences of robberies, while most of the high-rent non-TOD tracts have some robberies. Also, another confounding factor is that most TOD tracts are in the downtown area, which can often have the highest crime in a city. Overall, TOD tracts in Boston don’t appear to be safer than non-TOD tracts due to the complicated factors influencing crime. However, it appears that neighborhoods which have both high rent and TOD experience less robbery crime.

crime12 <- read.csv("E:/Upenn/CPLN508/TOD-Assignment/Crime12.csv", header = TRUE, sep = ",", quote = "\"",
                    dec = ".")
crime18 <- read.csv("E:/Upenn/CPLN508/TOD-Assignment/Crime18.csv", header = TRUE, sep = ",", quote = "\"",
                    dec = ".")

crime12.sf <- st_as_sf(crime12, coords = c("LON", "LAT"), crs = 4326, agr = "constant") %>% 
  st_transform(st_crs(allTracts.group))
crime18.sf <- st_as_sf(crime18, coords = c("LON", "LAT"), crs = 4326, agr = "constant")  %>% 
  st_transform(st_crs(allTracts.group))

################################Relate Crime to Census Tracts##############################

crime.sf <- rbind(crime12.sf, crime18.sf)

ggplot() +   
#  geom_sf(data=st_union(allTractsBos)) +
 scale_fill_manual(values = palette5,
                    labels = qBr(allTracts.group, "Rent.inf"),
                    name = "Rent\n(Quintile Breaks)") +
  geom_sf(data=allTracts.group, aes(fill = q5(Rent.inf))) +
  geom_sf(data=crime.sf, show.legend = "point", size= 0.1, color = "white") +
  labs(title="Crime, Transit Access, and Rent, 2012 and 2018", subtitle="Boston, MA", caption="Figure 5.1") +
  facet_wrap(~TOD) +
  mapTheme()

Conclusions

Overall, TOD don’t have significant influence in population growth and household income increase and is most influential in decreasing car ownership, but this influence is likely dependent on the density of subway lines and station in the TOD tracts. TOD area also seem to have a strong appeal to more civilized population compared to non-TOD area. As for downtown, its development is the most likely factor driving household income growth and increase of land value within the TOD area. Considering robbery, TOD area in Boston don’t appear to be safer than non-TOD area due to the complicated factors influencing crime. However, it’s noted that neighborhoods with high rent and also within TOD experience less robbery crime than those with high rent but outside the TOD area.

Based on the conclusions, it’s suggested that Boston increase the density of subways lines and add more stations in the fringe of the current TOD area. This action is likely to encourage more people choose more environmentally-friendly means of travel than driving. Also, more high-skilled jobs are preferred in the TOD area since more civilized population is moving there. This action will serve as a virtuous cycle to draw more educated people to Boston while boosting the local economics further more.

Despite all the analysis, it’s worth noting that since this brief is based on ACS tracts data, biases exist in data collection and the accountability of analysis is subject to the Modifiable Areal Unit Problem (MAUP) bias, and the possible ecological fallacy when we apply mean data.