Topic 8. Machine Learning I. Regression

We’ll start out our Machine Learning tutorial using a familiar dataset: the Opportunity Insights data on county-level inequality. We’ll predict inequality levels using the data. We’ll begin our journey into machine learning using a familiar method: linear regression. The ideas are largely the same as the previous tutorial, but this time we’ll start to think about regression from a machine learning perspective.

data <- read.csv("county_data.csv")

We’ll build a model predicting inequality at the county level, using several co-variates to make our prediction (in ML, covariates are often referred to as “features” in the model). This is considered “supervised machine learning” because we already know the actual inequality level in each county.

First, we’ll select our variables of interest. Then, we need to split our data into training and testing sets. We’ll use our training set to train the model, while our testing set will be used to validate our model.

#select the variables we'll use
library(dplyr)
model_data <- data %>% 
  select(c("inequality", "racial_segregation", "crime_rate", "unemp_rate", "percent_college_degree")) %>% 
  mutate(college = percent_college_degree/100) %>% 
  select(-"percent_college_degree")

#Check the mean of our dependent variable for reference
mean(model_data$inequality, na.rm=T)
sd(model_data$inequality, na.rm=T)

#split into training and testing data
#setting the seed allows us to recreate this analysis later
set.seed(8675309)

#we need to draw random samples for our training and testing data - this code does this for us
#The training data will be 75% of the data, and the remaining 25% will be the testing/validation data
split <- 0.75
rows  <- nrow(model_data)

train.entries <- sample(rows, rows*split)

model.train <- model_data[train.entries, ]
model.valid  <- model_data[-train.entries,  ]

model <- lm(inequality ~ racial_segregation, data=model.train)

summary(model)

Ok, now we have the model based on the test data. Let’s use this model on the validation data, and assess its fit.

model.valid <- model.valid %>%
  #use the model to predict inequality in the validation data
  #save the predicted coefficients from this model
    mutate(yhat = predict(model, newdata=model.valid)) %>%
  #calculate the residual so we can plot it
    mutate(residual = inequality - yhat)

#Save this for later, when we compare the two models
model.train <- model.train %>%
    mutate(yhat = predict(model, newdata=model.train)) %>%
    mutate(residual = inequality - yhat)

head(model.valid)

mean(model.valid$residual, na.rm=T)

Now, we’ll plot the residuals to take a look at the fit of the model:

library(ggplot2)
ggplot(model.valid) +
  geom_point(aes(x=yhat, y=residual)) +
  geom_hline(aes(yintercept=0), linetype="dashed", color="red") +
  xlab("Predicted Inequality") +
  ylab("Residual Inequality")+
  theme_minimal()+
  labs(title = "Residual Plot for Validation Model")

This residual is not perfect, but the spread around zero is relatively consistent, which is a good sign for us. The data does seem to be clustered around lower predicted levels of inequality though, rather than spread out across the distribution.

Is this an issue? Maybe not. This may indicate that one of our variables is skewed, and we can try transforming (e.g. taking the logarithm) the variable and see if that improves the model fit. This may also indicate that we’re missing a variable.

Let’s take a closer look at the model error - what kind of error are we looking at here?

ggplot(model.valid) +
  geom_histogram(aes(x=residual), bins=50, color="black", fill = "lightblue") +
  xlab("Residual (inequality)")+
  theme_minimal()+
  labs(title = "Distribution of Model Residuals")

There are some outliers, but most of the residual is relatively normally distributed with a mean of 0. We generally like to see residuals that are normally distributed. On the other hand, some of the variables have an error as high as 0.25! For a variable with a mean of 0.38, that is some substantial error.

Finally, we’ll check for overfitting today. To do this, we’ll compare the standard deviations of the training and test set residuals. If the test residuals sd is much larger than the sd for the training data, then we know that our model has been overfit.

#Check the training model
sd(model.train$residual, na.rm=T)

#now compare with the validation model
sd(model.valid$residual, na.rm=T)

These are very close! This tells me that my model is likely not overfitted. If my SD were, for example, 2x the size of the training model residuals, I’d be concerned.

Finally, the mean squared error (mse) and root mean squared error (rmse) are additional measures of model accuracy that are common in machine learning:

model.valid <- model.valid %>% 
  mutate(square_error = residual^2)

#mse
mean(model.valid$square_error, na.rm=T)

#calculate rmse - this is very similar to an SD for residuals. You'll often see this terminology used in ML instead of SD to assess model fit
sqrt(mean(model.valid$square_error, na.rm=T))

We’ll come back to this measure later!