set up and load libraries
setwd("D:/CP6570 Socioeconomic GIS/FinalProject")
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(vtable)
## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
library(jtools)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
read in data
# target data with census tract id
df <- read.csv("CensusTract_Regression.csv")
# join data: at the census tract level
UERate <- read.csv("UERate.csv") # unemployment rate variable, DP03_0009E
pubinc <- read.csv("pubinc.csv") # household with public assistance income (cash assistance and social insurance)
edu <- read.csv("EduAttain.csv") # percent population a Bachelor's degree
mg <- read.csv("hsmortg.csv") # percent household with a mortgage
# join data: at the zip code level
est <- read.csv("ct_est.csv") # number of business establishment
data cleaning
# keep 11-digit census tract id
UERate$GEOID <- substring(UERate$GEOID, 10)
pubinc$GEOID <- substring(pubinc$GEOID, 10)
edu$GEOID <- substring(edu$GEOID, 10)
mg$GEOID <- substring(mg$GEOID, 10)
# select interest columns only
est <- select(est, GEOID, zcta_est)
# convert variable to numeric
UERate$UERate <- as.numeric(UERate$UERate)
## Warning: NAs introduced by coercion
# Merge data (educational attainment, public assistance income, unemployment rate, mortgage)
df <- merge(df, edu, by = "GEOID")
df <- merge(df, pubinc, by = "GEOID")
df <- merge(df, UERate, by = "GEOID")
df <- merge(df, mg, by = "GEOID")
df <- merge(df, est, by = "GEOID")
# percent population with a bachelor's degree
df$B.pct <- (df$B / (df$LHSD + df$HSD + df$SC + df$B)) * 100
# percent owner-occupied housing units with a mortgage
df$mgpct <- (df$hsmort / df$totocc) * 100
# low-paying job access index
df$Norm_AC <- (df$AC - min(df$AC, na.rm=TRUE)) / (max(df$AC, na.rm = TRUE) - min(df$AC, na.rm = TRUE))
# percent household with public assistance income
df$PubincPct <- (df$pubinc / df$tot) * 100
# percent household without a vehicle (this variable is extracted from lab2)
df$NoCar <- df$NoCar * 100
# low-paying jobs to low income workers ratio
df$lowratio <- df$LowInc_job/df$LowInc_worker
# low-paying jobs dissimilarity
df$LowIncJob_DI <- abs((df$UE / sum(df$UE, na.rm = TRUE)) - (df$LowInc_job / sum(df$LowInc_job, na.rm = TRUE))) / 2 * 100
df$Norm_LPDI <- (df$LowIncJob_DI - min(df$LowIncJob_DI, na.rm = TRUE)) / (max(df$LowIncJob_DI, na.rm = TRUE) - min(df$LowIncJob_DI, na.rm = TRUE))
# initial dissimilarity
df$before_diff <- (df$UE / sum(df$UE, na.rm = TRUE)) - (df$LowInc_job / sum(df$LowInc_job, na.rm = TRUE))
# number of people need to relocated from each tract
#(positive value = number of pp need to move out from that area, vice versa)
df$pp_need_relocate <- round(sum(df$UE, na.rm=TRUE) * df$before_diff)
# number of people in each tract after relocation
df$after_pop <- df$UE - df$pp_need_relocate
# dissimilarity after people relocate
df$after_diff <- (df$after_pop / sum(df$after_pop, na.rm=TRUE)) -
(df$LowInc_job / sum(df$LowInc_job, na.rm = TRUE))
# before and after job-pop ratio
df$before_job_pop_ratio <- sum(df$LowInc_job, na.rm = TRUE) / sum(df$UE, na.rm = TRUE)
df$after_job_pop_ratio <- sum(df$LowInc_job, na.rm = TRUE) / sum(df$after_pop)
# ideal job-pop ratio
df$ideal <- sum(df$before_job) / sum(df$before_pop)
#test
test <- select(df, UE, LowInc_job, before_diff, pp_need_relocate, after_pop, after_diff)
imput missing values
# if the number of missing values in a interest variables are smaller than ten percent of total observations, then replace it with median value from that interest variable
if ((sum(is.na(df$UERate)) / nrow(df)) <= 0.1) {
a <- na.omit(df$UERate)
median_value <- median(a)
df$UERate[is.na(df$UERate)] <- median_value
}
if ((sum(is.na(df$B.pct)) / nrow(df)) <= 0.1) {
a <- na.omit(df$B.pct)
median_value <- median(a)
df$B.pct[is.na(df$B.pct)] <- median_value
}
if ((sum(is.na(df$mgpct)) / nrow(df)) <= 0.1) {
a <- na.omit(df$mgpct)
median_value <- median(a)
df$mgpct[is.na(df$mgpct)] <- median_value
}
if ((sum(is.na(df$NoCar)) / nrow(df)) <= 0.1) {
a <- na.omit(df$NoCar)
median_value <- median(a)
df$NoCar[is.na(df$NoCar)] <- median_value
}
if ((sum(is.na(df$PubincPct)) / nrow(df)) <= 0.1) {
a <- na.omit(df$PubincPct)
median_value <- median(a)
df$PubincPct[is.na(df$PubincPct)] <- median_value
}
if ((sum(is.na(df$Norm_AC)) / nrow(df)) <= 0.1) {
a <- na.omit(df$Norm_AC)
median_value <- median(a)
df$Norm_AC[is.na(df$Norm_AC)] <- median_value
}
if ((sum(is.na(df$LowIncJob_DI)) / nrow(df)) <= 0.1) {
a <- na.omit(df$LowIncJob_DI)
median_value <- median(a)
df$LowIncJob_DI[is.na(df$LowIncJob_DI)] <- median_value
}
if ((sum(is.na(df$lowratio)) / nrow(df)) <= 0.1) {
a <- na.omit(df$lowratio)
median_value <- median(a)
df$lowratio[is.na(df$lowratio)] <- median_value
}
if ((sum(is.na(df$zcta_est)) / nrow(df)) <= 0.1) {
a <- na.omit(df$zcta_est)
median_value <- median(a)
df$zcta_est[is.na(df$zcta_est)] <- median_value
}
if ((sum(is.na(df$Norm_LPDI)) / nrow(df)) <= 0.1) {
a <- na.omit(df$Norm_LPDI)
median_value <- median(a)
df$Norm_LPDI[is.na(df$Norm_LPDI)] <- median_value
}
if ((sum(is.na(df$pp_need_relocate)) / nrow(df)) <= 0.1) {
a <- na.omit(df$pp_need_relocate)
median_value <- median(a)
df$pp_need_relocate[is.na(df$pp_need_relocate)] <- median_value
}
rename variables
# rename field names
df <- df %>%
rename(c(BAPct=B.pct, MortgPct=mgpct, AI=Norm_AC, NoCarPct=NoCar, LPJD = LowIncJob_DI, LPJR=lowratio))
# select interest variables only
reg <- select(df, GEOID, UERate, BAPct, MortgPct, NoCarPct, PubincPct, AI, LPJD, LPJR, Norm_LPDI, pp_need_relocate)
build regression model
model <- lm(UERate ~
BAPct
+ MortgPct
+ NoCarPct
+ PubincPct
+ AI
+ LPJD
+ LPJR
+ pp_need_relocate
, data = reg)
summary(model)
##
## Call:
## lm(formula = UERate ~ BAPct + MortgPct + NoCarPct + PubincPct +
## AI + LPJD + LPJR + pp_need_relocate, data = reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.021 -1.907 -0.446 1.365 95.840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.229e+00 3.290e-01 12.854 < 2e-16 ***
## BAPct -2.301e-02 5.219e-03 -4.408 1.08e-05 ***
## MortgPct -3.910e-02 4.462e-03 -8.762 < 2e-16 ***
## NoCarPct 9.898e-02 1.009e-02 9.811 < 2e-16 ***
## PubincPct 4.774e-02 2.769e-02 1.724 0.084884 .
## AI 2.797e+00 8.100e-01 3.453 0.000563 ***
## LPJD 1.061e+02 3.436e+00 30.875 < 2e-16 ***
## LPJR 8.351e-02 2.031e-02 4.112 4.04e-05 ***
## pp_need_relocate 2.015e-02 4.764e-04 42.295 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.584 on 2787 degrees of freedom
## Multiple R-squared: 0.5096, Adjusted R-squared: 0.5082
## F-statistic: 362.1 on 8 and 2787 DF, p-value: < 2.2e-16
export data for mapping
#write.csv(reg, "clean_regression.csv")
descriptive statistics table
statstab <- reg
statstab$GEOID <- NULL
st(statstab)
Summary Statistics
|
Variable
|
N
|
Mean
|
Std. Dev.
|
Min
|
Pctl. 25
|
Pctl. 75
|
Max
|
|
UERate
|
2796
|
5.7
|
5.1
|
0
|
2.4
|
7.7
|
100
|
|
BAPct
|
2796
|
22
|
15
|
0.84
|
11
|
29
|
95
|
|
MortgPct
|
2796
|
40
|
20
|
0
|
25
|
56
|
94
|
|
NoCarPct
|
2796
|
6.6
|
8.3
|
0
|
1.2
|
8.8
|
59
|
|
PubincPct
|
2796
|
1.8
|
2.5
|
0
|
0
|
2.5
|
23
|
|
AI
|
2796
|
0.29
|
0.098
|
0
|
0.22
|
0.36
|
1
|
|
LPJD
|
2796
|
0.018
|
0.023
|
0
|
0.0047
|
0.023
|
0.32
|
|
LPJR
|
2796
|
1.2
|
4
|
0
|
0.23
|
1.3
|
164
|
|
Norm_LPDI
|
2796
|
0.057
|
0.071
|
0
|
0.015
|
0.071
|
1
|
|
pp_need_relocate
|
2796
|
-0.04
|
171
|
-1885
|
-50
|
76
|
723
|
data exploration: race-specific model
# explains very little variable about unemployment rate, thus, not including in the presentation slides
race <- read.csv("uer_by_race_cleaned.csv")
new <- merge(reg, race, by = "GEOID")
new$white <- as.numeric(new$white)
## Warning: NAs introduced by coercion
model_white <- lm(white ~
BAPct
+ MortgPct
+ NoCarPct
+ PubincPct
+ AI
+ LPJD
+ LPJR
, data = new)
summary(model_white)
##
## Call:
## lm(formula = white ~ BAPct + MortgPct + NoCarPct + PubincPct +
## AI + LPJD + LPJR, data = new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.869 -3.566 -1.267 1.638 95.906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.169842 0.626902 9.842 < 2e-16 ***
## BAPct -0.033301 0.010151 -3.281 0.00105 **
## MortgPct -0.031855 0.008555 -3.723 0.00020 ***
## NoCarPct 0.012889 0.019820 0.650 0.51556
## PubincPct 0.030995 0.052740 0.588 0.55679
## AI -0.748880 1.558657 -0.480 0.63094
## LPJD 70.869455 7.900186 8.971 < 2e-16 ***
## LPJR -0.645820 0.101902 -6.338 2.72e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.774 on 2733 degrees of freedom
## (55 observations deleted due to missingness)
## Multiple R-squared: 0.0505, Adjusted R-squared: 0.04807
## F-statistic: 20.77 on 7 and 2733 DF, p-value: < 2.2e-16
new$black <- as.numeric(new$black)
## Warning: NAs introduced by coercion
model_black <- lm(black ~
BAPct
+ MortgPct
+ NoCarPct
+ PubincPct
+ AI
+ LPJD
+ LPJR
, data = new)
summary(model_black)
##
## Call:
## lm(formula = black ~ BAPct + MortgPct + NoCarPct + PubincPct +
## AI + LPJD + LPJR, data = new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.485 -5.995 -3.399 2.851 94.874
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.01599 1.02048 6.875 7.70e-12 ***
## BAPct -0.02572 0.01641 -1.567 0.11722
## MortgPct -0.03527 0.01362 -2.590 0.00964 **
## NoCarPct 0.14061 0.03045 4.617 4.08e-06 ***
## PubincPct 0.01236 0.08425 0.147 0.88339
## AI 2.45775 2.55799 0.961 0.33674
## LPJD 121.41579 12.43001 9.768 < 2e-16 ***
## LPJR -1.28113 0.16145 -7.935 3.09e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.66 on 2635 degrees of freedom
## (153 observations deleted due to missingness)
## Multiple R-squared: 0.06356, Adjusted R-squared: 0.06108
## F-statistic: 25.55 on 7 and 2635 DF, p-value: < 2.2e-16
new$asian <- as.numeric(new$asian)
## Warning: NAs introduced by coercion
model_asian <- lm(asian ~
BAPct
+ MortgPct
+ NoCarPct
+ PubincPct
+ AI
+ LPJD
+ LPJR
, data = new)
summary(model_asian)
##
## Call:
## lm(formula = asian ~ BAPct + MortgPct + NoCarPct + PubincPct +
## AI + LPJD + LPJR, data = new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.002 -4.772 -3.841 -2.765 97.291
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.21932 1.91869 2.720 0.00659 **
## BAPct -0.02028 0.02567 -0.790 0.42957
## MortgPct -0.01274 0.02286 -0.557 0.57757
## NoCarPct 0.11990 0.06142 1.952 0.05108 .
## PubincPct 0.28177 0.15970 1.764 0.07785 .
## AI -2.75551 4.64259 -0.594 0.55291
## LPJD 4.69324 22.40626 0.209 0.83411
## LPJR -0.14085 0.28464 -0.495 0.62078
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.6 on 1698 degrees of freedom
## (1090 observations deleted due to missingness)
## Multiple R-squared: 0.00853, Adjusted R-squared: 0.004443
## F-statistic: 2.087 on 7 and 1698 DF, p-value: 0.04195
new$hispanic <- as.numeric(new$hispanic)
## Warning: NAs introduced by coercion
model_hispanic <- lm(hispanic ~
BAPct
+ MortgPct
+ NoCarPct
+ PubincPct
+ AI
+ LPJD
+ LPJR
, data = new)
summary(model_hispanic)
##
## Call:
## lm(formula = hispanic ~ BAPct + MortgPct + NoCarPct + PubincPct +
## AI + LPJD + LPJR, data = new)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.652 -4.812 -4.089 -2.805 96.455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.906648 1.334839 3.676 0.000242 ***
## BAPct 0.004005 0.020807 0.192 0.847376
## MortgPct -0.029062 0.017544 -1.657 0.097748 .
## NoCarPct -0.016177 0.041437 -0.390 0.696282
## PubincPct -0.026810 0.111134 -0.241 0.809389
## AI 2.455108 3.296496 0.745 0.456486
## LPJD 54.394997 16.011790 3.397 0.000692 ***
## LPJR -0.719729 0.205457 -3.503 0.000468 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.21 on 2450 degrees of freedom
## (338 observations deleted due to missingness)
## Multiple R-squared: 0.006991, Adjusted R-squared: 0.004154
## F-statistic: 2.464 on 7 and 2450 DF, p-value: 0.01612