What is statistical learning?

Statistical models attempt to summarize relationships between variables by reducing the dimensionality of the data. For example, here we have some simulated data on sales of Shamwow in 200 different markets. Our goal is to improve sales of the Shamwow. Since we cannot directly increase sales of the product (unless we go out and buy it ourselves), our only option is to increase advertising across three potential mediums: newspaper, radio, and TV.

In this example, the advertising budgets are our input variables, also called independent variables, features, or predictors. The sales of Shamwows is the output, also called the dependent variable or response.

By plotting the variables against one another using a scatterplot, we can see there is some sort of relationship between each medium’s advertising spending and Shamwow sales:

# get advertising data
advertising <- read_csv("data/Advertising.csv") %>%
  tbl_df() %>%
  select(-X1)

# plot separate facets for relationship between ad spending and sales
plot_ad <- advertising %>%
  gather(method, spend, -Sales) %>%
  ggplot(aes(spend, Sales)) +
  facet_wrap(~ method, scales = "free_x") +
  geom_point() +
  labs(x = "Spending (in thousands of dollars)")
plot_ad

But there seems to be a lot of noise in the data. How can we summarize this? We can do so by estimating a mathematical equation following the general form

\[Y = f(X) + \epsilon\]

where \(f\) is some fixed, unknown function of the relationship between the independent variable(s) \(X\) and the dependent variable \(Y\), with some random error \(\epsilon\).

Statistical learning refers to the set of approaches for estimating \(f\). There are many potential approaches to defining the functional form of \(f\). We already saw one method for doing so, ordinary least squares (OLS). Applied here, the results would look like:

plot_ad +
  geom_smooth(method = "lm", se = FALSE)

However statistical learning (and machine learning) allows us to use a wide range of functional forms beyond a simple linear model.

Why estimate \(f\)?

We’ve already discussed the two major goals of statistical modeling:

1. Prediction - use our knowledge of the relationship between \(X\) and \(Y\) to predict \(Y\) for given values of \(X\). Often the function \(f\) is treated as a black box - we don’t care what the function is, as long as it makes accurate predictions. If we are trying to boost sales of Shamwow, we may not care why specific factors drive an increase in sales - we just want to know how to adjust our advertising budgets to maximize sales.
2. Inference - use our knowledge of \(X\) and \(Y\) to understand the relationship between the variables. Here we are most interested in the explanation, not the prediction. So in the Shamwow example, we may not care about actual sales of the product - instead, we may be economists who wish to understand how advertising spending influences product sales. We don’t care about the actual product, we simply want to learn more about the process and generalize it to a wider range of settings.

How do we estimate \(f\)?

There are two major approaches to estimating \(f\): parametric and non-parametric methods.

method_model <- function(df) {
  lm(Sales ~ spend, data = df)
}

ad_pred <- advertising %>%
  gather(method, spend, -Sales) %>%
  group_by(method) %>%
  nest() %>%
  mutate(model = map(data, method_model),
         pred = map(model, broom::augment)) %>%
  unnest(pred)

plot_ad +
  geom_smooth(method = "lm", se = FALSE) +
  geom_linerange(data = ad_pred,
                 aes(ymin = Sales, ymax = .fitted),
                 color = "blue",
                 alpha = .5) 

Supervised vs. unsupervised learning

All the examples above implement supervised learning. That is, for each observation we have both the predictor measurements and the response measurements (i.e. an \(X\) and \(Y\)). We seek to fit a model that summarizes the relationship between the predictors and response.

In unsupervised learning, all we have is a set of measurements \(X\) for a series of observations, but no corresponding response \(Y\). Without an outcome measure, we cannot fit a linear regression model or employ a similar method. That does not mean we cannot use statistical learning to understand the data better. One example of unsupervised learning is cluster analysis. The goal is to determine whether the observations fall into distinct categories. Latent Direchlet allocation (LDA) is an example of cluster analysis applied to text data. In LDA, the individual words are the features or measurements we use to determine the best fitting clusters.

A linear regression approach

Let’s concentrate first on the relationship between age and survival. Using the methods we previously learned, we could estimate a linear regression model:

ggplot(titanic, aes(Age, Survived)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

Hmm. Not terrible, but you can immediately notice a couple of things. First, the only possible values for Survival are \(0\) and \(1\). Yet the linear regression model gives us predicted values such as \(.4\) and \(.25\). How should we interpret those?

One possibility is that these values are predicted probabilities. That is, the estimated probability a passenger will survive given their age. So someone with a predicted probability of \(.4\) has a 40% chance of surviving. Okay, but notice that because the line is linear and continuous, it extends infinitely in both directions of age.

ggplot(titanic, aes(Age, Survived)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, fullrange = TRUE) +
  xlim(0, 200)

What happens if a 200 year old person is on the Titanic? They would have a \(-.1\) probability of surviving. You cannot have a probability outside of the \([0,1]\) interval! Admittedly this is a trivial example, but in other circumstances this can become a more realistic scenario.

Or what if we didn’t want to predict survival, but instead predict the port from which an individual departed (Cherbourg, Queenstown, or Southampton). We could try and code this as a numeric response variable:

But why not instead code it:

Or even:

There is no inherent ordering to this variable. Any claimed linear relationship between a predictor and port of embarkation is completely dependent on how we convert the classes to numeric values.

Logistic regression

Rather than modeling the response \(Y\) directly, logistic regression instead models the probability that \(Y\) belongs to a particular category. In our first Titanic example, the probability of survival can be written as:

\[P(\text{survival} = \text{Yes} | \text{age})\]

The values of \(P(\text{survival} = \text{Yes} | \text{age})\) (or simply \(p(\text{survival})\) will range between 0 and 1. Given that predicted probability, we could predict anyone with for whom \(p(\text{survival}) > .5\) will survive the sinking, and anyone else will die.³

We can estimate the logistic regression model using the glm function.

survive_age <- glm(Survived ~ Age, data = titanic, family = binomial)
summary(survive_age)

## 
## Call:
## glm(formula = Survived ~ Age, family = binomial, data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1488  -1.0361  -0.9544   1.3159   1.5908  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -0.05672    0.17358  -0.327   0.7438  
## Age         -0.01096    0.00533  -2.057   0.0397 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 964.52  on 713  degrees of freedom
## Residual deviance: 960.23  on 712  degrees of freedom
##   (177 observations deleted due to missingness)
## AIC: 964.23
## 
## Number of Fisher Scoring iterations: 4

Which produces a line that looks like this:

ggplot(titanic, aes(Age, Survived)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"),
              se = FALSE)

It’s hard to tell, but the line is not perfectly linear. Let’s expand the range of the x-axis to prove this:

ggplot(titanic, aes(Age, Survived)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"),
              se = FALSE, fullrange = TRUE) +
  xlim(0,200)

No more predictions that a 200 year old has a \(-.1\) probability of surviving!

library(modelr)

titanic_age <- titanic %>%
  data_grid(Age)
titanic_age

## # A tibble: 88 × 1
##      Age
##    <dbl>
## 1   0.42
## 2   0.67
## 3   0.75
## 4   0.83
## 5   0.92
## 6   1.00
## 7   2.00
## 8   3.00
## 9   4.00
## 10  5.00
## # ... with 78 more rows

logit2prob <- function(x){
  exp(x) / (1 + exp(x))
}

library(modelr)

titanic_age <- titanic_age %>%
  add_predictions(survive_age) %>%
  mutate(pred = logit2prob(pred))
titanic_age

## # A tibble: 88 × 2
##      Age      pred
##    <dbl>     <dbl>
## 1   0.42 0.4846727
## 2   0.67 0.4839882
## 3   0.75 0.4837691
## 4   0.83 0.4835501
## 5   0.92 0.4833037
## 6   1.00 0.4830847
## 7   2.00 0.4803475
## 8   3.00 0.4776115
## 9   4.00 0.4748768
## 10  5.00 0.4721436
## # ... with 78 more rows

ggplot(titanic_age, aes(Age, pred)) +
  geom_line() +
  labs(title = "Relationship Between Age and Surviving the Titanic",
       y = "Predicted Probability of Survival")

Multiple predictors

But as the old principle of the sea goes, “women and children first”. What if age isn’t the only factor effecting survival? Fortunately logistic regression handles multiple predictors:

survive_age_woman <- glm(Survived ~ Age + Sex, data = titanic,
                         family = binomial)
summary(survive_age_woman)

## 
## Call:
## glm(formula = Survived ~ Age + Sex, family = binomial, data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7405  -0.6885  -0.6558   0.7533   1.8989  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.277273   0.230169   5.549 2.87e-08 ***
## Age         -0.005426   0.006310  -0.860     0.39    
## Sexmale     -2.465920   0.185384 -13.302  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 964.52  on 713  degrees of freedom
## Residual deviance: 749.96  on 711  degrees of freedom
##   (177 observations deleted due to missingness)
## AIC: 755.96
## 
## Number of Fisher Scoring iterations: 4

The coefficients essentially tell us the relationship between each individual predictor and the response, independent of other predictors. So this model tells us the relationship between age and survival, after controlling for the effects of gender. Likewise, it also tells us the relationship between gender and survival, after controlling for the effects of age. To get a better visualization of this, let’s use data_grid and add_predictions again:

titanic_age_sex <- titanic %>%
  data_grid(Age, Sex) %>%
  add_predictions(survive_age_woman) %>%
  mutate(pred = logit2prob(pred))
titanic_age_sex

## # A tibble: 176 × 3
##      Age    Sex      pred
##    <dbl>  <chr>     <dbl>
## 1   0.42 female 0.7815965
## 2   0.42   male 0.2330934
## 3   0.67 female 0.7813649
## 4   0.67   male 0.2328510
## 5   0.75 female 0.7812907
## 6   0.75   male 0.2327735
## 7   0.83 female 0.7812165
## 8   0.83   male 0.2326960
## 9   0.92 female 0.7811330
## 10  0.92   male 0.2326088
## # ... with 166 more rows

With these predicted probabilities, we can now plot the separate effects of age and gender:

ggplot(titanic_age_sex, aes(Age, pred, color = Sex)) +
  geom_line() +
  labs(title = "Probability of Surviving the Titanic",
       y = "Predicted Probability of Survival",
       color = "Sex")

This graph illustrates a key fact about surviving the sinking of the Titanic - age was not really a dominant factor. Instead, one’s gender was much more important. Females survived at much higher rates than males, regardless of age.

Quadratic terms

Logistic regression, like linear regression, assumes each predictor has an independent and linear relationship with the response. That is, it assumes the relationship takes the form \(y = \beta_0 + \beta_{1}x\) looks something like this:

sim_line <- tibble(x = runif(1000),
                   y = x * 1)

ggplot(sim_line, aes(x, y)) +
  geom_line()

But from algebra we know that variables can have non-linear relationships. Perhaps instead the relationship is parabolic like \(y = \beta_0 + \beta_{1}x + \beta_{2}x^2\):

sim_line <- tibble(x = runif(1000, -1, 1),
                   y = x^2 + x)

ggplot(sim_line, aes(x, y)) +
  geom_line()

Or a more general quadratic equation \(y = \beta_0 + \beta_{1}x + \beta_{2}x^2 + \beta_{3}x^3\):

sim_line <- tibble(x = runif(1000, -1, 1),
                   y = x^3 + x^2 + x)

ggplot(sim_line, aes(x, y)) +
  geom_line()

These can be accounted for in a logistic regression:

survive_age_square <- glm(Survived ~ Age + I(Age^2), data = titanic,
                          family = binomial)
summary(survive_age_square)

## 
## Call:
## glm(formula = Survived ~ Age + I(Age^2), family = binomial, data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2777  -1.0144  -0.9516   1.3421   1.4278  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.2688449  0.2722529   0.987   0.3234  
## Age         -0.0365193  0.0172749  -2.114   0.0345 *
## I(Age^2)     0.0003965  0.0002538   1.562   0.1183  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 964.52  on 713  degrees of freedom
## Residual deviance: 957.81  on 711  degrees of freedom
##   (177 observations deleted due to missingness)
## AIC: 963.81
## 
## Number of Fisher Scoring iterations: 4

titanic_age %>%
  add_predictions(survive_age) %>%
  mutate(pred = logit2prob(pred)) %>%
  ggplot(aes(Age, pred)) +
  geom_line() +
  labs(title = "Relationship Between Age and Surviving the Titanic",
       y = "Predicted Probability of Survival")

Here the parabolic term (\(\text{age}^2\)) is not actually meaningful. In other contexts it may be so.

Interactive terms

Another assumption of linear and logistic regression is that the relationships between predictors and responses are independent from one another. So for the age and gender example, we assume our function \(f\) looks something like:⁵

\[f = \beta_{0} + \beta_{1}\text{age} + \beta_{2}\text{gender}\]

However once again, that is an assumption. What if the relationship between age and the probability of survival is actually dependent on whether or not the individual is a female? This possibility would take the functional form:

\[f = \beta_{0} + \beta_{1}\text{age} + \beta_{2}\text{gender} + \beta_{3}(\text{age} \times \text{gender})\]

This is considered an interaction between age and gender. To estimate this in R, we simply specify Age * Sex in our formula for the glm function:⁶

survive_age_woman_x <- glm(Survived ~ Age * Sex, data = titanic,
                           family = binomial)
summary(survive_age_woman_x)

## 
## Call:
## glm(formula = Survived ~ Age * Sex, family = binomial, data = titanic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9401  -0.7136  -0.5883   0.7626   2.2455  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.59380    0.31032   1.913  0.05569 . 
## Age          0.01970    0.01057   1.863  0.06240 . 
## Sexmale     -1.31775    0.40842  -3.226  0.00125 **
## Age:Sexmale -0.04112    0.01355  -3.034  0.00241 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 964.52  on 713  degrees of freedom
## Residual deviance: 740.40  on 710  degrees of freedom
##   (177 observations deleted due to missingness)
## AIC: 748.4
## 
## Number of Fisher Scoring iterations: 4

As before, let’s estimate predicted probabilities and plot the interactive effects of age and gender.

titanic_age_sex_x <- titanic %>%
  data_grid(Age, Sex) %>%
  add_predictions(survive_age_woman_x) %>%
  mutate(pred = logit2prob(pred))
titanic_age_sex_x

## # A tibble: 176 × 3
##      Age    Sex      pred
##    <dbl>  <chr>     <dbl>
## 1   0.42 female 0.6461311
## 2   0.42   male 0.3245488
## 3   0.67 female 0.6472564
## 4   0.67   male 0.3233762
## 5   0.75 female 0.6476162
## 6   0.75   male 0.3230014
## 7   0.83 female 0.6479758
## 8   0.83   male 0.3226269
## 9   0.92 female 0.6483802
## 10  0.92   male 0.3222058
## # ... with 166 more rows

ggplot(titanic_age_sex_x, aes(Age, pred, color = Sex)) +
  geom_line() +
  labs(title = "Probability of Surviving the Titanic",
       y = "Predicted Probability of Survival",
       color = "Sex")

And now our minds are blown once again! For women, as age increases the probability of survival also increases. However for men, we see the opposite relationship: as age increases the probability of survival decreases. Again, the basic principle of saving women and children first can be seen empirically in the estimated probability of survival. Male children are treated similarly to female children, and their survival is prioritized. Even still, the odds of survival are always worse for men than women, but the regression function clearly shows a difference from our previous results.

You may think then that it makes sense to throw in interaction terms (and quadratic terms) willy-nilly to all your regression models since we never know for sure if the relationship is strictly linear and independent. You could do that, but once you start adding more predictors (3, 4, 5, etc.) that will get very difficult to keep track of (five-way interactions are extremely difficult to interpret - even three-way get to be problematic). The best advice is to use theory and your domain knowledge as your guide. Do you have a reason to believe the relationship should be interactive? If so, test for it. If not, don’t.

Evaluating model accuracy

age_accuracy <- titanic %>%
  add_predictions(survive_age) %>%
  mutate(pred = logit2prob(pred),
         pred = as.numeric(pred > .5))

mean(age_accuracy$Survived == age_accuracy$pred, na.rm = TRUE)

## [1] 0.5938375

x_accuracy <- titanic %>%
  add_predictions(survive_age_woman_x) %>%
  mutate(pred = logit2prob(pred),
         pred = as.numeric(pred > .5))

mean(x_accuracy$Survived == x_accuracy$pred, na.rm = TRUE)

## [1] 0.780112

titanic_split <- resample_partition(titanic, c(test = 0.3, train = 0.7))
map(titanic_split, dim)

## $test
## [1] 267  12
## 
## $train
## [1] 624  12

train_model <- glm(Survived ~ Age * Sex, data = titanic_split$train,
                   family = binomial)
summary(train_model)

## 
## Call:
## glm(formula = Survived ~ Age * Sex, family = binomial, data = titanic_split$train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7541  -0.6539  -0.5867   0.7727   2.1018  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.74926    0.36231   2.068 0.038638 *  
## Age          0.01095    0.01194   0.917 0.359097    
## Sexmale     -1.93347    0.50163  -3.854 0.000116 ***
## Age:Sexmale -0.02230    0.01586  -1.406 0.159580    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 674.37  on 502  degrees of freedom
## Residual deviance: 508.58  on 499  degrees of freedom
##   (121 observations deleted due to missingness)
## AIC: 516.58
## 
## Number of Fisher Scoring iterations: 4

x_test_accuracy <- titanic_split$test %>%
  tbl_df() %>%
  add_predictions(train_model) %>%
  mutate(pred = logit2prob(pred),
         pred = as.numeric(pred > .5))

mean(x_test_accuracy$Survived == x_test_accuracy$pred, na.rm = TRUE)

## [1] 0.7535545