STAT974 Project #1: Applying GAM to Homelessness in LA

By: Jeffrey Ellington

9:00am Mon-Wed Section

3/1/17

1. Introduction: Assessing Homelessness in LA

Public policy initiatives in Los Angeles use data on homelessness to properly direct funding to social programs. The information provides valuable insight into the counts and location of homeless people across LA. However, gathering data across over 2,000 census tracks in LA is burdensome and expensive. As a result, the city collected information from 509 of its over 2,000 census tracks to form the basis for a statistical analysis of homelessness in LA. This information can be analyzed to help allocate resources efficiently to the locations where homelessness is most prevalent.

2. Data Description

The data were generated from a random sample of 509 of the over 2,000 census tracts in Los Angeles. For our treatment, we will examine six variables:

  1. StreetTotal – The number of homeless individuals in a census tract 
  2. MedianIncome – Median household income in a census tract
  3. PropVacant – The proportion of vacant residential dwellings in a census tract
  4. PropMinority – The proportion of minority residents in a census tract
  5. PctCommercial – The percentage of a census tract zoned for commercial enterprises
  6. PctIndustrial – The percentage of a census tract zoned for industrial enterprises

StreetTotal is the response variable. It is expressed as a positive integer given it is count data. StreetTotal counts were collected in the 509 census tracts, where each individual count was performed over two consecutive days/evenings. The histogram for homelessness counts in LA can be seen below in Fig. 2.1. Note that the distribution is heavily skewed, with over 450 (92%) census tracts having counts of less than 100 homeless people observed during the sampling period. However, on the right extreme, 13 tracts had counts over 200. The highest count is 931, which is over 58x the median value of 16.

The other five variables will be examined in much greater detail once they’ve been properly cleaned. However, it’s worth stating up front that while these 5 variables seem important in an analysis of homelessness, they are far from complete. An expert on homelessness in Los Angeles could likely point to several other important variables, like crime or drug usage/access, not to mention how the data were actually collected and by whom. We also have no framework for existing policy in LA that may effect the analysis. Whatever our estimate of the mean, we can be assured it will be misspecified without greater information. With that caveat stated, there are still many avenues of analysis available.

hist(Homeless$StreetTotal, xlab = "Homeless Count", ylab = "Frequency", main = "Fig 2.1: Histogram of LA Homelessness Counts")

3. Data Cleaning

Before examining the relationships between StreetTotal and the other five variables, we will first need to clean and prepare the data. Reviewing the summary statistics for the 6 variables below, a few issues immediately pop out.

summary(Homeless)

  StreetTotal     MedianIncome      PropVacant       PropMinority       PctCommercial    PctIndustrial   
 Min.   :  0.0   Min.   :     0   Min.   :0.00000   Min.   :  0.05231   Min.   : 0.000   Min.   : 0.000  
 1st Qu.:  6.0   1st Qu.: 27315   1st Qu.:0.02504   1st Qu.:  0.36683   1st Qu.: 7.978   1st Qu.: 0.000  
 Median : 16.0   Median : 36800   Median :0.03701   Median :  0.56099   Median :14.614   Median : 0.000  
 Mean   : 37.6   Mean   : 41744   Mean   :0.04807   Mean   :  1.11793   Mean   :17.438   Mean   : 6.777  
 3rd Qu.: 36.0   3rd Qu.: 49853   3rd Qu.:0.05909   3rd Qu.:  0.69608   3rd Qu.:23.068   3rd Qu.: 7.503  
 Max.   :931.0   Max.   :200001   Max.   :0.43003   Max.   :100.00000   Max.   :99.050   Max.   :75.718  
                                  NA's   :4                                                              

First, PropVacant has 4 missing values. Looking at these four values more directly reveals several other quirks in the data. The MedianIncome for these 4 census tracts is zero, which is not believable. In addition, three of the four PropMinority values are 100.00, despite the acceptable range of values being 0.00 to 1.00. The information on these four census tracts is unreliable. Given the reasons above and that the StreetTotal counts are ordinary (range of 5-12), we will remove these four census tracts from the analysis. Luckily, there are no other cases of PropMinority or MedianIncome displaying these issues (MeidanIncome = 0 and PropMinority > 1.00).

Homeless[c(10,445,501,502),]

Let’s review the updated summary statistics for the new data set with the 4 faulty values removed to see if any other issues are immediately obvious. PropMinority’s max value of 1.00 and PropVacant/PctCommercial/PctIndustrial’s lower bounds of 0.00 are worth investigating. A quick glance at PropVacant/PctCommercial/PctIndustrial’s lower bounds of 0.00 did not reveal anything systematically divergent or suggest anything wrong in the data reporting. However, the census tract with PropMinority = 1.00 was a cause for skepticism.

summary(Homeless1)

  StreetTotal      MedianIncome      PropVacant       PropMinority     PctCommercial    PctIndustrial   
 Min.   :  0.00   Min.   :  3779   Min.   :0.00000   Min.   :0.07659   Min.   : 0.000   Min.   : 0.000  
 1st Qu.:  6.00   1st Qu.: 27654   1st Qu.:0.02504   1st Qu.:0.36683   1st Qu.: 7.978   1st Qu.: 0.000  
 Median : 16.00   Median : 37034   Median :0.03702   Median :0.56059   Median :14.614   Median : 0.000  
 Mean   : 37.84   Mean   : 42075   Mean   :0.04807   Mean   :0.53262   Mean   :17.471   Mean   : 6.657  
 3rd Qu.: 36.00   3rd Qu.: 49929   3rd Qu.:0.05909   3rd Qu.:0.69256   3rd Qu.:23.160   3rd Qu.: 7.263  
 Max.   :931.00   Max.   :200001   Max.   :0.43003   Max.   :1.00000   Max.   :99.050   Max.   :75.718  

The census tract with PropMinority = 1.00 proves interesting on multiple levels. It is also one of three census tracts with the highest median income (MedianIncome = $201,000). While it’s difficult to believe the census tract is literally 100% minority, it’s at least very directionally possible. There are 17 other census tracts with the proportion of minority residents greater than 90%, including a tract with 98.1%. What is less believable is the median income of $201,000, given a proportion of minority residents of 100%. Looking at the data, 304 of the 505 census tracts (60%) have a proportion of minority residents over 50%. Of this subset, the median MedianIncome is $30,830 and the mean MedianIncome is $33,560. Aside from the $201,000 value under examination, there is one other value greater than $100,000 ($105,652). For these reasons, it’s likely to believe that this value is either an extreme outlier that will skew the rest of the analysis -or- is misreported. In either case, it’s worth striking from the analysis.

One final adjustment we will make is to convert PropVacant and PropMinority to percentages (like PctCommercial and PctIndustrial). This will help center units in the same range (as opposed to 100x disparity), which will help for statistical reasons. Moving forward PropVacant and PropMinority will be seen as PctVacant and PctMinority.

At this point, after removing NA’s and looking for unreasonable values, the data seem clean and ready for analysis.

4. Analysis of Univariate Statistics

Quickly looking again with our response variable, StreetTotal, we see a skewed distribution. The median count of observed homeless people is 16. The 90th percentile is a count of 83. Yet, the top 5 are: 480, 595, 655, 666, 931. From a Level I perspective, we can state that based on the observed counts, homelessness is not evenly distributed across the 504 census tracts. Instead, homelessness clumps in specific areas as Fig 4.1 helps visualize below. Note: Fig 4.1 only shows the top 40 counts (8% of all census tracts).

hist(HomelessOver100$StreetTotal, xlab = "Homeless Count", ylab = "Frequency", main = "Fig 4.1: Histogram of LA Homelessness Counts Over 100")

Based on an intuitive first pass, at least for anyone who’s spent time in LA, we shouldn’t be too surprised by these findings. It’s fairly rare to see homeless people in affluent suburban neighborhoods and fairly common to see homeless people in either denser parts of the city where people cluster or in industrial parts of the city where homeless people are less likely to be bothered. Before we can generalize from this basic heuristic, we should look at the other variables to see how they are distributed. From there, we can build an understanding of how LA is constructed as a city in terms of racial, income, and property type distributions.

Looking at the other parameters, four of the five variables also form heavily skewed distributions as can be seen below in Fig 4.2a-d. Commonsense says that median income will follow a skewed or logarithmic distribution. It’s also not surprising that 92% of census tracts have residential vacancy less than 10% given LA is a desirable place to live. The distributions of commercial and industrial real estate percentages also reflect our general experiences: Most areas in a major metropolis are not industrial, but a few areas are very industrial. Likewise, commercial and residential tend to blend together more, given people like to shop near where they live for convenience. Again, nothing here is meant to imply causality; instead, we are simply stating that the distributions look reasonable to an intelligent person who has spent time in LA. That said, at this point, the analysis is trending toward Level II, given we are extrapolating from the specific observations within these 504 LA census tracts to inform a broader picture of LA as a city. Given the shape of the distributions, this feels reasonable.

par(mfrow=c(2,2))
hist(HomelessF$MedianIncome, xlab = "Median Income ($)", ylab = "Frequency", main = "Fig. 4.2a: LA Median Incomes")
hist(HomelessF$PctVacant, xlab = "Vacancy (%)", ylab = "Frequency", main = "Fig. 4.2b: LA % Vacant (Residential)")

hist(HomelessF$PctCommercial, xlab = "Commercial (%)", ylab = "Frequency", main = "Fig. 4.2c: LA % Commercial Real Esate")
hist(HomelessF$PctIndustrial, xlab = "Industrial (%)", ylab = "Frequency", main = "Fig. 4.2d: LA % Industrial Real Esate")

Looking at distribution of the percentage of residents that are minorities within a given census tract, we see a very different picture. The distribution is fairly normal, with a slightly fat tail built up on the left side (Fig 4.3). The mean and median are 53% and 56%, respectively, and reflect a diverse population. Looking at Wikipedia’s demographic information on LA for the 2010 census, Whites made up 29%, Latinos 49%, Blacks 10%, and Asians 11% (https://en.Wikipedia.org/wiki/Los_Angeles#Demographics). What our specific census data on the 504 tracts reflects is that most individual census tracts are comprised of a diverse mix of people. Again, on a Level I analysis, we see heterogeneity among and within the census tracts. Pushing to Level II, the Wikipedia information corroborates these findings for LA as a whole.

hist(HomelessF$PctMinority, xlab = "Minorities (%)", ylab = "Frequency", main = "Fig 4.3: LA Percentage Minorities") 

5. Analysis of Bivariate Statistics

Now that we have a firm understanding of the data, it’s time to explore how our response, StreetTotal, interacts with the other variables.

First, we will set a baseline using simple linear regression. For this, we will remain strictly in Level I analysis and not split our data in anyway. The goal here is merely to understand the direction and magnitude of effect each variable has on the response. In the code output below, we see the intercept and linear regression coefficient for each variable against StreetTotal (independently as single variable regressions).

co.Com

  (Intercept) PctCommercial 
   21.6521639     0.9289005 

co.Ind

  (Intercept) PctIndustrial 
    30.113596      1.169307 

co.MI

  (Intercept)  MedianIncome 
69.8774712945 -0.0007654172 

co.Min

(Intercept) PctMinority 
 10.3413383   0.5185553 

co.Vac

(Intercept)   PctVacant 
  10.981298    5.590858 

These coefficients state the following about the relationship with each variable and StreetTotal, for only these 504 census tracts. Again, to reiterate, this in not implying causality or statistical significance, merely a description of the data:

• PctCommercial: For every 1% point increase in percent of commercial land in a given census tract, the count of homeless observations increases by 0.93

• PctIndustrial: For every 1% point increase in percent of industrial land in a given census tract, the count of homeless observations increases by 1.17

• MedianIncome: For every $1 increase in the median income in a given census tract, the count of homeless observations decreases by 0.00076. Translated, for every increase of $1,000, the count of homeless observations decreases by 0.76

• PctMinority: For every 1% point increase in percent of the residents that are minorities in a given census tract, the count of homeless observations increases by 0.52

• PctVacant: For every 1% point increase in percent of residential property that is vacant in a given census tract, the count of homeless observations increases by 5.59

Reviewing these descriptions, it’s easy to see how some of these may generalize. On the whole, we’d intuitively expect to find higher counts of homeless people in places with industrial or vacant property, which the data reflect. Additionally, we’d expect to see homelessness counts drop in more affluent areas, which again the data reflect. The percent of minority residents has the weakest effect (among the four directly comparable percent values), which feels right given that we observed LA is a diverse community, both within and among census tracts. The slightly positive effect likely stems from the relationship we observed earlier between percent minority and median income. We look at this relationship again, plotted below as Fig 5.1. The relationship between Median Income and % Minority is weak, but noticeable. Again, all of the analysis directly above should be seen as a Level I description of the data (with the occasional sanity check), nothing more.

plot(MedianIncome~PctMinority, main = "Fig. 5.1: Median Income by % Minority")

Next, we look at an application of smoothing splines to the data to move beyond basic linear methods. The idea here is that we are looking at a highly skewed response distribution (as well as highly skewed prediction variables). Given this, we don’t anticipate a linear model will best fit the data. Instead of forcing a linear equation onto the data, we can instead using smoothing splines to identify knots (think of them as inflection points) in the data. The technical math is beyond the scope of this paper (and my abilities!), but the key concept is that for any given range of X, we will apply a given function, F(x) that isn’t constrained to being linear. When there is sufficient movement in the data away from F(x), we will introduce a knot, and then describe a new F(x) over a new range of x. As the number of knots increases, we become prone to over fitting the data. To combat this, we can segment the data into training, evaluation, and test samples. For the sake of speed, we will simply use the default method here of leave-one-out cross validation, which roughly approximates the results you’d typically see in a training, evaluation, and test split.

The three variables we will examine are: PctIndustrial (Fig 5.2a-d), PctVacant (Fig 5.3a-d), and MedianIncome (Fig 5.4a-d). We will examine four underlying concepts in univariate analysis. Each thought will be tied to a specific plot to help elucidate the point:

1. We are still in Level 1 analysis. Any attempt to interpret coefficients, degrees of freedom, or model fit is futile. There are too many statistical processes occurring under the hood to infer causality or statistical significant. The curves simply reflect that data as presented among the 504 random census tract samples. (Note: General, could be applied to all plots.)

2. When SPAR is small (0.5, in this case), the fits are extremely bumpy and heavily influenced by outliers. For example, in graph 5.2a: HL Count by % Indust.(spar=0.5), the sharp bump at PctIndustrial = 20% does not imply some magic increase in homelessness occurs then. It simply means that the large outliers in the count are pulling the function up at that moment. To be able to extend this to Level II analysis, we’d need a DEEP subject matter expert, who then could maybe infer that there are specific zoning laws/benefit programs tied to zoning laws in LA that have an artificial cliff at 20%. But for us to draw this type of conclusion from the data as presented would be irresponsible.

3. The distribution of the underlying variables being analyzed matters. To expand on this, let’s compare two plots: 5.2b: HL Count by % Indust.(spar=1.0) & 5.4b: HL Count by Med Income (spr=1.0). Even when SPAR = 1.0 (theoretically the max penalty for SPAR; essentially the algorithm is trying to add as few knots as possible), in 5.2b we see a large spike at Postindustrial = 65% before a sharp decline at PctIndustrial = 75%. Because there are so few observations at the tail end of this distribution (4 of 504 above 60%, i.e. <1%), the model will likely overfit. In contrast, in figure 5.4b, we see a very “linear” looking relationship all the way out to median income of $201,000. However, while homelessness does seem to decline with increasing median, the broader point is that the data in that tail is irrelevant given it has so little support. We could likely find a district in the Bay Area (my home city) that has a sky high median income and also a high level of homelessness: Downtown Palo Alto. The point is even with advanced models, do not place too much weight in any explanation that relies on weakly supported data. When dealing with income, this will almost always be the case given it’s often skewed underlying distribution.

4. Lastly, we discuss the Bias-Variance trade off by looking at high values of SPAR (SPAR = 1.5) as an approximation of a linear model. The simple conceptual point about the algorithm is: with enough penalty placed on each additional knot, the model will eventually approximate a linear model. In practice, we want to keep SPAR <= 1.0. However, to elucidate the underlying mechanisms of the tool, compare figures: 5.2c: HL Count by % Indust.(spar=1.5) & 5.2d: HL Count by % Indust.(linear). The lines are nearly indistinguishable. What is so important to generalize is: this highlights the larger and incredibly important overall Bias-Variance trade off. When SPAR is low, it better approximates the data as presented. When SPAR is sufficiently high, it approximates a linear model. In this case, we know the low value of SPAR is overfitting to the data in an attempt to capture all of the variance. When SPAR is high, it sacrifices trying to explain all the variance and instead biases its results toward the overall trend. To see this very clearly, look at the spike vs. non-spike at PctIndustrial = 65% in 5.2a vs. 5.2d. At SPAR = 1.5, the algorithm ignores the large outliers to preserve the overall trend.

par(mfrow=c(2,2))
out<-gam(StreetTotal~s(PctIndustrial, spar=.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.2a: HL Count by % Indust.(spar=0.5)")
out<-gam(StreetTotal~s(PctIndustrial, spar=1.0),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.2b: HL Count by % Indust.(spar=1.0)")

out<-gam(StreetTotal~s(PctIndustrial, spar=1.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.2c: HL Count by % Indust.(spar=1.5)")
plot(StreetTotal~PctIndustrial, ylab = "Count", main = "5.2d: HL Count by % Indust.(linear)")

abline(a=PctIndustrial, b=StreetTotal)
par(mfrow=c(2,2))

out<-gam(StreetTotal~s(PctVacant, spar=.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.3a: HL Count by % Vacant (spar=0.5)")
out<-gam(StreetTotal~s(PctVacant, spar=1.0),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.3b: HL Count by % Vacant (spar=1.0)")

out<-gam(StreetTotal~s(PctVacant, spar=1.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.3c: HL Count by % Vacant (spar=1.5)")
plot(StreetTotal~PctVacant, ylab = "Count", main = "5.3d: HL Count by % Vacant (linear)")

abline(a=PctVacant, b=StreetTotal)
par(mfrow=c(2,2))

out<-gam(StreetTotal~s(MedianIncome, spar=.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.4a: HL Count by Med Income (spr=0.5)")
out<-gam(StreetTotal~s(MedianIncome, spar=1.0),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.4b: HL Count by Med Income (spr=1.0)")

out<-gam(StreetTotal~s(MedianIncome, spar=1.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.4c: HL Count by Med Income (spr=1.5)")
plot(StreetTotal~MedianIncome, ylab = "Count", main = "5.4d: HL Count by Med Income (linear)")

abline(a=MedianIncome, b=StreetTotal)

6. Analysis of Multivariate Statistics

Finally, we get to the General Additive Model (GAM): our principle topic of inquiry. Before we start, it’s important to lay out a few statistical concepts. First, as noted in Section 2: Data Description, our estimate using GAM will lack many potentially explanatory variables. In short, we have the 5 predictors we’ve observed, but there is a lot more information required to specify the true mean function. Given that, we will still be OK moving ahead at Level I and potentially Level II analysis, but we will not be able to draw conclusions at Level III. Second, in Section 5. Analysis of Bivariate Statistics, we pointed out that it’s unwise to put significant weight on any explanation that is weakly supported by few observations. We saw that it took a SPAR = 1.5 to “ignore” the leverage points in 5.2c: HL Count by % Indust.(spar=1.5). Well, with GAM we are adding more variables to the models and hence creating hyper dimensional space. As a result, any observation in n-dimensional space will be less supported than it was in 2D space. Both response outliers and sparse patches of data will have to be examined cautiously. Finally, we will completely depart from the idea of regression coefficients when using GAM. The transformations involved in creating the model are sufficiently complex that any interpretation of the coefficients would depart from reality.

To evaluate our model, we will split the data into three buckets: training data (60%), evaluation data (20%), and test data (20%). This allows us to fit a model based on training data, tweak parameters via evaluation data (some snooping), and then test the model’s fit on the test data. This framework is a good way to prevent overfitting. One crucial point to emphasize: with our relatively small sample set (n=504) and existence of both massive outliers in the response and highly skewed underlying distributions for four of the five predictors, we will have to think critically about how to extrapolate results and extra beware of overfitting.

Based on my initial split, the data look reasonably stable. There are large outliers for StreetTotal among all three splits. In addition, the medians, mins, and maxs across the predictive variables are tolerably consistent. PctIndustrial’s 4 outliers above 60% may cause issues. Less noticeable, MedianIncome and PctCommercial may exhibit slight bias given outlier values. Overall, however, the split seems reasonable and in the spirit of a random split, we will move forward.

6.1: Test Data

After building the GAM model, we look at the fitted Homeless Counts on the training data for varying levels of SPAR (0.5, 0.75, 1.0, 1.5). Figure 6.1a-d displays these histograms below. From our thorough analysis before, we know that SPAR = 0.5 is highly sensitive and prone to overfit. In contrast, SPAR = 1.5 approaches a linear model, which is reflected in the (nearly) normal distribution produced in 6.1d. Looking at the underlying fits for SPAR = 1.5 on an individual element basis, the trend to approximate a linear model continues. SPAR = 0.75 and 1.0 are very hard to distinguish. On a more substantive level, how SPAR = 0.75 and 1.0 differ is interesting. SPAR = 0.75 (Fig. 6.1b) leaves more data in the tails (both at zero and the far right). From a policy level, 0 homeless vs. <5 homeless people is likely not the core concern. However, areas with massive homeless counts are likely very interesting to the policy makers. Given these factors, let’s error on the side of better fit and proceed with the evaluating the SPAR = 0.75 model on our evaluation data.

par(mfrow=c(2,2))
hist(fitted.values(out.5),breaks=20, main = "6.1a: Fitted HL Counts (SP=0.5)") 
hist(fitted.values(out.75),breaks=20, main = "6.1b: Fitted HL Counts (SP=0.75)") 

hist(fitted.values(out1),breaks=20, main = "6.1c: Fitted HL Counts (SP=1.0)") 
hist(fitted.values(out1.5),breaks=20, main = "6.1d: Fitted HL Counts (SP=1.5)") 

6.2: Evaluation Data

Interestingly, SPAR = 0.75 does a fairly poor job predicting on the evaluation data. To compare the MSE, which is the fit estimate for this type of model, we can see:

• SPAR 1.0 MSE = 9,235
• SPAR 0.75 MSE = 7,080
• SPAR 0.5 MSE = 2892

The bias-variance trade off favors the 0.5 model, but potentially too much at the expense of overfitting. Given this, let’s re-train our model to have a SPAR = 0.6, to see if we can get more “value” from at a SPAR value between the two.

After the re-fitting SPAR = 0.6 model on the training data, we see the histogram in Fig 6.2 below. On a descriptive basis, it’s hard to distinguish from either SPAR = 0.5 and SPAR = 0.75.

hist(fitted.values(out0.6),breaks=20, main = "6.2: Fitted HL Counts (SP=0.6) (Training Data)") 

Moving back to the evaluation data, we see the fit as shown in Fig 6.2.1a-b. The MSE is 4409. While we moved SPAR = 0.6 40% of the distance between SPAR = 0.5 and 0.75, we only lost 36% of the explanatory power. Most importantly, the SPAR = 0.6 model still seems open to outliers, having predicted a max over 400 compared to just 250 in the SPAR = 0.75 evaluation. While using rough heuristics at best, it seems the trade is worth it.

hist(preds.75, breaks=20, main = "6.2.1a: Evaluated HL Counts (SP=0.75) (Evaluation Data)")

hist(preds0.6, breaks=20, main = "6.2.1b: Evaluated HL Counts (SP=0.6) (Evaluation Data)")

6.3: Test Data

We now evaluate the SPAR = 0.6 model on test date. The model fits well, with an MSE of 1546 and the histogram as seen in Fig 6.3. What’s especially interesting is that the model predicts some values to be negative. The range is: [-29.73, 586.63]. From an interpretive standpoint, we’d just expect those areas to be “very low” homelessness zones, not to imply there can actually be negative homeless people. On a more substantive, nuanced level, it’s not irrational to assume just because an observation for homeless counts over two consecutive days was zero, that there not ever homeless people in the area. In that sense, with much more granular data and a much deeper subject matter expertise, we might be able to quantify/predict how “insulated” an area is from homelessness (i.e. it’s distance below zero). Finally, it’s important to note that the model showed two outliers above 420. Looking back at our entire n=504 sample, we had six observations above 420. Given the testing data comprises 20% of the data, we’d expect at random to see 1.2 observations above 420. Our observation of 2 seems totally plausible (and would happen 20% of the time).

hist(tpreds0.6, main = "6.3: Test HL Counts (SP=0.6) (Test Data)")

Section 7: Summary and Conclusions

As can be seen in Fig 7.1 below and in the MSE score of 1546, the GAM model with a SPAR = 0.6 does a good job of fitting both the training and test data.

par(mfrow=c(3,2))
plot(out0.6,residual=T,cex=1,pch=14,col="light blue", main = "Fig 7.1: SPAR = 0.6 Fit")

We can now summarize our findings and conclusions. There are three core findings I’d like to discuss:

1. From a Level I perspective, we can see that the data are complex with many features. Addressing only these 504 census tracts, we cannot approach the true mean specification function, however we can observe directional and magnitude relationships empirically for these observations. For example, low Median Income shows a rise in observed homelessness when Median Income is less than $50k. After that, with fewer data points the relationship a) becomes less stable and b) flattens out. Likewise, the percentage of commercial property is fairly bumpy (but flat overall) up until 50%. Beyond 50%, homeless count observations grow rapidly, before tapering back off near 80%. Again, towards the ends of this distribution there are both a) less observations and b) significant outliers.

2. From a Level II perspective, which I believe we can now just barely justify, we can likely extrapolate our GAM from the n=504 census tracts to the other n=1500+ census tracts in LA and get sensible results. I believe this for a two reasons. First, our GAM model for SPAR = 0.6 performed quite well on the testing hold-out. Assuming the n=504 census tracts were genuinely selected at random, it’s likely we’d see this trend continue in the other 75% of the sample hold out. Second, our data covered a seemingly broad sample of the potential values. From Minority rates of 98% and Vacancy rates up to 43%, to homeless counts of 931 and median incomes of $200k, we seem to have a stable slice of the overall population. If any of these distributions were severely contained, I’d be much more hesitant to extend to Level II. Of course, with that said, I’d forewarn of three potential error (among countless others) that might be seen extrapolating the data to the other 75% of census tracts. First, the underlying distributions of PctIndustrial, PctCommercial, PctVacant, and Median Income are all sensitive in the tail. In PctIndustrial and PctCommercial we see significant spikes where as Median Income and PctVacant do not. It’s not hard to imagine counter-factual in either direction. Second, PctMinority is a broad statistic. We have very little understanding of its underlying composition. It may be the case that when certain minority groups comprise the vast majority of the population in that census tract, the response variable changes. We simply have no way of knowing at this point given the summary data. Third, and finally, this information is temporally isolated to one moment in time. If it were to be extended, the actual factors describing where homeless people congregate may have changed in any given area. This is the issue with only having five predictive variables – we become more sensitive to any number of the many unobserved changes.

3. Finally, writing this paper was fun and made me think about a lot of the core underlying assumptions I carry into problems. When I decided to get my MBA at Wharton, the goal was to come statistically fluent. It’s a joy to have climbed up the ladder and gotten in-over-my-head, and now climb back down the ladder and unpack elements of what’s happening under the covers. Thank you!

Appendix: R Code

###### first glance at data
summary(Homeless)
hist(Homeless$StreetTotal, xlab = "Homeless Count", ylab = "Frequency", main = "Fig 2.1: Histogram of LA Homelessness Counts")

##### EDA
## Find NA's and remove if sensible
which(is.na(Homeless$PropVacant))
Homeless[c(10,445,501,502),]
Homeless1 <- Homeless[-c(10,445,501,502),]
Homeless1
## PropMinority == 1.00 ??? (removed)
sum(Homeless1$PropMinority > 0.00)
sum(Homeless1$MedianIncome == 0)
which.max(Homeless1$PropMinority)
Homeless1[496,]
sum(Homeless1$PropMinority > 0.90)
sum(Homeless1$PropMinority > 0.50)
hold1 <- which(Homeless1$PropMinority > 0.50)
hold2 <- Homeless1$MedianIncome[hold1]
sort(hold2)
summary(hold2)
## Final dataset with removed NA's/outliers (and Prop converted to %)
Homeless2 <- Homeless1[-496,]
Homeless2["PctVacant"] <- Homeless2$PropVacant*100
Homeless2["PctMinority"] <- Homeless2$PropMinority*100
summary(Homeless2)
HomelessF <- Homeless2[,-c(3,4)]
summary(HomelessF)

###### Univariate Analysis
## Setup / Basic Hist
summary(Homeless2)
sum(Homeless2$StreetTotal > 83)
sort(Homeless2$StreetTotal)
hist(Homeless2$StreetTotal)
## Top 100 counts & Histogram
HomelessOver100 <- subset(Homeless2, Homeless2$StreetTotal >= 100)
hist(HomelessOver100$StreetTotal, xlab = "Homeless Count", ylab = "Frequency", main = "Fig 4.1: Histogram of LA Homelessness Counts Over 100")
sum(Homeless2$StreetTotal >= 100)
## Reviewing other variable histograms
## skewed 
par(mfrow=c(2,2))
hist(Homeless2$MedianIncome, xlab = "Median Income", ylab = "Frequency", main = "Histogram of LA Median Incomes")
hist(Homeless2$PropVacant, xlab = "Proportion Vacant", ylab = "Frequency", main = "Histogram of LA Proportion Vacant (Residential)")
hist(Homeless2$PctCommercial, xlab = "Percent Commercial", ylab = "Frequency", main = "Histogram of LA % Commercial Real Esate")
hist(Homeless2$PctIndustrial, xlab = "Percent Industrial", ylab = "Frequency", main = "Histogram of LA % Industrial Real Esate")
sum(Homeless2$PropVacant < .1)/504
## normal-ish
hist(Homeless2$PropMinority, xlab = "Minorities (Proportion)", ylab = "Frequency", main = "Fig 4.3: LA Proportion Minorities") 
summary(Homeless2$PropMinority)


##### Bivariate Analysis (Regression --> SPAR)
## Normal GLM regression (extracting coefficients)
attach(HomelessF)
mod.MI <- glm(StreetTotal~MedianIncome)
summary(mod.MI)
co.MI <- coef(mod.MI)

mod.Vac <- glm(StreetTotal~PctVacant)
summary(mod.Vac)
co.Vac <- coef(mod.Vac)

mod.Com <- glm(StreetTotal~PctCommercial)
summary(mod.Com)
co.Com <- coef(mod.Com)

mod.Ind <- glm(StreetTotal~PctIndustrial)
summary(mod.Ind)
co.Ind <- coef(mod.Ind)

mod.Min <- glm(StreetTotal~PctMinority)
summary(mod.Min)
co.Min <- coef(mod.Min)

sum.stats <- c(co.Com, co.Ind, co.MI, co.Min, co.Vac)
sum.stats

mod.MIbyMin<- glm(MedianIncome~PctMinority)
summary(mod.MIbyMin)
co.MIbyMin <- coef(mod.MIbyMin)
plot(MedianIncome~PctMinority, main = "Median Income by % Minority")

## SPAR for single variable models (smoothing splines)
## Smoothing Splines
library(gam)
## % Ind
par(mfrow=c(2,2))
out<-gam(StreetTotal~s(PctIndustrial, spar=.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "Fig 5.2a: Homeless Count by % Industrial (spar=0.5)")
out<-gam(StreetTotal~s(PctIndustrial, spar=1.0),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "Fig 5.2b: Homeless Count by % Industrial (spar=1.0)")
out<-gam(StreetTotal~s(PctIndustrial, spar=1.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "Fig 5.2c: Homeless Count by % Industrial (spar=1.5)")
plot(StreetTotal~PctIndustrial, ylab = "Count", main = "Fig 5.2d: Homeless Count by % Industrial (linear)")
abline(a=PctIndustrial, b=StreetTotal)

## % Vacant
par(mfrow=c(2,2))
out<-gam(StreetTotal~s(PctVacant, spar=.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.3a: Homeless Count by % Vacant (spar=0.5)")
out<-gam(StreetTotal~s(PctVacant, spar=1.0),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.3b: Homeless Count by % Vacant (spar=1.0)")
out<-gam(StreetTotal~s(PctVacant, spar=1.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "5.3c: Homeless Count by % Vacant (spar=1.5)")
plot(StreetTotal~PctVacant, ylab = "Count", main = "5.3d: Homeless Count by % Vacant (linear)")
abline(a=PctVacant, b=StreetTotal)

## Med Inc
par(mfrow=c(2,2))
out<-gam(StreetTotal~s(MedianIncome, spar=.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "Fig 5.4a: Homeless Count by Median Inccome (spar=0.5)")
out<-gam(StreetTotal~s(MedianIncome, spar=1.0),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "Fig 5.4b: Homeless Count by Median Inccome (spar=1.0)")
out<-gam(StreetTotal~s(MedianIncome, spar=1.5),data=HomelessF)
plot(out, residuals=T, ylab = "Count", main = "Fig 5.4c: Homeless Count by Median Inccome (spar=1.5)")
plot(StreetTotal~MedianIncome, ylab = "Count", main = "Fig 5.4d: Homeless Count by Median Income (linear)")
abline(a=MedianIncome, b=StreetTotal)


##### Multivariate Analysis (GAM)
## Creating Training/Eval/Test Split 
## Set seed
set.seed(17)

## Split data
spec = c(train = .6, evaluation = .2, test = .2)

g = sample(cut(
  seq(nrow(HomelessF)), 
  nrow(HomelessF)*cumsum(c(0,spec)),
  labels = names(spec)
))

HomelessSplit = split(HomelessF, g)
sapply(HomelessSplit, nrow)/nrow(HomelessF)
## Look at summary statistics
summary(HomelessSplit$train)
hist(HomelessSplit$train$StreetTotal)
summary(HomelessSplit$evaluation)
hist(HomelessSplit$evaluation$StreetTotal)
summary(HomelessSplit$test)
hist(HomelessSplit$test$StreetTotal)

## GAM Model
## Gam library
library(gam)
## spar = 0.5
out.5<-gam(StreetTotal~s(MedianIncome,spar=.5)+s(PctCommercial,spar=.5)+
            s(PctIndustrial,spar=.5)+s(PctMinority,spar=.5)+s(PctVacant,spar=.5),
            family=gaussian) # Specify spar
summary(out.5)
par(mfrow=c(1,1))
plot(out.5,ask=T) # centered around the mean of the response
hist(fitted.values(out.5),breaks=20) 

## spar = 0.75
out.75<-gam(StreetTotal~s(MedianIncome,spar=.75)+s(PctCommercial,spar=.75)+
            s(PctIndustrial,spar=.75)+s(PctMinority,spar=.75)+s(PctVacant,spar=.75),
          family=gaussian) # Specify spar
summary(out.75)
par(mfrow=c(1,1))
plot(out.75,ask=T) # centered around the mean of the response
hist(fitted.values(out.75),breaks=20) 

## spar = 1.0
out1<-gam(StreetTotal~s(MedianIncome,spar=1)+s(PctCommercial,spar=1)+
            s(PctIndustrial,spar=1)+s(PctMinority,spar=1)+s(PctVacant,spar=1),
          family=gaussian) # Specify spar
summary(out1)
par(mfrow=c(1,1))
plot(out1,ask=T) # centered around the mean of the response
hist(fitted.values(out1),breaks=20) 

## spar = 1.5
out1.5<-gam(StreetTotal~s(MedianIncome,spar=1.5)+s(PctCommercial,spar=1.5)+
            s(PctIndustrial,spar=1.5)+s(PctMinority,spar=1.5)+s(PctVacant,spar=1.5),
          family=gaussian) # Specify spar
summary(out1.5)
par(mfrow=c(1,1))
plot(out1.5,ask=T) # centered around the mean of the response
hist(fitted.values(out1.5),breaks=20) 

par(mfrow=c(2,2))
hist(fitted.values(out.5),breaks=20) 
hist(fitted.values(out.75),breaks=20) 
hist(fitted.values(out1),breaks=20) 
hist(fitted.values(out1.5),breaks=20) 

## spar = 0.6
out0.6<-gam(StreetTotal~s(MedianIncome,spar=0.6)+s(PctCommercial,spar=0.6)+
              s(PctIndustrial,spar=0.6)+s(PctMinority,spar=0.6)+s(PctVacant,spar=0.6),
            family=gaussian) # Specify spar
summary(out0.6)
par(mfrow=c(1,1))
plot(out0.6,ask=T) # centered around the mean of the response
hist(fitted.values(out0.6),breaks=20) 

## Evaluation Data (SPAR = 1.0)
preds.1<-predict(out1, newdata=HomelessSplit$evaluation) # get fitted values
hist(preds.1) # what does the distribution look like?
resids.1<-HomelessSplit$evaluation$StreetTotal-preds.1 # get residuals
var(resids.1) # MSE for the residuals as a measure of fit

## Evaluation Data (SPAR = 0.75)
preds.75<-predict(out.75, newdata=HomelessSplit$evaluation) # get fitted values
hist(preds.75) # what does the distribution look like?
resids.75<-HomelessSplit$evaluation$StreetTotal-preds.75 # get residuals
var(resids.75) # MSE for the residuals as a measure of fit

## Evaluation Data (SPAR = 0.5)
preds.5<-predict(out.5, newdata=HomelessSplit$evaluation) # get fitted values
hist(preds.5) # what does the distribution look like?
resids.5<-HomelessSplit$evaluation$StreetTotal-preds.5 # get residuals
var(resids.5) # MSE for the residuals as a measure of fit

## Evaluation Data (SPAR = 0.6)
preds0.6<-predict(out0.6, newdata=HomelessSplit$evaluation) # get fitted values
hist(preds0.6) # what does the distribution look like?
resids0.6<-HomelessSplit$evaluation$StreetTotal-preds0.6 # get residuals
var(resids0.6) # MSE for the residuals as a measure of fit

## Test data
## Test Data (SPAR = 0.6)
tpreds0.6<-predict(out0.6, newdata=HomelessSplit$test) # get fitted values
hist(tpreds0.6) # what does the distribution look like?
tresids0.6<-HomelessSplit$test$StreetTotal-tpreds0.6 # get residuals
var(tresids0.6) # MSE for the residuals as a measure of fit

range(tpreds0.6)
sort(tpreds0.6)

##### Summary and Conclusion
par(mfrow=c(3,2))
plot(out0.6,residual=T,cex=1,pch=14,col="light blue", main = "Fig 7.1: SPAR = 0.6 Fit")