Epi 590 Final Project

Assessing the equity of the planned MARTA expansion in the City of Atlanta

1 Introduction

Access to transportation is an important social determinant of health (SDOH). People’s transit environment affects their ability to reach fresh food and healthcare institutions which can determine health on a daily basis. Transit access also determines what jobs people can reach which has important long-term socioeconomic consequences.

The purpose of my analysis will be to examine equity in the apportionment of funds from the $2.5 billion MARTA expansion tax approved in 2016 in the City of Atlanta. Disease mapping will be used to illustrate the distribution of funds at the census tract level and spatial regression will be used to determine if dollars received by a tract significantly covary with factors like median income, population, racial composition, etc. Necessary data on project budgets will be drawn from newspaper articles and the spatial location of new projects will be digitized into line files from project maps available on the MARTA website.

2 Data Import and Cleaning

Libraries

library(rio)
library(dplyr)
library(bit64)
library(sf)
library(tmap)
library(tigris)
library(lwgeom)
library(tidycensus)

City of Atlanta Census Tracts

City of Atlanta census tracts (prior to annexation of Emory) and associated attribute data were manually downloaded from the 500 cities project page on the CDC website. Next the geographies and attributes were imported and joined with the following code.

five_city <-import('C:/Users/dvern/Desktop/MARTA Project/500_cities.csv')
all_tracts <-st_read('C:/Users/dvern/Desktop/MARTA Project/tracts/500Cities_Tracts_Clip.shp')

atl_tracts <- all_tracts %>% filter(PlaceName == 'Atlanta')

atl_data <- five_city %>% filter(PlaceName == 'Atlanta')

atl_tracts$plctract10 <- as.character(atl_tracts$plctract10)

atl <- left_join(atl_tracts, atl_data, by= c('plctract10'='Place_TractID'))

#st_write(atl, 'C:/Users/dvern/Desktop/MARTA Project/atl_tracts.shp') 

#ESRI shapefile of atl census tracts saved to assist with georeferencing of project lines

Projection

The UTM Zone 16N projection was chosen for this project as it maintains distance relationships in small areas on the scale of cities and/or counties.

x <- 32616  #makes a variable to hold the EPSG code for UTM 16N
atl <-st_transform(atl, crs = x) %>% select(Population2010, TractFIPS)  #trim attributes I'm not concerned with for now
atl$TractFIPS <- as.character(atl$TractFIPS)
atl <- atl %>% rename(GEOID = TractFIPS) #change variable name to facilitate future join operation

Add in Emory/CDC

Emory and the CDC were not incorporated into the City of Atlanta at the time of the 500 cities project. This code will use the tidy census package to download the Dekalb county census tracts, filter down to Emory’s tract, then join it with the rest of Atlanta.

options(tigris_use_cache=TRUE) #prevent redownloading files
dek_tracts <- tracts('GA', county = 'Dekalb') %>% st_as_sf() #download census tracts to Dekalb county

emory <- dek_tracts %>% filter(NAME == '224.02') %>% st_transform(x) %>% select(GEOID) %>% 
  mutate(Population2010 = 1) # need to add placeholder population variable to enable rowbind with rest of ATL

atl_plus <- rbind(atl, emory)

Project Lines File

Digitization of the planned MARTA expansion project segments was conducted in ArcGIS. First the following project map was georeferenced to the atlanta census tracts file from earlier using directions found at this Brown university pdf.

Next the digitized segments were brought into the R environment.

dig_segs <-st_read('C:/Users/dvern/Desktop/MARTA Project/New segs/new_segs1.shp')
dig_segs.proj <-st_transform(dig_segs, crs=x)

dig_segs.proj$Name <- as.character(dig_segs.proj$Name) #assist future join

Finally the projected budget for each project was collated from this Atlanta Journal-Constitution article and imported as a table. This was the joined with the projects line file.

proj_vals <- import('C:/Users/dvern/Desktop/MARTA Project/proj_values.csv')

seg_vals <- left_join(dig_segs.proj, proj_vals) #associates each segment with the project budget for that segment

3 Geoprocessing Steps

Segment Length and budget ratio

The following steps calculate the individual segment lengths and the total length of each project (some comprise many segments, such as the addition of more frequent bus service along some routes). Then a cost per mile ratio was generated from total budget/total length. This ratio calculation took place in a separate data frame due to the use of groupby and the final step joins the ratio back to the original segments file.

#get individual segment length

seg_vals$length <- st_length(seg_vals)
seg_vals$length <- as.numeric(seg_vals$length) #st_length outputs in the units class which is incompatible with calculations
seg_vals <- seg_vals %>% mutate(km = length / 1000) %>% mutate(miles = km*0.621371)  #add lengths in km and miles


#forming the ratio variable for each project

ratio <- seg_vals %>% group_by(Name) %>% summarise(tot_len = sum(miles)) #sum lengths of segments that fall under the same project

rat1 <- left_join(ratio, proj_vals) #join in the proj values so that a $/mile ratio variable can be made
rat2 <- rat1 %>% mutate(ratio = Budg_mil / tot_len)

rat3 <- rat2 %>% st_drop_geometry() %>% select(ratio, Name) #select variables for joining back to main, ungrouped, segments file

#join ratio to original segments file

seg_vals1 <- left_join(seg_vals, rat3)

Joining segments to tracts, calculating total budgeted spending per tract

This step will join segments to tracts using st_intersection which will chop the lines into smaller sections based on what tract they intersect with. This will enable the calculation of the spending per section. Then, once the sections are joined with polygons, aggregation to the polygon will give us the total spending per polygon using the summarize step.

# intersection
int = st_intersection(seg_vals1, atl_plus)
# find out about the length of each line segment
int$len = st_length(int)
int1 <- int %>% mutate(km = len / 1000) %>% mutate(miles = km*0.621371)  #add lengths in km and miles
int2 <- int1 %>% mutate(spending = miles * ratio)

# spatial overlay
join_new = st_join(atl_plus, int2)
# use the ID of the polygon for the aggregation
out = group_by(join_new, GEOID.x) %>%
  summarize(tot_spend = sum(spending))

out[is.na(out)] <- 0 #replace NAs with zeros to assist future calculations

Next the major projects will be separated out by type to assist with future visulaization.

##separate out project types for mapping
LRT <- seg_vals %>% filter(Type == 'LRT')
FLBS <- seg_vals %>% filter(Type == 'FLBS')
BRT <- seg_vals %>% filter(Type == 'BRT')
ART <- seg_vals %>% filter(Type == 'ART' | Name == 'Cleveland AVE')

First Thematic Maps

The following code makes the first thematic maps of the spatial variation in budgeted MARTA spending. Major Project lines are shown for context. A version where the tracts with zero spending were shown as grey was experimented with.

  tm_shape(out) +
  tm_fill('tot_spend', 
          title = 'Budgeted Spending\n(Millions USD)',
          style = 'quantile',
          palette = 'BuPu') + 
  tm_borders() +
#tm_shape(out_none) + tm_fill('spend', col = 'grey', alpha = 1) + tm_borders() +
  tm_shape(LRT) + tm_lines(lwd = 2) +
  tm_shape(BRT) + tm_lines(lwd = 2, col = 'red') +
  tm_shape(ART) + tm_lines(lwd = 2, col = 'blue') + 
    #tm_shape(FLBS) + tm_lines(lwd = 1, col = 'navy') + 
  #tm_add_legend(type = 'fill', col = 'grey', labels = '0.00') +
  tm_add_legend(type = 'line', col = 'black', lwd = 3, labels = 'LRT') +
  tm_add_legend(type = 'line', col = 'blue', lwd = 3, labels = 'ART') +
  tm_add_legend(type = 'line', col = 'red', lwd = 3, labels = 'BRT') + 
  
  tm_credits('Source: Atlanta Journal-Constitution', fontface = 'bold', position = c('LEFT','BOTTOM')) +
  
  tm_layout(legend.title.fontface = 'bold')

*Note that this appears uncrowded in the normal viewer and upon export as a picture.

#make a dataset with just zero spending tracts
out$spend <- as.numeric(out$tot_spend)
out_none <- out %>% filter(spend == 0)
#map with greyed out zeros
  tm_shape(out) +
  tm_fill('tot_spend', 
          title = 'Budgeted Spending\n(Millions USD)',
          style = 'quantile',
          palette = 'BuPu') + 
  tm_borders() +
tm_shape(out_none) + tm_fill('spend', col = 'grey', alpha = 1) + tm_borders() +
  tm_shape(LRT) + tm_lines(lwd = 1) +
  tm_shape(BRT) + tm_lines(lwd = 1, col = 'red') +
  tm_shape(ART) + tm_lines(lwd = 1, col = 'blue') + 
    tm_shape(FLBS) + tm_lines(lwd = 1, col = 'navy') + 
  tm_add_legend(type = 'fill', col = 'grey', labels = '0.00') +
  tm_add_legend(type = 'line', col = 'black', lwd = 3, labels = 'LRT') +
  tm_add_legend(type = 'line', col = 'blue', lwd = 3, labels = 'ART') +
  tm_add_legend(type = 'line', col = 'red', lwd = 3, labels = 'BRT') + 
  
  tm_credits('Source: Atlanta Journal-Constitution', fontface = 'bold', position = c('LEFT','BOTTOM')) +
  
  tm_layout(legend.title.fontface = 'bold')

A close examination of this second map reveals that there are a few central tracts which probably shouldn’t be identified as zero spending. However, due to small inaccuracies in the georeferencing and digitization process, some segments and tracts don’t overlap when they probably should. I say probably because in most cases, the tract boundaries appear to follow major transportation routes.

Unfortunately, at this stage I’ve already spent a LOT of time geoprocessing and I need to get to the analysis. To improve this in the future, I would need to compare the tract boundaries in question and see if they indeed fall along major roads/ transit routes that are in the MARTA project list. From there, I’d need to redigitize segments that match up better.

On to analysis.

Downloading, cleaning, and joining other attribute data

Just kidding about being ready to do analysis. Now I need to download some covariate data (population and median income) from the ACS for my tracts. This took a fair amount of processing as Atlanta spreads across both Fulton and Dekalb counties.

#### attributes to tracts join
  
  library(tidycensus)
  
  all_vars <- load_variables(year = 2017, dataset = 'acs5', cache = T)
  
  # From my search, I found that I want the following variables: 'B25026_001','B21004_001'
  
  ## this part grabs the variables for the atl tracts in fulton county
  
  fult_pop <- get_acs(geography = 'tract',
                       variables = 'B25026_001',
                       county = 'Fulton',
                       state = 'Georgia',
                      year = 2017,
                       survey = 'acs5') %>% rename(population_2017= estimate)
  
  fult_inc <- get_acs(geography = 'tract',
                      variables = 'B21004_001',
                      county = 'Fulton',
                      state = 'Georgia',
                      year = 2017,
                      survey = 'acs5') %>% rename(med_income_17= estimate) %>% select(med_income_17, GEOID)
  
  
  atl_pop <-left_join(atl_plus, fult_pop)
  
  atl_pop_inc <- left_join(atl_pop, fult_inc)

  atl_atts_fult <- atl_pop_inc %>% select(GEOID, population_2017, Population2010, med_income_17)
  
  
  
  ### This part grabs the variables for the atl tracts in dekalb county
  
  dek_pop <- get_acs(geography = 'tract',
                      variables = 'B25026_001',
                      county = 'Dekalb',
                      state = 'Georgia',
                      year = 2017,
                      survey = 'acs5') %>% rename(population_2017= estimate)
  
  dek_inc <- get_acs(geography = 'tract',
                      variables = 'B21004_001',
                      county = 'Dekalb',
                      state = 'Georgia',
                      year = 2017,
                      survey = 'acs5') %>% rename(med_income_17= estimate) %>% select(med_income_17, GEOID)
  
  
  dek_pop <-left_join(atl_plus, dek_pop)
  
  dek_pop_inc <- left_join(dek_pop, dek_inc)
  
  atl_atts_dek <- dek_pop_inc %>% select(GEOID, population_2017, Population2010, med_income_17) %>% filter(is.na(population_2017) == F) %>% st_drop_geometry()
  
 
  
  ## now join fulton and dekalb polys into one atlanta attributes dataframe
  
  atl_atts <- left_join(atl_atts_fult, atl_atts_dek, by= 'GEOID')
  atl_atts1 <- atl_atts %>% 
    mutate(pop_17 = ifelse(is.na(population_2017.x),population_2017.y, population_2017.x)) %>%
    mutate(inc_17 = ifelse(is.na(med_income_17.x),med_income_17.y, med_income_17.x)) %>% select(GEOID, pop_17, inc_17) %>% st_drop_geometry()
    
  #above mutate code is necessary to consolidate the fulton / dekalb attributes into single columns
  
  
  # join with spending polygons
  
  
  atl_spend_atts <- left_join(out, atl_atts1, by = c('GEOID.x' = 'GEOID'))
  
  atl_spend_atts <- atl_spend_atts %>% filter(is.na(pop_17) == F) #remove some weird tracts that don't map and have no data

  

4 Analysis

4.1 Clusters

Looking at the map of budgeted spending by census tract it seems obvious that the constant risk assumption is violated. However, we’d like to know in what areas the spending is higher than expected based on the population of a tract. Thus, we will produce a map of poissson probabilities to test for and locate the local extremes.

library(spdep) #package for performing cluster analysis

atl_spend_atts$tot_spend <- as.numeric(atl_spend_atts$tot_spend) #conmvert from units class to numeric

atl_spend_atts <- atl_spend_atts %>% mutate(spend_ratio = ifelse( pop_17 > 0,(tot_spend * 1000000 / pop_17), NA)) %>%filter(is.na(spend_ratio) == F)

spend_probs <- probmap(n =atl_spend_atts$tot_spend,
                       x = atl_spend_atts$pop_17,
                       row.names = atl_spend_atts$GEOID.x)

atl_spend_atts$pmap <- spend_probs$pmap
#generate poisson probs

# This just summarizes (merges polygons) of all counties with significant HIGH risk
spend_hi <- atl_spend_atts %>% #form contiguous areas of sig high or low spending into polygons
  filter(pmap > 0.95) %>% summarise()
spend_lo <- atl_spend_atts %>%
  filter(pmap < 0.05) %>% summarise()
# make map of sig high or low spending areas
m1 <-tm_shape(atl_spend_atts) + 
  tm_fill('spend_ratio',
          palette = 'BuPu',
           style = 'quantile',
          title = 'Spending per\nPerson (USD)') +
  tm_borders() + 
  tm_layout(legend.position = c('RIGHT','TOP'),
            main.title = 'Clusters of MARTA\nExpansion Spending by Population') +
  
# Now plot only borders of high or low rate counties
  tm_shape(spend_lo) + 
  tm_borders(col = 'blue', lwd = 2) + 
  tm_shape(spend_hi) + 
  tm_borders(col = 'red', lwd = 2) +
# general layout stuff...  
tm_layout(legend.outside = F, legend.text.size = 0.55, 
          legend.title.size = 0.65, legend.title.fontface = 'bold',
          legend.position = c('LEFT', 'TOP'), inner.margins = c(0.02, .03, .02, 0),
         main.title.position = 'center', main.title.size = 1, main.title.fontface = 'bold') + 
   tm_add_legend(type = 'line',
              col = 'red',
              lwd = 2,
              labels = 'Significantly high') +
   tm_add_legend(type = 'line',
              col = 'blue',
              lwd = 2,
              labels = 'Significantly low')

Looking at the above map, almost all of the Atlanta census tracts fall into the category of significant high or low spending. Now, this was done with population count as the denominator, but what if using population density would give us a more clear picture?

The following code brings in the land area (km2) variable for each tract and calculates population density in persons per km2.

dek_tracts$ALAND <- as.numeric(dek_tracts$ALAND) #convert to numeric for calculations

dek_tracts_a <- dek_tracts %>% 
  select(GEOID, ALAND) %>% mutate(km2 = ALAND/1000000) %>% st_drop_geometry() 

#fulton tracts

ful_tracts <- dek_tracts <- tracts('GA', county = 'Fulton') %>% st_as_sf() %>% st_drop_geometry()

ful_tracts$ALAND <- as.numeric(ful_tracts$ALAND)

ful_tracts_a <- ful_tracts %>% 
  select(GEOID, ALAND) %>% mutate(km2 = ALAND/1000000)

                                                      
#joins and cleaning

fult_a <-left_join(atl_plus, ful_tracts_a)

dek_a <-left_join(atl_plus, dek_tracts_a)

all_a <- rbind(fult_a, dek_a) %>% filter(is.na(ALAND) == F) %>% select(GEOID, km2) %>% st_drop_geometry() #drop geometry to make a table that can be joined to main attributes file

area_atts <- left_join(atl_spend_atts, all_a, by=c('GEOID.x'='GEOID')) 
area_atts <-area_atts %>% mutate(pop_dens = pop_17 / km2) #create pop density variable in persons per km2

                               

Now we’ll run the test of local extremes again but with population density as the denominator. Finally, we’ll map the significant results to see how different it is than using pop counts.

Analysis:

spend_probs_dens <- probmap(n =area_atts$tot_spend,
                       x = area_atts$pop_dens,
                       row.names = area_atts$GEOID.x)
area_atts$pmap <- spend_probs_dens$pmap
#generate poisson probs
area_atts <- area_atts %>% mutate(spend_ratio_d = tot_spend *1000000/ pop_dens)

Map:

# This just summarizes (merges polygons) of all counties with significant HIGH risk
spend_hi_d <- area_atts %>% #form contiguous areas of sig high or low spending into polygons
  filter(pmap > 0.95) %>% summarise()
spend_lo_d <- area_atts %>%
  filter(pmap < 0.05) %>% summarise()
# make map of sig high or low spending areas
 m2 <-tm_shape(area_atts) + 
  tm_fill('spend_ratio_d',
          palette = 'BuPu',
           style = 'quantile',
          title = 'Spending per\nPerson Per Km2 (USD)') +
  tm_borders() + 
  tm_layout(legend.position = c('RIGHT','TOP'),
            main.title = 'Clusters of MARTA\nExpansion Spending by Pop. Density') +
  
# Now plot only borders of 'high-rate' counties
  tm_shape(spend_lo_d) + 
  tm_borders(col = 'blue', lwd = 2) + 
  tm_shape(spend_hi_d) + 
  tm_borders(col = 'red', lwd = 2) +
#general layout 
tm_layout(legend.outside = F, legend.text.size = 0.55, 
          legend.title.size = 0.65, legend.title.fontface = 'bold',
          legend.position = c('LEFT', 'TOP'), inner.margins = c(0.02, .03, .02, 0),
         main.title.position = 'center', main.title.size = 1, main.title.fontface = 'bold') + 
   tm_add_legend(type = 'line',
              col = 'red',
              lwd = 2,
              labels = 'Significantly high') +
   tm_add_legend(type = 'line',
              col = 'blue',
              lwd = 2,
              labels = 'Significantly low')

tmap_arrange(m1, m2)

The use of population density does change the distribution of significantly extreme tracts such that there are fewer significant tracts overall and some tracts switch from sig hi to low and vice versa. Either way you slice it, there are clear significant extreme clusters of spending. However, using pop density as the denominator seems to remove some of the noise caused by having similar populations in drastically differently sized tracts.

4.2 Regression

Now that I’ve established significant clustering of high and low spending exists,I’d like to see if this phenomenon can be explained by other underlying factors. In particular, I am interested to examine the correlation of median tract income with budgeted spending.

First I will fit my linear regression models.

#empty model

#before my tot_spend variable was in millions of dollars, I will convert it back to dollars for running my models

area_atts$tot_spend <- area_atts$tot_spend / 1000000

m0 <- lm(tot_spend ~ 1, data = area_atts)

#conditional model

m1 <- lm(tot_spend ~ pop_dens + inc_17, data = area_atts)

summary(m1)

The results of my conditional model show that, holding population density constant, 2017 median tract income was a significant negative predictor of budgeted MARTA spending. Based on the model, each dollar increase in median income is associated with a 796.8 dollar decrease in planned spending. However, it is important to note that since I have zero dollar values for 26 of 130 tracts, the residuals and residual standard error for this model are huge.

summary(area_atts$spend)

Now I’ll test the residuals from my model to see whether our factors explain any of the variation in spatial autocorrelation between tracts.

This will involve first choosing a neighbors definiton and then running the lm.morantest function. I have decided to use k-nearest neighbors of 3 due to the fact that emory is not touching any other atl tracts and several border tracts have only one contiguous neighbor.

##Creating nb object
#since we're using proximity based nb object, we need the tract centroid points
area_atts_sp <- as(area_atts, 'Spatial')
area_points  <- coordinates(area_atts_sp)
area_k3 <- knearneigh(area_points, k =3) #creates set of 3 nearest neighbors
area_k3_nb <- knn2nb(area_k3, row.names = area_atts_sp$GEOID.x, sym = T) #creates a symmetrical nb object
area_listw <- area_k3_nb %>% nb2listw() #creates listw object for Moran's I test
##performing lm.morantest()
#empty model residuals
lm.morantest(m0, listw = area_listw)

    Global Moran I for regression residuals

data:  
model: lm(formula = tot_spend ~ 1, data = area_atts)
weights: area_listw

Moran I statistic standard deviate = 4.5259, p-value = 3.007e-06
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.281758061     -0.007751938      0.004091894 

#conditional model residuals
#m1_1 <- lm(tot_spend ~ inc_17, data = area_atts) #I thought of to just test the effect of income since pop_dens wasn't significant but decided to include pop_dens after all
lm.morantest(m1, listw = area_listw)

    Global Moran I for regression residuals

data:  
model: lm(formula = tot_spend ~ pop_dens + inc_17, data = area_atts)
weights: area_listw

Moran I statistic standard deviate = 3.9708, p-value = 3.581e-05
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.233938193     -0.016975080      0.003992891 

From these tests we see there is relatively low spatial autocorrelation in model residuals for budgeted spending. Of that, only about 17% is explained by tract median income and population density. My conjecture is that much of the remaining autocorrelation is explained by the fact that many of the high-value projects are comprised of beltline rail. By definition, all these projects are contiguous and since the beltline follows old rail throughways in the city, their locations may be less entagled with other demographic factors.

For the final figures in this project, I’ll map the model residuals.

The maps illustrate that our factors slightly reduce spatial autocorrelation but there is still a fair amount of spatial structure present due to the periphery tracts with zero spending. Interestingly, some of the largest residuals are in the central portion of the city. This is where, as I mentioned previously, innacuracies in digitization may have missed assigning spending to some tracts.