Classification problems

The sinking of RMS Titanic provided the world with many things:

• A fundamental shock to the world as its faith in supposedly indestructible technology was shattered by a chunk of ice
• Perhaps the best romantic ballad of all time

• A tragic love story

  

  Titanic (1997)

Why did Jack have to die? Why couldn’t he have made it onto a lifeboat like Cal? We may never know the answer, but we can generalize the question a bit: why did some people survive the sinking of the Titanic while others did not?

In essence, we have a classification problem. The response is a binary variable, indicating whether a specific passenger survived. If we combine this with predictors that describe each passenger, we might be able to estimate a general model of survival.¹

Kaggle is an online platform for predictive modeling and analytics. They run regular competitions where they provide the public with a question and data, and anyone can estimate a predictive model to answer the question. They’ve run a popular contest based on a dataset of passengers from the Titanic. The datasets have been conveniently stored in a package called titanic. Let’s load the package and convert the desired data frame to a tibble.

library(titanic)
titanic <- titanic_train %>%
  as_tibble()

titanic %>%
  head() %>%
  knitr::kable()

The codebook contains the following information on the variables:

VARIABLE DESCRIPTIONS:
Survived        Survival
                (0 = No; 1 = Yes)
Pclass          Passenger Class
                (1 = 1st; 2 = 2nd; 3 = 3rd)
Name            Name
Sex             Sex
Age             Age
SibSp           Number of Siblings/Spouses Aboard
Parch           Number of Parents/Children Aboard
Ticket          Ticket Number
Fare            Passenger Fare
Cabin           Cabin
Embarked        Port of Embarkation
                (C = Cherbourg; Q = Queenstown; S = Southampton)

SPECIAL NOTES:
Pclass is a proxy for socio-economic status (SES)
 1st ~ Upper; 2nd ~ Middle; 3rd ~ Lower

Age is in Years; Fractional if Age less than One (1)
 If the Age is Estimated, it is in the form xx.5

With respect to the family relation variables (i.e. sibsp and parch)
some relations were ignored.  The following are the definitions used
for sibsp and parch.

Sibling:  Brother, Sister, Stepbrother, or Stepsister of Passenger Aboard Titanic
Spouse:   Husband or Wife of Passenger Aboard Titanic (Mistresses and Fiances Ignored)
Parent:   Mother or Father of Passenger Aboard Titanic
Child:    Son, Daughter, Stepson, or Stepdaughter of Passenger Aboard Titanic

Other family relatives excluded from this study include cousins,
nephews/nieces, aunts/uncles, and in-laws.  Some children travelled
only with a nanny, therefore parch=0 for them.  As well, some
travelled with very close friends or neighbors in a village, however,
the definitions do not support such relations.

So if this is our data, Survived is our response variable and the remaining variables are predictors. How can we determine who survives and who dies?

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

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

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!

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

## # A tibble: 88 x 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 x 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")

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

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

## # A tibble: 176 x 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

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")

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

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

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

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

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

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

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")

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

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 x 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")

Comparing models

How do we know if a logistic regression model is good or bad? One evalation criteria simply asks: how many prediction errors did the model make? For instance, how often did our basic model just using age perform? First we need to get the predicted probabilities for each individual in the original dataset, convert the probability to a prediction (I use a \(.5\) cut-point), then see what percentage of predictions were the same as the actual survivals?

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.4061625

\(40.6\%\) of the predictions based on age were incorrect. If we flipped a coin to make our predictions, we’d be right about 50% of the time. So this is a decent improvement. Of course, we also know that \(61.6\%\) of passengers died in the sinking, so if we just guessed that every passenger died we’d still make fewer mistakes than our predictive model. Maybe this model isn’t so great after all. What about our interactive age and gender model?

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.219888

This model is much better. Just by knowing an individual’s age and gender, our model is incorrect only 22% of the time.

Exercise: logistic regression with mental_health

Why do some people vote in elections while others do not? Typical explanations focus on a resource model of participation – individuals with greater resources, such as time, money, and civic skills, are more likely to participate in politics. An emerging theory assesses an individual’s mental health and its effect on political participation.⁶ Depression increases individuals’ feelings of hopelessness and political efficacy, so depressed individuals will have less desire to participate in politics. More importantly to our resource model of participation, individuals with depression suffer physical ailments such as a lack of energy, headaches, and muscle soreness which drain an individual’s energy and requires time and money to receive treatment. For these reasons, we should expect that individuals with depression are less likely to participate in election than those without symptoms of depression.

Use the mental_health data set in library(rcfss) and logistic regression to predict whether or not an individual voted in the 1996 presidental election.

library(rcfss)
mental_health

## # A tibble: 1,317 x 5
##    vote96   age  educ female mhealth
##     <dbl> <dbl> <dbl>  <dbl>   <dbl>
##  1      1    60    12      0       0
##  2      1    36    12      0       1
##  3      0    21    13      0       7
##  4      0    29    13      0       6
##  5      1    39    18      1       2
##  6      1    41    15      1       1
##  7      1    48    20      0       2
##  8      0    20    12      1       9
##  9      0    27    11      1       9
## 10      0    34     7      1       2
## # ... with 1,307 more rows

1. Estimate a logistic regression model of voter turnout with mhealth as the predictor. Estimate predicted probabilities and plot the logistic regression line using ggplot.

   Click for the solution

   # estimate model
   mh_model <- glm(vote96 ~ mhealth, data = mental_health,
                   family = binomial)
   tidy(mh_model)

   ##          term   estimate  std.error statistic      p.value
   ## 1 (Intercept)  1.2225348 0.08901013 13.734782 6.284476e-43
   ## 2     mhealth -0.1766612 0.02217096 -7.968135 1.610862e-15

   # estimate predicted probabilities
   mh_health <- mental_health %>%
     data_grid(mhealth) %>%
     add_predictions(mh_model) %>%
     mutate(pred = logit2prob(pred))
   mh_health

   ## # A tibble: 10 x 2
   ##    mhealth      pred
   ##      <dbl>     <dbl>
   ##  1       0 0.7725093
   ##  2       1 0.7399817
   ##  3       2 0.7045818
   ##  4       3 0.6665342
   ##  5       4 0.6261862
   ##  6       5 0.5840032
   ##  7       6 0.5405526
   ##  8       7 0.4964767
   ##  9       8 0.4524554
   ## 10       9 0.4091655

   # graph the line
   ggplot(mh_health, aes(mhealth, pred)) +
     geom_line() +
     labs(title = "Relationship Between Mental Health and Voter Turnout",
          y = "Predicted Probability of Voting") +
     scale_y_continuous(limits = c(0, 1))

   

2. Calculate the error rate of the model.

   Click for the solution

   mh_model_accuracy <- mental_health %>%
     as_tibble() %>%
     add_predictions(mh_model) %>%
     mutate(pred = logit2prob(pred),
            pred = as.numeric(pred > .5))
   
   (mh_model_err <- mean(mh_model_accuracy$vote96 != mh_model_accuracy$pred,
                         na.rm = TRUE))

   ## [1] 0.317388

3. Estimate a second logistic regression model of voter turnout using the training set and all the predictors. Calculate it’s error rate using the test set, and compare it to the original model. Which is better?

   Click for the solution

   # estimate model
   mh_model_all <- glm(vote96 ~ ., data = mental_health,
                   family = binomial)
   tidy(mh_model_all)

   ##          term    estimate   std.error  statistic      p.value
   ## 1 (Intercept) -4.25770320 0.469355720 -9.0713781 1.175211e-19
   ## 2         age  0.04460727 0.004459019 10.0038302 1.466142e-23
   ## 3        educ  0.25764980 0.026574296  9.6954513 3.152409e-22
   ## 4      female -0.03882557 0.130087258 -0.2984579 7.653537e-01
   ## 5     mhealth -0.11836876 0.024149285 -4.9015430 9.508685e-07

   # calculate error rate
   mh_model_all_accuracy <- mental_health %>%
     as_tibble() %>%
     add_predictions(mh_model_all) %>%
     mutate(pred = logit2prob(pred),
            pred = as.numeric(pred > .5))
   
   (mh_model_all_err <- mean(mh_model_all_accuracy$vote96 != mh_model_all_accuracy$pred,
                             na.rm = TRUE))

   ## [1] 0.2786636

   The model with all predictors has a 3.87% lower error rate than predictions based only on mental health.