Using OLS Regression to Predict Median House Values in Philadelphia

Version 3.0 | First Created Sep 26, 2024 | Updated Oct 14, 2024

Keywords: Ordinary Least Square Regression, K-fold Cross Validation, Sum of Square Errors, Root Mean Squared Error, Stepwise AIC

GitHub Repository: CPLN671-OLS-Regression

Introduction

Housing values in Philadelphia have seen a consistent increase in recent years, making it an important area of concern and research for urban planners and policymakers. This upward trend has raised issues around housing affordability, particularly for low and middle-income households. Many neighborhoods have experienced gentrification, leading to the displacement of long-term residents and changing the socio-economic landscape of the city. Examining housing values is critical for understanding these shifts and addressing concerns related to affordability, neighborhood stability, and equitable development across Philadelphia.

To better understand the potential factors influencing housing values in Philadelphia, we identified several key variables, including educational attainment, vacancy rates, the proportion of detached single-family homes, and poverty levels. These factors are closely linked to housing market trends and provide insights into the socio-economic conditions that may affect property values across the city.

Educational attainment is closely tied to economic outcomes and this includes property values. Studies suggest that areas with higher levels of education tend to have more robust local economies, higher incomes, and increased demand for housing, which pushes up home values Scott (2012). High vacancy rates are typically associated with declining neighborhoods and lower median house values. Vacant properties can be signs of economic distress, contributing to neighborhood decline through reduced upkeep, higher crime rates, and a negative perception of the area Mallach (2018). In Philadelphia, the issue of vacancy has been particularly problematic in certain neighborhoods, where a large number of abandoned properties has led to significant reductions in housing values Kromer (2002).

The proportion of single-family homes in an area can also influence housing values, depending on the local market. Single-family homes are generally more desirable in many U.S. housing markets, as they offer more space, privacy, and perceived stability compared to multi-family units Glaeser & Gyourko (2018). Research also consistently demonstrates a strong correlation between poverty levels and housing values. Areas with higher poverty rates often experience reduced demand for housing, as well as diminished investment in local infrastructure and services, which in turn can lead to lower property values Galster (2008).

In this study, we aim to explore the relationships between these socio-economic factors and median house values in Philadelphia using Ordinary Least Squares (OLS) regression. By analyzing these relationships, we hope to identify the key predictors of housing values in the city and provide insights that can inform policy decisions, urban planning strategies, and community development initiatives.

Methods

Data Cleaning

The original dataset used in this study consisted of 1816 census block groups in Philadelphia according to the 2000 Census, but to ensure the dataset was suitable for regression analysis, several data cleaning procedures were performed. First, every block group in the dataset with a population of less than 40 was removed. These block groups represent sparsely populated or unusually housed areas that may introduce noise into the analysis. Secondly, block groups without any housing units were also excluded due to their lack of significance in a study centered on housing attributes. Thirdly, block groups with a median housing value of less than $10,000 were not included. These block groups do not accurately reflect the reality of the typical housing market, and the unusually low home values in these block groups may bias the results of the study. Finally, a unique combination of values - median property values over $800,000 and typical household incomes less than $8,000 - led to the exclusion of one block group in North Philadelphia. As a clear outlier, it was excluded to avoid biasing the results. The final dataset includes 1,720 observations, each corresponding to a different Philadelphia block group, and has the following attributes:

• POLY_ID: Census block group ID
• MEDHVAL: Median value of all owner occupied housing units
• PCTBACHMOR: Proportion of residents in Block Group with at least a Bachelor’s degree
• PCTVACANT: Proportion of housing units that are vacant
• PCTSINGLES: Percent of housing units that are detached single family houses
• NBELPOV100: Number of households with incomes below 100% poverty level (i.e., number of households living in poverty)
• MEDHHINC: Median household income

Exploratory Data Analysis

We begin by examining the summary statistics and distributions of the dependent and independent variables. This will provide a better understanding of the data and help identify any potential outliers or unusual patterns that may need to be addressed before proceeding with regression analysis. In our case, the median home value (MEDHVAL) is the dependent variable. The independent variables are the proportion of the population with a bachelor’s degree or higher (PCBACHMORE), the number of households in poverty (NBELPOV100), the proportion of vacant homes (PCTVACANT), and the proportion of single family homes (PCTSINGLES). We computed the mean and standard deviation for each variable to understand the central tendency and variability of the data. If results show that the variable is non-normally distributed, we apply a log transformation to normalize the data.

This is followed by investigating the correlations between the predictor variables to avoid multicollinearity. Correlation is a statistical measure that describes the strength and direction of a relationship between variables. The correlation coefficient \(r\) quantifies the degree to which two variables are linearly related. The formula for the sample correlation coefficient \(r\) is as follows:

\[ r = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n} (x_i - \bar{x})^2 \sum_{i=1}^{n} (y_i - \bar{y})^2}} \] In a more concise way, this above formula is also equivalent to the following:

\[ r = \frac{1}{n-1} \sum_{i=1}^{n} \left( \frac{x_i - \bar{x}}{S_x} \right) \left( \frac{y_i - \bar{y}}{S_y} \right) \] where \(x_i\) and \(y_i\) are individual data points, \(\bar{x}\) and \(\bar{y}\) are the means of the variables \(x\) and \(y\), \(S_x\) and \(S_y\) are the standard deviations of the variables \(x\) and \(y\) respectively.

For each observation \(i\), the formula calculates the difference between the data point \(x_i\), \(y_i\) and the mean of the variables and then divides them by their standard deviation. Then, these standardized values are multiplied together to assess the extent to which the two variables deviate from their respective means in the same direction. The result of this multiplication is summed across all observations, which aggregates the individual contributions to the correlation. Finally, dividing this sum by \(n-1\) where \(n\) is the number of observations yields the sample correlation coefficient \(r\).

The correlation coefficient \(r\) can range from -1 to +1. A value of +1 indicates a perfect positive linear relationship, where increases in one variable correspond with increases in the other, while a value of -1 indicates a perfect negative linear relationship. A value of 0 suggests no linear relationship between the variables, meaning changes in one do not predict changes in the other. Understanding the correlations between the predictors is critical to identify whether any variables are highly correlated, which would suggest multicollinearity. Multicollinearity can distort the regression results by inflating the standard errors of the coefficients, which lead to less reliable interpretations.

Multiple Regression Analysis

After exploring the data, we conduct a multiple regression analysis to examine the relationship between the dependent variable (median house value) and the independent variables (educational attainment, poverty levels, vacancy rates, and single-family housing). Regression analysis is a statistical method used to examine the relationship between a dependent variable and one or more independent variables. It allows researchers to understand how the dependent variable changes as the independent variables vary and make inferences about the relationships among the variables. The model estimates coefficients for each independent variable, which quantify the expected change in the dependent variable for a one-unit increase in the independent variable, holding all other variables constant. For our problem, the equation for the multiple regression model is as follows:

\[ \text{LNMEDHVAL} = \beta_0 + \beta_1 \text{PCTVACANT} + \beta_2 \text{PCTSINGLES} + \beta_3 \text{PCTBACHMOR} + \beta_4 \text{LN(LNNBELPOV100)} + \epsilon \]

where \(\text{LNMEDHVAL}\) is the log-transformed median house value, \(\text{PCTVACANT}\) is the proportion of vacant housing units, \(\text{PCTSINGLES}\) is the proportion of single-family detached homes, \(\text{PCTBACHMOR}\) is the percentage of residents with a bachelor’s degree or higher, \(\text{LN(LNNBELPOV100)}\) is the log-transformed number of households living below the poverty line.

\(\beta_0\) is the intercept, \(\beta_1, \beta_2, \beta_3, \beta_4\) are the coefficients for each independent variable, and \(\epsilon\) is the residual error term. The coefficients \(\beta_1, \beta_2, \beta_3, \beta_4\) represent the expected change in \(\text{LNMEDHVAL}\) for a one-unit increase in the corresponding independent variable, holding all other variables constant. The error term \(\epsilon\) accounts for the variability in \(\text{LNMEDHVAL}\) that is not explained by the independent variables.

Regression Assumptions

There are several assumptions associated with regression analysis that needs to be checked before we further proceed. First, linearity asserts that there should be a linear relationship between the dependent variables and independent variable. To verify this assumption, we made scatterplots to visualize the relationship between the relationship between \(y\) and each of the predictor.

Second the assumption of independence of observations requires that the observations within the dataset be independent of one another, without any spatial, temporal, or other dependencies. In spatial contexts, we made choropleth maps to assess the presence of spatial autocorrelation in the residuals or dependent variable.

Third, the normality of residuals assumption states that the residuals should follow a normal distribution. While this assumption is less critical in large samples due to the Central Limit Theorem, it remains desirable for inference purposes. To evaluate this assumption, we made a histogram of the residuals - a bell-shaped curve would suggest that the residuals are normally distributed.

Another important assumption is homoscedasticity, which states that the variance of the residuals should remain constant across all levels of the independent variables. We assessed homoscedasticity by examining scatterplot of standardized residuals against predicted values. Ideally, these residuals should be evenly distributed around zero; any patterns may indicate the presence of heteroscedasticity.

The assumption of no multicollinearity asserts that the predictor variables should not be highly correlated with each other, as this can inflate the variance of the coefficient estimates and make them unstable. To check for multicollinearity, we created a correlation matrixfor the predictors to look for any correlations greater than 0.8 (or less than -0.8).

Finally, it is important to consider the ratio of observations to predictors in the model. A general guideline suggests having at least 10 to 15 observations for each predictor to ensure robust estimates. Since we have over 1700 observations in our dataset, this assumption is met.

Parameter Estimations

After verifying these assumptions, we proceed with the regression analysis. There are several parameters that we need to estimate here: \(\beta_0\), which is the intercept, \(\beta_1, \dots, \beta_k\), which are coefficients of the independent variables, as well as \(\sigma^2\), the variance of the error terms, which represents the variability in the dependent variable \(Y\) that is not explained by the predictors \(X_1, \dots, X_k\).

The least squares method estimates the coefficients by minimizing the sum of squared residuals. Given \(n\) observations on \(y\), and \(k\) predictors \(X_1, \dots, X_k\), the estimates \(\beta_1, \dots, \beta_k\) are chosen simultaneously to minimize the expression for the Error Sum of Squares (SSE), given by:

\[ \text{SSE} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 = \sum_{i=1}^{n} (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_{i1} - \dots - \hat{\beta}_k x_{ik})^2 \]

where \(y_i\) represents the actual values, \(\hat{\beta}_0\) is the estimated intercept, \(\beta_1\) is the coefficient for the predictor \(x_{1i}\), and \(x_{1i}\) is the predictor value for the \(i\)-th observation.

With the Error Sum of Squares (SSE), the equation of the estimated variance of the error term \(\sigma^2\) is given by:

\[ \sigma^2 = \frac{\text{SSE}}{n - (k+1)} \]

where \(n\) is the number of observations and \(k\) is the number of predictors in the model.

We evaluated the model’s goodness of fit using the coefficient of determination \(R^2\), which quantifies the proportion of variance in the dependent variable explained by the independent variables. The adjusted \(R^2\), which adjusts the \(R^2\) value based on the number of predictors in the model to provide a more accurate measure of model fit for multiple regression. In simple regression:

\[ SST = \sum_{i=1}^{n} (y_i - \bar{y})^2 \] where \(y_i\) is the observed value, \(\bar{y}\) is the mean of the observed values, and \(n\) is the number of observations. Then, \(R^2\) can be simply obtained by:

\[ R^2 = 1 - \frac{SSE}{SST} \] \(R^2\) is then typically adjusted as follows based on \(n\), the number of observations and \(k\), the number of predictors in the model.

\[ R^2_{\text{adj}} = 1 - \frac{(1 - R^2)(n - 1)}{n - k - 1} \]

Hypothesis Testing

In this analysis, we test several hypotheses to evaluate the significance of the model and its predictors. The overall significance of the regression model is assessed using the F-ratio. It compares the variance explained by the regression model to the variance that is not explained by the model, assessing whether the independent variables as a whole have a statistically significant effect on the dependent variable. A higher F-ratio indicates that the regression model explains a significant amount of variance in the dependent variable compared to the residual variance.

The null hypothesis (\(H_0\)) and alternative hypothesis (\(H_a\)) for the F-ratio are stated as follows:

• Null Hypothesis (\(H_0\)): The regression model does not explain a significant portion of the variance in the dependent variable (median house value). In our case, this means that the coefficients for all the predictors are equal to zero and none of them explain any variance in housing price. This can be formally stated as: \[ H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0 \]

• Alternative Hypothesis (\(H_a\)): At least one of the predictor variables has a significant effect on the dependent variable. In our case, this means that at least one of the coefficients for the predictors is not equal to zero, indicating that the predictor has a significant impact on housing prices. This can be stated as: \[ H_a: \text{At least one } \beta_i \neq 0 \]

For each individual predictor \(\beta_i\), we also test the following hypotheses:

• Null Hypothesis (\(H_{0i}\)): The coefficient for the predictor \(i\) is equal to zero, indicating that it does not have a significant effect on the dependent variable (median house value). This can be formally stated as: \[ H_{0i}: \beta_i = 0 \]

• Alternative Hypothesis (\(H_{ai}\)): The coefficient for the predictor \(i\) is not equal to zero, indicating that it has a significant effect on the dependent variable. This can be stated as: \[ H_{ai}: \beta_i \neq 0 \]

Each predictor’s significance is typically assessed using a t-test, which tests whether the estimated coefficient differs significantly from zero.

Additional Analyses

In addition, we conducted stepwise regression to identify the most important predictors of median house values in Philadelphia. Stepwise regression is a variable selection technique that automatically selects the best subset of predictors for the model. It involves adding or removing predictors based on their statistical significance, with the goal of creating a parsimonious model that explains the most variance in the dependent variable.

We acknowledge several limitations of stepwise regression. First, the final model selected by stepwise regression is not guaranteed to be optimal in any specific sense. There may be other models that are just as good, or even better, than the one identified through this process. Additionally, the method yields only a single final model, even though multiple models could provide equally valid or informative predictions. Moreover, stepwise regression does not take into account the researcher’s domain knowledge regarding the predictors. Important variables may be excluded unless manually forced into the model, which can lead to an incomplete or inaccurate understanding of the relationships between the variables. This method also runs the risk of Type I and Type II errors—meaning it may include unimportant variables or exclude important ones.

We also performed K-fold cross-validation to assess the predictive performance of the model. Cross-validation is a resampling technique used to evaluate the performance of a predictive model. In K-fold cross-validation, the dataset is divided into K subsets or folds. The model is trained on K-1 folds and tested on the remaining fold. This process is repeated K times, with each fold serving as the test set exactly once. The performance of the model is then averaged across the K folds to provide an overall estimate of predictive accuracy. In our analysis, we used K = 5, which is a common choice for cross-validation. The root mean squared error (RMSE) is used to compare the performance of different models. The RMSE is a measure of the differences between predicted and observed values, with lower values indicating better predictive accuracy. In previous analysis, we have calculated the sum of squared errors (SSE) as follows:

\[ \text{SSE} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \]

Then the mean square errors (MSE) is simply:

\[ \text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 \]

And the root mean square error (RMSE) is the square root of the MSE:

\[ \text{RMSE} = \sqrt{\text{MSE}} \]

Software and Pacakges

For our data analysis, we are utilizing the R programming language. The following libraries have been employed to facilitate our analysis:

• tidyverse: A collection of R packages designed for data science that includes tools for data manipulation, visualization, and analysis.
• here: A package that simplifies file path management, making it easier to work with project directories.
• kableExtra: A package for creating and customizing tables in R Markdown documents.
• ggplot2: A widely-used package for creating static graphics and visualizations based on the grammar of graphics.
• ggcorrplot: A package specifically designed for visualizing correlation matrices using ggplot2.
• sf: A package that provides support for spatial data in R, enabling the handling and analysis of geometric objects.
• patchwork: A package that allows for the combination of multiple ggplot2 plots into a single cohesive layout.
• MASS: A package that provides functions and datasets for various statistical methods, including linear and generalized linear models.
• caret: A package that streamlines the process of creating predictive models, including functions for data splitting, pre-processing, feature selection, and model tuning.

Results

Exploratory Results

We began by summarizing the statistical data for median house values and various socio-economic predictors. The median house value has an average of $66,287.73, with a standard deviation of $60,006.08. This is significantly lower than both the county and state averages and shows considerable variation in house values. For educational attainment, 16.08% of individuals hold a Bachelor’s degree or higher, much lower than the state average of 35.3%. The standard deviation of 17.77% suggests notable variability in this measure. The average number of households living in poverty is 189.77, with a wide range reflected by a standard deviation of 164.32. Vacant houses make up an average of 11.29%, slightly higher than the state average of 9.4%, with a standard deviation of 9.63%. Finally, single house units account for an average of 9.23%, with a standard deviation of 13.25%.

# examine the mean and standard deviation of dependent and independent variables. 
dist_results <- data.frame(Variable = character(), Mean = numeric(), SD = numeric(), stringsAsFactors = FALSE)
variables <- c("MEDHVAL", "PCTBACHMOR", "NBELPOV100", "PCTVACANT", "PCTSINGLES")
relabelled_variables <- c(
  "Median House Value",                                    # MEDHVAL
  "% of Individuals with Bachelor’s Degrees or Higher",    # PCTBACHMOR
  "# Households Living in Poverty",                        # NBELPOV100
  "% of Vacant Houses",                                    # PCTVACANT
  "% of Single House Units"                                # PCTSINGLES
)

for (i in seq_along(variables)) {
  
  mean_val <- round(mean(data[[variables[i]]], na.rm = TRUE), 3) 
  sd_val <- round(sd(data[[variables[i]]], na.rm = TRUE), 3) 
  
  # Store the relabeled variable names
  dist_results <- rbind(dist_results, data.frame(Variable = relabelled_variables[i], Mean = mean_val, SD = sd_val))
}


dist_results <- rbind(
  data.frame(Variable = "Dependent Variable", Mean = "", SD = ""),
  data.frame(Variable = "Median House Value", Mean = dist_results$Mean[1], SD = dist_results$SD[1]),
  data.frame(Variable = "Predictors", Mean = "", SD = ""),
  dist_results[-1, ] # Remove the first row because it has already been used above
)

dist_results %>%  
  kable(row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% 
  row_spec(0, bold = TRUE) %>%  # Bold the header row
  row_spec(1, bold = TRUE) %>%  # Bold the 'Dependent Variable' row
  row_spec(3, bold = TRUE)

Variable	Mean	SD
Dependent Variable
Median House Value	66287.733	60006.076
Predictors
% of Individuals with Bachelor’s Degrees or Higher	16.081	17.77
# Households Living in Poverty	189.771	164.318
% of Vacant Houses	11.289	9.628
% of Single House Units	9.226	13.249

According to the histograms, none of the distributions of the dependent and independent variables were normal. To address this, we applied a logarithmic transformation to the variables.

# histogram

data %>%
  pivot_longer(cols = c("MEDHVAL", "PCTBACHMOR", "NBELPOV100", "PCTVACANT", "PCTSINGLES"),
               names_to = "Variable",
               values_to = "Value") %>% 
  ggplot(aes(x = Value)) +
  geom_histogram(aes(y = ..count..), fill = "#283d3b", alpha = 0.7) +  
  facet_wrap(~Variable, scales = "free", ncol = 3, labeller = as_labeller(c(
    "MEDHVAL" = "Median House Value",
    "PCTBACHMOR" = "% with Bachelor’s Degrees or Higher",
    "NBELPOV100" = "# Households Living in Poverty",
    "PCTVACANT" = "% of Vacant Houses",
    "PCTSINGLES" = "% of Single House Units"
  ))) +  
  labs(x = "Value", y = "Count", title = "Histograms of Dependent and Predictor Variables") +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=6),
        axis.text.y=element_text(size=6), 
        axis.title=element_text(size=8))

This transformation has helped to normalize some of the distributions, particularly for variables like median house value, number of households in poverty, and percentage of vacant houses. These variables now display more symmetric, bell-shaped distributions. However, for variables like the percentage of individuals with a bachelor’s degree and percentage of single house units, the transformation was less effective, and they still deviate from a normal distribution.

Moving forward, we will use the log-transformed results for our analysis. The remaining regression assumptions will be examined in detail in a separate section titled Regression Assumption Checks.

# histograms of the logged variables
reg_data %>%
  pivot_longer(cols = c("LNMEDHVAL", "LNPCTBACHMOR", "LNNBELPOV100", "LNPCTVACANT", "LNPCTSINGLES"),
               names_to = "Variable",
               values_to = "Value") %>% 
  ggplot(aes(x = Value)) +
  geom_histogram(aes(y = ..count..), fill = "#283d3b", alpha = 0.7) +  
  facet_wrap(~Variable, scales = "free", ncol = 3, labeller = as_labeller(c(
    "LNMEDHVAL" = "Median House Value",
    "LNPCTBACHMOR" = "% with Bachelor’s Degrees or Higher",
    "LNNBELPOV100" = "# Households Living in Poverty",
    "LNPCTVACANT" = "% of Vacant Houses",
    "LNPCTSINGLES" = "% of Single House Units"
  ))) +  
  labs(x = "Value", y = "Count", title = "Histograms with Logged Transform Variables") +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=6),
        axis.text.y=element_text(size=6), 
        axis.title=element_text(size=8))

The choropleth maps of the dependent variable and the predictor variables reveal several notable patterns. In particular, the map of median house values shows higher values in the northwest Philadelphia and center city while lower values in the south and west of Philadelphia. This pattern is consistent with the distribution of single house units and the percentage of individuals with a bachelor’s degree or higher, which also exhibit higher concentrations in the northwest and center city. In contrast, the maps of households living in poverty and vacant housing units show higher concentrations in the north, west, and parts of the south of Philadelphia. These patterns suggest an inverse relationship between poverty and vacant housing with educational attainment and single-family homeownership. It also appears that there might some degree of multicolinearity between percentage of population with a bachelor’s degree or higher and percentage of single family units, given that they share similar spatial patterns.

# choropleth maps of predictor and dependent variables
ggplot(philly) +
  geom_sf(aes(fill = LNMEDHVAL), color = "transparent") +
  scale_fill_gradientn(colors = c("#FAF9F6", "#c44536"), 
                       name = "LNMEDHVAL", 
                       na.value = "transparent") +  # Handle NAs
  theme(legend.text = element_text(size = 9),
        legend.title = element_text(size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill = NA, size = 0.8)) +
  labs(title = "Log Transformed Median House Value")

map_pctvacant <- ggplot(philly) +
  geom_sf(aes(fill = PCTVACANT), color = "transparent") +
  scale_fill_gradientn(colors = c("#FAF9F6", "#c44536"), 
                       name = "PCTVACANT", 
                       na.value = "transparent") +  # Handle NAs
  theme(legend.text = element_text(size = 9),
        legend.title = element_text(size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill = NA, size = 0.8)) +
  labs(title = "% of Vacant Houses")

map_pctsingles <- ggplot(philly) +
  geom_sf(aes(fill = PCTSINGLES), color = "transparent") +
  scale_fill_gradientn(colors = c("#FAF9F6", "#c44536"), 
                       name = "PCTSINGLES", 
                       na.value = "transparent") +  # Handle NAs
  theme(legend.text = element_text(size = 9),
        legend.title = element_text(size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill = NA, size = 0.8)) +
  labs(title = "% of Single House Units")

map_pcbachmore <- ggplot(philly) +
  geom_sf(aes(fill = PCTBACHMOR), color = "transparent") +
  scale_fill_gradientn(colors = c("#FAF9F6", "#c44536"), 
                       name = "PCTBACHMOR", 
                       na.value = "transparent") +  # Handle NAs
  theme(legend.text = element_text(size = 9),
        legend.title = element_text(size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill = NA, size = 0.8)) +
  labs(title = "% Bachelor's Degree or Higher")


map_lnbelpov100 <- ggplot(philly) +
  geom_sf(aes(fill = LNNBELPOV), color = "transparent") +
  scale_fill_gradientn(colors = c("#FAF9F6", "#c44536"), 
                       name = "LNNBELPOV", 
                       na.value = "transparent") +  # Handle NAs
  theme(legend.text = element_text(size = 9),
        legend.title = element_text(size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill = NA, size = 0.8)) +
  labs(title = "Logged Transformed Poverty")

map_pctvacant + map_pctsingles + map_pcbachmore + map_lnbelpov100 + 
  plot_layout(ncol = 2)

The correlation matrix supports our observations from the map. There’s a positive linear relationship between the percentage of people with at least a Bachelor’s Degree or higher and the percentage of single-family units, with a correlation coefficient of 0.2. The percentage of households living below the poverty line also positively correlated with the percentage of vacant houses, with a correlation coefficient of 0.25. These correlations suggest that the predictors are related to each other, but since none of these correlations exceed 0.8, we can proceed with the regression analysis without concerns about multicollinearity.

custom_labels <- c(
  "% of Individuals with Bachelor’s Degrees or Higher" = "PCTBACHMOR",
  "% of Vacant Houses" = "PCTVACANT",
  "% of Single House Units" = "PCTSINGLES",
  "# Households Living in Poverty" = "LNNBELPOV100"
)

cor_matrix <- cor(reg_data %>% dplyr::select(PCTBACHMOR, PCTVACANT, PCTSINGLES, LNNBELPOV100))
rownames(cor_matrix) <- names(custom_labels)
colnames(cor_matrix) <- names(custom_labels)

ggcorrplot(cor_matrix, 
           type = "lower", 
           lab = TRUE, 
           lab_size = 3, 
           colors = c("#283d3b", "white", "#c44536")) +
  labs(title = "Correlation Matrix for all Predictor Variables") +
  theme(plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x = element_text(size = 7),
        axis.text.y = element_text(size = 7), 
        axis.title = element_text(size = 8))

Regression Results

Parameter Estimates

Four independent variables, the proportion of vacant housing units (PCTVACANT), the proportion of single-family detached homes (PCTSINGLES), the proportion of residents with at least a Bachelor’s degree (PCTBACHMOR), and the number of log-transformed households living below the poverty line (LNNBELPOV100), were regressed to predict the median log-transformed home values (LNMEDHVAL) for the Philadelphia census block group.

The regression output tells us that the median home values are highly significant and are associated with the proportion of single-family detached homes(PCTSINGLES), the proportion of residents with at least a Bachelor’s degree (PCTBACHMOR), the proportion of vacant housing units (PCTVACANT), and the number of log-transformed households living below the poverty line (LNNBELPOV100), with \(p < 0.001\) for all predictors.

As the proportions of vacant housing units (PCTVACANT) goes up by one unit (i.e.percentage point), the expected change in median home value is:

\[(e^{\beta} - 1) \times 100\% = (e^{-0.0191569} - 1) \times 100\% \approx -1.897\% \]

That is, as percentage of population living in poverty goes up by 1 unit (1%), median home value goes down by approximately 1.897%, holding other predictors constant. The p-value of of less than 0.0001 for PCTVACANT tells us that if there is actually no relationship between PCTVACANT and the dependent variable (i.e., if the null hypothesis that \(\beta_1 = 0\) is actually true), then the probability of getting a \(\beta_1\) coefficient estimate of -0.0191569 is less than 0.0001.

Similarly, as the proportion of single-family detached homes (PCTSINGLES) goes up by one unit (i.e.percentage point), the expected change in median home value is:

\[(e^{\beta} - 1) \times 100\% = (e^{0.0029769} - 1) \times 100\% \approx 0.298\% \]

That is, as percentage of population living in poverty goes up by 1 unit (1%), median home value goes up by approximately 0.298%, holding other predictors constant. The p-value of of less than 0.0001 for PCTSINGLES tells us that if there is actually no relationship between PCTSINGLES and the dependent variable (i.e., if the null hypothesis that \(\beta_2 = 0\) is actually true), then the probability of getting a \(\beta_2\) coefficient estimate of 0.0029769 is less than 0.0001.

Moreover, as the proportion of people with a Bachelor’s degree or higher (PCTBACHMOR) goes up by one unit (i.e.percentage point), the expected change in median home value is:

\[(e^{\beta} - 1) \times 100\% = (e^{0.0209098} - 1) \times 100\% \approx 2.11\% \]

That is, as percentage of population with a Bachelor’s degree or higher goes up by 1 unit (1%), median home value goes up by approximately 2.1%, holding other predictors constant. The p-value of of less than 0.0001 for PCTBACHMOR tells us that if there is actually no relationship between PCTBACHMOR and the dependent variable (i.e., if the null hypothesis that \(\beta_3 = 0\) is actually true), then the probability of getting a \(\beta_3\) coefficient estimate of 0.0209098 is less than 0.0001.

And finally, as the number of log-transformed households living below the poverty line (LNNBELPOV100) changes by 1%, the expected value of median home value changes by

\[(1.01^{\beta} - 1) \times 100\% = (1.01^{-0.0789054} - 1) \times 100\% \approx -0.078\%\]

That is, as the number of households living below the poverty line goes up 1%, median home value changes by -0.078%, holding % in other predictors constant. The p-value of of less than 0.0001 for LNNBELPOV100 tells us that if there is actually no relationship between LNNBELPOV100 and the dependent variable (i.e., if the null hypothesis that \(\beta_4 = 0\) is actually true), then the probability of getting a \(\beta_4\) coefficient estimate of -0.0789054 is less than 0.0001.

All of these low probabilities indicate that we can safely reject \(H_0: \beta_1 = 0\) for \(H_a: \beta_1 \neq 0\), \(H_0: \beta_2 = 0\) for \(H_a: \beta_2 \neq 0\), \(H_0: \beta_3 = 0\) for \(H_a: \beta_3 \neq 0\), and \(H_0: \beta_4 = 0\) for \(H_a: \beta_4 \neq 0\) (at most reasonable levels of \(\alpha = P(\text{Type I error})\)).

Model Fit

The residuals from the model show a fairly reasonable distribution of errors, with a minimum residual value of -2.26, a first quartile of -0.20, a median of 0.04, a third quartile of 0.22, and a maximum value of 2.24. This distribution suggests that the model captures changes in housing values reasonably well, although some of the extreme residuals suggest that outliers may exist.

The model had a multiple R-squared value of 0.6623 and an adjusted R-squared value of 0.6615, indicating that approximately 66% of the variance in LNMEDHVAL is explained by the four predictors. The F-statistic for this model was 840.9 (p < 0.001), further confirming the statistical significance of the model as a whole.

fit <- lm(LNMEDHVAL ~ PCTVACANT + PCTSINGLES + PCTBACHMOR + LNNBELPOV100, data=reg_data)
summary(fit)

## 
## Call:
## lm(formula = LNMEDHVAL ~ PCTVACANT + PCTSINGLES + PCTBACHMOR + 
##     LNNBELPOV100, data = reg_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.25825 -0.20391  0.03822  0.21744  2.24347 
## 
## Coefficients:
##                Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  11.1137661  0.0465330 238.836 < 0.0000000000000002 ***
## PCTVACANT    -0.0191569  0.0009779 -19.590 < 0.0000000000000002 ***
## PCTSINGLES    0.0029769  0.0007032   4.234            0.0000242 ***
## PCTBACHMOR    0.0209098  0.0005432  38.494 < 0.0000000000000002 ***
## LNNBELPOV100 -0.0789054  0.0084569  -9.330 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3665 on 1715 degrees of freedom
## Multiple R-squared:  0.6623, Adjusted R-squared:  0.6615 
## F-statistic: 840.9 on 4 and 1715 DF,  p-value: < 0.00000000000000022

The analysis of variance (ANOVA) table provides additional insight into the significance of each predictor in explaining the variation in LNMEDHVAL. For example, the sum of squares for PCTVACANT is 180.392. This value represents the variation in LNMEDHVAL that can be attributed to changes in the percentage of vacant housing units. A higher sum of squares indicates a stronger relationship between the predictor and the response variable. The sum of squares for PCTSINGLES is 24.543, showing that the percentage of single-parent households also contributes to explaining the variation in LNMEDHVAL, though to a lesser extent than PCTVACANT. The sum of squares for PCTBACHMOR is 235.118, indicating that the percentage of residents with a bachelor’s degree or higher explains a significant portion of the variance in LNMEDHVAL, more than the previous two predictors. This reinforces the importance of educational attainment in driving housing values. Lastly, the sum of squares for LNNBELPOV100 is 11.692, suggesting that while this predictor does contribute to explaining LNMEDHVAL, its impact is comparatively smaller than the others.

The total sum of squares, which represents the unexplained variation in the model, is 230.34, with a mean square error of 0.134.Overall, the ANOVA results, combined with the OLS regression results, confirmed that all four predictors-PCTVACANT, PCTSINGLES, PCTBACHMOR, and LNNBELPOV100-contribute significantly to explaining the changes in LNMEDHVAL were all significant contributors. These findings provide strong evidence of the robustness of the model and the relevance of neighbourhood characteristics in predicting housing values.

anova_table <- anova(fit)
anova_table

## Analysis of Variance Table
## 
## Response: LNMEDHVAL
##                Df  Sum Sq Mean Sq  F value                Pr(>F)    
## PCTVACANT       1 180.392 180.392 1343.087 < 0.00000000000000022 ***
## PCTSINGLES      1  24.543  24.543  182.734 < 0.00000000000000022 ***
## PCTBACHMOR      1 235.118 235.118 1750.551 < 0.00000000000000022 ***
## LNNBELPOV100    1  11.692  11.692   87.054 < 0.00000000000000022 ***
## Residuals    1715 230.344   0.134                                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Regression Assumption Check

In this section, we conducted a variety of analyses to check if the assumptions of the linear regression model were met. Earlier in this report, we have already checked for variable distribution (see the stats table above) and multicolinearity (see the corr diagram above). Thus, we began here by examining the linearity assumption of the model. We created scatter plots of the dependent variable (LNMEDHVAL) against each of the predictors (PCTVACANT, PCTSINGLES, PCTBACHMOR, and LNNBELPOV100).

The scatter plot for PCTVACANT shows a general negative relationship between the proportion of vacant housing units and LNMEDHVAL. In the graph, there are clusters of block groups where high vacancy rates correspond to lower housing values, particularly on the left side of the plot where the values of LNMEDHVAL are higher. However, the plot is not perfectly linear. Instead, it reveals a more complex relationship with significant variation in housing values across different vacancy rates. In the middle and lower ranges of vacancy rates, the decline in housing values becomes less pronounced. This suggests that while vacant housing is associated with lower property values, the impact diminishes or becomes less predictable as the vacancy rate changes.

In the scatter plot for PCTSINGLES, a weak positive trend can be observed, though the points are scattered widely around the line, indicating a weak linear relationship. As the percentage of single-family homes increases, the values of LNMEDHVAL tend to increase, but this pattern is not consistent across all block groups. The scatter of points shows significant variability in housing values at various levels of single-family housing, with some block groups exhibiting high values even at lower percentages of single-family homes. This variability suggests that while there is a positive association between single-family housing and higher values, the relationship is more nuanced, and other factors may be influencing house values at the same time.

The scatter plot for PCTBACHMOR, however, shows a much clearer linear relationship. As the proportion of residents with at least a bachelor’s degree increases, the LNMEDHVAL values increase in a more consistent and linear fashion. The points on the graph form a relatively tight band along a positive trend line, indicating that higher education levels within a neighborhood are strongly associated with higher housing values. The consistent upward trend in this graph suggests that educational attainment is a key driver of housing prices, and the linear relationship implies that this variable is well-suited for inclusion in the linear regression model.

The scatter plot for LNNBELPOV100 reveals a predominantly negative relationship, but like PCTVACANT, the association is not strictly linear. In the graph, as the number of households living below the poverty line increases, the LNMEDHVAL values generally decrease. However, the points are scattered widely, particularly at lower levels of poverty, where housing values vary significantly. This indicates that while poverty levels have a negative impact on housing values, the effect is not uniform across all block groups. There is substantial variation in the housing values even in neighborhoods with similar poverty rates, suggesting that other factors may be at play. The non-linear pattern observed in this graph highlights the complexity of poverty’s influence on the housing market.

In summary, the scatterplot shows that while PCTBACHMOR and LNMEDHVAL exhibit a strong and clear linear relationship, the other predictors - in particular PCTVACANT and LNNBELPOV100 - exhibit more complex non-linear patterns. The scatterplots provide important insights into the relationship between the independent variables and housing values, suggesting that while the model is assumed to be linear, it may not fully capture the underlying dynamics. For variables such as PCTVACANT and LNNBELPOV100, the nonlinearity observed in the plots suggests that more complex nonlinear models may provide a better fit and more accurate predictions.

# scatterplot

reg_data %>%
  pivot_longer(cols = c("PCTBACHMOR", "LNNBELPOV100", "PCTVACANT", "PCTSINGLES"),
               names_to = "Variable",
               values_to = "Value") %>% 
ggplot(aes(x = Value, y = LNMEDHVAL)) +
  geom_point(color = "#283d3b", alpha = 0.7, size = 0.4) +
  geom_smooth(method = "lm", color = "#c44536", se = FALSE) +  # Add a linear regression line
  facet_wrap(~ Variable, scales = "free", labeller = as_labeller(c(
    "PCTBACHMOR" = "% with Bachelor’s Degrees or Higher",
    "LNNBELPOV100" = "Households Living in Poverty",
    "PCTVACANT" = "% of Vacant Houses",
    "PCTSINGLES" = "% of Single House Units"
  )))  +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=6),
        axis.text.y=element_text(size=6), 
        axis.title=element_text(size=8)) +
  labs(title = "Scatter Plots of Dependent Variable vs. Predictors", 
       x = "Predictor Value", 
       y = "Log of Median House Value")

Next, we examined the assumption of normality for residuals of the regression model. The histogram of standardized residuals provides a detailed view of the distribution of the residuals from the OLS regression model. This graph maps the frequency of residual values along the x-axis, with the residuals centered around zero. The shape of the histogram approximates a bell curve, which is characteristic of a normal distribution—a key assumption for OLS regression.

The majority of the residuals are concentrated between -2 and 2 on the x-axis, with the highest frequency occurring near zero. This concentration suggests that most of the predicted values are close to the observed values, leading to residuals that are small and distributed symmetrically around zero. This is a positive sign, as it indicates that the model does not systematically overpredict or underpredict the dependent variable (LNMEDHVAL).

At the center of the distribution, the residuals form a tall, narrow peak, representing a large number of residuals close to zero. This suggests that the model fits most observations well. As we move further away from the center, the frequency of residuals decreases symmetrically, with fewer residuals having large positive or negative values. This tapering off on both sides of the histogram suggests that extreme residuals are relatively rare, which is typical of a normal distribution.

However, there are some noticeable deviations from perfect normality at the tails of the distribution. On the left side of the graph, there are a few residuals with values smaller than -3, and on the right side, there are a few residuals larger than 3. These residuals represent outliers—cases where the model’s predictions deviate significantly from the observed values. While these extreme values do not appear to dominate the distribution, they are worth further examination, as they could indicate areas where the model’s assumptions are being violated or where influential points might be affecting the overall model fit.

In summary, the histogram of standardized residuals indicates that the residuals are approximately normally distributed, which supports the assumption of normality for the regression model. Although there are some outliers present, their overall impact on the normality of the residuals seems minimal. The general shape of the histogram suggests that the OLS regression model fits the data reasonably well, with residuals mostly falling within an acceptable range around zero.

ggplot(reg_data, aes(x = Standardized_Residuals)) +
  geom_histogram(bins = 30, fill = "#283d3b", alpha = 0.7) +
  labs(title = "Histogram of Standardized Residuals", 
       x = "Standardized Residuals", 
       y = "Frequency") +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=6),
        axis.text.y=element_text(size=6), 
        axis.title=element_text(size=8))

Following the normality check, we examined the assumption of homoscedasticity, which requires that the residuals have constant variance across all levels of the fitted values. The scatter plot of standardized residuals versus fitted values is a valuable diagnostic for testing essential regression assumptions, including homoscedasticity and the presence of outliers. The fitted values that the model predicts are plotted on the x-axis, and the standardized residuals—that is, residuals adjusted by their standard deviation—are shown on the y-axis. Because standardized residuals give an indication of the standard deviation of each residual’s divergence from the expected value, they make it simple to identify outliers and significant data points. Relatives are generally regarded as normal when they fall between -2 and +2, and anything outside of these bounds can lead to anomalous points that could have an impact on the model’s performance.

Homoscedasticity, or the need that residuals show constant variance across all levels of fitted values, is one of the fundamental presumptions of regression. The residuals in this scatter plot show no discernible trend of rising or falling variance with increasing fitted values; instead, they seem to be pretty uniformly distributed around the horizontal line at zero. The fact that there is no obvious funnel shape or systematic variation trend suggests that the assumption of homoscedasticity is reasonably satisfied. The notion of constant variance is supported by the residuals’ generally constant spread across the fitted value range, which indicates that the variance of errors is still stable.The symmetry of the residuals around the zero line is another positive aspect of the plot. This symmetry means that the model does not systematically over- or under-predict the dependent variable over the range of fitted values (LNMEDHVAL).

However, the plot does show some potential outliers, particularly those points that are outside the range of the -3 and +3 standardized residuals. These extreme residuals are located above and below the main residual clusters and may indicate that the model’s predicted values deviate significantly from the actual values. Although the number of outliers is relatively small, their presence indicates that the model may not fit all data points well.

In addition, although the residuals are generally well clustered around zero, some slight dispersion occurs as the fitted values increase, suggesting that the accuracy of the model predictions may decrease slightly at higher values of LNMEDHVAL. The higher the fitted value, the greater the residual variance, which may indicate the presence of slight heteroskedasticity, but not of sufficient severity to raise significant concerns about the validity of the model. However, this pattern could be explored further to ensure that the model performs well over the entire range of predicted values.

ggplot(reg_data, aes(x = Fitted, y = Standardized_Residuals)) +
  geom_point(color = "#283d3b", alpha = 0.6) +    # Scatter plot points
  geom_hline(yintercept = 0, linetype = "dashed", color = "#c44536") +   # Add a horizontal line at y = 0
  labs(
    title = "Scatter Plot of Standardized Residuals vs Fitted Values",
    x = "Fitted Values",
    y = "Standardized Residuals"
  ) +
  theme_light() +   
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=6),
        axis.text.y=element_text(size=6), 
        axis.title=element_text(size=8))

Another assumption we test here is independence of observations, especially spatial autocorrelations. In the map of LNMEDHVAL (log-transformed median house values), there appears to be a spatial pattern, with higher values concentrated in certain areas, such as Center City and other affluent neighborhoods, and lower values clustered in more economically distressed areas. This suggests that the housing market is not spatially random and that nearby block groups tend to have similar house values.

Similarly, the map of PCTVACANT (proportion of vacant housing units) shows clear geographic clusters, with higher vacancy rates concentrated in certain parts of the city, particularly in areas experiencing economic decline. For PCTSINGLES (proportion of single-family homes), the maps reveal spatial patterns as well, with block groups characterised by a high proportion of single-family homes, often located in residential areas on the outskirts of a suburb or city. The PCTBACHMOR (proportion of Residents with at least a Bachelor’s Degree) map shows a strong spatial pattern, with higher levels of educational attainment concentrated in the more affluent central areas of the city, while lower levels of educational attainment are concentrated in the peripheral areas, which have fewer economic resources. Unsurprisingly, The LNNBELPOV100 (log-transformed number of households living below the poverty line) map shows that in some parts of the city, particularly in economically distressed neighborhoods, there is a clear concentration of households with higher poverty rates. All of these suggests that observations (block groups) are not independent of each other and it is important to account for such dependencies in the model through geographic weighted regression techniques.

philly %>%
  left_join(reg_data %>% dplyr::select(c(POLY_ID, Standardized_Residuals)), by = "POLY_ID") %>% 
  ggplot(.) +
  geom_sf(aes(fill = Standardized_Residuals), color = "transparent") +  # Use geom_sf to map the polygons
  scale_fill_gradientn(colors = c("#FAF9F6", "#c44536"), 
                       name = "Std Residuals", 
                       na.value = "transparent") +  # Choose a color palette, invert direction if needed
  labs(title = "Choropleth Map of Standardized Residuals") +
  theme(legend.text = element_text(size = 9),
        legend.title = element_text(size = 10),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(size = 9, face = "italic"),
        plot.title = element_text(size = 12, face = "bold"),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill = NA, size = 0.8))

The choropleth map of standardized regression residuals provides a detailed geographic breakdown of how well the OLS regression model fits the data across different areas of the city. The map uses a color gradient to represent the magnitude of the residuals, with darker red shades indicating large positive residuals (where the model underestimates the actual values), and lighter shades (toward white) indicating large negative residuals (where the model overestimates the actual values).

In the northwestern part of the city, particularly along some of the more suburban and affluent areas, dark red shading indicates large positive residuals. This suggests that the model underestimates housing values in these neighborhoods. These areas are often characterized by higher-income households and larger single-family homes, suggesting that the model may not fully capture the higher market value trends in these regions, possibly due to unaccounted factors like proximity to high-demand amenities or specific neighborhood features that influence property values.

In contrast, in parts of South and Southwest Philadelphia, as well as parts of North Philadelphia, lighter shaded areas indicate large negative residuals. Here, the model tends to overestimate housing values, meaning that the predicted values are higher than the actual observed values. These neighborhoods may be facing economic challenges, such as higher vacancy rates, lower household incomes, or poverty, which the model has difficulty fully capturing through the selected predictors. The residuals in these areas suggest that housing markets in these neighborhoods are more complex and influenced by localized factors not fully represented in the regression model, such as higher levels of economic distress or neighborhood disinvestment.

Some neighborhoods near Center City display more neutral residuals, where the model’s predictions are relatively accurate, with values closer to zero. These areas appear in lighter red or pink shades, indicating that the model is performing reasonably well in these parts of the city. This may be because the central areas are more homogeneous in terms of housing values and more likely to align with the trends captured by the model’s predictors, such as education levels and vacancy rates.

The overall pattern observed in the map suggests spatial clustering, as block groups with similar residual values tend to be geographically close to each other. The model tends to underestimate housing values in more affluent suburbs and overestimate housing values in economically distressed communities. This highlights the need for further improvement to address spatial auto-correlations.

Additional Analyses

We performed stepwise regression to determine the most relevant predictors for the model. The initial model included all four predictors, with a residual sum of squares (RSS) of 230.34 and an Akaike Information Criterion (AIC) of -3448.07, where a lower AIC indicates a better balance between goodness of fit and model complexity.

During the stepwise procedure, the process examined whether removing each predictor would reduce the AIC and improve the model. However, removing any predictor resulted in a higher AIC and an increase in RSS, indicating a poorer model fit. For instance, removing PCTSINGLES increased the RSS to 232.75 and the AIC to -3432.2, while removing LNNBELPOV100 increased the RSS to 242.04 and the AIC to -3364.9. Removing PCTVACANT had a more substantial impact, increasing the RSS to 281.89 and the AIC to -3102.7. Finally, removing PCTBACHMOR caused the RSS to jump to 429.36 and the AIC to -2379.0, indicating a severe decline in model performance.

Given these results, none of the predictors were removed, as eliminating any of them worsened the fit of the model. Therefore, the final model remained unchanged from the initial one, including all four predictors. The final model had a residual deviance of 230.34 and an AIC of -3448.07, indicating that the model successfully balanced complexity with fit.

stepwise_model <-  stepAIC(fit, direction = "both")

## Start:  AIC=-3448.07
## LNMEDHVAL ~ PCTVACANT + PCTSINGLES + PCTBACHMOR + LNNBELPOV100
## 
##                Df Sum of Sq    RSS     AIC
## <none>                      230.34 -3448.1
## - PCTSINGLES    1     2.407 232.75 -3432.2
## - LNNBELPOV100  1    11.692 242.04 -3364.9
## - PCTVACANT     1    51.546 281.89 -3102.7
## - PCTBACHMOR    1   199.020 429.36 -2379.0

stepwise_model$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## LNMEDHVAL ~ PCTVACANT + PCTSINGLES + PCTBACHMOR + LNNBELPOV100
## 
## Final Model:
## LNMEDHVAL ~ PCTVACANT + PCTSINGLES + PCTBACHMOR + LNNBELPOV100
## 
## 
##   Step Df Deviance Resid. Df Resid. Dev       AIC
## 1                       1715   230.3435 -3448.073

We performed cross-validation to assess the model’s predictive accuracy and determine whether a reduced model with fewer predictors could provide similar performance. The original model, which includes all four predictors, has a lower RMSE (0.3676) compared to the reduced model with only two predictors (0.4424). This indicates that the original model provides better predictive accuracy. The increase in RMSE when removing the additional predictors (PCTSINGLES and LNNBELPOV100) suggests that excluding these variables results in a less accurate model. Overall, the original model outperforms the reduced model, demonstrating that all four predictors contribute to improved model accuracy.

lm <-  trainControl(method = "cv", number = 5)

cvlm_model <- train(LNMEDHVAL ~ PCTVACANT + PCTSINGLES + PCTBACHMOR + LNNBELPOV100, data=reg_data, method = "lm", trControl = lm)

print(cvlm_model)

## Linear Regression 
## 
## 1720 samples
##    4 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 1376, 1376, 1376, 1376, 1376 
## Resampling results:
## 
##   RMSE       Rsquared  MAE      
##   0.3664127  0.658518  0.2726616
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

cvlm_model_reduced = train(LNMEDHVAL ~ PCTVACANT + MEDHHINC, data = reg_data, method = "lm", trControl = lm)

print(cvlm_model_reduced)

## Linear Regression 
## 
## 1720 samples
##    2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 1376, 1376, 1376, 1376, 1376 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.4427044  0.5050842  0.3174302
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Discussion and Limitations

In summary, in this report, we conducted a multiple regression to examine the relationship between median house value (the dependent variable) and four independent variables: vacancy rates (PCTVACANT), percentage of single-family housing (PCTSINGLES), educational attainment (PCTBACHMOR), and poverty levels (LNNBELPOV100).

Our findings show that all four predictors were statistically significant. Specifically, vacancy rates and poverty levels had a negative impact on median house value, while single-family housing and educational attainment had a positive influence. The fact that educational attainment was highly significant and positively related to house value is not surprising, as areas with higher education levels tend to attract higher-income households. Similarly, the negative effect of vacancy rates and poverty levels aligns with expectations, as these factors typically lower property values.

The overall quality of the model was good, with an \(R^2\) value of 0.6623, indicating that approximately 66% of the variance in median house value was explained by the predictors. Additionally, the F-ratio test, with its extremely low p-value, confirms that the model as a whole is statistically significant and that at least one predictor has a significant impact on the outcome.

Through stepwise regression results, all four predictors were retained in the final model. This indicates that each of these variables contributes meaningfully to explaining the variation in median house value, and none were deemed unnecessary based on the stepwise selection process.When examining the cross-validation results using RMSE (Root Mean Square Error), the four-predictor model generally showed a better fit compared to a simpler model with fewer predictors. This suggests that including all four predictors leads to more accurate predictions and better overall model performance.

While the four-predictor model shows better predictive accuracy, there are likely other variables that could enhance the model’s explanatory power if included. For instance, factors like crime rates, transportation accessibility, and proximity to amenities (such as parks, schools, and shopping centers) are commonly associated with housing values. Spatial autocorrelation, as suggested by the choropleth maps, indicates that house prices tend to cluster geographically—nearby properties often have similar values. This spatial dependency suggests that including neighborhood or location-based variables, or even applying spatial lag regression, could better capture the clustered nature of house prices. Housing values also depend on the characteristics of the house itself, such as age, size, storage options, and the presence of amenities like parking or green space. Additionally, local environmental factors, such as proximity to high-crime areas or noise pollution, can impact property values. Yet, many of these variables are often correlated with each other (e.g., transportation access and crime rates), so care must be taken in selecting which predictors to include to avoid multicollinearity.

In addition to this limitation of not including some potentially important predictors, we have to acknowledge there are some other limitations in terms of meeting the regression assumptions. One key assumption of OLS regression is that linear relationships exist between the dependent and independent variables. While the scatter plots showed that PCTBACHMOR had a clear linear relationship with median house value, the relationships for PCTVACANT, PCTSINGLES, and LNNBELPOV100 were more complex and not entirely linear. This suggests that the model may not fully capture the underlying dynamics of these predictors and their impact on house values.

Another key assumption is the independence of observations, which is violated due to spatial autocorrelation. This suggests that the model may not fully account for the spatial dependencies in the data, which could lead to biased parameter estimates and inaccurate predictions. To address this limitation, future analyses could explore spatial regression techniques, such as spatial lag or spatial error models, to account for spatial autocorrelation and improve model performance.

On top of that, for our current model, we are using the raw number of households living in poverty (LNNBELPOV100) as a predictor. The raw count does not account for differences in population size across neighborhoods, which could lead to misleading results. Larger neighborhoods with more households may naturally have more households in poverty, even if the poverty rate is relatively low. A more meaningful measure would be to use the percentage of households in poverty, which normalizes the data across neighborhoods, providing a better sense of poverty’s prevalence relative to the overall population.

Possibility of Ridge or LASSO Regression

Ridge and LASSO regression are regularization techniques used to address multicollinearity and overfitting in models with many predictors. Ridge regression penalizes the size of the coefficients, shrinking them towards zero but keeping all predictors in the model, which is useful when all predictors are believed to contribute to the outcome. LASSO (Least Absolute Shrinkage and Selection Operator) goes a step further by forcing some coefficients to be exactly zero, effectively performing variable selection.

In this case, since the stepwise regression suggests that all four predictors significantly contribute to the model, there’s no strong indication of multicollinearity or overfitting. Therefore, Ridge or LASSO might not be necessary here unless there is evidence of high multicollinearity or many more predictors. If the model were more complex or had a larger number of predictors, LASSO could help in selecting the most important ones, while Ridge could help shrink coefficients to prevent overfitting.

References

Galster, G. (2008). Neighbourhood dynamics and housing markets (pp. 153–171). https://doi.org/10.1002/9780470690680.ch9

Glaeser, E., & Gyourko, J. (2018). The economic implications of housing supply. Journal of Economic Perspectives, 32(1), 3–30.

Kromer, J. (2002). Vacant-property policy and practice: Baltimore and philadelphia. Brookings Institution Center on Urban; Metropolitan Policy Washington, DC.

Mallach, A. (2018). The empty house next door. Cambridge, MA: Lincoln Institute of Land Policy.

Scott, A. (2012). Triumph of the city: How our greatest invention makes us richer, smarter, greener, healthier, and happier – by edward glaeser. Economic Geography, 88, 97–100. https://doi.org/10.2307/41474071