# The Paradox of Housing Supply Growth: Conservative Sprawl and Liberal Infill in the 21st Century United States
# Manuscript published in Housing Policy Debate

# Analysis file
# Yonah Freemark, yfreemark@urban.org
# May 2026

## setwd("/Users/yfreemark/Library/CloudStorage/Box-Box/Blog Posts/2025 - Ideology and housing/")

#### Libraries ####

library(foreign)
require(lmtest)
require(corrplot)
require(RColorBrewer)
require(sandwich)
require(ggplot2)
require(hrbrthemes)
require(scales)
require(psych)
library(stringi)
library(tidyverse)
require(ggeffects)
library(urbnthemes)
library(shiny)
library(sf)
library(urbnmapr)
library(data.table)
library(ggrepel)
# remotes::install_github("UrbanInstitute/urbnthemes")

set_urbn_defaults(style = "print")

#### BEGIN HERE -- DATA LOAD ####

censusCounties <- read.csv("census_counties_full_data.csv", header = TRUE)
cityDatabase <- read.csv("census_cities_full_data.csv", header = TRUE)
censusTracts <- read.csv("census_tracts_full_data.csv", header = TRUE)

# convert MRP data by reversing direction and having all data begin at 0 to avoid negative problem
censusCounties$mrp_ideology_04_11_invert_base_0 <- censusCounties$mrp_ideology_04_11*-1
censusCounties$mrp_ideology_04_11_invert_base_0 <- (censusCounties$mrp_ideology_04_11_invert_base_0+(-1*min(censusCounties$mrp_ideology_04_11_invert_base_0,na.rm=T)))
censusCounties$mrp_ideology_12_16_invert_base_0 <- censusCounties$mrp_ideology_12_16*-1
censusCounties$mrp_ideology_12_16_invert_base_0 <- (censusCounties$mrp_ideology_12_16_invert_base_0+(-1*min(censusCounties$mrp_ideology_12_16_invert_base_0,na.rm=T)))
censusCounties$mrp_ideology_17_21_invert_base_0 <- censusCounties$mrp_ideology_17_21*-1
censusCounties$mrp_ideology_17_21_invert_base_0 <- (censusCounties$mrp_ideology_17_21_invert_base_0+(-1*min(censusCounties$mrp_ideology_17_21_invert_base_0,na.rm=T)))

cityDatabase$mrp_ideology_04_11_invert_base_0 <- cityDatabase$mrp_ideology_04_11*-1
cityDatabase$mrp_ideology_04_11_invert_base_0 <- (cityDatabase$mrp_ideology_04_11_invert_base_0+(-1*min(cityDatabase$mrp_ideology_04_11_invert_base_0,na.rm=T)))
cityDatabase$mrp_ideology_12_16_invert_base_0 <- cityDatabase$mrp_ideology_12_16*-1
cityDatabase$mrp_ideology_12_16_invert_base_0 <- (cityDatabase$mrp_ideology_12_16_invert_base_0+(-1*min(cityDatabase$mrp_ideology_12_16_invert_base_0,na.rm=T)))
cityDatabase$mrp_ideology_17_21_invert_base_0 <- cityDatabase$mrp_ideology_17_21*-1
cityDatabase$mrp_ideology_17_21_invert_base_0 <- (cityDatabase$mrp_ideology_17_21_invert_base_0+(-1*min(cityDatabase$mrp_ideology_17_21_invert_base_0,na.rm=T)))

countyLarge <- subset(censusCounties, censusCounties$A41AA2000>50000)
cityLarge <- subset(cityDatabase, cityDatabase$units_2000>25000)

#### APPENDIX TABLE 1 (MSAs BY SHARE TRACTS IN EACH GROUP) ####

msaList <- unique(subset(censusCounties$NAME, !is.na(censusCounties$NAME)))

msaShare <- data.frame(matrix(ncol=9, nrow=length(msaList)))
colnames(msaShare) <- c("MSA","Units 2000","Units 2020","Share Units in Tracts <500","500 to 999","1000 to 1999","2000 to 2999","3000 to 4999","5000 Plus")
msaShare[,1] <- msaList

for(i in 1:nrow(msaShare)){
  group <- subset(censusCounties, censusCounties$NAME==msaShare[i,1])
  msaShare[i,2] <- sum(group$A41AA2000, na.rm=T)
  msaShare[i,3] <- sum(group$A41AA2020, na.rm=T)
  msaShare[i,4] <- sum(group$hsg_2000_in_tracts_below_500, na.rm=T)/msaShare[i,2]
  msaShare[i,5] <- (sum(group$hsg_2000_in_tracts_above_500, na.rm=T)-sum(group$hsg_2000_in_tracts_above_1000, na.rm=T))/msaShare[i,2]
  msaShare[i,6] <- (sum(group$hsg_2000_in_tracts_above_1000, na.rm=T)-sum(group$hsg_2000_in_tracts_above_2000, na.rm=T))/msaShare[i,2]
  msaShare[i,7] <- (sum(group$hsg_2000_in_tracts_above_2000, na.rm=T)-sum(group$hsg_2000_in_tracts_above_3000, na.rm=T))/msaShare[i,2]
  msaShare[i,8] <- (sum(group$hsg_2000_in_tracts_above_3000, na.rm=T)-sum(group$hsg_2000_in_tracts_above_5000, na.rm=T))/msaShare[i,2]
  msaShare[i,9] <- (sum(group$hsg_2000_in_tracts_above_5000, na.rm=T))/msaShare[i,2]
}

msaShare <- msaShare[order(-msaShare$`Units 2000`),]
write.csv(msaShare[1:21,],"appendix_table_1.csv",row.names = T)

#### APPENDIX TABLE 2 (MSAs BY SHARE AREA IN EACH GROUP) ####

msaList <- unique(subset(censusCounties$NAME, !is.na(censusCounties$NAME)))

msaShare <- data.frame(matrix(ncol=9, nrow=length(msaList)))
colnames(msaShare) <- c("MSA","Area","Share Area in Tracts <500","500 to 999","1000 to 1999","2000 to 2999","3000 to 4999","5000 Plus","Units 2000")
msaShare[,1] <- msaList

for(i in 1:nrow(msaShare)){
  group <- subset(censusCounties, censusCounties$NAME==msaShare[i,1])
  msaShare[i,2] <- sum(group$area_sq_mi, na.rm=T)
  msaShare[i,3] <- sum(group$area_in_tracts_below_500, na.rm=T)/msaShare[i,2]
  msaShare[i,4] <- (sum(group$area_in_tracts_above_500, na.rm=T)-sum(group$area_in_tracts_above_1000, na.rm=T))/msaShare[i,2]
  msaShare[i,5] <- (sum(group$area_in_tracts_above_1000, na.rm=T)-sum(group$area_in_tracts_above_2000, na.rm=T))/msaShare[i,2]
  msaShare[i,6] <- (sum(group$area_in_tracts_above_2000, na.rm=T)-sum(group$area_in_tracts_above_3000, na.rm=T))/msaShare[i,2]
  msaShare[i,7] <- (sum(group$area_in_tracts_above_3000, na.rm=T)-sum(group$area_in_tracts_above_5000, na.rm=T))/msaShare[i,2]
  msaShare[i,8] <- (sum(group$area_in_tracts_above_5000, na.rm=T))/msaShare[i,2]
  msaShare[i,9] <- sum(group$A41AA2000, na.rm=T)
}

msaShare <- msaShare[order(-msaShare$`Units 2000`),]
write.csv(msaShare[1:21,],"appendix_table_2.csv",row.names = T)

#### DISSIMILARITY INDEX - COUNTIES ####

disMat <- data.frame(matrix(ncol=20, nrow=length(unique(countyLarge$GISJOIN))))
colnames(disMat) <- c("GISJOIN","State","County","mrp_ideology_04_11_invert_base_0","hsg_density_2000","race_wht_alone","median_value","dissimilarity","StateNum","CountyNum","unit_change_2000_20","units_2000","mrp_ideology_17_21_invert_base_0","hsg_density_2010","race_wht_alone_2010","median_value_2010","dissim0010","dissim1020","state_legislative0010","state_legislative1020")
disMat[,1] <- unique(countyLarge$GISJOIN)

for(i in 1:nrow(disMat)){
  group <- subset(countyLarge, countyLarge$GISJOIN==disMat[i,1])
  disMat[i,2] <- group$STATE
  disMat[i,9] <- group$STATEFP
  disMat[i,3] <- group$COUNTY
  disMat[i,10] <- group$COUNTYFP
  disMat[i,11] <- group$chg_00_20
  disMat[i,12] <- group$A41AA2000
  disMat[i,4] <- group$mrp_ideology_04_11_invert_base_0
  disMat[i,5] <- group$hsg_density_2000
  disMat[i,6] <- group$race_wht_alone
  disMat[i,7] <- group$median_value
  disMat[i,13] <- group$mrp_ideology_17_21_invert_base_0
  disMat[i,14] <- group$A41AA2010/group$area_sq_mi
  disMat[i,15] <- group$race_wht_alone_2010
  disMat[i,16] <- group$median_value_2010 
  disMat[i,19] <- as.character(group$partisan_leg_2000s)
  disMat[i,20] <- as.character(group$partisan_leg_2010s)
  dissim <- 0
  dissim0010 <- 0
  dissim1020 <- 0
  groupB <- subset(censusTracts, censusTracts$COUNTYA==disMat[i,10] & censusTracts$STATEA==disMat[i,9])
  if(nrow(groupB)>0){
    groupB$newUnits <- pmax(groupB$units_2020-groupB$units_2000,0) # Not allowing non-negative
    groupB$newUnits0010 <- pmax(groupB$units_2010-groupB$units_2000,0) # Not allowing non-negative
    groupB$newUnits1020 <- pmax(groupB$units_2020-groupB$units_2010,0) # Not allowing non-negative
    for(j in 1:nrow(groupB)){
      dissim <- sum(dissim, abs((groupB$newUnits[j]/sum(groupB$newUnits,na.rm=T)) - (groupB$units_2000[j]/sum(groupB$units_2000,na.rm=T))),na.rm=T)
      dissim0010 <- sum(dissim0010, abs((groupB$newUnits0010[j]/sum(groupB$newUnits0010,na.rm=T)) - (groupB$units_2000[j]/sum(groupB$units_2000,na.rm=T))),na.rm=T)
      dissim1020 <- sum(dissim1020, abs((groupB$newUnits1020[j]/sum(groupB$newUnits1020,na.rm=T)) - (groupB$units_2010[j]/sum(groupB$units_2010,na.rm=T))),na.rm=T)
    }
    disMat[i,8] <- 0.5 * dissim # population-weighted dissimilarity index
    disMat[i,17] <- 0.5 * dissim0010 # population-weighted dissimilarity index
    disMat[i,18] <- 0.5 * dissim1020 # population-weighted dissimilarity index
  }}

# If dissimilarity = 0 : New units are added in perfect proportion to existing units
# If dissimilarity = 1 : New units are added only in tracts with no existing units (or vice versa).

### DISSIMMILARITY INDEX: PANEL - COUNTIES ####

dissPanel <- data.frame(matrix(ncol=8,nrow=2*nrow(disMat)))
colnames(dissPanel) <- c("GISJOIN","Period","Density","Ideology","StatePartisan","White_Share","Median_Value","Dissimilarity")
dissPanel[,1] <- rep(disMat$GISJOIN,2)
dissPanel[,2] <- c(rep("Period0010",nrow(disMat)),rep("Period1020",nrow(disMat)))
for(i in 1:nrow(dissPanel)){
  group <- subset(disMat, disMat$GISJOIN==dissPanel[i,1])
  ifelse(dissPanel[i,2]=="Period0010",{
    dissPanel$Density[i] <- group$hsg_density_2000
    dissPanel$Ideology[i] <- group$mrp_ideology_04_11_invert_base_0
    dissPanel$StatePartisan[i] <- group$state_legislative0010
    dissPanel$White_Share[i] <- group$race_wht_alone
    dissPanel$Median_Value[i] <- group$median_value*1.28 # inflation adjusted to 2010
    dissPanel$Dissimilarity[i] <- group$dissim0010
  },{
    dissPanel$Density[i] <- group$hsg_density_2010
    dissPanel$Ideology[i] <- group$mrp_ideology_17_21_invert_base_0
    dissPanel$StatePartisan[i] <- group$state_legislative1020
    dissPanel$White_Share[i] <- group$race_wht_alone_2010
    dissPanel$Median_Value[i] <- group$median_value_2010 # inflation
    dissPanel$Dissimilarity[i] <- group$dissim1020
  })
}

#### TABLE 2, MODEL V (DISSIMILARITY) - COUNTIES ####

# Using - interaction
regI <- (lm(Dissimilarity ~ I((Density/1000) * Ideology) + I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(Median_Value) + as.factor(Period), data=dissPanel)) #
summary(regI)
nrow(dissPanel)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = dissPanel)) # Robust SEs, clustered on counties

#### DISSIMILARITY INDEX - CITIES ####

disMatCities <- data.frame(matrix(ncol=20, nrow=nrow(cityLarge)))
colnames(disMatCities) <- c("GISJOIN","State","City","mrp_ideology_04_11_invert_base_0","hsg_density_2000","race_wht_alone","median_value","dissimilarity","StateNum","CityNum","unit_change_2000_20","units_2000","mrp_ideology_17_21_invert_base_0","hsg_density_2010","race_wht_alone_2010","median_value_2010","dissim0010","dissim1020","state_legislative0010","state_legislative1020")
disMatCities[,1] <- unique(cityLarge$GISJOIN)

for(i in 1:nrow(disMatCities)){
  group <- subset(cityLarge, cityLarge$GISJOIN==disMatCities[i,1])
  disMatCities[i,2] <- group$STATE
  disMatCities[i,9] <- group$STATEFP10
  disMatCities[i,3] <- group$NAME10
  disMatCities[i,10] <- group$GEOID10
  disMatCities[i,11] <- group$chg_00_20
  disMatCities[i,12] <- group$units_2000
  disMatCities[i,4] <- group$mrp_ideology_04_11_invert_base_0
  disMatCities[i,5] <- group$hsg_density_2000
  disMatCities[i,6] <- group$race_wht_alone
  disMatCities[i,7] <- group$median_value
  disMatCities[i,13] <- group$mrp_ideology_17_21_invert_base_0
  disMatCities[i,14] <- group$units_2010/(group$ALAND10/2589988.11)
  disMatCities[i,15] <- group$race_wht_alone_2010
  disMatCities[i,16] <- group$median_value_2010 
  disMatCities[i,19] <- as.character(group$partisan_leg_2000s)
  disMatCities[i,20] <- as.character(group$partisan_leg_2010s)
  dissim <- 0
  dissim0010 <- 0
  dissim1020 <- 0
  groupB <- subset(censusTracts, censusTracts$PLC_GEOID==disMatCities[i,10])
  if(nrow(groupB)>0){
    groupB$newUnits <- pmax(groupB$units_2020-groupB$units_2000,0) # Not allowing non-negative
    groupB$newUnits0010 <- pmax(groupB$units_2010-groupB$units_2000,0) # Not allowing non-negative
    groupB$newUnits1020 <- pmax(groupB$units_2020-groupB$units_2010,0) # Not allowing non-negative
    for(j in 1:nrow(groupB)){
      dissim <- sum(dissim, abs((groupB$newUnits[j]/sum(groupB$newUnits,na.rm=T)) - (groupB$units_2000[j]/sum(groupB$units_2000,na.rm=T))),na.rm=T)
      dissim0010 <- sum(dissim0010, abs((groupB$newUnits0010[j]/sum(groupB$newUnits0010,na.rm=T)) - (groupB$units_2000[j]/sum(groupB$units_2000,na.rm=T))),na.rm=T)
      dissim1020 <- sum(dissim1020, abs((groupB$newUnits1020[j]/sum(groupB$newUnits1020,na.rm=T)) - (groupB$units_2010[j]/sum(groupB$units_2010,na.rm=T))),na.rm=T)
    }
    disMatCities[i,8] <- 0.5 * dissim # population-weighted dissimilarity index
    disMatCities[i,17] <- 0.5 * dissim0010 # population-weighted dissimilarity index
    disMatCities[i,18] <- 0.5 * dissim1020 # population-weighted dissimilarity index
  }}

# If dissimilarity = 0 : New units are added in perfect proportion to existing units
# If dissimilarity = 1 : New units are added only in tracts with no existing units (or vice versa).

### DISSIMMILARITY INDEX: PANEL - CITIES ####

dissPanel <- data.frame(matrix(ncol=8,nrow=2*nrow(disMatCities)))
colnames(dissPanel) <- c("GISJOIN","Period","Density","Ideology","StatePartisan","White_Share","Median_Value","Dissimilarity")
dissPanel[,1] <- rep(disMatCities$GISJOIN,2)
dissPanel[,2] <- c(rep("Period0010",nrow(disMatCities)),rep("Period1020",nrow(disMatCities)))
for(i in 1:nrow(dissPanel)){
  group <- subset(disMatCities, disMatCities$GISJOIN==dissPanel[i,1])
  ifelse(dissPanel[i,2]=="Period0010",{
    dissPanel$Density[i] <- group$hsg_density_2000
    dissPanel$Ideology[i] <- group$mrp_ideology_04_11_invert_base_0
    dissPanel$StatePartisan[i] <- group$state_legislative0010
    dissPanel$White_Share[i] <- group$race_wht_alone
    dissPanel$Median_Value[i] <- group$median_value*1.28 # inflation adjusted to 2010
    dissPanel$Dissimilarity[i] <- group$dissim0010
  },{
    dissPanel$Density[i] <- group$hsg_density_2010
    dissPanel$Ideology[i] <- group$mrp_ideology_17_21_invert_base_0
    dissPanel$StatePartisan[i] <- group$state_legislative1020
    dissPanel$White_Share[i] <- group$race_wht_alone_2010
    dissPanel$Median_Value[i] <- group$median_value_2010 # inflation
    dissPanel$Dissimilarity[i] <- group$dissim1020
  })
}

#### TABLE 2, MODEL VI (DISSIMILARITY) - CITIES ####

# Using - interactions
regI <- (lm(Dissimilarity ~ I((Density/1000) * Ideology) + I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(as.numeric(Median_Value)) + as.factor(Period), data=dissPanel)) #
summary(regI)
nrow(dissPanel)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = dissPanel)) # Robust SEs, clustered on cities

#### COUNTY PANEL ####

# Define all combinations of GISJOIN, Period, and Change_Type
unique_GISJOIN <- unique(countyLarge$GISJOIN)
periods <- c("period_2000_2010", "period_2010_2020")
change_types <- c(
  "change_overall", "change_below_1000", "change_1000_1999", "change_1000_plus",
  "change_1000_2999", "change_2000_2999", "change_3000_4999", "change_3000_plus", "change_5000_plus"
)

# Create all combinations
altMat <- expand.grid(GISJOIN = unique_GISJOIN, Period = periods, Change_Type = change_types)

# Merge with countyLarge
altMat <- altMat %>%
  left_join(
    countyLarge %>%
      select(
        GISJOIN, race_wht_alone, race_blk_alone, median_value, median_income, median_rent,
        built_1990_2000, BA_plus, mrp_ideology_04_11_invert_base_0, partisan_leg_2000s, partisan_overall_2000s,
        hsg_density_2000, A41AA2000, A41AA2010, A41AA2020, area_sq_mi,
        race_wht_alone_2010, race_blk_alone_2010, median_value_2010, median_income_2010, median_rent_2010,
        built_2000_2010, BA_plus_2010, mrp_ideology_17_21_invert_base_0, partisan_leg_2010s, partisan_overall_2010s,
        chg_00_10_in_tracts_below_1000_pct, chg_00_10_in_tracts_above_1000_pct, chg_00_10_in_tracts_1000_1999_pct,
        chg_00_10_in_tracts_1000_2999_pct, chg_00_10_in_tracts_2000_2999_pct, chg_00_10_in_tracts_3000_4999_pct,
        chg_00_10_in_tracts_above_3000_pct, chg_00_10_in_tracts_above_5000_pct,
        chg_10_20_in_tracts_below_1000_pct, chg_10_20_in_tracts_above_1000_pct, chg_10_20_in_tracts_1000_1999_pct,
        chg_10_20_in_tracts_1000_2999_pct, chg_10_20_in_tracts_2000_2999_pct, chg_10_20_in_tracts_3000_4999_pct,
        chg_10_20_in_tracts_above_3000_pct, chg_10_20_in_tracts_above_5000_pct
      ),
    by = "GISJOIN"
  )

# Calculate Change_Value and other variables
altMat <- altMat %>%
  mutate(
    # 2000-2010 period
    race_wht_alone = ifelse(Period == "period_2000_2010", race_wht_alone, race_wht_alone_2010),
    race_blk_alone = ifelse(Period == "period_2000_2010", race_blk_alone, race_blk_alone_2010),
    median_value = ifelse(Period == "period_2000_2010", median_value * 1.28, median_value_2010),
    median_income = ifelse(Period == "period_2000_2010", median_income * 1.28, median_income_2010),
    median_rent = ifelse(Period == "period_2000_2010", median_rent * 1.28, median_rent_2010),
    built_last_decade = ifelse(Period == "period_2000_2010", built_1990_2000, built_2000_2010),
    BA_plus = ifelse(Period == "period_2000_2010", BA_plus, BA_plus_2010),
    mrp_ideology = ifelse(Period == "period_2000_2010", mrp_ideology_04_11_invert_base_0, mrp_ideology_17_21_invert_base_0),
    hsg_density = ifelse(Period == "period_2000_2010", hsg_density_2000, A41AA2010 / area_sq_mi),
    partisan_leg = ifelse(Period == "period_2000_2010", as.character(partisan_leg_2000s), as.character(partisan_leg_2010s)),
    partisan_overall = ifelse(Period == "period_2000_2010", as.character(partisan_overall_2000s), as.character(partisan_overall_2010s)),
    
    # Change_Value logic
    Change_Value = case_when(
      Period == "period_2000_2010" & Change_Type == "change_overall" ~ (A41AA2010 - A41AA2000) / A41AA2000,
      Period == "period_2000_2010" & Change_Type == "change_below_1000" ~ chg_00_10_in_tracts_below_1000_pct,
      Period == "period_2000_2010" & Change_Type == "change_1000_plus" ~ chg_00_10_in_tracts_above_1000_pct,
      Period == "period_2000_2010" & Change_Type == "change_1000_1999" ~ chg_00_10_in_tracts_1000_1999_pct,
      Period == "period_2000_2010" & Change_Type == "change_1000_2999" ~ chg_00_10_in_tracts_1000_2999_pct,
      Period == "period_2000_2010" & Change_Type == "change_2000_2999" ~ chg_00_10_in_tracts_2000_2999_pct,
      Period == "period_2000_2010" & Change_Type == "change_3000_4999" ~ chg_00_10_in_tracts_3000_4999_pct,
      Period == "period_2000_2010" & Change_Type == "change_3000_plus" ~ chg_00_10_in_tracts_above_3000_pct,
      Period == "period_2000_2010" & Change_Type == "change_5000_plus" ~ chg_00_10_in_tracts_above_5000_pct,
      Period == "period_2010_2020" & Change_Type == "change_overall" ~ (A41AA2020 - A41AA2010) / A41AA2010,
      Period == "period_2010_2020" & Change_Type == "change_below_1000" ~ chg_10_20_in_tracts_below_1000_pct,
      Period == "period_2010_2020" & Change_Type == "change_1000_plus" ~ chg_10_20_in_tracts_above_1000_pct,
      Period == "period_2010_2020" & Change_Type == "change_1000_1999" ~ chg_10_20_in_tracts_1000_1999_pct,
      Period == "period_2010_2020" & Change_Type == "change_1000_2999" ~ chg_10_20_in_tracts_1000_2999_pct,
      Period == "period_2010_2020" & Change_Type == "change_2000_2999" ~ chg_10_20_in_tracts_2000_2999_pct,
      Period == "period_2010_2020" & Change_Type == "change_3000_4999" ~ chg_10_20_in_tracts_3000_4999_pct,
      Period == "period_2010_2020" & Change_Type == "change_3000_plus" ~ chg_10_20_in_tracts_above_3000_pct,
      Period == "period_2010_2020" & Change_Type == "change_5000_plus" ~ chg_10_20_in_tracts_above_5000_pct,
      TRUE ~ NA_real_
    )
  ) %>%
  select(
    GISJOIN, Period, Change_Type, Change_Value, race_wht_alone, race_blk_alone,
    median_value, median_income, median_rent, built_last_decade, BA_plus,
    hsg_density, mrp_ideology, partisan_leg, partisan_overall
  )

#### APPENDIX TABLE 4 — COUNTIES ####

# Model 1 -- all
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_overall")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_overall"))) # Robust SEs
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_overall"))
nrow(subset(altMat, altMat$Change_Type=="change_overall")) # Note -- this doesn't actually exclude any tracts! Nothing has this problem.

# Model 2 -- <1000
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_below_1000")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_below_1000")))
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_below_1000"))

# Model 3 -- 1000-1999
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-2) & altMat$Change_Value < (2) & altMat$Change_Type=="change_1000_1999")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-2) & altMat$Change_Value < (2) & altMat$Change_Type=="change_1000_1999")))
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_1000_1999"))

# Model 4 -- 2000-2999
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_2000_2999")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_2000_2999")))
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_2000_2999"))

# Model 5 -- 3000-4999
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_3000_4999")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_3000_4999")))
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_3000_4999"))

# Model 6 -- 5000 plus
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_5000_plus")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_5000_plus")))
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_5000_plus"))

#### TABLE 1 — COUNTIES ####

# Model 1 -- all, intersected
regI <- (lm(Change_Value ~ I((hsg_density/1000) * mrp_ideology) + I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_overall")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_overall")))
nrow(subset(altMat, !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2) & altMat$Change_Type=="change_overall"))

#### CITY PANEL ####

# Define all combinations of GISJOIN, Period, and Change_Type
unique_GISJOIN <- unique(cityLarge$GISJOIN)
periods <- c("period_2000_2010", "period_2010_2020")
change_types <- c(
  "change_overall", "change_below_1000", "change_1000_1999", "change_1000_plus",
  "change_1000_2999", "change_2000_2999", "change_3000_4999", "change_3000_plus", "change_5000_plus"
)

# Create all combinations
altMatCity <- expand.grid(GISJOIN = unique_GISJOIN, Period = periods, Change_Type = change_types)

# Merge with cityLarge
altMatCity <- altMatCity %>%
  left_join(
    cityLarge %>%
      select(
        GISJOIN, race_wht_alone, race_blk_alone, median_value, median_income, median_rent,
        built_1990_2000, BA_plus, mrp_ideology_04_11_invert_base_0, dem_share_council_2000_2009, dem_share_council_2010_2019,
        partisan_leg_2000s, partisan_overall_2000s,
        hsg_density_2000, hsg_density_2010, units_2000, units_2010, units_2020,
        race_wht_alone_2010, race_blk_alone_2010, median_value_2010, median_income_2010, median_rent_2010,
        built_2000_2010, BA_plus_2010, mrp_ideology_17_21_invert_base_0, partisan_leg_2010s, partisan_overall_2010s,
        chg_00_10_in_tracts_below_1000_pct, chg_00_10_in_tracts_above_1000_pct, chg_00_10_in_tracts_1000_1999_pct,
        chg_00_10_in_tracts_1000_2999_pct, chg_00_10_in_tracts_2000_2999_pct, chg_00_10_in_tracts_3000_4999_pct,
        chg_00_10_in_tracts_above_3000_pct, chg_00_10_in_tracts_above_5000_pct,
        chg_10_20_in_tracts_below_1000_pct, chg_10_20_in_tracts_above_1000_pct, chg_10_20_in_tracts_1000_1999_pct,
        chg_10_20_in_tracts_1000_2999_pct, chg_10_20_in_tracts_2000_2999_pct, chg_10_20_in_tracts_3000_4999_pct,
        chg_10_20_in_tracts_above_3000_pct, chg_10_20_in_tracts_above_5000_pct
      ),
    by = "GISJOIN"
  )

# Calculate Change_Value and other variables
altMatCity <- altMatCity %>%
  mutate(
    # 2000-2010 period
    race_wht_alone = ifelse(Period == "period_2000_2010", race_wht_alone, race_wht_alone_2010),
    race_blk_alone = ifelse(Period == "period_2000_2010", race_blk_alone, race_blk_alone_2010),
    median_value = ifelse(Period == "period_2000_2010", median_value * 1.28, median_value_2010),
    median_income = ifelse(Period == "period_2000_2010", median_income * 1.28, median_income_2010),
    median_rent = ifelse(Period == "period_2000_2010", median_rent * 1.28, median_rent_2010),
    built_last_decade = ifelse(Period == "period_2000_2010", built_1990_2000, built_2000_2010),
    BA_plus = ifelse(Period == "period_2000_2010", BA_plus, BA_plus_2010),
    mrp_ideology = ifelse(Period == "period_2000_2010", mrp_ideology_04_11_invert_base_0, mrp_ideology_17_21_invert_base_0),
    hsg_density = ifelse(Period == "period_2000_2010", hsg_density_2000, hsg_density_2010),
    dem_share_council = ifelse(Period == "period_2000_2010", dem_share_council_2000_2009, dem_share_council_2010_2019),
    partisan_leg = ifelse(Period == "period_2000_2010", as.character(partisan_leg_2000s), as.character(partisan_leg_2010s)),
    partisan_overall = ifelse(Period == "period_2000_2010", as.character(partisan_overall_2000s), as.character(partisan_overall_2010s)),
    
    # Change_Value logic
    Change_Value = case_when(
      Period == "period_2000_2010" & Change_Type == "change_overall" ~ (units_2010 - units_2000) / units_2000,
      Period == "period_2000_2010" & Change_Type == "change_below_1000" ~ chg_00_10_in_tracts_below_1000_pct,
      Period == "period_2000_2010" & Change_Type == "change_1000_plus" ~ chg_00_10_in_tracts_above_1000_pct,
      Period == "period_2000_2010" & Change_Type == "change_1000_1999" ~ chg_00_10_in_tracts_1000_1999_pct,
      Period == "period_2000_2010" & Change_Type == "change_1000_2999" ~ chg_00_10_in_tracts_1000_2999_pct,
      Period == "period_2000_2010" & Change_Type == "change_2000_2999" ~ chg_00_10_in_tracts_2000_2999_pct,
      Period == "period_2000_2010" & Change_Type == "change_3000_4999" ~ chg_00_10_in_tracts_3000_4999_pct,
      Period == "period_2000_2010" & Change_Type == "change_3000_plus" ~ chg_00_10_in_tracts_above_3000_pct,
      Period == "period_2000_2010" & Change_Type == "change_5000_plus" ~ chg_00_10_in_tracts_above_5000_pct,
      Period == "period_2010_2020" & Change_Type == "change_overall" ~ (units_2020 - units_2010) / units_2010,
      Period == "period_2010_2020" & Change_Type == "change_below_1000" ~ chg_10_20_in_tracts_below_1000_pct,
      Period == "period_2010_2020" & Change_Type == "change_1000_plus" ~ chg_10_20_in_tracts_above_1000_pct,
      Period == "period_2010_2020" & Change_Type == "change_1000_1999" ~ chg_10_20_in_tracts_1000_1999_pct,
      Period == "period_2010_2020" & Change_Type == "change_1000_2999" ~ chg_10_20_in_tracts_1000_2999_pct,
      Period == "period_2010_2020" & Change_Type == "change_2000_2999" ~ chg_10_20_in_tracts_2000_2999_pct,
      Period == "period_2010_2020" & Change_Type == "change_3000_4999" ~ chg_10_20_in_tracts_3000_4999_pct,
      Period == "period_2010_2020" & Change_Type == "change_3000_plus" ~ chg_10_20_in_tracts_above_3000_pct,
      Period == "period_2010_2020" & Change_Type == "change_5000_plus" ~ chg_10_20_in_tracts_above_5000_pct,
      TRUE ~ NA_real_
    )
  ) %>%
  select(
    GISJOIN, Period, Change_Type, Change_Value, race_wht_alone, race_blk_alone,
    median_value, median_income, median_rent, built_last_decade, BA_plus,
    hsg_density, mrp_ideology, partisan_leg, partisan_overall,
    dem_share_council
  )

#### APPENDIX TABLE 3 - CITIES ####

# Model 1 -- all
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall"))-3

# Model 2 -- <1000
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_below_1000")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_below_1000")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_below_1000"))-5

# Model 3 -- 1000-1999
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-2) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_1000_1999")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-2) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_1000_1999")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_1000_1999"))-5

# Model 4 -- 2000-2999
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_2000_2999")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_2000_2999")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_2000_2999"))-2

# Model 5 -- 3000-4999
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_3000_4999")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_3000_4999")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_3000_4999"))

# Model 6 -- 5000 plus
regI <- (lm(Change_Value ~ I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_5000_plus")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_5000_plus")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_5000_plus"))

#### TABLE 1 - CITIES ####

# Model 2 -- intersected
regI <- (lm(Change_Value ~ I((hsg_density/1000) * mrp_ideology) + I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall"))-3

# (for background -- not using -- excluding cities that had a negative change in housing supply)
regI <- (lm(Change_Value ~ I((hsg_density/1000) * mrp_ideology) + I(hsg_density/1000) + mrp_ideology + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value >= (0) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value >= 0 & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value >= (0) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall"))-2

# Model 3 -- all, intersected, with council council instead of ideology
regI <- (lm(Change_Value ~ I((hsg_density/1000) * dem_share_council) + I(hsg_density/1000) + dem_share_council + as.factor(partisan_leg) + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall")))
nrow(subset(altMatCity, !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2) & altMatCity$Change_Type=="change_overall"))-488

#### TRACT PANEL - In cities with >25,000 housing units ####

# Add tract-level data

tract2010Adjust <- read.csv("Additional Political Data/Tracts 2000 to 2010/Tracts_2000_2010_data.csv", header = T)
tract2010Adjust$tractpopWhtAlone2000shr <- as.numeric(tract2010Adjust$popWhtAlone2000)/as.numeric(tract2010Adjust$popTotal2000)
tract2010Adjust$tractpopWhtAlone2010shr <- as.numeric(tract2010Adjust$popWhtAlone2010)/as.numeric(tract2010Adjust$popTotal2010)
tract2010Adjust$tracteduBA2010shr <- as.numeric(tract2010Adjust$eduBA2010)/as.numeric(tract2010Adjust$eduTotal2010)
tract2010Adjust$tracteduBA2000shr <- as.numeric(tract2010Adjust$eduBA2000)/as.numeric(tract2010Adjust$eduTotal2000)
tract2010Adjust$tractmedHHinc2010 <- as.numeric(tract2010Adjust$medHHinc2010)
tract2010Adjust$tractmedHHinc2000 <- as.numeric(tract2010Adjust$medHHinc2000)
tract2010Adjust$tractmedRent2010 <- as.numeric(tract2010Adjust$medRent2010)
tract2010Adjust$tractmedRent2000 <- as.numeric(tract2010Adjust$medRent2000)
tract2010Adjust$tractmedValue2000 <- as.numeric(tract2010Adjust$medValue2000)
tract2010Adjust$tractmedValue2010 <- as.numeric(tract2010Adjust$medValue2010)
tract2010Adjust$tractbuilt20002010shr <- as.numeric(tract2010Adjust$built20002010)/as.numeric(tract2010Adjust$builtTotal2010)
tract2010Adjust$tractbuilt19902000shr <- as.numeric(tract2010Adjust$built19902000)/as.numeric(tract2010Adjust$builtTotal2000)

# Add tract-level data to panel

censusTractsLargeCity <- subset(censusTracts, censusTracts$PLC_GJ %in% subset(cityDatabase$PLC_GJ, cityDatabase$units_2000>25000))

# Pre-compute unique GISJOINs and periods
unique_GISJOIN <- unique(censusTractsLargeCity$GISJOIN)
periods <- c("period_2000_2010", "period_2010_2020")

# Create a data frame with all combinations of GISJOIN and period
altMatTract <- expand.grid(GISJOIN = unique_GISJOIN, Period = periods)

# Merge with cityDatabase and tract2010Adjust
altMatTract <- altMatTract %>%
  left_join(
    censusTractsLargeCity %>%
      select(GISJOIN, PLC_GJ, units_2000, units_2010, units_2020, ALAND10),
    by = "GISJOIN"
  ) %>%
  left_join(
    cityDatabase %>%
      select(PLC_GJ, race_wht_alone, race_blk_alone, median_value, median_income, median_rent,
             built_1990_2000, BA_plus, mrp_ideology_04_11_invert_base_0, partisan_leg_2000s, partisan_overall_2000s,
             dem_share_council_2000_2009, dem_share_council_2010_2019,
             hsg_density_2000, hsg_density_2010,
             race_wht_alone_2010, race_blk_alone_2010, median_value_2010, median_income_2010, median_rent_2010,
             built_2000_2010, BA_plus_2010, mrp_ideology_17_21_invert_base_0, partisan_leg_2010s, partisan_overall_2010s),
    by = c("PLC_GJ" = "PLC_GJ")
  ) %>%
  left_join(
    tract2010Adjust,
    by = "GISJOIN"
  )

# Calculate Change_Value and hsg_density_tract
altMatTract <- altMatTract %>%
  mutate(
    Change_Value = case_when(
      Period == "period_2000_2010" ~ (units_2010 - units_2000) / units_2000,
      Period == "period_2010_2020" ~ (units_2020 - units_2010) / units_2010,
      TRUE ~ NA_real_
    ),
    Change_Value_Nominal = case_when(
      Period == "period_2000_2010" ~ (units_2010 - units_2000),
      Period == "period_2010_2020" ~ (units_2020 - units_2010),
      TRUE ~ NA_real_
    ),
    hsg_density_tract = case_when(
      Period == "period_2000_2010" ~ units_2000 / (ALAND10 / 2589988.11),
      Period == "period_2010_2020" ~ units_2010 / (ALAND10 / 2589988.11),
      TRUE ~ NA_real_
    )
  )

# Adjust for inflation and select columns
altMatTract <- altMatTract %>%
  mutate(
    # City-level variables
    housing_density_city = ifelse(Period == "period_2000_2010", hsg_density_2000, hsg_density_2010),
    race_wht_alone_city = ifelse(Period == "period_2000_2010", race_wht_alone, race_wht_alone_2010),
    race_blk_alone_city = ifelse(Period == "period_2000_2010", race_blk_alone, race_blk_alone_2010),
    median_value_city = ifelse(Period == "period_2000_2010", median_value * 1.28, median_value_2010),
    median_income_city = ifelse(Period == "period_2000_2010", median_income * 1.28, median_income_2010),
    median_rent_city = ifelse(Period == "period_2000_2010", median_rent * 1.28, median_rent_2010),
    built_last_decade_city = ifelse(Period == "period_2000_2010", built_1990_2000, built_2000_2010),
    BA_plus_city = ifelse(Period == "period_2000_2010", BA_plus, BA_plus_2010),
    mrp_ideology_city = ifelse(Period == "period_2000_2010", mrp_ideology_04_11_invert_base_0, mrp_ideology_17_21_invert_base_0),
    dem_share_council_city = ifelse(Period == "period_2000_2010", dem_share_council_2000_2009, dem_share_council_2010_2019),
    partisan_leg_state = ifelse(Period == "period_2000_2010", as.character(partisan_leg_2000s), as.character(partisan_leg_2010s)),
    partisan_overall_state = ifelse(Period == "period_2000_2010", as.character(partisan_overall_2000s), as.character(partisan_overall_2010s)),
    
    # Tract-level variables
    units_tract_begin = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", units_2000, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", units_2010, NA)),
    units_tract_end = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", units_2010, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", units_2020, NA)),
    race_wht_alone_tract = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", tractpopWhtAlone2000shr, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", tractpopWhtAlone2010shr, NA)),
    median_value_tract = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", tractmedValue2000 * 1.28, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", tractmedValue2010, NA)),
    median_income_tract = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", tractmedHHinc2000 * 1.28, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", tractmedHHinc2010, NA)),
    median_rent_tract = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", tractmedRent2000 * 1.28, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", tractmedRent2000, NA)),
    built_last_decade_tract = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", tractbuilt19902000shr, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", tractbuilt20002010shr, NA)),
    BA_plus_tract = ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2000_2010", tracteduBA2000shr, ifelse(nrow(tract2010Adjust) > 0 & Period == "period_2010_2020", tracteduBA2010shr, NA))
  ) %>%
  select(
    GISJOIN, Period, City = PLC_GJ, Change_Value, Change_Value_Nominal, race_wht_alone_city, race_blk_alone_city,
    median_value_city, median_income_city, median_rent_city, built_last_decade_city, BA_plus_city,
    hsg_density_tract, mrp_ideology_city, partisan_leg_state, partisan_overall_state,
    housing_density_city,
    race_wht_alone_tract, median_value_tract, median_income_tract, median_rent_tract,
    built_last_decade_tract, BA_plus_tract, dem_share_council_city,
    units_tract_begin, units_tract_end
  )

#### FIGURE 7 / APPENDIX FIGURE 8 / APPENDIX FIGURE 9 - TRACTS BY CITY REGRESSIONS - BINS AND PLOT ####

# Convert to data.table for faster operations
setDT(altMatTract)

# Pre-allocate binTab
binTab <- data.table(
  min = c(1, 1, 500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 15000),
  max = c(500000, 500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 15000, 20000),
  coef = NA_real_,
  se = NA_real_,
  p_val = NA_real_,
  num_cities = NA_real_
)

# Iterate over bins
for (i in 1:nrow(binTab)) {
  # Subset once per bin
  group <- altMatTract[
    hsg_density_tract > binTab$min[i] & hsg_density_tract < binTab$max[i]
  ]
  
  # Get unique cities
  cityGrp <- unique(group$City)
  
  # Pre-allocate cityTab
  cityTab <- data.table(
    Period = rep(c("period_2000_2010", "period_2010_2020"), each = length(cityGrp)),
    City = rep(cityGrp, 2),
    mrp_ideology_city = NA_real_,
    partisan_leg_state = NA_real_,
    race_wht_alone_city = NA_real_,
    median_value_city = NA_real_,
    housing_density_city = NA_real_,
    Change_Value = NA_real_,
    Change_Value_Median = NA_real_,
    Change_Value_Mean = NA_real_
  )
  
  # Fill cityTab using vectorized operations
  for (j in 1:nrow(cityTab)) {
    groupCity <- group[
      City == cityTab$City[j] & Period == cityTab$Period[j]
    ]
    if (nrow(groupCity) > 0) {
      cityTab$housing_density_city[j] <- groupCity$housing_density_city[1]
      cityTab$mrp_ideology_city[j] <- groupCity$mrp_ideology_city[1]
      cityTab$partisan_leg_state[j] <- groupCity$partisan_leg_state[1]
      cityTab$race_wht_alone_city[j] <- groupCity$race_wht_alone_city[1]
      cityTab$median_value_city[j] <- groupCity$median_value_city[1]
      cityTab$Change_Value[j] <- (sum(groupCity$units_tract_end, na.rm = TRUE)-sum(groupCity$units_tract_begin, na.rm = TRUE))/sum(groupCity$units_tract_begin, na.rm = TRUE)
      groupCityMedMean <- subset(groupCity, groupCity$units_tract_end>0 & groupCity$units_tract_begin>0) # Eliminates situations where the tract had zero units at the beginning or end
      cityTab$Change_Value_Median[j] <- median((groupCityMedMean$units_tract_end-groupCityMedMean$units_tract_begin)/groupCityMedMean$units_tract_begin, na.rm = TRUE)
      cityTab$Change_Value_Mean[j] <- mean((groupCityMedMean$units_tract_end-groupCityMedMean$units_tract_begin)/groupCityMedMean$units_tract_begin, na.rm = TRUE)
    }
  }
  
  cityTab$Change_Value_Use <- cityTab$Change_Value # Can change to Change_Value_Median or Change_Value_Mean or Change_Value
  
  # Filter and run regression
  regData <- cityTab[ !is.na(Change_Value_Use) & Change_Value_Use < 5 & Change_Value_Use > (-0.9) # Note — excluding major exceptions
  ]
  if (nrow(regData) > 0) {
    regI <- lm(
      Change_Value_Use ~ I(housing_density_city / 1000) + mrp_ideology_city +
        as.factor(partisan_leg_state) + race_wht_alone_city +
        log(as.numeric(median_value_city)) + as.factor(Period),
      data = regData
    )
    # robust <- coeftest(regI, vcov = sandwich)
    robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~City, data = regData))
    binTab$coef[i] <- robust[3, 1]
    binTab$se[i] <- robust[3, 2]
    binTab$p_val[i] <- robust[3, 4]
    binTab$num_cities[i] <- length(unique(regData$City))
  }
}

for(i in 1:nrow(binTab)){
  binTab$p_val_stars[i] <- ifelse(binTab$p_val[i]<0.001,"***",
                                  ifelse(binTab$p_val[i]<0.01,"**",
                                         ifelse(binTab$p_val[i]<0.05,"**","")))
}

binTab$min_max <- paste(binTab$min, binTab$max, sep=" to ")
binTab$min_max <- paste(binTab$min_max, binTab$p_val_stars, sep=" ")
binTab$min_max[1] <- paste("All tracts",binTab$p_val_stars[1],sep=" ")

binTab$min_max <- factor(binTab$min_max, levels=c(binTab$min_max))

# Create the plot
p <- ggplot(binTab, aes(x = min_max, y = coef)) +
  geom_errorbar(aes(ymin = coef - se, ymax = coef + se), width = 0.2, color = "#1696d2") +  # Error bars
  geom_point(size = 2, color = "black") +  # Points for coefficients
  labs(
    #title = "Coefficient Plot with Standard Errors",
    x = "Housing density in units per square mile among tracts in bin, by city",
    y = NULL #"Effect of more liberal ideology on percent change in housing units" # represents the coefficient in the regression
  ) +
  #theme_minimal() +
  coord_cartesian(ylim=c(-0.8,0.5)) + 
  geom_hline(yintercept = 0, linetype = "dotted") +
  scale_y_continuous(breaks = c(-0.75,-0.5,-0.25,0,0.25,0.5)) +
  theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
        axis.text.y = element_text(size = 10),
        axis.title.x = element_text(size = 12))

# Print the plots
ggsave(filename = "figure_7.png", plot = p, width = 8, height = 5)
ggsave(filename = "appendix_figure_8.png", plot = p, width = 8, height = 5)
ggsave(filename = "appendix_figure_9.png", plot = p, width = 8, height = 5)

#### TABLE 1 - TRACTS ####

# The below should have already been activated

# Model 4 -- all, intersected
regI <- (lm(Change_Value ~ I((hsg_density_tract/1000) * mrp_ideology_city) + I(hsg_density_tract/1000) + mrp_ideology_city + as.factor(partisan_leg_state) + race_wht_alone_tract + log(as.numeric(median_value_tract)) + as.factor(Period), data = subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2))))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~City, data = subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2))))
nrow(subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2)))-114

# Model 5 -- all, intersected, with fixed effects
regI <- (lm(Change_Value ~ I((hsg_density_tract/1000) * mrp_ideology_city) + I(hsg_density_tract/1000) + mrp_ideology_city + as.factor(partisan_leg_state) + race_wht_alone_tract + log(as.numeric(median_value_tract)) + as.factor(Period) + as.factor(City), data = subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2))))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~City, data = subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2))))
nrow(subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2)))-114

# Model 6 --  with council control, no fixed effects
regI <- (lm(Change_Value ~ I(((hsg_density_tract/1000)) * dem_share_council_city)+ I((hsg_density_tract/1000)) + dem_share_council_city + as.factor(partisan_leg_state) + race_wht_alone_tract + log(as.numeric(median_value_tract)) + as.factor(Period), data = subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2)))) # altMatTract$hsg_density_tract>0 & 
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~City, data = subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2)))) #
nrow(subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2)))-16198

#### APPENDIX FIGURE 7 - GRAPH SHOWING EFFECTS OF HOUSING GROWTH ON DENSITY ####

altMatUse <- subset(altMatTract, !is.na(altMatTract$Change_Value) & altMatTract$hsg_density_tract>0 & altMatTract$median_value_tract>0 & altMatTract$Change_Value > (-0.5) & altMatTract$Change_Value < (2)) #

altMatUse$hsg_density_tract_1000 <- altMatUse$hsg_density_tract/1000
altMatUse$partisan_leg_state <- as.factor(altMatUse$partisan_leg_state)
altMatUse$Period <- as.factor(altMatUse$Period)

regII <- lm(Change_Value ~ hsg_density_tract_1000 * mrp_ideology_city + 
  partisan_leg_state + race_wht_alone_tract + log(as.numeric(median_value_tract)) + 
  Period, data = altMatUse)

# Define the values you want
ideology_vals <- c(0.63, 0.76, 0.91) # from most conservative to most liberal
density_vals <- seq(min(altMatUse$hsg_density_tract_1000, na.rm = TRUE),
                    max(altMatUse$hsg_density_tract_1000, na.rm = TRUE),
                    length.out = 100)

# Create a data frame for prediction
pred_df <- expand.grid(
  hsg_density_tract_1000 = density_vals,
  mrp_ideology_city = ideology_vals
)
# Add other variables at their mean (or typical) values
pred_df$partisan_leg_state <- levels(altMatUse$partisan_leg_state)[1] # or use mode
pred_df$race_wht_alone_tract <- mean(altMatUse$race_wht_alone_tract, na.rm = TRUE)
pred_df$median_value_tract <- mean(altMatUse$median_value_tract, na.rm = TRUE)
pred_df$Period <- levels(altMatUse$Period)[1] # or use mode

pred_df$predicted <- predict(regII, newdata = pred_df)

pred_df$mrp_ideology_city_factor <- factor(pred_df$mrp_ideology_city, levels =  c(0.63, 0.76, 0.91), labels = c("1st quartile city (more conservative)","Median city","3rd quartile city (more liberal)"))

p <- ggplot(pred_df, aes(x = hsg_density_tract_1000, y = predicted,
                    color = (mrp_ideology_city_factor))) +
  geom_line() +
  labs(#title = "Effect of Ideology on Housing Growth by Density",
       x = "Housing density (units/square mile)",
       y = NULL, # Predicted housing growth
       color = "Ideology value") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1)), labels = scales::percent) +
  scale_x_continuous(breaks = c(0,2.5,5,7.5,10), labels = c(0, "2,500", "5,000", "7,500", "10,000")) +
  coord_cartesian(xlim=c(0,10), ylim=c(-0.1, 0.2)) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        legend.text = element_text(size = 12),
        axis.title.x = element_text(size = 12))

ggsave(filename = "appendix_figure_7.png", plot = p, width = 8, height = 5)

### APPENDIX TABLE 6 ####

# Model 1 (new) -- Intersections, ideology, county-level
regI <- (lm(Change_Value ~ I((hsg_density/1000) * mrp_ideology) + I(hsg_density/1000) + mrp_ideology + race_wht_alone + log(median_value) + as.factor(Period), data = subset(altMat, altMat$Change_Type=="change_overall" & !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2)))) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMat, altMat$Change_Type=="change_overall" & !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2)))) 
nrow(subset(altMat, altMat$Change_Type=="change_overall" & !is.na(altMat$Change_Value) & altMat$Change_Value > (-0.5) & altMat$Change_Value < (2)))

# Model 2 (new) -- Intersections, ideology, city-level
regI <- (lm(Change_Value ~ I((hsg_density/1000) * mrp_ideology) + I(hsg_density/1000) + mrp_ideology + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, altMatCity$Change_Type=="change_overall" & !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2)))) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(altMatCity, altMatCity$Change_Type=="change_overall" & !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2)))) #
nrow(subset(altMatCity, altMatCity$Change_Type=="change_overall" & !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2)))-3

# Model 3 (new) -- Intersections, city council, city-level
regI <- (lm(Change_Value ~ I((hsg_density/1000) * dem_share_council) + I(hsg_density/1000) + dem_share_council + race_wht_alone + log(as.numeric(median_value)) + as.factor(Period), data = subset(altMatCity, altMatCity$Change_Type=="change_overall" & !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2)))) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = ubset(altMatCity, altMatCity$Change_Type=="change_overall" & !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2)))) #
nrow(subset(altMatCity, altMatCity$Change_Type=="change_overall" & !is.na(altMatCity$Change_Value) & altMatCity$Change_Value > (-0.5) & altMatCity$Change_Value < (2)))-488

#### APPENDIX TABLE 7 ####

# Model 1 (new): County, intersected
regI <- (lm(I(chg_00_20/A41AA2000) ~ I((hsg_density_2000/1000) * mrp_ideology_04_11_invert_base_0) + I(hsg_density_2000/1000) + mrp_ideology_04_11_invert_base_0 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20/countyLarge$A41AA2000<2 & countyLarge$chg_00_20/countyLarge$A41AA2000>-0.5)))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(countyLarge, countyLarge$chg_00_20/countyLarge$A41AA2000<2 & countyLarge$chg_00_20/countyLarge$A41AA2000>-0.5)))
nrow(subset(countyLarge, countyLarge$chg_00_20/countyLarge$A41AA2000<2 & countyLarge$chg_00_20/countyLarge$A41AA2000>-0.5))

# Model 2 (new): City, intersected
regI <- (lm(I((units_2020-units_2000)/units_2000) ~ I((hsg_density_2000/1000) * mrp_ideology_04_11_invert_base_0) + I(hsg_density_2000/1000) + mrp_ideology_04_11_invert_base_0 + race_wht_alone + log(median_value), data=subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000<2 & (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000>-0.5)))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000<2 & (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000>-0.5)))
nrow(subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000<2 & (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000>-0.5))-1

# Model 3 (new): City, intersected, council
regI <- (lm(I((units_2020-units_2000)/units_2000) ~ I((hsg_density_2000/1000) * dem_share_council_2000_2009) + I(hsg_density_2000/1000) + dem_share_council_2000_2009 + race_wht_alone + log(median_value), data=subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000<2 & (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000>-0.5)))
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000<2 & (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000>-0.5)))
nrow(subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000<2 & (cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000>-0.5))-301

#### ROBUSTNESS TESTS (NOT PRESENTED IN PAPER) ####

# Using alternative ideology specification
regI <- (lm(chg_00_20_in_tracts_below_1000_pct ~ I(hsg_density_2000/1000) + mrp_ideology_17_21_invert_base_0 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_1000_1999_pct ~ I(hsg_density_2000/1000) + mrp_ideology_17_21_invert_base_0 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_2000_2999_pct ~ I(hsg_density_2000/1000) + mrp_ideology_17_21_invert_base_0 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_3000_4999_pct ~ I(hsg_density_2000/1000) + mrp_ideology_17_21_invert_base_0 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_above_5000_pct ~ I(hsg_density_2000/1000) + mrp_ideology_17_21_invert_base_0 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

## Using presidential specification (2008)
regI <- (lm(chg_00_20_in_tracts_below_1000_pct ~ I(hsg_density_2000/1000) + demshare_pres_08 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_1000_1999_pct ~ I(hsg_density_2000/1000) + demshare_pres_08 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_2000_2999_pct ~ I(hsg_density_2000/1000) + demshare_pres_08 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_3000_4999_pct ~ I(hsg_density_2000/1000) + demshare_pres_08 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_above_5000_pct ~ I(hsg_density_2000/1000) + demshare_pres_08 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

## Using presidential specification (2020)
regI <- (lm(chg_00_20_in_tracts_below_1000_pct ~ I(hsg_density_2000/1000) + demshare_pres_20 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_1000_1999_pct ~ I(hsg_density_2000/1000) + demshare_pres_20 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_2000_2999_pct ~ I(hsg_density_2000/1000) + demshare_pres_20 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_3000_4999_pct ~ I(hsg_density_2000/1000) + demshare_pres_20 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

regI <- (lm(chg_00_20_in_tracts_above_5000_pct ~ I(hsg_density_2000/1000) + demshare_pres_20 + race_wht_alone + log(median_value), data=subset(countyLarge, countyLarge$chg_00_20_in_tracts_below_1000_pct<2 & countyLarge$chg_00_20_in_tracts_below_1000_pct>-0.5)))
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

#### FIGURE 6 (IDEOLOGY VS HOUSING GROWTH BY BIN) ####

# City Level — Change within the respective tracts -- using this in paper
plotChange <- cityLarge %>%
  ggplot(mapping = aes(x = mrp_ideology_04_11_invert_base_0)) +
  geom_smooth(aes(y = chg_00_20_in_tracts_above_3000/hsg_2000_in_tracts_above_3000, color = "Tracts with at least 3,000 units/sq. mi."), show.legend = TRUE) +
  geom_smooth(aes(y = (chg_00_20_in_tracts_above_1000-chg_00_20_in_tracts_above_3000)/(hsg_2000_in_tracts_above_1000-hsg_2000_in_tracts_above_3000), color = "Tracts with 1,000 to 2,999 units/sq. mi."), show.legend = TRUE) +
  geom_smooth(aes(y = (units_2020-units_2000)/(units_2000), color = "All tracts"), show.legend = TRUE) +
  scale_color_manual(values = c(
    "Tracts with at least 3,000 units/sq. mi." = "#ec008b",
    "Tracts with 1,000 to 2,999 units/sq. mi." = "#008bec",
    "All tracts" = "#000"
  )) +
  scale_x_continuous(expand = expansion(mult = c(0.002, 0)),
                     labels = scales::number) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.002)),
                     labels = scales::percent) +
  coord_cartesian(ylim=c(-0.05,0.4), xlim=c(quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.025,na.rm=T),quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.975,na.rm=T))) +
  labs(x = "Ideology (conservative to liberal)",
       y = NULL, #"Percent change in housing, 2000 to 2020, by city",
       color = "Tract Type") +
  theme(axis.text.x = element_text(size = 10),
        legend.text = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10)) +
  scatter_grid()

ggsave(filename = "figure_6.png", plot = plotChange, width = 8, height = 5)

#### Data comparison -- ideologies among counties and cities (footnote 4) ####

# Most liberal cities
nrow(subset(cityLarge, cityLarge$mrp_ideology_17_21_invert_base_0>(0.95))) # 74
# Most conservative cities
nrow(subset(cityLarge, cityLarge$mrp_ideology_17_21_invert_base_0<(0.6))) # 63

# Overall change
mean(subset((cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000, cityLarge$mrp_ideology_17_21_invert_base_0>(0.95)),na.rm=T) # 18.1
mean(subset((cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000, cityLarge$mrp_ideology_17_21_invert_base_0<(0.6)),na.rm=T) # 30.0
t.test(subset((cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000, cityLarge$mrp_ideology_17_21_invert_base_0>(0.95)),
       subset((cityLarge$units_2020-cityLarge$units_2000)/cityLarge$units_2000, cityLarge$mrp_ideology_17_21_invert_base_0<(0.6)),na.rm=T) # **

# tracts 1,000 to 2,999
mean(subset(((cityLarge$hsg_2020_in_tracts_above_1000_in_2000-cityLarge$hsg_2020_in_tracts_above_3000_in_2000)-(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000))/(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000), cityLarge$mrp_ideology_17_21_invert_base_0>(0.95)),na.rm=T) # 16.7%
mean(subset(((cityLarge$hsg_2020_in_tracts_above_1000_in_2000-cityLarge$hsg_2020_in_tracts_above_3000_in_2000)-(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000))/(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000), cityLarge$mrp_ideology_17_21_invert_base_0<(0.6)),na.rm=T) # 6.3%
t.test(subset(((cityLarge$hsg_2020_in_tracts_above_1000_in_2000-cityLarge$hsg_2020_in_tracts_above_3000_in_2000)-(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000))/(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000), cityLarge$mrp_ideology_17_21_invert_base_0>(0.95)),subset(((cityLarge$hsg_2020_in_tracts_above_1000_in_2000-cityLarge$hsg_2020_in_tracts_above_3000_in_2000)-(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000))/(cityLarge$hsg_2000_in_tracts_above_1000-cityLarge$hsg_2000_in_tracts_above_3000), cityLarge$mrp_ideology_17_21_invert_base_0<(0.6)),na.rm=T) # ***

# tracts >= 3,000
mean(subset((cityLarge$hsg_2020_in_tracts_above_3000_in_2000-cityLarge$hsg_2000_in_tracts_above_3000)/cityLarge$hsg_2000_in_tracts_above_3000, cityLarge$mrp_ideology_17_21_invert_base_0>(0.95)),na.rm=T) # 13.2%
mean(subset((cityLarge$hsg_2020_in_tracts_above_3000_in_2000-cityLarge$hsg_2000_in_tracts_above_3000)/cityLarge$hsg_2000_in_tracts_above_3000, cityLarge$mrp_ideology_17_21_invert_base_0<(0.6)),na.rm=T) # -0.00%
t.test(subset((cityLarge$hsg_2020_in_tracts_above_3000_in_2000-cityLarge$hsg_2000_in_tracts_above_3000)/cityLarge$hsg_2000_in_tracts_above_3000, cityLarge$mrp_ideology_17_21_invert_base_0>(0.95)),subset((cityLarge$hsg_2020_in_tracts_above_3000_in_2000-cityLarge$hsg_2000_in_tracts_above_3000)/cityLarge$hsg_2000_in_tracts_above_3000, cityLarge$mrp_ideology_17_21_invert_base_0<(0.6)),na.rm=T) # ***

#### APPENDIX FIGURE 5 - (HOUSING VALUES VS HOUSING GROWTH BY BIN) ####

# City data -- using
plotChange <- cityLarge %>%
  ggplot(mapping = aes(x = median_value)) +
  geom_smooth(aes(y = chg_00_20_in_tracts_above_3000_pct, color = "Tracts with at least 3,000 units/sq. mi."), show.legend = TRUE) +
  geom_smooth(aes(y = (chg_00_20_in_tracts_1000_2999_pct), color = "Tracts with 1,000 to 2,999 units/sq. mi."), show.legend = TRUE) +
  geom_smooth(aes(y = (chg_00_20)/(A41AA2000), color = "All tracts"), show.legend = TRUE) +
  scale_color_manual(values = c(
    "Tracts with at least 3,000 units/sq. mi." = "#ec008b",
    "Tracts with 1,000 to 2,999 units/sq. mi." = "#008bec",
    "All tracts" = "#000"
  )) +
  scale_x_continuous(expand = expansion(mult = c(0.002, 0)),
                     labels = scales::dollar,
                     breaks = c(100000,150000,200000,250000,300000,350000)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.002)),
                     labels = scales::percent) +
  coord_cartesian(ylim=c(-0.1,0.4), xlim=c(quantile(cityLarge$median_value,0.025,na.rm=T),quantile(cityLarge$median_value,0.975,na.rm=T))) +
  labs(x = "Median housing value",
       y = NULL, #"Percent change in housing, 2000 to 2020, by city",
       color = "Tract Type") +
  theme(axis.text.x = element_text(size = 10),
        legend.text = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10)) +
  scatter_grid()

ggsave(filename = "appendix_figure_5.png", plot = plotChange, width = 8, height = 5)

#### APPENDIX FIGURE 6 - (SHARE WHITE VS HOUSING GROWTH BY BIN) ####

# City data -- using
plotChange <- cityLarge %>%
  ggplot(mapping = aes(x = race_wht_alone)) +
  geom_smooth(aes(y = chg_00_20_in_tracts_above_3000_pct, color = "Tracts with at least 3,000 units/sq. mi."), show.legend = TRUE) +
  geom_smooth(aes(y = (chg_00_20_in_tracts_1000_2999_pct), color = "Tracts with 1,000 to 2,999 units/sq. mi."), show.legend = TRUE) +
  geom_smooth(aes(y = (chg_00_20)/(A41AA2000), color = "All tracts"), show.legend = TRUE) +
  scale_color_manual(values = c(
    "Tracts with at least 3,000 units/sq. mi." = "#ec008b",
    "Tracts with 1,000 to 2,999 units/sq. mi." = "#008bec",
    "All tracts" = "#000"
  )) +
  scale_x_continuous(expand = expansion(mult = c(0.002, 0)),
                     labels = scales::percent,
                     breaks = c(0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.002)),
                     labels = scales::percent) +
  coord_cartesian(ylim=c(-0.05,0.3), xlim=c(quantile(cityLarge$race_wht_alone,0.025,na.rm=T),quantile(cityLarge$race_wht_alone,0.975,na.rm=T))) +
  labs(x = "Share population white",
       y = NULL, #"Percent change in housing, 2000 to 2020, by city",
       color = "Tract Type") +
  theme(axis.text.x = element_text(size = 10),
        legend.text = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10)) +
  scatter_grid()

ggsave(filename = "appendix_figure_6.png", plot = plotChange, width = 8, height = 5)

#### APPENDIX FIGURE 4 (DEM PRESIDENTIAL SHARE VS HOUSING CHANGE) ####

# For cities (using)

plotChange <- cityLarge %>%
  ggplot(mapping = aes(x = demshare_pres_20)) +
  geom_smooth(aes(y = chg_00_20_in_tracts_above_3000/hsg_2000_in_tracts_above_3000, color = "Tracts with at least 3,000 units/sq. mi."), show.legend = TRUE, se=FALSE) +
  geom_smooth(aes(y = (chg_00_20_in_tracts_above_1000)/(hsg_2000_in_tracts_above_1000), color = "Tracts with at least 1,000 units/sq. mi."), show.legend = TRUE, se=FALSE) +
  geom_smooth(aes(x = demshare_pres_08, y = (chg_00_20_in_tracts_above_1000)/(hsg_2000_in_tracts_above_1000)), color = "#008bec", lty="dashed", show.legend = FALSE, se=FALSE) +
  geom_smooth(aes(x = demshare_pres_08, y = (chg_00_20_in_tracts_above_3000)/(hsg_2000_in_tracts_above_3000)), color = "#ec008b", lty="dashed", show.legend = FALSE, se=FALSE) +
  geom_smooth(aes(y = (chg_00_20)/(A41AA2000), color = "All tracts"), show.legend = TRUE, se=FALSE) +
  geom_smooth(aes(x = demshare_pres_08, y = (chg_00_20)/(A41AA2000)), color = "#000", lty="dashed", show.legend = FALSE, se=FALSE) +
  scale_color_manual(values = c(
    "Tracts with at least 3,000 units/sq. mi." = "#ec008b",
    "Tracts with at least 1,000 units/sq. mi." = "#008bec",
    "All tracts" = "#000"
  )) +
  scale_x_continuous(expand = expansion(mult = c(0.002, 0)),
                     labels = scales::percent) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.002)),
                     labels = scales::percent) +
  labs(x = "Democratic presidential vote share, 2008 (dashed) and 2020 (solid)",
       # title = "Percent change in housing units, 2000 to 2020, by county",
       y = NULL,
       color = "Tract Type") +
  theme(axis.text.x = element_text(size = 10),
        legend.text = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10)) +
  scatter_grid()

ggsave(filename = "appendix_figure_4.png", plot = plotChange, width = 8, height = 5)

#### MSA -- LAND USE REGULATIONS VERSUS CONSTRUCTION -- NOT USING IN PAPER -- BUT NEEDED TO IMPLEMENT THE BELOW ####

# Linking MSAs to cities
msaLinkTab <- data.frame(matrix(nrow=21,ncol=3))
colnames(msaLinkTab) <- c("msa_name","city_name","PLC_GJ")
msaLinkTab[,2] <- sort(subset(cityLarge$NAME10, cityLarge$NAME10 %in% c("New York","Los Angeles","Chicago","Washington","San Francisco","Philadelphia","Boston","Detroit","Dallas","Houston","Miami","Atlanta","Seattle","Phoenix","Cleveland","Minneapolis","Tampa","St. Louis","Denver","Pittsburgh","San Diego")))
msaLinkTab[,1] <- c("Atlanta, GA MSA","Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA","Chicago--Gary--Kenosha, IL--IN--WI CMSA","Cleveland--Akron, OH CMSA","Dallas--Fort Worth, TX CMSA","Denver--Boulder--Greeley, CO CMSA","Detroit--Ann Arbor--Flint, MI CMSA",
                    "Houston--Galveston--Brazoria, TX CMSA","Los Angeles--Riverside--Orange County, CA CMSA","Miami--Fort Lauderdale, FL CMSA","Minneapolis--St. Paul, MN--WI MSA","New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA","Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA","Phoenix--Mesa, AZ MSA","Pittsburgh, PA MSA","San Diego, CA MSA","San Francisco--Oakland--San Jose, CA CMSA","Seattle--Tacoma--Bremerton, WA CMSA","St. Louis, MO--IL MSA","Tampa--St. Petersburg--Clearwater, FL MSA","Washington--Baltimore, DC--MD--VA--WV CMSA")
for(i in 1:nrow(msaLinkTab)){
  msaLinkTab$PLC_GJ[i] <- subset(cityLarge$GISJOIN, cityLarge$NAME10==msaLinkTab[i,2])
}

nzlud <- read.csv("nzlud_data.csv", header = TRUE)

msaMat <- data.frame(matrix(nrow=length(unique(censusCounties$NAME)),ncol=25))
colnames(msaMat) <- c("MSA","Hsg_2000","Change_00_20","Change_in_tracts_0_1000","Change_in_tracts_1000_2000","Change_in_tracts_2000_3000","Change_in_tracts_3000_5000","Change_in_tracts_5000_plus","Units_tracts_0_1000","Units_tracts_1000_2000","Units_tracts_2000_3000","Units_tracts_3000_5000","Units_tracts_5000_plus","nzlud_zri_full_st","nTracts_0_1000","nTracts_1000_2000","nTracts_2000_3000","nTracts_3000_5000","nTracts_5000_plus",
                      "Change_in_tracts_0_500","Change_in_tracts_500_1000",
                      "Units_tracts_0_500","Units_tracts_500_1000",
                      "nTracts_0_500","nTracts_500_1000")
msaMat[,1] <- unique(censusCounties$NAME)
for(i in 1:nrow(msaMat)){
  group <- subset(censusCounties, censusCounties$NAME==msaMat[i,1])
  msaMat[i,2] <- sum(group$A41AA2000, na.rm=T)
  msaMat[i,3] <- sum(group$chg_00_20, na.rm=T)
  msaMat[i,20] <- sum(group$chg_00_20_in_tracts_below_500, na.rm=T) # 0 to 500
  msaMat[i,21] <- sum(group$chg_00_20_in_tracts_500_999, na.rm=T) # 500 to 1000
  msaMat[i,4] <- msaMat[i,3]-sum(group$chg_00_20_in_tracts_above_1000, na.rm=T) # 0 to 1000
  msaMat[i,5] <- msaMat[i,3]-msaMat[i,4]-sum(group$chg_00_20_in_tracts_above_2000, na.rm=T) # 1000 to 2000
  msaMat[i,6] <- msaMat[i,3]-msaMat[i,4]-msaMat[i,5]-sum(group$chg_00_20_in_tracts_above_3000, na.rm=T) # 2000 to 3000
  msaMat[i,7] <- sum(group$chg_00_20_in_tracts_above_3000, na.rm=T)-sum(group$chg_00_20_in_tracts_above_5000, na.rm=T) # 3000 to 5000
  msaMat[i,8] <- sum(group$chg_00_20_in_tracts_above_5000, na.rm=T) # 5000 plus
  
  msaMat[i,22] <- sum(group$hsg_2000_in_tracts_below_500, na.rm=T) # 0 to 500
  msaMat[i,23] <- sum(group$hsg_2000_in_tracts_above_1000, na.rm=T)-sum(group$hsg_2000_in_tracts_below_500, na.rm=T) # 500 to 1000
  msaMat[i,9] <- sum(group$hsg_2000_in_tracts_below_1000, na.rm=T) # 0 to 1000
  msaMat[i,10] <- sum(group$hsg_2000_in_tracts_above_1000, na.rm=T)-sum(group$hsg_2000_in_tracts_above_2000, na.rm=T) # 1000 to 2000
  msaMat[i,11] <- sum(group$hsg_2000_in_tracts_above_2000, na.rm=T)-sum(group$hsg_2000_in_tracts_above_3000, na.rm=T) # 2000 to 3000
  msaMat[i,12] <- sum(group$hsg_2000_in_tracts_above_3000, na.rm=T)-sum(group$hsg_2000_in_tracts_above_5000, na.rm=T) # 3000 to 5000
  msaMat[i,13] <- sum(group$hsg_2000_in_tracts_above_5000, na.rm=T) # 5000 plus
  groupB <- subset(nzlud, nzlud$name_2000==msaMat[i,1])
  if(nrow(groupB)>0){
    msaMat[i,14] <- groupB$zri_full_st
  }
  groupC <- subset(censusTracts, censusTracts$state_county %in% group$state_county)
  msaMat[i,24] <- nrow(subset(groupC,groupC$density_00 < 500)) # 0 to 500
  msaMat[i,25] <- nrow(subset(groupC,groupC$density_00 < 1000 & groupC$density_00 >= 500)) # 500 to 1000
  msaMat[i,15] <- nrow(subset(groupC,groupC$density_00 < 1000)) # 0 to 1000
  msaMat[i,16] <- nrow(subset(groupC,groupC$density_00 < 2000 & groupC$density_00 >= 1000)) # 1000 to 2000
  msaMat[i,17] <- nrow(subset(groupC,groupC$density_00 < 3000 & groupC$density_00 >= 2000)) # 2000 to 3000
  msaMat[i,18] <- nrow(subset(groupC,groupC$density_00 < 5000 & groupC$density_00 >= 3000)) # 3000 to 5000
  msaMat[i,19] <- nrow(subset(groupC,groupC$density_00 >= 5000)) # >5000
}

msaMat$Change_in_tracts_above_5000_pct <- msaMat$Change_in_tracts_5000_plus/msaMat$Units_tracts_5000_plus
msaMat$Change_in_tracts_3000_5000_pct <- (msaMat$Change_in_tracts_3000_5000)/(msaMat$Units_tracts_3000_5000)
msaMat$Change_in_tracts_1000_2000_pct <- (msaMat$Change_in_tracts_1000_2000)/(msaMat$Units_tracts_1000_2000)
msaMat$Change_in_tracts_2000_3000_pct <- (msaMat$Change_in_tracts_2000_3000)/(msaMat$Units_tracts_2000_3000)
msaMat$Change_in_tracts_1000_3000_pct <- (msaMat$Change_in_tracts_1000_2000+msaMat$Change_in_tracts_2000_3000)/(msaMat$Units_tracts_1000_2000+msaMat$Units_tracts_2000_3000)
msaMat$Change_in_tracts_below_1000_pct <- msaMat$Change_in_tracts_0_1000/msaMat$Units_tracts_0_1000
msaMat$Change_in_tracts_below_500_pct <- msaMat$Change_in_tracts_0_500/msaMat$Units_tracts_0_500
msaMat$Change_in_tracts_500_1000_pct <- (msaMat$Change_in_tracts_500_1000)/(msaMat$Units_tracts_500_1000)

msaMat$Change_in_tracts_below_500_pct_contribution <- (msaMat$Change_in_tracts_0_500/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)
msaMat$Change_in_tracts_500_1000_pct_contribution <- (msaMat$Change_in_tracts_500_1000/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)
msaMat$Change_in_tracts_below_1000_pct_contribution <- (msaMat$Change_in_tracts_0_1000/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)
msaMat$Change_in_tracts_1000_2000_pct_contribution <- (msaMat$Change_in_tracts_1000_2000/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)
msaMat$Change_in_tracts_2000_3000_pct_contribution <- (msaMat$Change_in_tracts_2000_3000/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)
msaMat$Change_in_tracts_3000_5000_pct_contribution <- (msaMat$Change_in_tracts_3000_5000/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)
msaMat$Change_in_tracts_above_5000_pct_contribution <- (msaMat$Change_in_tracts_5000_plus/msaMat$Change_00_20) * (msaMat$Change_00_20/msaMat$Hsg_2000)

summary(lm(Change_00_20/Hsg_2000 ~ nzlud_zri_full_st, data=msaMat))

summary(lm(Change_in_tracts_0_1000/Units_tracts_0_1000 ~ nzlud_zri_full_st, data=msaMat))
# No significant association

summary(lm((Change_in_tracts_1000_2000+Change_in_tracts_2000_3000+Change_in_tracts_3000_5000+Change_in_tracts_5000_plus)/(Units_tracts_1000_2000+Units_tracts_2000_3000+Units_tracts_3000_5000+Units_tracts_5000_plus) ~ nzlud_zri_full_st, data=msaMat))
# This suggests that a higher regulation is associated with more housing production at higher densities

summary(lm((Change_in_tracts_2000_3000+Change_in_tracts_3000_5000+Change_in_tracts_5000_plus)/(Units_tracts_2000_3000+Units_tracts_3000_5000+Units_tracts_5000_plus) ~ nzlud_zri_full_st, data=msaMat))
# Same, stronger

summary(lm((Change_in_tracts_3000_5000+Change_in_tracts_5000_plus)/(Units_tracts_3000_5000+Units_tracts_5000_plus) ~ nzlud_zri_full_st, data=msaMat))
# Not significant

#### FIGURE 2, PANEL A (MSA GROWTH OVERALL) ####

msaMat <- msaMat %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

msaMatUse <- subset(msaMat, msaMat$Hsg_2000>1000000)

# Add in central-city data
for(i in 1:nrow(msaMatUse)){
  cityUse <- subset(msaLinkTab, msaLinkTab$msa_name==msaMatUse$MSA[i])
  group <- subset(censusTractsLargeCity, censusTractsLargeCity$PLC_GJ==cityUse$PLC_GJ)
  msaMatUse$central_city_change_2000_2020[i] <- (sum(group$units_2020,na.rm=T)-sum(group$units_2000,na.rm=T))/sum(group$units_2000,na.rm=T)
}

msaMatUseGraph <- data.frame(matrix(nrow=nrow(msaMatUse)*2,ncol=3))
colnames(msaMatUseGraph) <- c("MSA_new","type","change_2000_2020")
msaMatUseGraph[,1] <- rep(msaMatUse$MSA_new,2)
msaMatUseGraph[,2] <- c(rep("Metropolitan area",length(msaMatUse$MSA_new)),rep("Central city",length(msaMatUse$MSA_new)))
for(i in 1:nrow(msaMatUseGraph)){
  group <- subset(msaMatUse, msaMatUse$MSA_new==msaMatUseGraph[i,1])
  msaMatUseGraph[i,3] <- ifelse(msaMatUseGraph[i,2]=="Metropolitan area", group$Change_00_20/group$Hsg_2000, group$central_city_change_2000_2020)
}

plotOverall <- msaMatUseGraph %>%
  ggplot(mapping = aes(x = MSA_new, y = change_2000_2020, fill = factor(type))) +
  geom_col(position = "dodge") +
  scale_color_manual(values = c("white", "black")) +
  scale_y_continuous(labels = scales::percent, 
                     limits = c(-0.18,0.54),
                     breaks = c(-0.10,0,0.1,0.2,0.3,0.4,0.5),
                     expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL,
       y = NULL,
       title = "A. Percent Increase in Housing Units"
  ) +
  guides(color = "none",fill = guide_legend(ncol = 2)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10),
        legend.position = "none") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

ggsave(filename = "figure_2_panel_a.png", plot = plotOverall, width = 8, height = 3)

#### FIGURE 3 (MSAs BY NUMERICAL CHANGE IN HOUSING UNITS) ####

msaLong <- msaMat %>% # Change_in_tracts_0_1000
  pivot_longer(
    cols = c(Change_in_tracts_0_500,Change_in_tracts_500_1000,Change_in_tracts_1000_2000,Change_in_tracts_2000_3000, Change_in_tracts_3000_5000, Change_in_tracts_5000_plus),
    names_to = "Variable",
    values_to = "Change"
  )

msaLong <- msaLong %>%
  mutate(Variable = recode(Variable,
                           Change_in_tracts_0_500 = "Tracts with less than 500 units/sq. mi.",
                           Change_in_tracts_500_1000 = "Tracts with 500 to 999 units/sq. mi.",
                           Change_in_tracts_1000_2000 = "Tracts with 1,000 to 1,999 units/sq. mi.",
                           Change_in_tracts_2000_3000 = "Tracts with 2,000 to 2,999 units/sq. mi.",
                           Change_in_tracts_3000_5000 = "Tracts with 3,000 to 4,999 units/sq. mi.",
                           Change_in_tracts_5000_plus = "Tracts with at least 5,000 units/sq. mi."
  ))

msaLong$Variable <- factor(msaLong$Variable, levels=c("Tracts with less than 500 units/sq. mi.","Tracts with 500 to 999 units/sq. mi.","Tracts with 1,000 to 1,999 units/sq. mi.","Tracts with 2,000 to 2,999 units/sq. mi.","Tracts with 3,000 to 4,999 units/sq. mi.","Tracts with at least 5,000 units/sq. mi."))

msaLong <- subset(msaLong, msaLong$Hsg_2000>1000000)

msaLong <- msaLong %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

plotMSA <- msaLong %>%
  ggplot() +
  geom_col(mapping = aes(x = MSA_new, y = Change, fill = Variable)) +
  scale_color_manual(values = c("white", "black")) +
  scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL,
       y = NULL) +
  guides(color = "none",fill = guide_legend(ncol = 2)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=10),
        axis.text.y = element_text(size=10),
        legend.text = element_text(size=12)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) +
  scale_fill_manual(values = c("#55b748", "#1696d2", "#fdbf11", "#000000", "#d2d2d2", "#ec008b"))

ggsave(filename = "figure_3.png", plot = plotMSA, width = 8, height = 6)

## RATIOS COMPARED TO CURRENT HOUSING

msaMatLong <- subset(msaMat, msaMat$Hsg_2000>1000000)
msaMatLong <- msaMat

msaMatLong$share_units_overall_2000_in_0_1000 <- msaMatLong$Units_tracts_0_1000/msaMatLong$Hsg_2000
msaMatLong$share_units_growth_2000_in_0_1000 <- msaMatLong$Change_in_tracts_0_1000/msaMatLong$Change_00_20

msaMatLong$ratio <- msaMatLong$share_units_growth_2000_in_0_1000/msaMatLong$share_units_overall_2000_in_0_1000

subset(msaMatLong[,c('MSA','ratio')], msaMatLong$Change_00_20>0 & msaMatLong$ratio<1) # only Binghamton and Chico
subset(msaMatLong[,c('MSA','ratio')], msaMatLong$Change_00_20>0 & msaMatLong$Hsg_2000>1000000) # 

msaMatLong$share_units_overall_2000_in_3000 <- (msaMatLong$Units_tracts_5000_plus+msaMatLong$Units_tracts_3000_5000)/msaMatLong$Hsg_2000
msaMatLong$share_units_growth_2000_in_3000 <- (msaMatLong$Change_in_tracts_5000_plus+msaMatLong$Change_in_tracts_3000_5000)/msaMatLong$Change_00_20

msaMatLong$ratioB <- msaMatLong$share_units_growth_2000_in_3000/msaMatLong$share_units_overall_2000_in_3000

subset(msaMatLong[,c('MSA','ratioB')], !is.na(msaMatLong$ratioB) & msaMatLong$Change_00_20>0 & msaMatLong$Hsg_2000>1000000)
subset(msaMatLong[,c('MSA','ratioB')], !is.na(msaMatLong$ratioB) & msaMatLong$ratioB>1 & msaMatLong$Change_00_20>0)

nrow(msaMatLong) # 273

#### FIGURE 4 (MSA HOUSING UNIT GROWTH, BY DENSITY BIN, PERCENT) ####

msaMat <- msaMat %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

msaLong <- msaMat %>%
  pivot_longer(
    cols = c(Change_in_tracts_above_5000_pct,Change_in_tracts_3000_5000_pct,Change_in_tracts_2000_3000_pct,Change_in_tracts_1000_2000_pct),
    names_to = "Variable",
    values_to = "Change"
  )

msaLong <- msaLong %>%
  mutate(Variable = recode(Variable,
                                "Change_in_tracts_above_5000_pct" = "Tracts with at least 5,000 units/sq. mi.",
                                "Change_in_tracts_3000_5000_pct" = "Tracts with 3,000 to 4,999 units/sq. mi.",
                                "Change_in_tracts_2000_3000_pct" = "Tracts with 2,000 to 2,999 units/sq. mi.",
                                "Change_in_tracts_1000_2000_pct" = "Tracts with 1,000 to 1,999 units/sq. mi."))

msaLong <- subset(msaLong, msaLong$Hsg_2000>1000000)

plotMSAbyTractDensity <- msaLong %>%
  ggplot(mapping = aes(x = MSA_new, y = Change, fill = factor(Variable))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.1))) +
  labs(
    x = NULL,
    y = NULL
  ) +
  guides(color = "none", fill = guide_legend(ncol = 2)) +
  theme(
    legend.text = element_text(size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text.y = element_text(size = 10)
  ) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) +
  scale_fill_manual(
    values = c(
      "Tracts with 1,000 to 1,999 units/sq. mi." = "#fdbf11",
      "Tracts with 2,000 to 2,999 units/sq. mi." = "#000000",
      "Tracts with 3,000 to 4,999 units/sq. mi." = "#d2d2d2",
      "Tracts with at least 5,000 units/sq. mi." = "#ec008b"
    )
  )

ggsave(filename = "figure_4.png", plot = plotMSAbyTractDensity, width = 8, height = 6)

#### APPENDIX FIGURE 2 (MSA HOUSING UNIT GROWTH, BY DENSITY BIN, NUMERICAL) ####

msaMat <- msaMat %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

msaMat$Change_in_tracts_5000_plus_per_tract <- msaMat$Change_in_tracts_5000_plus/msaMat$nTracts_5000_plus
msaMat$Change_in_tracts_3000_5000_per_tract <- msaMat$Change_in_tracts_3000_5000/msaMat$nTracts_3000_5000
msaMat$Change_in_tracts_2000_3000_per_tract <- msaMat$Change_in_tracts_2000_3000/msaMat$nTracts_2000_3000
msaMat$Change_in_tracts_1000_2000_per_tract <- msaMat$Change_in_tracts_1000_2000/msaMat$nTracts_1000_2000

msaLong <- msaMat %>%
  pivot_longer(
    cols = c(Change_in_tracts_5000_plus_per_tract,Change_in_tracts_3000_5000_per_tract,Change_in_tracts_2000_3000_per_tract,Change_in_tracts_1000_2000_per_tract),
    names_to = "Variable",
    values_to = "Change"
  )

msaLong <- msaLong %>%
  mutate(Variable = recode(Variable,
                           "Change_in_tracts_5000_plus_per_tract" = "Tracts with at least 5,000 units/sq. mi.",
                           "Change_in_tracts_3000_5000_per_tract" = "Tracts with 3,000 to 4,999 units/sq. mi.",
                           "Change_in_tracts_2000_3000_per_tract" = "Tracts with 2,000 to 2,999 units/sq. mi.",
                           "Change_in_tracts_1000_2000_per_tract" = "Tracts with 1,000 to 1,999 units/sq. mi."))

msaLong <- subset(msaLong, msaLong$Hsg_2000>1000000)

plotMSAbyTractDensity <- msaLong %>%
  ggplot(mapping = aes(x = MSA_new, y = Change, fill = factor(Variable))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.1))) +
  labs(
    x = NULL,
    y = NULL
  ) +
  guides(color = "none", fill = guide_legend(ncol = 2)) +
  theme(
    legend.text = element_text(size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text.y = element_text(size = 10)
  ) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) +
  scale_fill_manual(
    values = c(
      "Tracts with 1,000 to 1,999 units/sq. mi." = "#fdbf11",
      "Tracts with 2,000 to 2,999 units/sq. mi." = "#000000",
      "Tracts with 3,000 to 4,999 units/sq. mi." = "#d2d2d2",
      "Tracts with at least 5,000 units/sq. mi." = "#ec008b"
    )
  )

ggsave(filename = "appendix_figure_2.png", plot = plotMSAbyTractDensity, width = 8, height = 6)

#### APPENDIX FIGURE 1 (GROWTH IN LOW-DENSITY TRACTS) ####

msaMatBig <- subset(msaMat, msaMat$Hsg_2000>1000000)

plotmsaMatBig <- msaMatBig %>%
  ggplot(mapping = aes(x=MSA_new, y = Change_in_tracts_below_1000_pct)) +
  geom_col() +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.1))) +
  labs(#title = "Change in housing units, 2000 to 2020, by MSA (2000) and by tract density per square mile (2000)", 
       x = NULL,
       y = NULL) +
  guides(color = "none") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=10),
        axis.text.y = element_text(size=10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

ggsave(filename = "appendix_figure_1.png", plot = plotmsaMatBig, width = 8, height = 5)

#### TABLE 2, MODEL III (AVERAGE TRACT DENSITY) ####

# Create panel
densityPanel <- data.frame(matrix(ncol=9,nrow=2*nrow(censusCounties)))
colnames(densityPanel) <- c("GISJOIN","Period","Density","Ideology","StatePartisan","White_Share","Median_Value","Avg_Density","Units_2000")
densityPanel[,1] <- rep(censusCounties$GISJOIN,2)
densityPanel[,2] <- c(rep("Period0010",nrow(censusCounties)),rep("Period1020",nrow(censusCounties)))
for(i in 1:nrow(densityPanel)){
  group <- subset(censusCounties, censusCounties$GISJOIN==densityPanel[i,1])
  ifelse(densityPanel[i,2]=="Period0010",{
    densityPanel$Density[i] <- group$hsg_density_2000
    densityPanel$Ideology[i] <- group$mrp_ideology_04_11_invert_base_0
    densityPanel$StatePartisan[i] <- as.character(group$partisan_leg_2000s)
    densityPanel$White_Share[i] <- group$race_wht_alone
    densityPanel$Median_Value[i] <- group$median_value*1.28 # inflation adjusted to 2010
    densityPanel$Avg_Density[i] <- group$average_tract_density_new_units_00_10
    densityPanel$Units_2000[i] <- group$A41AA2000
  },{
    densityPanel$Density[i] <- group$A41AA2010/group$area_sq_mi
    densityPanel$Ideology[i] <- group$mrp_ideology_17_21_invert_base_0
    densityPanel$StatePartisan[i] <- as.character(group$partisan_leg_2010s)
    densityPanel$White_Share[i] <- group$race_wht_alone_2010
    densityPanel$Median_Value[i] <- group$median_value_2010 # inflation
    densityPanel$Avg_Density[i] <- group$average_tract_density_new_units_10_20
    densityPanel$Units_2000[i] <- group$A41AA2000
  })
}

densityPanelLarge <- subset(densityPanel, densityPanel$Units_2000>50000) # using just the larger counties

# Interactive - using

regI <- (lm(Avg_Density ~ I((Density/1000) * Ideology) + I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(Median_Value) + as.factor(Period), data=densityPanelLarge)) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = densityPanelLarge))
nrow(densityPanelLarge)-20

# For robustness: Without New York City
nycCounties <- c("New York County","Kings County","Queens County","Bronx County","Richmond County")
nycCountiesGJ <- subset(censusCounties$GISJOIN, censusCounties$COUNTY %in% nycCounties & censusCounties$STATE=="New York")

regI <- (lm(Avg_Density ~ I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(Median_Value) + as.factor(Period), data=subset(densityPanelLarge, !(densityPanelLarge$GISJOIN %in% nycCountiesGJ)))) #
summary(regI)
robust <- coeftest(regI, vcov=sandwich) # Robust SEs

##### City version, Model IV #####

cityDatabase$median_value <- as.numeric(cityDatabase$median_value)
cityDatabase$median_value_2010 <- as.numeric(cityDatabase$median_value_2010)

# Create panel
# Reshape the data in a vectorized way
densityPanelCities <- cityDatabase %>%
  # Create a row for each period
  bind_rows(
    # Period 0010
    cityDatabase %>%
      mutate(
        Period = "Period0010",
        Median_Value = median_value * 1.28  # inflation adjusted to 2010
      ) %>%
      rename(
        Density = hsg_density_2000,
        Ideology = mrp_ideology_04_11_invert_base_0,
        StatePartisan = partisan_leg_2000s,
        White_Share = race_wht_alone,
        Avg_Density = average_tract_density_new_units_00_10,
        Units_2000 = units_2000,
        Dem_City_Council = dem_share_council_2000_2009
      ),
    # Period 1020
    cityDatabase %>%
      mutate(Period = "Period1020") %>%
      rename(
        Density = hsg_density_2010,
        Ideology = mrp_ideology_17_21_invert_base_0,
        StatePartisan = partisan_leg_2010s,
        White_Share = race_wht_alone_2010,
        Median_Value = median_value_2010,  # no adjustment needed
        Avg_Density = average_tract_density_new_units_10_20,
        Units_2000 = units_2010,
        Dem_City_Council = dem_share_council_2010_2019
      )
  ) %>%
  select(GISJOIN, Period, Density, Ideology, StatePartisan, White_Share, Median_Value, Avg_Density, Units_2000, Dem_City_Council)

densityPanelLargeCities <- subset(densityPanelCities, densityPanelCities$Units_2000>25000)

# Interaction — using
regI <- (lm(Avg_Density ~ I((Density/1000) * Ideology) + I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(Median_Value) + as.factor(Period), data=densityPanelLargeCities)) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = densityPanelLargeCities))
nrow(densityPanelLargeCities)-12

#### FIGURE 2, PANEL C (AVERAGE TRACT DENSITY OF NEW UNITS -- MSAS) ####

msaMatDensity <- data.frame(matrix(nrow=length(unique(censusCounties$NAME)),ncol=3))
colnames(msaMatDensity) <- c("MSA","Average_tract_density_of_new_units","Housing_2000")
msaMatDensity[,1] <- unique(censusCounties$NAME)
for(i in 1:nrow(msaMatDensity)){
  group <- subset(censusCounties, censusCounties$NAME==msaMatDensity[i,1])
  msaMatDensity[i,2] <- sum(group$average_tract_density_new_units*group$A41AA2000,na.rm=T)/sum(group$A41AA2000,na.rm=T)
  msaMatDensity[i,3] <- sum(group$A41AA2000,na.rm=T)
}

msaMatBig <- subset(msaMatDensity, msaMatDensity$Housing_2000>1000000)

# Add in central-city data
for(i in 1:nrow(msaMatBig)){
  cityUse <- subset(msaLinkTab, msaLinkTab$msa_name==msaMatBig$MSA[i])
  msaMatBig$Average_tract_density_of_new_units_city[i] <- subset(cityDatabase$average_tract_density_new_units, cityDatabase$PLC_GJ==cityUse$PLC_GJ)
}

msaMatBig <- msaMatBig %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

# Just MSAs
plotAvgDens <- msaMatBig %>%
  ggplot(mapping = aes(x=MSA_new, y = Average_tract_density_of_new_units)) +
  geom_col() +
  scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.1))) +
  labs(, x = NULL,
       y = NULL,
       title = "C. Average Tract Density Where a Unit Was Added") +
  guides(color = "none") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

ggsave(filename = "figure_2_panel_c_old.png", plot = plotAvgDens, width = 8, height = 2.5)

# To add city data

msaMatUseGraph <- data.frame(matrix(nrow=nrow(msaMatBig)*2,ncol=3))
colnames(msaMatUseGraph) <- c("MSA_new","type","average_tract_density_new_units")
msaMatUseGraph[,1] <- rep(msaMatBig$MSA_new,2)
msaMatUseGraph[,2] <- c(rep("Metropolitan area",length(msaMatBig$MSA_new)),rep("Central city",length(msaMatBig$MSA_new)))
for(i in 1:nrow(msaMatUseGraph)){
  group <- subset(msaMatBig, msaMatBig$MSA_new==msaMatUseGraph[i,1])
  msaMatUseGraph[i,3] <- ifelse(msaMatUseGraph[i,2]=="Metropolitan area", group$Average_tract_density_of_new_units, group$Average_tract_density_of_new_units_city)
}

# With cities

plotAvgDens <- msaMatUseGraph %>%
  ggplot(mapping = aes(x = MSA_new, y = average_tract_density_new_units, fill = factor(type))) +
  geom_col(position = "dodge") +
  scale_color_manual(values = c("white", "black")) +
  scale_y_continuous(labels = scales::comma, 
                     limits = c(0,20000),
                     #breaks = c(0,0.15,0.3,0.45),
                     expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL,
       y = NULL,
       title = "C. Average Tract Density Where a Unit Was Added"
  ) +
  guides(color = "none",fill = guide_legend(ncol = 2)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10),
        legend.position = "none") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

ggsave(filename = "figure_2_panel_c.png", plot = plotAvgDens, width = 8, height = 3)

#### FOR AVERAGE PERSON, HOW MUCH HOUSING IS ADDED (NOT USING IN PAPER, BUT NEEDED FOR THE FOLLOWING) ####

MSAs <- unique(censusCounties$NAME)
MSAs <- MSAs[!is.na(MSAs)]

typicalMat <- data.frame(matrix(ncol=8, nrow=length(MSAs)))
colnames(typicalMat) <- c("MSA","typical","unit_change_2000_20","units_2000","typical_in_less_1000","typical_in_1_3000","typical_in_3_5000","typical_in_5000_plus")
typicalMat[,1] <- MSAs

for(i in 1:nrow(typicalMat)){
  # Subset counties for the current MSA/group
  groupCounties <- subset(censusCounties$GISJOIN, censusCounties$NAME == typicalMat[i, 1])
  group <- subset(censusCounties, censusCounties$GISJOIN %in% groupCounties)
  group$state_county <- paste(group$STATEFP, group$COUNTYFP, sep = "-")
  
  # Store total new and existing units for the group
  typicalMat[i, 3] <- sum(group$chg_00_20, na.rm = TRUE)
  typicalMat[i, 4] <- sum(group$A41AA2000, na.rm = TRUE)
  
  # Subset tracts for all counties in the group
  groupB <- subset(censusTracts, censusTracts$state_county %in% group$state_county)
  groupB$newUnits_no_zero <- pmax(groupB$CM7AA2020 - groupB$CM7AA2000, 0)
  
  typicalMat[i, 2] <- sum(groupB$newUnits_no_zero * groupB$CM7AA2000, na.rm=T)/sum(groupB$CM7AA2000, na.rm=T)
  
  groupC <- subset(censusTracts, censusTracts$density_00 < 1000 & censusTracts$state_county %in% group$state_county)
  groupC$newUnits_no_zero <- pmax(groupC$CM7AA2020 - groupC$CM7AA2000, 0)
  typicalMat[i, 5] <- sum(groupC$newUnits_no_zero * groupC$CM7AA2000, na.rm=T)/sum(groupC$CM7AA2000, na.rm=T)
  groupD <- subset(censusTracts, censusTracts$density_00 >= 1000 & censusTracts$density_00 < 3000 & censusTracts$state_county %in% group$state_county)
  groupD$newUnits_no_zero <- pmax(groupD$CM7AA2020 - groupD$CM7AA2000, 0)
  typicalMat[i, 6] <- sum(groupD$newUnits_no_zero * groupD$CM7AA2000, na.rm=T)/sum(groupD$CM7AA2000, na.rm=T)
  groupE <- subset(censusTracts, censusTracts$density_00 >= 3000 & censusTracts$density_00 < 5000 & censusTracts$state_county %in% group$state_county)
  groupE$newUnits_no_zero <- pmax(groupE$CM7AA2020 - groupE$CM7AA2000, 0)
  typicalMat[i, 7] <- sum(groupE$newUnits_no_zero * groupE$CM7AA2000, na.rm=T)/sum(groupE$CM7AA2000, na.rm=T)
  groupF <- subset(censusTracts, censusTracts$density_00 >= 5000 & censusTracts$state_county %in% group$state_county)
  groupF$newUnits_no_zero <- pmax(groupF$CM7AA2020 - groupF$CM7AA2000, 0)
  typicalMat[i, 8] <- sum(groupF$newUnits_no_zero * groupF$CM7AA2000, na.rm=T)/sum(groupF$CM7AA2000, na.rm=T)
  
}

typicalMatUse <- subset(typicalMat, typicalMat$units_2000>1000000)

typicalLong <- typicalMatUse %>%
  pivot_longer(
    cols = c(typical_in_less_1000,typical_in_1_3000,typical_in_3_5000,typical_in_5000_plus),
    names_to = "Variable",
    values_to = "Change"
  )

typicalLong <- typicalLong %>%
  mutate(Variable = recode(Variable,
                           typical_in_less_1000 = "Tracts with 0 to 1000 units/sq. mi.",
                           typical_in_1_3000 = "Tracts with 1000 to 3000 units/sq. mi.",
                           typical_in_3_5000 = "Tracts with 3000 to 4000 units/sq. mi.",
                           typical_in_5000_plus = "Tracts with greater than 5000 units/sq. mi."
  ))

typicalLong <- typicalLong %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

plotTypical <- typicalLong %>%
  ggplot(mapping = aes(x=MSA_new, y = Change, fill = factor(Variable))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL,
       y = "Increase in housing units, weighted by residents") +
  guides(color = "none",fill = guide_legend(ncol = 2)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

#### FIGURE 2, PANEL D (DISSIMILARITY BY MSA) ####

# Dissimilarity index.

MSAs <- unique(censusCounties$NAME)
MSAs <- MSAs[!is.na(MSAs)]

disMat <- data.frame(matrix(ncol=4, nrow=length(MSAs)))
colnames(disMat) <- c("MSA","dissimilarity","unit_change_2000_20","units_2000")
disMat[,1] <- MSAs

for(i in 1:nrow(disMat)){
  # Subset counties for the current MSA/group
  #groupCounties <- subset(censusCountiesMSAs2000$GISJOIN, censusCountiesMSAs2000$NAME == disMat[i, 1])
  groupCounties <- subset(censusCounties$GISJOIN, censusCounties$NAME == typicalMat[i, 1])
  group <- subset(censusCounties, censusCounties$GISJOIN %in% groupCounties)
  group$state_county <- paste(group$STATEFP, group$COUNTYFP, sep = "-")
  
  # Store total new and existing units for the group
  disMat[i, 3] <- sum(group$chg_00_20, na.rm = TRUE)
  disMat[i, 4] <- sum(group$A41AA2000, na.rm = TRUE)
  
  # Subset tracts for all counties in the group
  groupB <- subset(censusTracts, censusTracts$state_county %in% group$state_county)
  groupB$newUnits <- pmax(groupB$CM7AA2020 - groupB$CM7AA2000, 0)
  
  dissim <- 0
  
  if(nrow(groupB) > 0){
    # Calculate total new and existing units ONCE
    total_new <- sum(groupB$newUnits, na.rm = TRUE)
    total_existing <- sum(groupB$CM7AA2000, na.rm = TRUE)
    
    # Avoid division by zero
    if(total_new > 0 && total_existing > 0){
      # Calculate shares for all tracts at once
      groupB$share_new <- groupB$newUnits / total_new
      groupB$share_existing <- groupB$CM7AA2000 / total_existing
      
      # Sum absolute differences between shares
      dissim <- sum(abs(groupB$share_new - groupB$share_existing), na.rm = TRUE)
    }
  }
  
  # Calculate the Dissimilarity Index
  disMat[i, 2] <- 0.5 * dissim
}

# If dissimilarity = 0 : New units are added in perfect proportion to existing units
# If dissimilarity = 1 : New units are added only in tracts with no existing units (or vice versa).

disMatUse <- subset(disMat, disMat$units_2000>1000000)

disMatUse <- disMatUse %>%
  mutate(
    MSA_new = case_when(
      MSA == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      MSA == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      MSA == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      MSA == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      MSA == "San Diego, CA MSA" ~ "San Diego",
      MSA == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      MSA == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      MSA == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      MSA == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      MSA == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      MSA == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      MSA == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      MSA == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      MSA == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      MSA == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      MSA == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      MSA == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      MSA == "St. Louis, MO--IL MSA" ~ "St. Louis",
      MSA == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      MSA == "Atlanta, GA MSA" ~ "Atlanta",
      MSA == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ "" # Default value if no match
    )
  )

# Add in central-city data
for(i in 1:nrow(disMatUse)){
  cityUse <- subset(msaLinkTab, msaLinkTab$msa_name==disMatUse$MSA[i])
  disMatUse$dissimilarity_cities[i] <- subset(disMatCities$dissimilarity, disMatCities$GISJOIN==cityUse$PLC_GJ)
}

disMatUseGraph <- data.frame(matrix(nrow=nrow(disMatUse)*2,ncol=3))
colnames(disMatUseGraph) <- c("MSA_new","type","dissimilarity")
disMatUseGraph[,1] <- rep(disMatUse$MSA_new,2)
disMatUseGraph[,2] <- c(rep("Metropolitan area",length(disMatUse$MSA_new)),rep("Central city",length(disMatUse$MSA_new)))
for(i in 1:nrow(disMatUseGraph)){
  group <- subset(disMatUse, disMatUse$MSA_new==disMatUseGraph[i,1])
  disMatUseGraph[i,3] <- ifelse(disMatUseGraph[i,2]=="Metropolitan area", group$dissimilarity, group$dissimilarity_cities)
}

plotDissim <- disMatUseGraph %>%
  ggplot(mapping = aes(x = MSA_new, y = dissimilarity, fill = factor(type))) +
  geom_col(position = "dodge") +
  scale_color_manual(values = c("white", "black")) +
  scale_y_continuous(labels = scales::comma, 
                     limits = c(0,0.9),
                     #breaks = c(0,0.15,0.3,0.45),
                     expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL,
       y = NULL,
       title = "D. Dissimilarity"
  ) +
  guides(color = "none",fill = guide_legend(ncol = 2)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10),
        legend.position = "none") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

ggsave(filename = "figure_2_panel_d.png", plot = plotDissim, width = 8, height = 2.5)

#### TABLE 2, MODEL I (COUNTY WEIGHTED POINTS NEARBY) ####

pointsNearPanel <- data.frame(matrix(ncol=9,nrow=2*nrow(countyLarge)))
colnames(pointsNearPanel) <- c("GISJOIN","Period","Density","Ideology","StatePartisan","White_Share","Median_Value","NearNumber","Units_2000")
pointsNearPanel[,1] <- rep(countyLarge$GISJOIN,2)
pointsNearPanel[,2] <- c(rep("Period0010",nrow(countyLarge)),rep("Period1020",nrow(countyLarge)))
for(i in 1:nrow(pointsNearPanel)){
  group <- subset(countyLarge, countyLarge$GISJOIN==pointsNearPanel[i,1])
  ifelse(pointsNearPanel[i,2]=="Period0010",{
    pointsNearPanel$Density[i] <- group$hsg_density_2000
    pointsNearPanel$Ideology[i] <- group$mrp_ideology_04_11_invert_base_0
    pointsNearPanel$StatePartisan[i] <- as.character(group$partisan_leg_2000s)
    pointsNearPanel$White_Share[i] <- group$race_wht_alone
    pointsNearPanel$Median_Value[i] <- group$median_value*1.28 # inflation adjusted to 2010
    pointsNearPanel$NearNumber[i] <- group$Average_NUMPOINTS_weighted_00_10
    pointsNearPanel$Units_2000[i] <- group$A41AA2000
  },{
    pointsNearPanel$Density[i] <- group$A41AA2010/group$area_sq_mi
    pointsNearPanel$Ideology[i] <- group$mrp_ideology_17_21_invert_base_0
    pointsNearPanel$StatePartisan[i] <- as.character(group$partisan_leg_2010s)
    pointsNearPanel$White_Share[i] <- group$race_wht_alone_2010
    pointsNearPanel$Median_Value[i] <- group$median_value_2010 # inflation
    pointsNearPanel$NearNumber[i] <- group$Average_NUMPOINTS_weighted_10_20
    pointsNearPanel$Units_2000[i] <- group$A41AA2000
  })
}

# Interactive -- using
regI <- (lm(NearNumber ~ I((Density/1000) * Ideology) + I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(Median_Value) + as.factor(Period), data=pointsNearPanel)) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = pointsNearPanel))
nrow(pointsNearPanel)

##### CITY ANALYSIS, Model II #####

pointsNearPanelCities <- data.frame(matrix(ncol=10,nrow=2*nrow(cityLarge)))
colnames(pointsNearPanelCities) <- c("GISJOIN","Period","Density","Ideology","StatePartisan","White_Share","Median_Value","NearNumber","Units_2000","DemCouncil")
pointsNearPanelCities[,1] <- rep(cityLarge$GISJOIN,2)
pointsNearPanelCities[,2] <- c(rep("Period0010",nrow(cityLarge)),rep("Period1020",nrow(cityLarge)))
for(i in 1:nrow(pointsNearPanelCities)){
  group <- subset(cityLarge, cityLarge$GISJOIN==pointsNearPanelCities[i,1])
  ifelse(pointsNearPanelCities[i,2]=="Period0010",{
    pointsNearPanelCities$Density[i] <- group$hsg_density_2000
    pointsNearPanelCities$Ideology[i] <- group$mrp_ideology_04_11_invert_base_0
    pointsNearPanelCities$StatePartisan[i] <- as.character(group$partisan_leg_2000s)
    pointsNearPanelCities$White_Share[i] <- group$race_wht_alone
    pointsNearPanelCities$Median_Value[i] <- group$median_value*1.28 # inflation adjusted to 2010
    pointsNearPanelCities$NearNumber[i] <- group$Average_NUMPOINTS_weighted_00_10
    pointsNearPanelCities$Units_2000[i] <- group$units_2000
    pointsNearPanelCities$DemCouncil[i] <- group$dem_share_council_2000_2009
  },{
    pointsNearPanelCities$Density[i] <- group$hsg_density_2010
    pointsNearPanelCities$Ideology[i] <- group$mrp_ideology_17_21_invert_base_0
    pointsNearPanelCities$StatePartisan[i] <- as.character(group$partisan_leg_2010s)
    pointsNearPanelCities$White_Share[i] <- group$race_wht_alone_2010
    pointsNearPanelCities$Median_Value[i] <- group$median_value_2010 # inflation
    pointsNearPanelCities$NearNumber[i] <- group$Average_NUMPOINTS_weighted_10_20
    pointsNearPanelCities$Units_2000[i] <- group$units_2010
    pointsNearPanelCities$DemCouncil[i] <- group$dem_share_council_2010_2019
  })
}

# Interactions -- using
regI <- (lm(NearNumber ~ I((Density/1000) * Ideology) + I(Density/1000) + Ideology + as.factor(StatePartisan) + White_Share + log(as.numeric(Median_Value)) + as.factor(Period), data=pointsNearPanelCities)) #
summary(regI)
robust <- coeftest(regI, vcov=vcovCL(regI, cluster = ~GISJOIN, data = pointsNearPanelCities))
nrow(pointsNearPanelCities) - 5

#### FIGURE 2, PANEL B (POINTS NEARBY BY METRO AREA) ####

censusTracts <- censusTracts %>%
  left_join(
    select(censusCounties, state_county, MSA_GJ, NAME),
    by = "state_county"
  )

# Calculate weighted average and median for each metro area
weighted_stats <- censusTracts %>%
  group_by(NAME) %>%
  summarise(
    Average_NUMPOINTS_weighted_00_10 = weighted.mean(Average_NUMPOINTS_00_10, w = units_2000, na.rm = TRUE),
    Average_NUMPOINTS_weighted_10_20 = weighted.mean(Average_NUMPOINTS_10_20, w = units_2010, na.rm = TRUE),
    Average_NUMPOINTS_weighted_00_20 = weighted.mean(Average_NUMPOINTS_00_20, w = units_2000, na.rm = TRUE)
  )

metroTab <- data.frame(matrix(ncol=6,nrow=length(unique(censusCounties$MSA_GJ))))
colnames(metroTab) <- c("MSA_GJ","NAME","FL5001","Average_NUMPOINTS_weighted_00_10","Average_NUMPOINTS_weighted_10_20","Average_NUMPOINTS_weighted_00_20")
metroTab[,1] <- unique(censusCounties$MSA_GJ)

metroTab <- censusCounties %>%
  # Keep only unique MSA_GJ rows
  distinct(MSA_GJ, .keep_all = TRUE) %>%
  # Select and rename columns
  select(MSA_GJ, NAME, FL5001) %>%
  # Join with weighted_stats to add the weighted averages and medians
  left_join(
    weighted_stats %>%
      select(NAME, Average_NUMPOINTS_weighted_00_10, Average_NUMPOINTS_weighted_10_20, Average_NUMPOINTS_weighted_00_20),
    by = "NAME"
  )

metroTab <- metroTab %>%
  mutate(
    MSA_new = case_when(
      NAME == "Los Angeles--Riverside--Orange County, CA CMSA" ~ "Los Angeles",
      NAME == "Chicago--Gary--Kenosha, IL--IN--WI CMSA" ~ "Chicago",
      NAME == "Houston--Galveston--Brazoria, TX CMSA" ~ "Houston",
      NAME == "Phoenix--Mesa, AZ MSA" ~ "Phoenix",
      NAME == "San Diego, CA MSA" ~ "San Diego",
      NAME == "New York--Northern New Jersey--Long Island, NY--NJ--CT--PA CMSA" ~ "New York",
      NAME == "Dallas--Fort Worth, TX CMSA" ~ "Dallas",
      NAME == "Miami--Fort Lauderdale, FL CMSA" ~ "Miami",
      NAME == "Detroit--Ann Arbor--Flint, MI CMSA" ~ "Detroit",
      NAME == "Seattle--Tacoma--Bremerton, WA CMSA" ~ "Seattle",
      NAME == "Philadelphia--Wilmington--Atlantic City, PA--NJ--DE--MD CMSA" ~ "Philadelphia",
      NAME == "Cleveland--Akron, OH CMSA" ~ "Cleveland",
      NAME == "Pittsburgh, PA MSA" ~ "Pittsburgh",
      NAME == "San Francisco--Oakland--San Jose, CA CMSA" ~ "San Francisco",
      NAME == "Boston--Worcester--Lawrence, MA--NH--ME--CT CMSA" ~ "Boston",
      NAME == "Tampa--St. Petersburg--Clearwater, FL MSA" ~ "Tampa",
      NAME == "Minneapolis--St. Paul, MN--WI MSA" ~ "Minneapolis",
      NAME == "St. Louis, MO--IL MSA" ~ "St. Louis",
      NAME == "Washington--Baltimore, DC--MD--VA--WV CMSA" ~ "Washington",
      NAME == "Atlanta, GA MSA" ~ "Atlanta",
      NAME == "Denver--Boulder--Greeley, CO CMSA" ~ "Denver",
      TRUE ~ NA # Default value if no match
    )
  )

metroTabUse <- subset(metroTab, !is.na(metroTab$MSA_new))

# Add in central-city data
for(i in 1:nrow(metroTabUse)){
  cityUse <- subset(msaLinkTab, msaLinkTab$msa_name==metroTabUse$NAME[i])
  metroTabUse$Average_NUMPOINTS_weighted_00_20_city[i] <- subset(cityLarge$Average_NUMPOINTS_weighted_00_20, cityLarge$GISJOIN==cityUse$PLC_GJ)
}

msaMatUseGraph <- data.frame(matrix(nrow=nrow(metroTabUse)*2,ncol=3))
colnames(msaMatUseGraph) <- c("MSA_new","type","Average_NUMPOINTS_weighted_00_20_city")
msaMatUseGraph[,1] <- rep(metroTabUse$MSA_new,2)
msaMatUseGraph[,2] <- c(rep("Metropolitan area",length(metroTabUse$MSA_new)),rep("Central city",length(metroTabUse$MSA_new)))
for(i in 1:nrow(msaMatUseGraph)){
  group <- subset(metroTabUse, metroTabUse$MSA_new==msaMatUseGraph[i,1])
  msaMatUseGraph[i,3] <- ifelse(msaMatUseGraph[i,2]=="Metropolitan area", group$Average_NUMPOINTS_weighted_00_20, group$Average_NUMPOINTS_weighted_00_20_city)
}

plotOverall <- msaMatUseGraph %>%
  ggplot(mapping = aes(x = MSA_new, y = Average_NUMPOINTS_weighted_00_20_city, fill = factor(type))) +
  geom_col(position = "dodge") +
  scale_color_manual(values = c("white", "black")) +
  scale_y_continuous(labels = scales::comma, 
                     limits = c(0,2000),
                     #breaks = c(0,0.15,0.3,0.45),
                     expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL,
       y = NULL,
       title = "B. Weighted Average Number of Housing Units Added Within Half-Mile"
  ) +
  guides(color = "none",fill = guide_legend(ncol = 2)) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
        axis.text.y = element_text(size = 10),
        legend.position = "none") +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20))

ggsave(filename = "figure_2_panel_b.png", plot = plotOverall, width = 8, height = 3)

#### FIGURE 5, PANEL A (MSA BY US REGION) ####

msaList <- unique(censusCounties$NAME)
msaList <- subset(msaList, !is.na(msaList))

# Bureau of Economic Analysis Regions
regionFarWest <- c("Alaska","Hawaii","California","Nevada","Oregon","Washington")
regionRockyMountain <- c("Colorado","Utah","Wyoming","Idaho","Montana")
regionSouthwest <- c("Arizona","New Mexico","Oklahoma","Texas")
regionPlains <- c("North Dakota","South Dakota","Nebraska","Kansas","Minnesota","Iowa","Missouri")
regionGreatLakes <- c("Illinois","Wisconsin","Indiana","Michigan","Ohio")
regionNewEngland <- c("Maine","New Hampshire","Vermont","Rhode Island","Massachusetts","Connecticut")
regionMideast <- c("New York","New Jersey","Pennsylvania","Delaware","District of Columbia","Maryland")
regionSoutheast <- c("Arkansas","Louisiana","Mississippi","Alabama","Tennessee","Kentucky","West Virginia","Virginia","North Carolina","South Carolina","Georgia","Florida")

msaTabLarge <- data.frame(matrix(ncol=19, nrow=length(msaList)))
colnames(msaTabLarge) <- c("MSA","Units_2000","Units_2020","BEA_Region","Housing_Value",
                           "Base_less_1000","Base_1000_2000","Base_2000_3000","Base_3000_5000","Base_5000_plus","Added_less_1000","Added_1000_2000","Added_2000_3000","Added_3000_5000","Added_5000_plus",
                           "Base_less_500","Base_500_1000", # 16, 17
                           "Added_less_500","Added_500_1000")
msaTabLarge[,1] <- msaList
for(i in 1:nrow(msaTabLarge)){
  group <- subset(censusCounties, censusCounties$NAME==msaTabLarge[i,1])
  msaTabLarge[i,2] <- sum(group$A41AA2000,na.rm=T)
  msaTabLarge[i,3] <- sum(group$A41AA2020,na.rm=T)
  msaTabLarge[i,5] <- sum(group$A41AA2000*group$median_value,na.rm=T)/sum(group$A41AA2000,na.rm=T)
  
  msaTabLarge[i,16] <- sum(group$hsg_2000_in_tracts_below_500,na.rm=T)
  msaTabLarge[i,17] <- sum(group$hsg_2000_in_tracts_above_500,na.rm=T)-sum(group$hsg_2000_in_tracts_above_1000,na.rm=T)
  
  msaTabLarge[i,6] <- sum(group$hsg_2000_in_tracts_below_1000,na.rm=T)
  msaTabLarge[i,7] <- sum(group$hsg_2000_in_tracts_above_1000,na.rm=T)-sum(group$hsg_2000_in_tracts_above_2000,na.rm=T)
  msaTabLarge[i,8] <- sum(group$hsg_2000_in_tracts_above_2000,na.rm=T)-sum(group$hsg_2000_in_tracts_above_3000,na.rm=T)
  msaTabLarge[i,9] <- sum(group$hsg_2000_in_tracts_above_3000,na.rm=T)-sum(group$hsg_2000_in_tracts_above_5000,na.rm=T)
  msaTabLarge[i,10] <- sum(group$hsg_2000_in_tracts_above_5000,na.rm=T)
  
  msaTabLarge[i,18] <- sum(group$hsg_2020_in_tracts_below_500_in_2000,na.rm=T)-msaTabLarge[i,16]
  msaTabLarge[i,19] <- sum(group$hsg_2020_in_tracts_above_500_in_2000,na.rm=T)-sum(group$hsg_2020_in_tracts_above_1000_in_2000,na.rm=T)-msaTabLarge[i,17]
  
  msaTabLarge[i,11] <- sum(group$hsg_2020_in_tracts_below_1000_in_2000,na.rm=T)-msaTabLarge[i,6]
  msaTabLarge[i,12] <- sum(group$hsg_2020_in_tracts_above_1000_in_2000,na.rm=T)-sum(group$hsg_2020_in_tracts_above_2000_in_2000,na.rm=T)-msaTabLarge[i,7]
  msaTabLarge[i,13] <- sum(group$hsg_2020_in_tracts_above_2000_in_2000,na.rm=T)-sum(group$hsg_2020_in_tracts_above_3000_in_2000,na.rm=T)-msaTabLarge[i,8]
  msaTabLarge[i,14] <- sum(group$hsg_2020_in_tracts_above_3000_in_2000,na.rm=T)-sum(group$hsg_2020_in_tracts_above_5000_in_2000,na.rm=T)-msaTabLarge[i,9]
  msaTabLarge[i,15] <- sum(group$hsg_2020_in_tracts_above_5000_in_2000,na.rm=T)-msaTabLarge[i,10]
  
  biggestCountyState <- group[order(-group$A41AA2000),]$STATE[1]
  msaTabLarge[i,4] <- ifelse(biggestCountyState %in% regionFarWest,"Far West",
                      ifelse(biggestCountyState %in% regionRockyMountain,"Rocky Mountains",
                      ifelse(biggestCountyState %in% regionSouthwest,"Southwest",
                      ifelse(biggestCountyState %in% regionPlains,"Plains",
                      ifelse(biggestCountyState %in% regionGreatLakes,"Great Lakes",
                      ifelse(biggestCountyState %in% regionNewEngland,"New England",
                             ifelse(biggestCountyState %in% regionMideast,"Mideast",
                                    ifelse(biggestCountyState %in% regionSoutheast,"Southeast",""))))))))
}

msaTabSummary <- data.frame(matrix(ncol=33,nrow=length(unique(msaTabLarge$BEA_Region))))
colnames(msaTabSummary) <- c("Region","Units_2000","Base_less_1000","Base_1000_2000","Base_2000_3000","Base_3000_5000","Base_5000_plus","Added_all","Added_less_1000","Added_1000_2000","Added_2000_3000","Added_3000_5000","Added_5000_plus","Pct_change_all","Pct_change_less_1000","Pct_change_1000_2000","Pct_change_2000_3000","Pct_change_3000_5000","Pct_change_5000_plus","Pct_change_all_med","Pct_change_less_1000_med","Pct_change_1000_2000_med","Pct_change_2000_3000_med","Pct_change_3000_5000_med","Pct_change_5000_plus_med",
                             "Base_less_500","Base_500_1000", # 26, 27
                             "Added_less_500","Added_500_1000", # 28, 29
                             "Pct_change_less_500","Pct_change_500_1000", # 30, 31
                             "Pct_change_less_500_med","Pct_change_500_1000_med" # 32, 33
                             )

msaTabSummary[,1] <- unique(msaTabLarge$BEA_Region)
for(i in 1:nrow(msaTabSummary)){
  group <- subset(msaTabLarge, msaTabLarge$BEA_Region==msaTabSummary[i,1])
  msaTabSummary[i,2] <- sum(group$Units_2000,na.rm=T)
  
  msaTabSummary[i,26] <- sum(group$Base_less_500,na.rm=T)
  msaTabSummary[i,27] <- sum(group$Base_500_1000,na.rm=T)
  
  msaTabSummary[i,3] <- sum(group$Base_less_1000,na.rm=T)
  msaTabSummary[i,4] <- sum(group$Base_1000_2000,na.rm=T)
  msaTabSummary[i,5] <- sum(group$Base_2000_3000,na.rm=T)
  msaTabSummary[i,6] <- sum(group$Base_3000_5000,na.rm=T)
  msaTabSummary[i,7] <- sum(group$Base_5000_plus,na.rm=T)
  
  msaTabSummary[i,8] <- sum(group$Units_2020,na.rm=T)-msaTabSummary[i,2]
  
  msaTabSummary[i,28] <- sum(group$Added_less_500,na.rm=T)
  msaTabSummary[i,29] <- sum(group$Added_500_1000,na.rm=T)
  
  msaTabSummary[i,9] <- sum(group$Added_less_1000,na.rm=T)
  msaTabSummary[i,10] <- sum(group$Added_1000_2000,na.rm=T)
  msaTabSummary[i,11] <- sum(group$Added_2000_3000,na.rm=T)
  msaTabSummary[i,12] <- sum(group$Added_3000_5000,na.rm=T)
  msaTabSummary[i,13] <- sum(group$Added_5000_plus,na.rm=T)
  
  msaTabSummary[i,14] <- msaTabSummary[i,8]/msaTabSummary[i,2]
  
  msaTabSummary[i,30] <- msaTabSummary[i,28]/msaTabSummary[i,26]
  msaTabSummary[i,31] <- msaTabSummary[i,29]/msaTabSummary[i,27]
  
  msaTabSummary[i,15] <- msaTabSummary[i,9]/msaTabSummary[i,3]
  msaTabSummary[i,16] <- msaTabSummary[i,10]/msaTabSummary[i,4]
  msaTabSummary[i,17] <- msaTabSummary[i,11]/msaTabSummary[i,5]
  msaTabSummary[i,18] <- msaTabSummary[i,12]/msaTabSummary[i,6]
  msaTabSummary[i,19] <- msaTabSummary[i,13]/msaTabSummary[i,7]
  
  msaTabSummary[i,20] <- median((group$Units_2020-group$Units_2000)/group$Units_2000,na.rm=T)
  
  msaTabSummary[i,32] <- median(group$Added_less_500/group$Base_less_500,na.rm=T)
  msaTabSummary[i,33] <- median(group$Added_500_1000/group$Base_500_1000,na.rm=T)
  
  msaTabSummary[i,21] <- median(group$Added_less_1000/group$Base_less_1000,na.rm=T)
  msaTabSummary[i,22] <- median(group$Added_1000_2000/group$Base_1000_2000,na.rm=T)
  msaTabSummary[i,23] <- median(group$Added_2000_3000/group$Base_2000_3000,na.rm=T)
  msaTabSummary[i,24] <- median(group$Added_3000_5000/group$Base_3000_5000,na.rm=T)
  msaTabSummary[i,25] <- median(group$Added_5000_plus/group$Base_5000_plus,na.rm=T)
}

write.csv(msaTabSummary,"Msa_summary_by_region.csv", row.names = T)

msaTabLong <- msaTabSummary %>%
  pivot_longer(
    cols = c(Pct_change_less_500, Pct_change_500_1000, Pct_change_1000_2000, Pct_change_2000_3000,
             Pct_change_3000_5000, Pct_change_5000_plus),
    names_to = "Unit_Category",
    values_to = "Pct_Change"
  )

msaTabLong$Unit_Category <- factor(msaTabLong$Unit_Category,
                                   levels = c("Pct_change_less_500","Pct_change_500_1000", "Pct_change_1000_2000",
                                              "Pct_change_2000_3000", "Pct_change_3000_5000",
                                              "Pct_change_5000_plus"),
   labels = c("Tracts with fewer than 500 units/sq. mi.","Tracts with 500 to 999 units/sq. mi.",
              "Tracts with 1,000 to 1,999 units/sq. mi.",
   "Tracts with 2,000 to 2,999 units/sq. mi.", "Tracts with 3,000 to 4,999 units/sq. mi.",
   "Tracts with at least 5,000 units/sq. mi."))

plotRegionMSA <- msaTabLong %>%
  ggplot(mapping = aes(x=Region, y = Pct_Change, fill = factor(Unit_Category))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::percent, breaks=c(0,0.25,0.5,0.75,1,1.25), expand = expansion(mult = c(0, 0.1))) +
  labs(title = "A. Change by Region", 
    x = NULL,
    y = NULL) +
  guides(fill = "none") +
  theme(legend.position = "none",
        axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) +
  scale_fill_manual(values = c("#55b748", "#1696d2", "#fdbf11", "#000000", "#d2d2d2", "#ec008b"))

ggsave(filename = "figure_5_panel_a.png", plot = plotRegionMSA, width = 8, height = 3)

#### APPENDIX FIGURE 3, PANEL A, (MSA GROWTH BY REGION -- MEDIAN) ####

msaTabLong <- msaTabSummary %>%
  pivot_longer(
    cols = c(Pct_change_less_500_med, Pct_change_500_1000_med, Pct_change_1000_2000_med, Pct_change_2000_3000_med,
             Pct_change_3000_5000_med, Pct_change_5000_plus_med),
    names_to = "Unit_Category",
    values_to = "Pct_Change"
  )

msaTabLong$Unit_Category <- factor(msaTabLong$Unit_Category,
                                   levels = c("Pct_change_less_500_med",
                                              "Pct_change_500_1000_med", "Pct_change_1000_2000_med",
                                              "Pct_change_2000_3000_med", "Pct_change_3000_5000_med",
                                              "Pct_change_5000_plus_med"),
             labels = c("Tracts with less than 500 units/sq. mi.",
                        "Tracts with 500 to 999 units/sq. mi.", "Tracts with 1,000 to 1,999 units/sq. mi.",
             "Tracts with 2,000 to 2,999 units/sq. mi.", "Tracts with 3,000 to 4,999 units/sq. mi.",
             "Tracts with at least 5,000 units/sq. mi."))

plotRegionMSA <- msaTabLong %>%
  ggplot(mapping = aes(x=Region, y = Pct_Change, fill = factor(Unit_Category))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.1))) +
  labs(title = "A. Change by Region, Median Metropolitan Area in Group", 
       x = NULL,
       y = NULL) +
  guides(fill = "none") +
  theme(legend.position = "none",
        axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=10)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 20)) +
  scale_fill_manual(values = c("#55b748", "#1696d2", "#fdbf11", "#000000", "#d2d2d2", "#ec008b"))

ggsave(filename = "appendix_figure_3_panel_a.png", plot = plotRegionMSA, width = 8, height = 3)

#### FIGURE 5, PANEL B (MSA GROWTH BY HOUSING VALUE) ####

msaTabLarge$Housing_Value_Group <- ifelse(msaTabLarge$Housing_Value>=150000,"≥ $150,000",
                      ifelse(msaTabLarge$Housing_Value>=125000,"$125,000 to $149,999",
                      ifelse(msaTabLarge$Housing_Value>=110000,"$110,000 to $124,999",
                      ifelse(msaTabLarge$Housing_Value>=100000,"$100,000 to $109,999",
                      ifelse(msaTabLarge$Housing_Value>=90000,"$90,000 to $99,999",
                      ifelse(msaTabLarge$Housing_Value>=80000,"$80,000 to $89,999","< $80,000"))))))

msaTabSummary <- data.frame(matrix(ncol=33,nrow=length(unique(msaTabLarge$Housing_Value_Group))))
colnames(msaTabSummary) <- c("Housing_value","Units_2000","Base_less_1000","Base_1000_2000","Base_2000_3000","Base_3000_5000","Base_5000_plus","Added_all","Added_less_1000","Added_1000_2000","Added_2000_3000","Added_3000_5000","Added_5000_plus","Pct_change_all","Pct_change_less_1000","Pct_change_1000_2000","Pct_change_2000_3000","Pct_change_3000_5000","Pct_change_5000_plus","Pct_change_all_med","Pct_change_less_1000_med","Pct_change_1000_2000_med","Pct_change_2000_3000_med","Pct_change_3000_5000_med","Pct_change_5000_plus_med",
                             "Base_less_500","Base_500_1000", # 26, 27
                             "Added_less_500","Added_500_1000", # 28, 29
                             "Pct_change_less_500","Pct_change_500_1000", # 30, 31
                             "Pct_change_less_500_med","Pct_change_500_1000_med" # 32, 33
                             )
msaTabSummary[,1] <- unique(msaTabLarge$Housing_Value_Group)
for(i in 1:nrow(msaTabSummary)){
  group <- subset(msaTabLarge, msaTabLarge$Housing_Value_Group==msaTabSummary[i,1])
  msaTabSummary[i,2] <- sum(group$Units_2000,na.rm=T)
  
  msaTabSummary[i,26] <- sum(group$Base_less_500,na.rm=T)
  msaTabSummary[i,27] <- sum(group$Base_500_1000,na.rm=T)
  
  msaTabSummary[i,3] <- sum(group$Base_less_1000,na.rm=T)
  msaTabSummary[i,4] <- sum(group$Base_1000_2000,na.rm=T)
  msaTabSummary[i,5] <- sum(group$Base_2000_3000,na.rm=T)
  msaTabSummary[i,6] <- sum(group$Base_3000_5000,na.rm=T)
  msaTabSummary[i,7] <- sum(group$Base_5000_plus,na.rm=T)
  
  msaTabSummary[i,8] <- sum(group$Units_2020,na.rm=T)-msaTabSummary[i,2]
  
  msaTabSummary[i,28] <- sum(group$Added_less_500,na.rm=T)
  msaTabSummary[i,29] <- sum(group$Added_500_1000,na.rm=T)
  
  msaTabSummary[i,9] <- sum(group$Added_less_1000,na.rm=T)
  msaTabSummary[i,10] <- sum(group$Added_1000_2000,na.rm=T)
  msaTabSummary[i,11] <- sum(group$Added_2000_3000,na.rm=T)
  msaTabSummary[i,12] <- sum(group$Added_3000_5000,na.rm=T)
  msaTabSummary[i,13] <- sum(group$Added_5000_plus,na.rm=T)
  
  msaTabSummary[i,14] <- msaTabSummary[i,8]/msaTabSummary[i,2]
  
  msaTabSummary[i,30] <- msaTabSummary[i,28]/msaTabSummary[i,26]
  msaTabSummary[i,31] <- msaTabSummary[i,29]/msaTabSummary[i,27]
  
  msaTabSummary[i,15] <- msaTabSummary[i,9]/msaTabSummary[i,3]
  msaTabSummary[i,16] <- msaTabSummary[i,10]/msaTabSummary[i,4]
  msaTabSummary[i,17] <- msaTabSummary[i,11]/msaTabSummary[i,5]
  msaTabSummary[i,18] <- msaTabSummary[i,12]/msaTabSummary[i,6]
  msaTabSummary[i,19] <- msaTabSummary[i,13]/msaTabSummary[i,7]
  
  msaTabSummary[i,20] <- median((group$Units_2020-group$Units_2000)/group$Units_2000,na.rm=T)
  
  msaTabSummary[i,32] <- median(group$Added_less_500/group$Base_less_500,na.rm=T)
  msaTabSummary[i,33] <- median(group$Added_500_1000/group$Base_500_1000,na.rm=T)
  
  msaTabSummary[i,21] <- median(group$Added_less_1000/group$Base_less_1000,na.rm=T)
  msaTabSummary[i,22] <- median(group$Added_1000_2000/group$Base_1000_2000,na.rm=T)
  msaTabSummary[i,23] <- median(group$Added_2000_3000/group$Base_2000_3000,na.rm=T)
  msaTabSummary[i,24] <- median(group$Added_3000_5000/group$Base_3000_5000,na.rm=T)
  msaTabSummary[i,25] <- median(group$Added_5000_plus/group$Base_5000_plus,na.rm=T)
}

write.csv(msaTabSummary,"Msa_summary_by_region_value.csv", row.names = T)

msaTabLong <- msaTabSummary %>%
  pivot_longer(
    cols = c(Pct_change_less_500, Pct_change_500_1000, Pct_change_1000_2000, Pct_change_2000_3000,
             Pct_change_3000_5000, Pct_change_5000_plus),
    names_to = "Unit_Category",
    values_to = "Pct_Change"
  )

msaTabLong$Unit_Category <- factor(msaTabLong$Unit_Category,
                                   levels = c("Pct_change_less_500","Pct_change_500_1000", "Pct_change_1000_2000",
                                              "Pct_change_2000_3000", "Pct_change_3000_5000",
                                              "Pct_change_5000_plus"),
                  labels = c("Tracts with fewer than 500 units/sq. mi.",
                             "Tracts with 500 to 999 units/sq. mi.", "Tracts with 1,000 to 1,999 units/sq. mi.",
                  "Tracts with 2,000 to 2,999 units/sq. mi.", "Tracts with 3,000 to 4,999 units/sq. mi.",
                  "Tracts with at least 5,000 units/sq. mi."))

msaTabLong$Housing_value <- factor(msaTabLong$Housing_value, levels=c("< $80,000","$80,000 to $89,999","$90,000 to $99,999","$100,000 to $109,999","$110,000 to $124,999","$125,000 to $149,999","≥ $150,000"))

plotValueMSA <- msaTabLong %>%
  ggplot(mapping = aes(x = Housing_value, y = Pct_Change, fill = factor(Unit_Category))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.1))) +
  labs(title = "B. Change by Housing Value",
       x = NULL,
       y = NULL) +
  guides(color = "none", fill = guide_legend(ncol = 2)) +
  theme(
    legend.position = "bottom",  # Place legend below the x-axis
    legend.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10)
  ) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c("#55b748", "#1696d2", "#fdbf11", "#000000", "#d2d2d2", "#ec008b"))

ggsave(filename = "figure_5_panel_b.png", plot = plotValueMSA, width = 8, height = 4.5)

#### APPENDIX FIGURE 3, PANEL B (MSA GRAPH BY HOUSING VALUE MEDIANS) ####

msaTabLong <- msaTabSummary %>%
  pivot_longer(
    cols = c(Pct_change_less_500_med,Pct_change_500_1000_med,Pct_change_1000_2000_med,Pct_change_2000_3000_med,Pct_change_3000_5000_med,Pct_change_5000_plus_med),
    names_to = "Unit_Category",
    values_to = "Pct_Change"
  )

msaTabLong$Unit_Category <- factor(msaTabLong$Unit_Category,
                                   levels = c("Pct_change_less_500_med","Pct_change_500_1000_med", "Pct_change_1000_2000_med",
                                              "Pct_change_2000_3000_med", "Pct_change_3000_5000_med",
                                              "Pct_change_5000_plus_med"),
   labels = c("Tracts with fewer than 500 units/sq. mi.",
              "Tracts with 500 to 999 units/sq. mi.", "Tracts with 1,000 to 1,999 units/sq. mi.",
   "Tracts with 2,000 to 2,999 units/sq. mi.", "Tracts with 3,000 to 4,999 units/sq. mi.",
   "Tracts with at least 5,000 units/sq. mi."))

msaTabLong$Housing_value <- factor(msaTabLong$Housing_value, levels=c("< $80,000","$80,000 to $89,999","$90,000 to $99,999","$100,000 to $109,999","$110,000 to $124,999","$125,000 to $149,999","≥ $150,000"))

plotValueMSA <- msaTabLong %>%
  ggplot(mapping = aes(x = Housing_value, y = Pct_Change, fill = factor(Unit_Category))) +
  geom_col(position = "dodge") +
  scale_y_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.1))) +
  labs(title = "B. Change by Housing Value, Median Metropolitan Area in Group",
       x = NULL,
       y = NULL) +
  guides(color = "none", fill = guide_legend(ncol = 2)) +
  theme(
    legend.position = "bottom",  # Place legend below the x-axis
    legend.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10)
  ) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c("#55b748", "#1696d2", "#fdbf11", "#000000", "#d2d2d2", "#ec008b"))

ggsave(filename = "appendix_figure_3_panel_b.png", plot = plotValueMSA, width = 8, height = 4.5)

# Example MSAs
groupMSA <- subset(msaTabLarge, msaTabLarge$Housing_Value_Group=="> $150,000")
nrow(groupMSA)
groupMSA[order(-groupMSA$Units_2000),1][1]

#### COUNTY CORRELATIONS -- BACKGROUND INFORMATION ####

corGrp <- subset(countyLarge, !is.na(countyLarge$hsg_density_2000))

# correlations
cor(countyLarge$race_wht_alone,countyLarge$race_blk_alone) # too close
cor(countyLarge$race_wht_alone,countyLarge$median_value) # OK!
cor(countyLarge$race_wht_alone,countyLarge$median_income) # OK!
cor(countyLarge$race_wht_alone,countyLarge$median_rent) # OK!
cor(countyLarge$race_wht_alone,countyLarge$built_1990_2000) # OK!
cor(countyLarge$race_wht_alone,countyLarge$BA_plus) # OK!
cor(countyLarge$race_wht_alone,countyLarge$mrp_ideology_04_11) # OK!
cor(corGrp$race_wht_alone,corGrp$hsg_density_2000) # OK!
cor(countyLarge$race_blk_alone,countyLarge$median_value) # OK!
cor(countyLarge$race_blk_alone,countyLarge$median_income) # OK!
cor(countyLarge$race_blk_alone,countyLarge$median_rent) # OK!
cor(countyLarge$race_blk_alone,countyLarge$built_1990_2000) # OK!
cor(countyLarge$race_blk_alone,countyLarge$BA_plus) # OK!
cor(countyLarge$race_blk_alone,countyLarge$mrp_ideology_04_11) # OK!
cor(corGrp$race_blk_alone,corGrp$hsg_density_2000) # OK!
cor(countyLarge$median_value,countyLarge$median_income) # too close
cor(countyLarge$median_value,countyLarge$median_rent) # too close
cor(countyLarge$median_value,countyLarge$built_1990_2000) # OK!
cor(countyLarge$median_value,countyLarge$BA_plus) # too close
cor(countyLarge$median_value,countyLarge$mrp_ideology_04_11) # OK!
cor(log(corGrp$median_value),corGrp$hsg_density_2000) # OK!
cor(countyLarge$median_income,countyLarge$median_rent) # too close
cor(countyLarge$median_income,countyLarge$built_1990_2000) # OK!
cor(countyLarge$median_income,countyLarge$BA_plus) # too close
cor(countyLarge$median_income,countyLarge$mrp_ideology_04_11) # OK!
cor(corGrp$median_income,corGrp$hsg_density_2000) # OK!
cor(countyLarge$median_rent,countyLarge$built_1990_2000) # OK!
cor(countyLarge$median_rent,countyLarge$BA_plus) # too close
cor(countyLarge$median_rent,countyLarge$mrp_ideology_04_11) # OK!
cor(corGrp$median_rent,corGrp$hsg_density_2000) # OK!
cor(countyLarge$built_1990_2000,countyLarge$BA_plus) # OK!
cor(countyLarge$built_1990_2000,countyLarge$mrp_ideology_04_11) # OK! -- VERY CLOSE
cor(corGrp$built_1990_2000,corGrp$hsg_density_2000) # OK!
cor(countyLarge$BA_plus,countyLarge$mrp_ideology_04_11) # OK!
cor(corGrp$BA_plus,corGrp$hsg_density_2000) # OK!
cor(corGrp$mrp_ideology_04_11,corGrp$hsg_density_2000) # OK!
# so, we can do: race_wht_alone, median_value, built_1990_2000, hsg_density_2000 to avoid multicollinearity

#### STATES AND REGIONS -- FOR BACKGROUND, and footnote 3 ####

# Bureau of Economic Analysis Regions
regionFarWest <- c("Alaska","Hawaii","California","Nevada","Oregon","Washington")
regionRockyMountain <- c("Colorado","Utah","Wyoming","Idaho","Montana")
regionSouthwest <- c("Arizona","New Mexico","Oklahoma","Texas")
regionPlains <- c("North Dakota","South Dakota","Nebraska","Kansas","Minnesota","Iowa","Missouri")
regionGreatLakes <- c("Illinois","Wisconsin","Indiana","Michigan","Ohio")
regionNewEngland <- c("Maine","New Hampshire","Vermont","Rhode Island","Massachusetts","Connecticut")
regionMideast <- c("New York","New Jersey","Pennsylvania","Delaware","District of Columbia","Maryland")
regionSoutheast <- c("Arkansas","Louisiana","Mississippi","Alabama","Tennessee","Kentucky","West Virginia","Virginia","North Carolina","South Carolina","Georgia","Florida")

for(i in 1:nrow(censusCounties)){
  
censusCounties$region[i] <- ifelse(censusCounties$STATE[i] %in% regionFarWest,"Far West",
                             ifelse(censusCounties$STATE[i] %in% regionRockyMountain,"Rocky Mountains",
                                    ifelse(censusCounties$STATE[i] %in% regionSouthwest,"Southwest",
                                           ifelse(censusCounties$STATE[i] %in% regionPlains,"Plains",
                                                  ifelse(censusCounties$STATE[i] %in% regionGreatLakes,"Great Lakes",
                                                         ifelse(censusCounties$STATE[i] %in% regionNewEngland,"New England",
                                                                ifelse(censusCounties$STATE[i] %in% regionMideast,"Mideast",
                                                                       ifelse(censusCounties$STATE[i] %in% regionSoutheast,"Southeast",""))))))))
}

quantile(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$A41AA2020>50000), 0.05, na.rm=T)
quantile(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$A41AA2020>50000), 0.95, na.rm=T)

median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Far West" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Great Lakes" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Mideast" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="New England" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Plains" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Rocky Mountains" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Southeast" & censusCounties$A41AA2020>50000),na.rm=T)
median(subset(censusCounties$mrp_ideology_17_21_invert_base_0, censusCounties$region=="Southwest" & censusCounties$A41AA2020>50000),na.rm=T)

#### APPENDIX TABLE 5 - INTERPRETATION (CITY VERSION) ####

# Using the regression results from table 1, model 2

# Table of housing growth from 2000 to 2010
interpretationMat <- data.frame(matrix(nrow=12,ncol=3))
rownames(interpretationMat) <- c("GOP State (overall)","Mixed State (overall)","Democratic State (overall)","GOP State (80% pct. density)","Mixed State (80% pct. density)","Democratic State (80% pct. density)","GOP State (95% pct. density)","Mixed State (95% pct. density)","Democratic State (95% pct. density)","GOP State (20% pct. density)","Mixed State (20% pct. density)","Democratic State (20% pct. density)")
colnames(interpretationMat) <- c("10% Most Liberal","Baseline ideology","10% Most Conservative")

densityScore <- mean(cityLarge$hsg_density_2000) # 1594
whiteScore <- mean(cityLarge$race_wht_alone) # 0.62
medValue <- mean(cityLarge$median_value,na.rm=T) # 145303.5
secondPeriod <- 0

coefInter <- 0.12
coefDensity <- -0.13
coefIdeo <- -0.26
coefGOP <- 0.06
coefMix <- 0.02
coefWhite <- -0.04
coefVal <- 0.07
coefPeriod <- -0.07
coefIntersect <- -0.49

# GOP state, baseline ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- mean(cityLarge$mrp_ideology_04_11_invert_base_0,na.rm=T) # 0.64
interpretationMat[1,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue)) + (coefPeriod * secondPeriod) + coefIntersect 
# 0.07 with GOP = 1, Mix = 0, and baseline data
# Mixed state, baseline ideology, overall
GOPScore <- 0
MixScore <- 1
interpretationMat[2,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, baseline ideology, overall
GOPScore <- 0
MixScore <- 0
interpretationMat[3,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, liberal ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[1,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, liberal ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[2,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, liberal ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[3,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, conservative ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[1,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, conservative ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[2,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, conservative ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[3,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# Higher density cities

densityScore <- quantile(cityLarge$hsg_density_2000,0.8) # 2326

# GOP state, baseline ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- mean(cityLarge$mrp_ideology_04_11_invert_base_0,na.rm=T) # 0.64
interpretationMat[4,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# 0.07 with GOP = 1, Mix = 0, and baseline data
# Mixed state, baseline ideology, overall
GOPScore <- 0
MixScore <- 1
interpretationMat[5,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, baseline ideology, overall
GOPScore <- 0
MixScore <- 0
interpretationMat[6,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, liberal ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[4,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, liberal ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[5,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, liberal ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[6,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, conservative ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[4,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, conservative ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[5,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, conservative ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[6,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# Highest density cities

densityScore <- quantile(cityLarge$hsg_density_2000,0.95) # 4166

# GOP state, baseline ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- mean(cityLarge$mrp_ideology_04_11_invert_base_0,na.rm=T) # 0.64
interpretationMat[7,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# 0.07 with GOP = 1, Mix = 0, and baseline data
# Mixed state, baseline ideology, overall
GOPScore <- 0
MixScore <- 1
interpretationMat[8,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, baseline ideology, overall
GOPScore <- 0
MixScore <- 0
interpretationMat[9,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, liberal ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[7,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, liberal ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[8,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, liberal ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[9,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, conservative ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[7,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, conservative ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[8,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, conservative ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[9,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# Lowest density cities

densityScore <- quantile(cityLarge$hsg_density_2000,0.2) # 687

# GOP state, baseline ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- mean(cityLarge$mrp_ideology_04_11_invert_base_0,na.rm=T) # 0.64
interpretationMat[10,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# 0.07 with GOP = 1, Mix = 0, and baseline data
# Mixed state, baseline ideology, overall
GOPScore <- 0
MixScore <- 1
interpretationMat[11,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, baseline ideology, overall
GOPScore <- 0
MixScore <- 0
interpretationMat[12,2] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, liberal ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[10,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, liberal ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[11,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, liberal ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.9,na.rm=T) # 0.86
interpretationMat[12,1] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

# GOP state, conservative ideology, overall
GOPScore <- 1
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[10,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Mixed state, conservative ideology, overall
GOPScore <- 0
MixScore <- 1
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[11,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 
# Dem state, conservative ideology, overall
GOPScore <- 0
MixScore <- 0
ideologyScore <- quantile(cityLarge$mrp_ideology_04_11_invert_base_0,0.1,na.rm=T) # 0.45
interpretationMat[12,3] <- (coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect) 

interpretationMat

#### For background (in text) -- overall medians, and just cities with non-negative housing growth over the 2000 to 2010 period ####

# totals
nrow(cityLarge) # 486
nrow(subset(cityLarge, (cityLarge$units_2010-cityLarge$units_2000)>=0)) # 440
nrow(subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2010)>=0)) # 441
nrow(subset(cityLarge, (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 451
summary(subset(cityLarge$chg_00_20, (cityLarge$units_2020-cityLarge$units_2000)<0)) # -2875 units for median city
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, (cityLarge$units_2020-cityLarge$units_2000)<0)) # -5.1% for median city
summary(subset(cityLarge$units_2000, (cityLarge$units_2010-cityLarge$units_2000)<0)) # 40,293
(subset(c(cityLarge$NAME10, cityLarge$STATE), (cityLarge$units_2020-cityLarge$units_2000)<0))

densityScore <- mean(subset(cityLarge$hsg_density_2000, (cityLarge$units_2010-cityLarge$units_2000)>=0)) # 1549 # mean density
densityScore <- quantile(subset(cityLarge$hsg_density_2000, (cityLarge$units_2010-cityLarge$units_2000)>=0),0.95) # 4180 # high density
whiteScore <- mean(subset(cityLarge$race_wht_alone, (cityLarge$units_2010-cityLarge$units_2000)>=0)) # 0.62
medValue <- mean(subset(cityLarge$median_value, (cityLarge$units_2010-cityLarge$units_2000)>=0), na.rm=T) # 152360
secondPeriod <- 0

coefInter <- 0.122
coefDensity <- -0.135
coefIdeo <- -0.213
coefGOP <- 0.052
coefMix <- 0.018
coefWhite <- -0.046
coefVal <- 0.045
coefPeriod <- -0.075
coefIntersect <- -0.18

GOPScore <- 1
MixScore <- 0
ideologyScore <- mean(subset(cityLarge$mrp_ideology_04_11_invert_base_0,(cityLarge$units_2010-cityLarge$units_2000)>=0),na.rm=T) # 0.63

(coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect)
# mean ideology, GOP state, mean density = 0.16
# mean ideology, GOP state, high density = 0.004

ideologyScore <- quantile(subset(cityLarge$mrp_ideology_04_11_invert_base_0,(cityLarge$units_2010-cityLarge$units_2000)>=0),0.9,na.rm=T) # 0.86
(coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect)
# liberal ideology, GOP state, mean density = 0.15
# liberal ideology, GOP state, high density = 0.07

ideologyScore <- quantile(subset(cityLarge$mrp_ideology_04_11_invert_base_0,(cityLarge$units_2010-cityLarge$units_2000)>=0),0.1,na.rm=T) # 0.45
(coefInter * (densityScore/1000) * ideologyScore) + (coefDensity * densityScore/1000) + (coefIdeo * ideologyScore) + (coefGOP * GOPScore) + (coefMix * MixScore) + (coefWhite * whiteScore) + (coefVal * log(medValue) + (coefPeriod * secondPeriod) + coefIntersect)
# conservative ideology, GOP state, mean density = 0.16
# conservative ideology, GOP state, high density = -0.05

summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 16.5% median
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$hsg_density_2000 > 4166 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 10.8% median
nrow(subset(cityLarge, cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 100 cities
summary(subset(cityLarge$mrp_ideology_04_11_invert_base_0, cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0))
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 8.9% median

## Liberal cities
# Overall, 25 most liberal
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.931)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.931)) # 12.2% median
# Overall, 25 most liberal, with no housing loss
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.92 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.92 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 12.6% median
# High densities, 25 most liberal
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.9175 & cityLarge$hsg_density_2000 > 2000)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.9175 & cityLarge$hsg_density_2000 > 2000)) # 11.7% median
# High densities, 25 most liberal, with no housing loss
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.88 & cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 > 0.88 & cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 11.7% median

## Conservative cities
# Overall, 25 most conservative
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.391)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.391)) # 25.5% median
# Overall, 25 most conservative, with no housing loss
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.391 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.391 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 25.5% median
# High densities, 25 most conservative
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.653 & cityLarge$hsg_density_2000 > 2000)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.653 & cityLarge$hsg_density_2000 > 2000)) # 4.4% median
# High densities, 25 most conservative, with no housing loss
nrow(subset(cityLarge, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.653 & cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 25 cities
summary(subset(cityLarge$chg_00_20/cityLarge$units_2000, cityLarge$mrp_ideology_04_11_invert_base_0 < 0.653 & cityLarge$hsg_density_2000 > 2000 & (cityLarge$units_2020-cityLarge$units_2000)>=0)) # 4.4% median

#### Ideology by Metro - footnote 7 ####

metroIdeology <- data.frame(matrix(ncol=22, nrow=length(unique(censusCounties$NAME))))
colnames(metroIdeology) <- c("metro","units_2000","units_2020","weighted_ideology_04_11","weighted_ideology_17_21","units_2000_tracts_above_1000","units_2020_tracts_above_1000","units_2000_tracts_above_2000","units_2020_tracts_above_2000","units_2000_tracts_above_3000","units_2020_tracts_above_3000","units_2000_tracts_above_5000","units_2020_tracts_above_5000","units_2000_tracts_above_500","units_2020_tracts_above_500",
                             "n_tracts_total","n_tracts_0_500","n_tracts_500_1000","n_tracts_1000_2000",
                             "n_tracts_2000_3000","n_tracts_3000_5000","n_tracts_5000_plus")
metroIdeology[,1] <- unique(censusCounties$NAME)
for(i in 1:nrow(metroIdeology)){
  group <- subset(censusCounties, censusCounties$NAME==metroIdeology[i,1])
  metroIdeology[i,2] <- sum(group$A41AA2000,na.rm=T)
  metroIdeology[i,3] <- sum(group$A41AA2020,na.rm=T)
  metroIdeology[i,4] <- sum(group$mrp_ideology_04_11_invert_base_0*group$A41AA2000,na.rm=T)/sum(group$A41AA2000,na.rm=T)
  metroIdeology[i,5] <- sum(group$mrp_ideology_17_21_invert_base_0*group$A41AA2020,na.rm=T)/sum(group$A41AA2020,na.rm=T)
  metroIdeology[i,6] <- sum(group$hsg_2000_in_tracts_above_1000,na.rm=T)
  metroIdeology[i,7] <- sum(group$hsg_2020_in_tracts_above_1000_in_2000,na.rm=T)
  metroIdeology[i,8] <- sum(group$hsg_2000_in_tracts_above_2000,na.rm=T)
  metroIdeology[i,9] <- sum(group$hsg_2020_in_tracts_above_2000_in_2000,na.rm=T)
  metroIdeology[i,10] <- sum(group$hsg_2000_in_tracts_above_3000,na.rm=T)
  metroIdeology[i,11] <- sum(group$hsg_2020_in_tracts_above_3000_in_2000,na.rm=T)
  metroIdeology[i,12] <- sum(group$hsg_2000_in_tracts_above_5000,na.rm=T)
  metroIdeology[i,13] <- sum(group$hsg_2020_in_tracts_above_5000_in_2000,na.rm=T)
  metroIdeology[i,14] <- sum(group$hsg_2000_in_tracts_above_500,na.rm=T)
  metroIdeology[i,15] <- sum(group$hsg_2020_in_tracts_above_500_in_2000,na.rm=T)
  groupA <- subset(censusTracts, paste(censusTracts$STATEA, censusTracts$COUNTYA,sep="-") %in% group$state_county)
  metroIdeology[i,16] <- nrow(groupA)
  metroIdeology[i,17] <- nrow(subset(groupA, groupA$density_00<500))
  metroIdeology[i,18] <- nrow(subset(groupA, groupA$density_00>=500 & groupA$density_00<1000))
  metroIdeology[i,19] <- nrow(subset(groupA, groupA$density_00>=1000 & groupA$density_00<2000))
  metroIdeology[i,20] <- nrow(subset(groupA, groupA$density_00>=2000 & groupA$density_00<3000))
  metroIdeology[i,21] <- nrow(subset(groupA, groupA$density_00>=3000 & groupA$density_00<5000))
  metroIdeology[i,22] <- nrow(subset(groupA, groupA$density_00>=5000))
}

metroIdeology$growth_overall <- (metroIdeology$units_2020-metroIdeology$units_2000)/metroIdeology$units_2000
metroIdeology$growth_0_500 <- ((metroIdeology$units_2020-metroIdeology$units_2020_tracts_above_500)-(metroIdeology$units_2000-metroIdeology$units_2000_tracts_above_500))/(metroIdeology$units_2000-metroIdeology$units_2000_tracts_above_500)
metroIdeology$growth_500_1000 <- ((metroIdeology$units_2020_tracts_above_500-metroIdeology$units_2020_tracts_above_1000)-(metroIdeology$units_2000_tracts_above_500-metroIdeology$units_2000_tracts_above_1000))/(metroIdeology$units_2000_tracts_above_500-metroIdeology$units_2000_tracts_above_1000)
metroIdeology$growth_1000_2000 <- ((metroIdeology$units_2020_tracts_above_1000-metroIdeology$units_2020_tracts_above_2000)-(metroIdeology$units_2000_tracts_above_1000-metroIdeology$units_2000_tracts_above_2000))/(metroIdeology$units_2000_tracts_above_1000-metroIdeology$units_2000_tracts_above_2000)
metroIdeology$growth_2000_3000 <- ((metroIdeology$units_2020_tracts_above_2000-metroIdeology$units_2020_tracts_above_3000)-(metroIdeology$units_2000_tracts_above_2000-metroIdeology$units_2000_tracts_above_3000))/(metroIdeology$units_2000_tracts_above_2000-metroIdeology$units_2000_tracts_above_3000)
metroIdeology$growth_3000_5000 <- ((metroIdeology$units_2020_tracts_above_3000-metroIdeology$units_2020_tracts_above_5000)-(metroIdeology$units_2000_tracts_above_3000-metroIdeology$units_2000_tracts_above_5000))/(metroIdeology$units_2000_tracts_above_3000-metroIdeology$units_2000_tracts_above_5000)
metroIdeology$growth_5000_plus <- (metroIdeology$units_2020_tracts_above_5000-metroIdeology$units_2000_tracts_above_5000)/metroIdeology$units_2000_tracts_above_5000

# Nominal growth per tract
metroIdeology$growth_overall_nom <- (metroIdeology$units_2020-metroIdeology$units_2000)/metroIdeology$n_tracts_total
metroIdeology$growth_0_500_nom <- ((metroIdeology$units_2020-metroIdeology$units_2020_tracts_above_500)-(metroIdeology$units_2000-metroIdeology$units_2000_tracts_above_500))/metroIdeology$n_tracts_0_500
metroIdeology$growth_500_1000_nom <- ((metroIdeology$units_2020_tracts_above_500-metroIdeology$units_2020_tracts_above_1000)-(metroIdeology$units_2000_tracts_above_500-metroIdeology$units_2000_tracts_above_1000))/metroIdeology$n_tracts_500_1000
metroIdeology$growth_1000_2000_nom <- ((metroIdeology$units_2020_tracts_above_1000-metroIdeology$units_2020_tracts_above_2000)-(metroIdeology$units_2000_tracts_above_1000-metroIdeology$units_2000_tracts_above_2000))/metroIdeology$n_tracts_1000_2000
metroIdeology$growth_2000_3000_nom <- ((metroIdeology$units_2020_tracts_above_2000-metroIdeology$units_2020_tracts_above_3000)-(metroIdeology$units_2000_tracts_above_2000-metroIdeology$units_2000_tracts_above_3000))/metroIdeology$n_tracts_2000_3000
metroIdeology$growth_3000_5000_nom <- ((metroIdeology$units_2020_tracts_above_3000-metroIdeology$units_2000_tracts_above_5000)-(metroIdeology$units_2000_tracts_above_3000-metroIdeology$units_2000_tracts_above_5000))/metroIdeology$n_tracts_3000_5000
metroIdeology$growth_5000_plus_nom <- (metroIdeology$units_2020_tracts_above_5000-metroIdeology$units_2000_tracts_above_5000)/metroIdeology$n_tracts_5000_plus

metroIdeology <- metroIdeology[order(-metroIdeology$units_2000),]

metroIdeology[1:20,]

nrow(subset(metroIdeology,metroIdeology$units_2000>400000)) #51

metroIdeology %>% filter(units_2000>400000) %>%
  ggplot(mapping = aes(x = weighted_ideology_04_11, y = (units_2020-units_2000)/units_2000)) +
  geom_point(alpha = 0.2) +
  geom_smooth(color =  "#ec008b") +
  scale_x_continuous(expand = expansion(mult = c(0.002, 0))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.002)),
                     #limits = c(0, 20000), 
                     labels = scales::percent) +  
  labs(x = "Ideology (2004-2011)",
       y = "Change in housing units (2000-2020)") +
  scatter_grid()

metroIdeology %>% filter(units_2000>400000) %>%
  ggplot(mapping = aes(x = weighted_ideology_04_11, y = (units_2020_tracts_above_3000-units_2000_tracts_above_3000)/units_2000_tracts_above_3000)) +
  geom_point(alpha = 0.2) +
  geom_smooth(color =  "#ec008b") +
  scale_x_continuous(expand = expansion(mult = c(0.002, 0))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.002)),
                     #limits = c(0, 20000), 
                     labels = scales::percent) +  
  labs(x = "Ideology (2004-2011)",
       y = "Change in housing units (2000-2020)") +
  scatter_grid()

liberalQuant <- quantile(subset(metroIdeology$weighted_ideology_04_11,metroIdeology$units_2000>100000),0.75)
modQuant <- quantile(subset(metroIdeology$weighted_ideology_04_11,metroIdeology$units_2000>100000),0.5)
conservativeQuant <- quantile(subset(metroIdeology$weighted_ideology_04_11,metroIdeology$units_2000>100000),0.25)

# Summary, among metros with at least 100,000 units in 2000
summaryMetTab <- data.frame(matrix(ncol=3,nrow=28))
colnames(summaryMetTab) <- c("Group","Value","Portion")
summaryMetTab[,1] <- c(rep("Liberal",7),rep("Liberal-Moderate",7),rep("Moderate-Conservative",7),rep("Conservative",7))
summaryMetTab[,3] <- c("Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus",
                       "Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus",
                       "Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus",
                             "Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus")

summaryMetTab[1,2] <- median(subset(metroIdeology$growth_overall, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T) # 17.9%
summaryMetTab[8,2] <- median(subset(metroIdeology$growth_overall, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>modQuant & metroIdeology$weighted_ideology_04_11<=liberalQuant),na.rm=T)
summaryMetTab[15,2] <- median(subset(metroIdeology$growth_overall, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=conservativeQuant & metroIdeology$weighted_ideology_04_11<=modQuant),na.rm=T)
summaryMetTab[22,2] <- median(subset(metroIdeology$growth_overall, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<conservativeQuant),na.rm=T) # 24.8%

summaryMetTab[2,2] <- median(subset(metroIdeology$growth_0_500, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T) #29.7%
summaryMetTab[9,2] <- median(subset(metroIdeology$growth_0_500, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>modQuant & metroIdeology$weighted_ideology_04_11<=liberalQuant),na.rm=T)
summaryMetTab[16,2] <- median(subset(metroIdeology$growth_0_500, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=conservativeQuant & metroIdeology$weighted_ideology_04_11<=modQuant),na.rm=T)
summaryMetTab[23,2] <- median(subset(metroIdeology$growth_0_500, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<conservativeQuant),na.rm=T)

summaryMetTab[3,2] <- median(subset(metroIdeology$growth_500_1000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T)
summaryMetTab[10,2] <- median(subset(metroIdeology$growth_500_1000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=liberalQuant & metroIdeology$weighted_ideology_04_11>modQuant),na.rm=T)
summaryMetTab[17,2] <- median(subset(metroIdeology$growth_500_1000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=modQuant & metroIdeology$weighted_ideology_04_11>=conservativeQuant),na.rm=T)
summaryMetTab[24,2] <- median(subset(metroIdeology$growth_500_1000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<conservativeQuant),na.rm=T)

summaryMetTab[4,2] <- median(subset(metroIdeology$growth_1000_2000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T)
summaryMetTab[11,2] <- median(subset(metroIdeology$growth_1000_2000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>modQuant & metroIdeology$weighted_ideology_04_11<=liberalQuant),na.rm=T)
summaryMetTab[18,2] <- median(subset(metroIdeology$growth_1000_2000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=conservativeQuant & metroIdeology$weighted_ideology_04_11<=modQuant),na.rm=T)
summaryMetTab[25,2] <- median(subset(metroIdeology$growth_1000_2000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<conservativeQuant),na.rm=T)

summaryMetTab[5,2] <- median(subset(metroIdeology$growth_2000_3000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T)
summaryMetTab[12,2] <- median(subset(metroIdeology$growth_2000_3000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>modQuant & metroIdeology$weighted_ideology_04_11<=liberalQuant),na.rm=T)
summaryMetTab[19,2] <- median(subset(metroIdeology$growth_2000_3000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=conservativeQuant & metroIdeology$weighted_ideology_04_11<=modQuant),na.rm=T)
summaryMetTab[26,2] <- median(subset(metroIdeology$growth_2000_3000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<conservativeQuant),na.rm=T)

summaryMetTab[6,2] <- median(subset(metroIdeology$growth_3000_5000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T)
summaryMetTab[13,2] <- median(subset(metroIdeology$growth_3000_5000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>modQuant & metroIdeology$weighted_ideology_04_11<=liberalQuant),na.rm=T)
summaryMetTab[20,2] <- median(subset(metroIdeology$growth_3000_5000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=conservativeQuant & metroIdeology$weighted_ideology_04_11<=modQuant),na.rm=T)
summaryMetTab[27,2] <- median(subset(metroIdeology$growth_3000_5000, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<conservativeQuant),na.rm=T)

summaryMetTab[7,2] <- median(subset(metroIdeology$growth_5000_plus, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>liberalQuant),na.rm=T)
summaryMetTab[14,2] <- median(subset(metroIdeology$growth_5000_plus, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>modQuant & metroIdeology$weighted_ideology_04_11<=liberalQuant),na.rm=T)
summaryMetTab[21,2] <- median(subset(metroIdeology$growth_5000_plus, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=conservativeQuant & metroIdeology$weighted_ideology_04_11<=modQuant),na.rm=T)
summaryMetTab[28,2] <- median(subset(metroIdeology$growth_5000_plus, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab$Portion <- factor(summaryMetTab$Portion, levels=c(summaryMetTab$Portion[1:7]))
summaryMetTab$Group <- factor(summaryMetTab$Group, levels=c("Liberal","Liberal-Moderate","Moderate-Conservative","Conservative"))

summaryMetTab %>%
  ggplot(mapping = aes(x=Portion, y = Value, fill = factor(Group))) +
  geom_col(position = "dodge") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = "Tract density",
       y = NULL)

summary(subset(metroIdeology$weighted_ideology_04_11,metroIdeology$units_2000>100000))
length(subset(metroIdeology$weighted_ideology_04_11,metroIdeology$units_2000>100000)) # 152

#### NOMINAL GROWTH PER TRACT ####

# Summary, among metros with at least 100,000 units in 2000
summaryMetTab <- data.frame(matrix(ncol=3,nrow=28))
colnames(summaryMetTab) <- c("Group","Value","Portion")
summaryMetTab[,1] <- c(rep("Liberal",7),rep("Liberal-Moderate",7),rep("Moderate-Conservative",7),rep("Conservative",7))
summaryMetTab[,3] <- c("Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus",
                       "Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus",
                       "Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus",
                       "Overall","< 500","500-1000","1000-2000","2000-3000","3000-5000","5000 plus")

summaryMetTab[1,2] <- median(subset(metroIdeology$growth_overall_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[8,2] <- median(subset(metroIdeology$growth_overall_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<modQuant & metroIdeology$weighted_ideology_04_11>=liberalQuant),na.rm=T)
summaryMetTab[15,2] <- median(subset(metroIdeology$growth_overall_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=conservativeQuant & metroIdeology$weighted_ideology_04_11>=modQuant),na.rm=T)
summaryMetTab[22,2] <- median(subset(metroIdeology$growth_overall_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab[2,2] <- median(subset(metroIdeology$growth_0_500_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[9,2] <- median(subset(metroIdeology$growth_0_500_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<modQuant & metroIdeology$weighted_ideology_04_11>=liberalQuant),na.rm=T)
summaryMetTab[16,2] <- median(subset(metroIdeology$growth_0_500_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=conservativeQuant & metroIdeology$weighted_ideology_04_11>=modQuant),na.rm=T)
summaryMetTab[23,2] <- median(subset(metroIdeology$growth_0_500_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab[3,2] <- median(subset(metroIdeology$growth_500_1000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[10,2] <- median(subset(metroIdeology$growth_500_1000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=liberalQuant & metroIdeology$weighted_ideology_04_11<modQuant),na.rm=T)
summaryMetTab[17,2] <- median(subset(metroIdeology$growth_500_1000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>=modQuant & metroIdeology$weighted_ideology_04_11<=conservativeQuant),na.rm=T)
summaryMetTab[24,2] <- median(subset(metroIdeology$growth_500_1000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab[4,2] <- median(subset(metroIdeology$growth_1000_2000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[11,2] <- median(subset(metroIdeology$growth_1000_2000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<modQuant & metroIdeology$weighted_ideology_04_11>=liberalQuant),na.rm=T)
summaryMetTab[18,2] <- median(subset(metroIdeology$growth_1000_2000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=conservativeQuant & metroIdeology$weighted_ideology_04_11>=modQuant),na.rm=T)
summaryMetTab[25,2] <- median(subset(metroIdeology$growth_1000_2000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab[5,2] <- median(subset(metroIdeology$growth_2000_3000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[12,2] <- median(subset(metroIdeology$growth_2000_3000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<modQuant & metroIdeology$weighted_ideology_04_11>=liberalQuant),na.rm=T)
summaryMetTab[19,2] <- median(subset(metroIdeology$growth_2000_3000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=conservativeQuant & metroIdeology$weighted_ideology_04_11>=modQuant),na.rm=T)
summaryMetTab[26,2] <- median(subset(metroIdeology$growth_2000_3000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab[6,2] <- median(subset(metroIdeology$growth_3000_5000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[13,2] <- median(subset(metroIdeology$growth_3000_5000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<modQuant & metroIdeology$weighted_ideology_04_11>=liberalQuant),na.rm=T)
summaryMetTab[20,2] <- median(subset(metroIdeology$growth_3000_5000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=conservativeQuant & metroIdeology$weighted_ideology_04_11>=modQuant),na.rm=T)
summaryMetTab[27,2] <- median(subset(metroIdeology$growth_3000_5000_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab[7,2] <- median(subset(metroIdeology$growth_5000_plus_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<liberalQuant),na.rm=T)
summaryMetTab[14,2] <- median(subset(metroIdeology$growth_5000_plus_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<modQuant & metroIdeology$weighted_ideology_04_11>=liberalQuant),na.rm=T)
summaryMetTab[21,2] <- median(subset(metroIdeology$growth_5000_plus_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11<=conservativeQuant & metroIdeology$weighted_ideology_04_11>=modQuant),na.rm=T)
summaryMetTab[28,2] <- median(subset(metroIdeology$growth_5000_plus_nom, metroIdeology$units_2000>100000 & metroIdeology$weighted_ideology_04_11>conservativeQuant),na.rm=T)

summaryMetTab$Portion <- factor(summaryMetTab$Portion, levels=c(summaryMetTab$Portion[1:7]))
summaryMetTab$Group <- factor(summaryMetTab$Group, levels=c("Liberal","Liberal-Moderate","Moderate-Conservative","Conservative"))

summaryMetTab %>%
  ggplot(mapping = aes(x=Portion, y = Value, fill = factor(Group))) +
  geom_col(position = "dodge") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = "Tract density",
       y = NULL)

median(subset(metroIdeology$n_tracts_5000_plus, metroIdeology$units_2000>100000 & metroIdeology$n_tracts_5000_plus>0 & metroIdeology$weighted_ideology_04_11>conservativeQuant)) # 1, with growth of 39 units/tract
median(subset(metroIdeology$n_tracts_5000_plus, metroIdeology$units_2000>100000 & metroIdeology$n_tracts_5000_plus>0 & metroIdeology$weighted_ideology_04_11<liberalQuant)) # 24, with growth of 135 units/tract

median(subset(metroIdeology$n_tracts_0_500, metroIdeology$units_2000>100000 & metroIdeology$n_tracts_0_500>0 & metroIdeology$weighted_ideology_04_11>conservativeQuant)) # 50, with growth of 691 units/tract
median(subset(metroIdeology$n_tracts_0_500, metroIdeology$units_2000>100000 & metroIdeology$n_tracts_0_500>0 & metroIdeology$weighted_ideology_04_11<liberalQuant)) # 90, with growth of 514 units/trat

subset(metroIdeology$metro, metroIdeology$units_2000>500000 & metroIdeology$weighted_ideology_04_11>conservativeQuant) # Nashville, OK City, Salt Lake
subset(metroIdeology$metro, metroIdeology$units_2000>500000 & metroIdeology$weighted_ideology_04_11<liberalQuant) # NYC, LA, Chicago, etc.