#1 Load packages
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.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(pastecs)
##
## Attaching package: 'pastecs'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## The following object is masked from 'package:tidyr':
##
## extract
library(MASS)
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
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
#2 Load data and construct analysis variables
ahs.household.data <- read_csv("household.csv")
## Rows: 55669 Columns: 1180
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (635): CONTROL, JACPRIMARY, JACSECNDRY, JADEQUACY, JAIRRATE, JBATHEXCLU,...
## dbl (545): TOTROOMS, PERPOVLVL, OUTAGEFRQ, RENT, DINING, LAUNDY, RATINGHS, R...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ahs.household.data <- ahs.household.data %>% mutate(tot.cost.burden = ((TOTHCAMT * 12) + MAINTAMT) / HINCP)
clean.research.data <- ahs.household.data %>%
filter(is.finite(tot.cost.burden),tot.cost.burden > 0,TENURE == "'1'",MARKETVAL > 0,HINCP > 0, MAINTAMT >= 0,PERPOVLVL >= 0,RATINGHS >= 1 & RATINGHS <= 10,YRBUILT > 0 )
clean.research.data <- clean.research.data %>%
mutate(log.cost = log(tot.cost.burden))
clean.research.data <- clean.research.data %>% mutate(housing.age = 2023 - YRBUILT)
clean.research.data <- clean.research.data %>%
mutate(age.group = case_when(
housing.age < 15 ~ "Under 15 years",
housing.age >= 15 & housing.age <= 30 ~ "15-30 years",
housing.age >= 31 & housing.age <= 50 ~ "31-50 years",
housing.age >= 51 ~ "51+ years"))
clean.research.data$age.group <- factor(
clean.research.data$age.group,
levels = c("Under 15 years", "15-30 years", "31-50 years", "51+ years"))
clean.research.data <- clean.research.data %>%
mutate(
income.group = case_when(
PERPOVLVL < 200 ~ "Lower income",
PERPOVLVL >= 200 & PERPOVLVL < 400 ~ "Moderate income",
PERPOVLVL >= 400 ~ "Higher income" ))
clean.research.data$income.group <- factor(
clean.research.data$income.group,
levels = c("Lower income",
"Moderate income",
"Higher income"))
#3 Descriptive statistics
# Continuous variables used in core models
vars <- clean.research.data %>%
dplyr::select(
tot.cost.burden,
housing.age,
MARKETVAL,
PERPOVLVL,
RATINGHS)
stat.desc(vars)
## tot.cost.burden housing.age MARKETVAL PERPOVLVL
## nbr.val 2.468900e+04 2.468900e+04 2.468900e+04 2.468900e+04
## nbr.null 0.000000e+00 1.900000e+01 0.000000e+00 0.000000e+00
## nbr.na 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## min 4.049340e-03 0.000000e+00 1.000000e+03 1.000000e+00
## max 2.900400e+04 1.040000e+02 9.999998e+06 5.010000e+02
## range 2.900400e+04 1.040000e+02 9.998998e+06 5.000000e+02
## sum 7.458311e+04 1.217191e+06 1.303739e+10 9.054553e+06
## median 2.039832e-01 4.300000e+01 3.922040e+05 4.340000e+02
## mean 3.020904e+00 4.930094e+01 5.280649e+05 3.667444e+02
## SE.mean 1.223885e+00 1.720404e-01 4.004727e+03 9.788923e-01
## CI.mean.0.95 2.398888e+00 3.372095e-01 7.849505e+03 1.918688e+00
## var 3.698152e+04 7.307423e+02 3.959582e+11 2.365774e+04
## std.dev 1.923058e+02 2.703225e+01 6.292521e+05 1.538107e+02
## coef.var 6.365835e+01 5.483109e-01 1.191619e+00 4.193949e-01
## RATINGHS
## nbr.val 2.468900e+04
## nbr.null 0.000000e+00
## nbr.na 0.000000e+00
## min 1.000000e+00
## max 1.000000e+01
## range 9.000000e+00
## sum 2.160780e+05
## median 9.000000e+00
## mean 8.751995e+00
## SE.mean 8.556223e-03
## CI.mean.0.95 1.677071e-02
## var 1.807456e+00
## std.dev 1.344416e+00
## coef.var 1.536126e-01
# Group frequencies for ANOVA categories
table(clean.research.data$age.group)
##
## Under 15 years 15-30 years 31-50 years 51+ years
## 2396 3693 6672 11928
table(clean.research.data$income.group)
##
## Lower income Moderate income Higher income
## 4709 6574 13406
#4 Baseline linear model using logged cost burden
affordable.housing.model <- lm(log.cost ~ housing.age + MARKETVAL + PERPOVLVL + RATINGHS, data = clean.research.data)
summary(affordable.housing.model)
##
## Call:
## lm(formula = log.cost ~ housing.age + MARKETVAL + PERPOVLVL +
## RATINGHS, data = clean.research.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3973 -0.4714 0.0030 0.4416 9.7602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.005e-02 3.818e-02 1.573 0.116
## housing.age -1.502e-03 1.910e-04 -7.863 3.91e-15 ***
## MARKETVAL 2.551e-07 8.282e-09 30.807 < 2e-16 ***
## PERPOVLVL -4.456e-03 3.391e-05 -131.410 < 2e-16 ***
## RATINGHS 2.552e-03 3.843e-03 0.664 0.507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8001 on 24684 degrees of freedom
## Multiple R-squared: 0.4124, Adjusted R-squared: 0.4123
## F-statistic: 4330 on 4 and 24684 DF, p-value: < 2.2e-16
#5 Regression assumption checks
# a) Linearity (plot and raintest)
raintest(affordable.housing.model)
##
## Rainbow test
##
## data: affordable.housing.model
## Rain = 0.97456, df1 = 12345, df2 = 12339, p-value = 0.9238
plot(affordable.housing.model, which = 1)
# b) Independence of errors (durbin-watson)
durbinWatsonTest(affordable.housing.model)
## lag Autocorrelation D-W Statistic p-value
## 1 0.03018269 1.939538 0
## Alternative hypothesis: rho != 0
# c) Homoscedasticity (plot, bptest)
plot(affordable.housing.model,which=3)
bptest(affordable.housing.model)
##
## studentized Breusch-Pagan test
##
## data: affordable.housing.model
## BP = 963.12, df = 4, p-value < 2.2e-16
# d) Normality of residuals (QQ plot, KS test)
plot(affordable.housing.model,which=2)
ks.test(affordable.housing.model$residuals, "pnorm", mean = mean(affordable.housing.model$residuals), sd = sd(affordable.housing.model$residuals))
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: affordable.housing.model$residuals
## D = 0.04736, p-value < 2.2e-16
## alternative hypothesis: two-sided
# e) No multicolinarity (VIF, cor)
vif(affordable.housing.model)
## housing.age MARKETVAL PERPOVLVL RATINGHS
## 1.028008 1.047420 1.049293 1.029629
kitchen_sink_vars <- clean.research.data %>% dplyr::select(housing.age, MARKETVAL, PERPOVLVL, RATINGHS)
cor(kitchen_sink_vars)
## housing.age MARKETVAL PERPOVLVL RATINGHS
## housing.age 1.00000000 -0.02412773 -0.09321871 -0.14118329
## MARKETVAL -0.02412773 1.00000000 0.19605153 0.09351354
## PERPOVLVL -0.09321871 0.19605153 1.00000000 0.05683212
## RATINGHS -0.14118329 0.09351354 0.05683212 1.00000000
#6 Two-way ANOVA: age group and income group
anova.model <- aov(
tot.cost.burden ~ age.group * income.group,
data = clean.research.data
)
summary(anova.model)
## Df Sum Sq Mean Sq F value Pr(>F)
## age.group 3 77232 25744 0.697 0.554
## income.group 2 779052 389526 10.540 2.66e-05 ***
## age.group:income.group 6 189442 31574 0.854 0.528
## Residuals 24677 911953954 36956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#7 Median burden by age and income group
interaction.plot(
clean.research.data$age.group,
clean.research.data$income.group,
clean.research.data$tot.cost.burden,
fun = median,
ylab = "Median Housing Cost Burden",
xlab = "Housing Age Group",
trace.label = "Income Group")
#8 Supplemental repair mechanism models
# Does housing age predict resident housing quality ratings?
condition.model <- lm(
RATINGHS ~ housing.age,
data = clean.research.data
)
summary(condition.model)
##
## Call:
## lm(formula = RATINGHS ~ housing.age, data = clean.research.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.0631 -0.7962 0.2038 1.1335 1.6321
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.0981655 0.0176189 516.39 <2e-16 ***
## housing.age -0.0070216 0.0003134 -22.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.331 on 24687 degrees of freedom
## Multiple R-squared: 0.01993, Adjusted R-squared: 0.01989
## F-statistic: 502.1 on 1 and 24687 DF, p-value: < 2.2e-16
# Does housing age predict maintenance spending?
maintenance.model <- lm(
MAINTAMT ~ housing.age + PERPOVLVL + MARKETVAL,
data = clean.research.data
)
summary(maintenance.model)
##
## Call:
## lm(formula = MAINTAMT ~ housing.age + PERPOVLVL + MARKETVAL,
## data = clean.research.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9169 -1127 -716 120 98351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.354e+02 6.905e+01 6.305 2.93e-10 ***
## housing.age 5.498e+00 7.714e-01 7.127 1.05e-12 ***
## PERPOVLVL 9.421e-01 1.382e-01 6.816 9.57e-12 ***
## MARKETVAL 7.860e-04 3.365e-05 23.358 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3262 on 24685 degrees of freedom
## Multiple R-squared: 0.02815, Adjusted R-squared: 0.02803
## F-statistic: 238.4 on 3 and 24685 DF, p-value: < 2.2e-16
#9 Gamma GLM robustness check
glm.data <- clean.research.data %>%
filter(tot.cost.burden < quantile(tot.cost.burden,.99))
burden.glm <- glm(
tot.cost.burden ~ housing.age + MARKETVAL + PERPOVLVL + RATINGHS,
data = glm.data,
family = Gamma(link=log)
)
summary(burden.glm)
##
## Call:
## glm(formula = tot.cost.burden ~ housing.age + MARKETVAL + PERPOVLVL +
## RATINGHS, family = Gamma(link = log), data = glm.data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.498e-01 4.955e-02 5.041 4.66e-07 ***
## housing.age -1.377e-03 2.467e-04 -5.581 2.42e-08 ***
## MARKETVAL 3.059e-07 1.084e-08 28.231 < 2e-16 ***
## PERPOVLVL -4.341e-03 4.489e-05 -96.699 < 2e-16 ***
## RATINGHS -1.970e-03 4.977e-03 -0.396 0.692
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.058201)
##
## Null deviance: 26417 on 24441 degrees of freedom
## Residual deviance: 12674 on 24437 degrees of freedom
## AIC: -19875
##
## Number of Fisher Scoring iterations: 9
#10 Final visualization: logged burden by income group
ggplot(clean.research.data,
aes(x=income.group,y=log.cost)) +
geom_boxplot() +
labs(
title="Housing Cost Burden by Income Group",
x="Income Group",
y="Logged Housing Cost Burden")
#test to include in final maybe
#cuz it shows something i saw on the anova more clearly which is the slope for housing age vs burden by income groups
age.by.income <- clean.research.data %>%
group_by(housing.age, income.group) %>%
summarize(median.burden = median(tot.cost.burden, na.rm = TRUE),.groups = "drop")
ggplot(age.by.income,
aes(housing.age, median.burden, color=income.group)) +
geom_line() +
coord_cartesian(xlim=c(0,50), ylim=c(0,2))
# Examine extremely low household incomes
clean.research.data %>%
summarise(
total.households = n(),
income.under.100 = sum(HINCP < 100),
income.under.1000 = sum(HINCP < 1000),
income.under.5000 = sum(HINCP < 5000),
income.under.10000 = sum(HINCP < 10000),
low.income.high.value = sum(HINCP < 1000 & MARKETVAL > 500000),
low.income.million.value = sum(HINCP < 1000 & MARKETVAL > 1000000)
)
## # A tibble: 1 × 7
## total.households income.under.100 income.under.1000 income.under.5000
## <int> <int> <int> <int>
## 1 24689 28 99 636
## # ℹ 3 more variables: income.under.10000 <int>, low.income.high.value <int>,
## # low.income.million.value <int>
# View the 15 households with the lowest reported income
clean.research.data %>%
arrange(HINCP) %>%
head(15) %>%
dplyr::select(HINCP, MARKETVAL, TOTHCAMT, MAINTAMT,
tot.cost.burden, housing.age, PERPOVLVL)
## # A tibble: 15 × 7
## HINCP MARKETVAL TOTHCAMT MAINTAMT tot.cost.burden housing.age PERPOVLVL
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 4 623472 181 1 543. 23 1
## 2 4 241651 1441 743 4509. 63 1
## 3 4 70192 424 0 1272 53 1
## 4 4 372592 741 206 2274. 33 1
## 5 4 2150337 9668 0 29004 73 1
## 6 8 395868 484 0 726 33 1
## 7 8 364111 2953 5114 5069. 53 1
## 8 8 163910 609 0 914. 53 1
## 9 8 935418 609 1 914. 43 1
## 10 8 1211079 1315 600 2048. 43 1
## 11 10 468339 955 978 1244. 63 1
## 12 10 99715 400 810 561 73 1
## 13 10 282194 1002 0 1202. 53 1
## 14 12 633853 1176 192 1192 43 1
## 15 20 59074 1022 2 613. 104 1