Overview

Local intelligence is a strong predictor of housing prices in Philadelphia that can optimize Zillow’s popular and consequential Zestimate tool. Real estate professionals as well as home buyers and sellers alike trust Zestimate to be an accurate prediction of local real estate values. Investing in the tool’s reliability stands to not only benefit Zillow’s customers, but the company’s dominant market position, too.

The following study explores a series of local intelligence variables that range from tree-lined street to local median income. A series of statistical processes reveal a handful of key variables that are then processed by prediction models. The best prediction model is then tested on 2023 sale price data for accuracy. The study is mindful of creating models that are generalizable across Philadelphia, and hopefully generalizable across major metro areas throughout the country.

Select and Load Data

Firstly, sales data is a tremendous resource in terms of the internal characteristics of each property which can certainly influence price. This study leverages data such as year built, total square footage, number of bathrooms and bedrooms to holistically understand what influences price.

However, the property itself does not tell the full story. Where it is located and the characteristics of that area are key, too. There are myriad local intelligence variables to chose from between a number of key open data sources (2021-2023), including: U.S. Census American Community Survey, Pennsylvania Spatial Data Access (PASDA), Philadelphia Open Data Portal, or from the Delaware Valley Regional Planning Commission (DVRPC). This study priorizited complete datasets that affect the city at large. For example, litter is a systemic issue whereas airports are very localized. The following table breaks down viable sources and from which we selected local intelligence variables that we deemed to be generalizable for the purposes of this study:

# Create a data frame with your data
data_frame <- data.frame(
  "Dataset" = c(
    "Floodplain Data",
    "Sanitation Centers",
    "Litter Index",
    "Demographics and Socioeconomic Data",
    "Housing Typology",
    "Trees",
    "Neighborhood",
    "Proximity to Highways",
    "School Catchment Performance",
    "Walkscore", 
    "1937 Redlining Score"
  ),
  "Generalizability" = c(
    "",   # Floodplain Data 
    "",   # Sanitation Centers 
    "✔",  # Litter Index 
    "✔",  # Demographics and Socioeconomic Data 
    "",   # Housing Typology 
    "✔",  # Trees 
    "✔",  # Neighborhood
    "✔",  # Proximity to Highways 
    "",   # School Catchment Performance 
    "",   # Walkscore 
    "✔"  # 1937 Redlining Score
  )
)


dt_url <- c("https://data-phl.opendata.arcgis.com/datasets/phl::fema-100-flood-plain/explore",
            "https://opendataphilly.org/datasets/sanitation-areas/",
            "https://opendataphilly.org/datasets/litter-index/",
            "https://www.census.gov/",
            "https://www.census.gov/",
            "https://opendataphilly.org/datasets/phillytreemap/",
           "https://www.pasda.psu.edu/uci/DataSummary.aspx?dataset=7047",
           "https://dvrpc-dvrpcgis.opendata.arcgis.com/datasets/dvrpcgis::freight-highway-network/explore",
           "https://www.philasd.org/performance/programsservices/open-data/school-information/",
           "http://walkscore.com/",
           "https://ncrc.org/holc/")

data_frame %>%
  mutate(Dataset = cell_spec((Dataset), "html", link = dt_url)) %>%
  kable(format = "html", escape = FALSE) %>%
  kable_styling(bootstrap_options = c("hover", "condensced"), full_width = FALSE, 
                position = "left") %>%
  column_spec(column = 2, width = "20%", extra_css = "text-align:center;"
               )
Dataset Generalizability
Floodplain Data
Sanitation Centers
Litter Index ✔
Demographics and Socioeconomic Data ✔
Housing Typology
Trees ✔
Neighborhood ✔
Proximity to Highways ✔
School Catchment Performance
Walkscore
1937 Redlining Score ✔

Boundaries

From highways to rivers, man-made and natural boundaries alike delineate neighborhoods where certain people and types of businesses congregate. Neighborhood boundaries are arguably more organic than census block groups (which are optimized for consistent populations) or highway. Additionally, because residential sales are generalized by residential neighborhoods not by census block groups, neighborhoods (as defined by the Philly Planning Dept) are likely significant on local intelligence.You can see where areas of high income and low income cluster across Philadelphia, and how these boundaries are more recognizable at a neighbrhood level than block group level.

# plot sales data using median cut off  

sales_plot_mod <- sales2023 %>%
  filter(toPredict == "MODELLING")

median_sale_price <- median(sales_plot_mod$sale_price)
custom_colors <- ifelse(sales_plot_mod$sale_price <= median_sale_price, "#A5C4D4", "#4B7740")

sales_plot <- ggplot() +
  geom_sf(data = blocks21.sf, color = "gray70") +
  geom_sf(data = sales_plot_mod, aes(color = custom_colors), size = 0.5, alpha = 0.5) +
  scale_color_manual(
    values = c("#A5C4D4" = "#A5C4D4", "#4B7740" = "#4B7740"),
    labels = c("#4B7740" = "Above Median", "#A5C4D4" = "Below Median")) +
  labs(title = "Block Groups", color = "") +
  theme(legend.position = "none") + 
  guides(color = guide_legend(override.aes = list(size = 3))) +
  labs(subtitle = "2021 Philadelphia Median \nSale Price is $229,900") +
  mapTheme()

#Bring in neighborhoods

hoods <- st_read("/Users/LauraFrances_1/Library/CloudStorage/GoogleDrive-lauramuellersoppart@gmail.com/My Drive/Penn/Courses/MUSA508_PPA/PPA_Midterm/PPA_Midterm_Git/data/PhillyPlanning_Neighborhoods/PhillyPlanning_Neighborhoods.shp") %>%
  select(NAME) %>%
  rename(Neighborhood = NAME)

sales2023 <- st_join(sales2023, hoods, join = st_within)

hood_plot <- ggplot() +
  geom_sf(data = blocks21.sf, fill = "gray90", color = "gray90") +
  geom_sf(data = hoods, color = "gray60") +
  geom_sf(data = sales_plot_mod, aes(color = custom_colors), size = 0.5, alpha = 0.5) +
  scale_color_manual(
    values = c("#A5C4D4" = "#A5C4D4", "#4B7740" = "#4B7740"),
    labels = c("#4B7740" = "Above Median", "#A5C4D4" = "Below Median")) +
  labs(title = "Neighborhoods") +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  labs(subtitle = "Neighborhoods as defined by Philadelphia \nPlanning Dept") +
  mapTheme()

#Bring in highways and apply buffer of 1/8 mile (660')

highways <- st_read("/Users/LauraFrances_1/Library/CloudStorage/GoogleDrive-lauramuellersoppart@gmail.com/My Drive/Penn/Courses/MUSA508_PPA/PPA_Midterm/PPA_Midterm_Git/data/Freight_Highway_Network.geojson")

highways <- highways %>%
  filter(state == "PA") %>%
  st_transform('EPSG:2272')

highways_phl <-  st_intersection(highways, blocks21.sf) %>%
  select(geometry)

highway_buffer <- st_buffer(highways_phl, 660)
highway_buffer$hwy_grade <- 1

# Make a dummy variable where 1 is a house within highway buffer and 0 is outside of buffer. 

saleshighways <- st_intersection(sales2023, highway_buffer)

homes_near_hwys <- saleshighways$objectid  #creates a list using unique ids from saleshighways

sales2023 <- sales2023 %>%
  mutate(hwy_grade = ifelse(objectid %in% homes_near_hwys, 1, 0))

hwy_plot <- ggplot() + 
  geom_sf(data = blocks21.sf, color = "gray70") +
   geom_sf(data = sales_plot_mod, aes(color = custom_colors), size = 0.5, alpha = 0.5) +
  scale_color_manual(
    values = c("#A5C4D4" = "#A5C4D4", "#4B7740" = "#4B7740"),
    labels = c("#4B7740" = "Above Median", "#A5C4D4" = "Below Median")) +
  geom_sf(data = highway_buffer, aes(fill = "#36151E", color = "#36151E")) +
  scale_fill_identity(name = "Highway", labels = "#36151E") +
  scale_color_identity(name = "Highway", labels = "#36151E") +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  labs(title = "Highway Proximity") +
  labs(subtitle = "The highway buffer is 660 feet, \nabout 2 Philadelphia blocks.") +
  mapTheme()

#plot all three together

plot_grid(sales_plot, hood_plot, hwy_plot, ncol = 3, align = "v")

It is important to be aware of how the sales are distributed around the city. There are many block groups without any housing sales, whereas there are fewer neighborhoods where this is true because neighborhoods are residential in nature. Geographically, because neighborhoods are generally divided along residential vs industrial zoning lines, it is simpler to comprehend how it is generalizable across the city.

Interestingly, when considering the impact of highways in Philadelphia, it is reasonable to assess that it is not generalizable. However, given how consistently home prices drop below the median within proximity to the highway regardless of the sale is in Center City or Port Richmond, we decided that it is a variable worth further investigation.

Redlining and Price Clustering

In the 1930s, a federal agency, the Home Owners’ Loan Corporation (HOLC), created maps of mortgage risk which denied entire neighborhoods access to public lending benefits and overall economic opportunity. The institutionalize of marginalization of poor, Black communities ingrained persistent residential segregation that is still evidenced by income inequality across the city.

# Get and wrangle redlining data

redlining.sf <- st_read("/Users/LauraFrances_1/Library/CloudStorage/GoogleDrive-lauramuellersoppart@gmail.com/My Drive/Penn/Courses/MUSA508_PPA/PPA_Midterm/PPA_Midterm_Git/data/Redlining.geojson") %>% 
  st_transform('EPSG:2272') %>% 
  filter(state == "PA", city == "Philadelphia") 

redlining.sf <- redlining.sf %>%
  dplyr::select(-c(holc_id, state, city, neighborhood_id, area_description_data, name)) %>%
  drop_na() 

# Join each house with historic redlining grade 

sales2023 <- st_join(sales2023, redlining.sf, join = st_within) 

sales2023 <- sales2023 %>% 
  dplyr::mutate(holc_grade = ifelse(is.na(holc_grade) == TRUE, "NA", holc_grade))

# Plotting Redlining Maps 

redlining_colors <- c("NA" = "gray85", "A" = "#197991", "B" = "#4B7740", "C" = "#E3D26F", "D" = "#D52F4D")
redlining_palette <- c("#197991", "#4B7740", "#E3D26F", "#D52F4D", "#A9233B")

redlining.sf.phl <- st_intersection(redlining.sf, blocks21.sf) #clip redlining to PHL boundaries
  
Redlining_plot <- ggplot() +
  geom_sf(data = blocks21.sf, color = "gray70") +
  geom_sf(data = redlining.sf.phl, aes(fill = holc_grade, text = holc_grade), alpha = 0.8) +
  scale_fill_manual(values = redlining_colors, name = "HOLC Grade") +
  labs(title = "1937 HOLC Redlining", caption = "by 2020 Block Group") +
  theme(legend.position = "right") +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  mapTheme() 

Income_plot <- ggplot() +
  geom_sf(data = blocks21.sf, color = "grey70") +
  geom_sf(data = sales2023, aes(colour = q5(MedHHInc)), 
          show.legend = "point", size = .5, alpha = 0.5) +
  scale_colour_manual(values = rev(redlining_palette),
                   labels=qBr(sales2023,"MedHHInc"),
                   name="Income\nBrackets") +
  labs(title="2023 Median Incomes", caption = "by 2020 Block Group") +
  theme(legend.position = "right") +
  guides(color = guide_legend(override.aes = list(size = 3))) +
  mapTheme()

plot_grid(Redlining_plot, Income_plot, ncol = 2, align = "v")

NCRC states that, “There is significantly greater economic inequality in cities where more of the HOLC graded high-risk or “Hazardous” areas are currently minority neighborhoods.” (NCRC) The impact of HOLC is generalizable. One exception to note is Center City where it has many high income earners despite its high-risk grade likely due to post-war downtown revitalization.

When considering the patterns of redlining and neighborhoods, it is likely that neighborhoods such as Center City will account for divergences in the redlining data, and neighborhoods such as Kensington will account for consistencies in the redlining data when we develop our definition of local intelligence.

## K NEAREST NEIGHBOR: PRICE

sales2023 <- sales2023 %>% 
  mutate(nb = st_knn(geometry, k = 5),
         wt = st_weights(nb),
         price_lag = st_lag(sale_price, nb, wt))

Price Lag : To systemically investigate price clustering beyond redlining and in today’s data, a function called “k-nearest neighbor” calculates a weighted average of the 5 nearest properties’ sales price. This price lag variable accounts Philadelphia’s ‘block-by-block’ dynamics borne from its density and long history where many uses from mansions to power plants have converged within mere feet of one another - but can drastically impact sale price potentially.

Environmental Conditions

Proximity to amenities and distance from disamenities certainly has a colloquial effect on housing prices. Broker descriptions love to tout, “2 Blocks from Park!” or “Tucked away, quiet street!” To test the effect of environmental data, the study inclues tree density and litter conditions.

## TREES by BLOCK GROUP

trees.sf <- st_read("https://opendata.arcgis.com/api/v3/datasets/30ef36e9e880468fa74e2d5b18da4cfb_0/downloads/data?format=geojson&spatialRefId=4326") %>%
  st_transform('EPSG:2272') 

trees.sf <- st_join(trees.sf, blocks21.sf, join = st_within)

trees.sf <- trees.sf %>%
  select(c(OBJECTID, TREE_NAME, YEAR, GEOID))

trees.sf <- trees.sf %>%
  group_by(GEOID) %>%  # census block group
  summarize(trees = n()) 

blocks21.sf.trees <- st_join(blocks21.sf, trees.sf, by = "GEOID")

blocks21.sf.trees <- mutate(blocks21.sf.trees, trees = ifelse(is.na(trees), 0, trees)) %>% 
  select(trees, GEOID.y) 

sales2023 <- sales2023 %>%
  st_join(blocks21.sf.trees, join = st_intersects)

trees_plot <- ggplot() +
  geom_sf(data = blocks21.sf, fill = "gray20", color = "gray40") +
  stat_density2d(data = data.frame(st_coordinates(trees.sf)),
                 aes(X, Y, fill = ..level.., alpha = ..level..),
                 size = 0.01, bins = 30, geom = 'polygon') +
  scale_fill_gradient(
    high = "#4B7740", 
    low = "#ABCDA2", 
    name = "",
   breaks = c(0.0000000001, 0.000000001),  # Set the breaks for low and high
    labels = c("low", "high")) +  # Set the corresponding labels
  scale_alpha(range = c(0, .7), guide = "none") +
  labs(title = "Tree Density", color = "") +
  theme(legend.position = "bottom") + 
  guides(color = guide_legend(override.aes = list(size = 3))) +
  labs(caption = "Tree Density is measured by \n number of trees in the block group") +
  mapTheme()

# NEAREST LITTER TO HOUSE

litterindex <- st_read("https://opendata.arcgis.com/datasets/04fa63e09b284dbfbde1983eab367319_0.geojson") %>%
  st_transform('EPSG:2272') %>% 
  rename(litterscore = HUNDRED_BLOCK_SCORE) %>% 
  rename(littergrade = SCORE_COLOR)

nearest_litter <- st_nearest_feature(sales2023, litterindex)
sales2023$litterscore <- litterindex$litterscore[nearest_litter]
sales2023$littergrade <- litterindex$littergrade[nearest_litter]

sales2023 <- sales2023 %>%
  mutate(litterscore = round(litterscore, 2))

litter_plot <- ggplot() +
  geom_sf(data = blocks21.sf, fill = "gray20", color = "gray40") +
  stat_density2d(data = data.frame(st_coordinates(litterindex)),
                 aes(X, Y, fill = ..level.., alpha = ..level..),
                 size = 0.01, bins = 30, geom = 'polygon') +
  scale_fill_gradient(low = "#E3D26F", high = "#36151E", name = "",
                      breaks = c(0.0000000002, 0.0000000006),
                      labels = c("low", "high")) +
  scale_alpha(range = c(0, .7), guide = "none") +
  labs(title = "Litter Score", color = "") +
  theme(legend.position = "bottom") + 
  guides(color = guide_legend(override.aes = list(size = 3))) +
  labs(caption = "Litter Score is based on the Litter Index \ndeveloped by the City of Philadelphia") +
  mapTheme()

# plot two maps
plot_grid(trees_plot, litter_plot, ncol = 2, align = "v")

Interestingly, these two environmental data points may send mixed signals. Where there are more trees, there is more litter. Since tree density and litter conditions send strong visual cues to potential buyers of the level of care of the neighborhood, the impact of these variables on the model may strengthen or dilute the results.

Build Prediction Models

“All models are wrong, but some are useful.” - George Box

Zestimate seeks to predict housing prices. Predicting city trends is fraught. Oftentimes unpredictable political decisions such as post-war state-sponsored suburbanization (aka white flight) buck trend predictions. There is no model than account for every unforeseeable scenario. However as previously discussed, useful models are generalizable. Building a useful prediction model requires feature selection and feature engineering. These features are then used to train a model to recognize patterns associated with sale price. The predictive model will use those patterns to produce a Zestimate. While the nearest prices of nearby sales are likely a fine starting point to build Zestimate, including local intelligence will certainly strengthen the product.

Categorize Variables

This study uses the hedonic model (reference) theoretical framework that deconstructs housing prices into three categories:

  1. Internal Characteristics
  2. Public Services and Disamenities
  3. Spatial Process of Prices
# SELECTING VARIABLES TO INCLUDE IN MODEL

## PREPPING SALES DATA INTO MODEL V CHALLENGE

sales2023.mod <- sales2023 %>%
  filter(toPredict == "MODELLING")

sales2023.chal <- sales2023 %>%
  filter(toPredict == "CHALLENGE")

all_var <- c("objectid", "census_tract","GEOID", "toPredict","sale_price", "Neighborhood", "fireplaces", "garage_spaces", "number_of_bathrooms", "number_of_bedrooms", "number_stories", "total_area", "total_livable_area", "year_built","historic", "holc_grade", "price_lag", "trees", "litterscore", "TotalPop", "MedHHInc", "pctdetached", "pctmultifam_small", "hwy_grade")

sales2023.select <- sales2023.mod %>% 
  dplyr::select(all_of(all_var)) %>%
  na.omit()

sales2023.summary <- sales2023.select %>%
  st_drop_geometry() %>% 
  select(-c(objectid, census_tract, GEOID, toPredict))


# table of summary statistics with variable descriptions. Sort these variables by their category (internal characteristics, amenities/public services or spatial structure). Check out the `stargazer` package for this.

categories_to_include <- c("Amenities / Public Services", "Internal Characteristic", "Spatial Process")

#subsetting categories

#internal characteristics
InternalCharacteristics <- sales2023 %>%
  dplyr::select("fireplaces", "garage_spaces", 
          "number_of_bathrooms", "number_of_bedrooms", "number_stories", 
           "total_area", "total_livable_area", "year_built", "historic", "sale_price")

stargazer(as.data.frame(InternalCharacteristics), type="text", digits=1, title = "Internal Characteristics", out = "Training_PHLInternal.txt", omit.summary.stat = "N")

Internal Characteristics
=====================================================
Statistic             Mean    St. Dev.  Min    Max   
-----------------------------------------------------
fireplaces            0.05       0.3     0      5    
garage_spaces          0.4       0.7     0     72    
number_of_bathrooms    1.1       0.7     0      8    
number_of_bedrooms     2.6       1.3     0     31    
number_stories         2.0       0.6     1      5    
total_area           1,828.9   3,818.9   0   226,512 
total_livable_area   1,333.6    548.8    0   15,246  
year_built           1,936.7    51.1     0    2,024  
historic               1.9       0.4     1      3    
sale_price          271,846.8 240,083.4  0  5,990,000
-----------------------------------------------------
#amenities/public services

Amenities_Public_Services <- sales2023 %>%
  dplyr::select("Neighborhood", "trees", "TotalPop", "MedHHInc", "pctdetached", "pctmultifam_small", "litterscore")

stargazer(as.data.frame(Amenities_Public_Services), type="text", digits=1, title = "Public Services and Disamenities", out = "Training_PHLInternal.txt", omit.summary.stat = "N")

Public Services and Disamenities
===============================================
Statistic           Mean   St. Dev. Min   Max  
-----------------------------------------------
trees              116.7    112.9    2   2,136 
TotalPop          1,445.2   708.6    0   4,518 
MedHHInc          56,282.7 39,156.1  0  250,001
pctdetached         9.0      12.7   0.0  100.0 
pctmultifam_small   12.4     12.2   0.0  75.0  
litterscore         1.8      0.6    0.0   4.0  
-----------------------------------------------
#spatial process

Spatial <- sales2023 %>%
  dplyr::select("price_lag")

stargazer(as.data.frame(Spatial), type="text", digits=1, title = "Spatial Process of Prices", out = "Training_PHLInternal.txt", omit.summary.stat = "N")

Spatial Process of Prices
==================================================
Statistic   Mean    St. Dev.    Min        Max    
--------------------------------------------------
price_lag 267,985.6 194,111.3 20,000.0 3,292,800.0
--------------------------------------------------

As we feature select and engineer, here are some important takeaways from the above summary tables:

  • “price_lag” is the weighted average of nearby sales. These values are closely mirrored by actual sale price. Nearby prices are likely very influential so it is important to select features that can tell another side of the story.
  • Sales price ranges considerably in Philadelphia, and identifying other data that contains the same variety may be influential to help inform spatial price variation.
  • The variation of trees is quite high, with a minimum of only 2 in some areas, and a maximum of nearly 2150. This may be error or it may be important disparity to account for.
  • Total Population is based on census block groups which are engineered to normalize to 600 to 3,000 people.

Test for Correlation

The goal is to select variables that most significantly correlate to sale price to include in the predictive model. Correlation is a type of association test. Are sale prices more closely associated to square footage or to number of bedrooms? These are the types of subtle but important distinction we want to seek out.

# CORRELATION MATRIX OF VARIABLES

numericVars <- sales2023.select %>%
  st_drop_geometry() %>%
  select(-toPredict, -GEOID, -objectid, -census_tract) %>%
  select(-holc_grade, -hwy_grade, -historic, -Neighborhood) #remove categorical variables for corrtest

price.corr <- corr.test(numericVars)

if (sum(is.na(numericVars)) > 0) {
  # Handle missing values (e.g., impute or remove)
  numericVars <- na.omit(numericVars)  # Remove rows with missing values
}

correlation_matrix <- round(cor(numericVars), 1)
p_values <- cor_pmat(numericVars)

ggcorrplot(
  correlation_matrix,
  p.mat = p_values,
  colors = c("#36151E", "white", "#197991"),
  type = "lower",
  insig = "blank"
) +
labs(title = "Degrees of Correlation to Sale Price") +
  annotate(
  geom = "rect",
  xmin = .5, xmax = 16, ymin = .5, ymax = 1.5,
  fill = "transparent", color = "red", alpha = 0.5
)

The correlation matrix revealed strong positive and negative correlations. These are the primary variables to note for further engineering:

  • Fireplaces
  • Number of bathrooms
  • Total livable area (sf)
  • Tree density
  • Median Household Income
  • Litter conditions
  • Price lag

It is important to note where variables correlate with each other as well. Price lag and median income had very strong positive correlation with one another, as did number of bathrooms and number of bedrooms. The number of bedrooms and the number of bathrooms are highly correlated to one another. Therefore to avoid multicollinearity, or two variables telling the same story, we will focus on number of bathrooms since it is also more correlated to sale price.

# str(sales2023.select)

# Create 4 home price correlation scatterplots that you think are of interest (using open data variables)

#MedHHInc, trees, price_lag, litterscore

st_drop_geometry(sales2023.select) %>% 
  dplyr::select(sale_price, MedHHInc, trees, price_lag, litterscore) %>%
  gather(Variable, Value, -sale_price) %>% 
   ggplot(aes(Value, sale_price)) +
     geom_point(size = .5) + geom_smooth(method = "lm", se=F, colour = "#D52F4D") +
  facet_wrap(~Variable, ncol = 2, scales = "free", labeller = labeller(Variable = c(price_lag = "Price Lag", MedHHInc = "Median Income", trees = "Tree Density", litterscore = "Litter Conditions"))) +
     labs(title = "Price as a function of select variables") +
     plotTheme()

The relationship between sale price and median income and tree density are similarly positive, where as litter is negative. The correlation between price lag and sale price is nearly 1, which indicates a very tight fit (as supported by the values in the categorizing summary table.)

Examine Spatially

# Create 3 maps of the 3 most interesting variables.

plot_MedHHInc <- ggplot() +
  geom_sf(data = blocks21.sf, color = "gray70") +
  geom_sf(data = sales2023, aes(colour = q5(MedHHInc)), 
          show.legend = "point", size = .5) +
  scale_colour_manual(values = rev(palette),
                   labels=qBr(sales2023,"MedHHInc"),
                   name="Median\nIncome") +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  mapTheme()

plot_price_lag <- ggplot() +
  geom_sf(data = blocks21.sf, color = "gray70") +
  geom_sf(data = sales2023, aes(colour = q5(price_lag)), 
          show.legend = "point", size = .5) +
  scale_colour_manual(values = rev(palette),
                   labels=qBr(sales2023,"price_lag"),
                   name="Nearby\nPrice") +
  guides(color = guide_legend(override.aes = list(size = 5))) +
  mapTheme()

trees_plot <- ggplot() +
  geom_sf(data = blocks21.sf, color = "gray70") +
  stat_density2d(data = data.frame(st_coordinates(trees.sf)),
                 aes(X, Y, fill = ..level.., alpha = ..level..),
                 size = 0.01, bins = 50, geom = 'polygon') +
  scale_fill_gradient(
    low = "#ABCDA2",
    high = "#4B7740", 
    name = "Tree Density") +
  scale_alpha(range = c(0, .5), guide = "none") +
mapTheme()

plot_grid(plot_MedHHInc, plot_price_lag, trees_plot, align = "v", ncol = 2)