Logistic Regression

STAT 20: Introduction to Probability and Statistics

Concept Questions

A logistic regression model was fit in an attempt to predict the sex of a penguin "male" or "female" based on its body mass (grams).

Assuming that no change to the penguins dataset was made, will the model be predicting the probability of the penguin being male or the probability of the penguin being female?

01:00

m1 <- glm(sex ~ body_mass_g, data = penguins, family = "binomial")

 (Intercept)  body_mass_g 
-5.162541644  0.001239819 

Which of the expressions given in the poll (math or code) will correctly calculate the predicted probability that a penguin that weighs 4000 g is a female? Select all that apply

01:00

To answer this one correctly, they’ll need to know (or remember from the RQ) that the reference level is female, so this model is predicting the probability male.

You can show that they can either use R as a calculator or use the predict function.

1/(1 + exp(-(-5.15 + .00124 * 4000))) predict(m1, newdata = data.frame(body_mass_g = 4000), type = “response”)

The latter is more precise, which is useful here since p-hat (prob of male) is .449 therefore prob of female is .55.

As a bonus, try sketching this function on a scatterplot!

When sketching this, the positive slope means the s-curve goes up and to the right. The negative intercept shifts the whole curve a tiny bit to the left, reflecting the base rate of a couple more males than females (168 vs 165).

What is the misclassification rate of this model?

# A tibble: 4 × 3
# Groups:   sex [2]
  sex    y_hat      n
  <fct>  <chr>  <int>
1 female female   109
2 female male      56
3 male   female    74
4 male   male      94

01:00

(56+74)/(109+56+74+94) = 130/333

m2 <- glm(sex ~ body_mass_g + bill_length_mm, 
          data = penguins, family = "binomial")

   (Intercept)    body_mass_g bill_length_mm 
   -6.91208086     0.00101530     0.06112808 

Open up RStudio and fit the model here in the slides. What are the predicted sexes of these two penguins?

1. body mass = 3900 g, bill length = 50
2. body mass = 4100 g, bill length = 35

01:00

You can have students work on their laptops for this question.

This one extends the task of the previous one by having them include another covariate and also to do the thresholding procedure to go from p-hat to y-hat. This doesn’t specify the threshold, so using .5 is probably a good default.

predict(m2, newdata = data.frame(body_mass_g = 3900, bill_length_mm = 50), type = “response”) # male predict(m2, newdata = data.frame(body_mass_g = 4100, bill_length_mm = 35), type = “response”) # female

Or they can do it “by hand”, in which case, the first one uses \(z = -6.91208086 + 0.00101530 \times 3900 + 0.06112808\times 50 = 0.1039931\) so \(\hat{p} = \frac{1}{1+exp(-0.1039931)} = 0.5259749\) and we would predict male. For the second, we have that \(z = -6.91208086 + 0.00101530 \times 4100 + 0.06112808 \times 35 = -0.6098681\). Plugging this into the logistic function gives \(\hat{p} = \frac{1}{1+exp(0.6098681)} = 0.3520893\), so we would predict female.

Misclassification: Appendix

Building a predictive model

1. Decide on the mathematical form of the model: logistic linear regression

2. Select a metric that defines the “best” fit: the coefficients in logistic regression are the ones that minimize not the RSS function but a function called log-loss (which we don’t have time to cover)

3. Estimating the coefficients of the model that are best using the training data: we know how to do this: test + train + glm()!

4. Evaluating predictive accuracy using a test data set:\(R^2\) isn’t relevant here. We need a new metric!

Example: penguins

set.seed(132)

# randomly sample train/test set split
set_type <- sample(x = c('train', 'test'), 
                   size = nrow(penguins), 
                   replace = TRUE, 
                   prob = c(0.8, 0.2))

Example: penguins

set.seed(132)

# randomly sample train/test set split
set_type <- sample(x = c('train', 'test'), 
                   size = nrow(penguins), 
                   replace = TRUE, 
                   prob = c(0.8, 0.2))

train <- penguins |>
  filter(set_type == "train")

test <- penguins |>
  filter(set_type == "test")

Predicting into test set

m2 <- glm(sex ~ body_mass_g + bill_length_mm,
          data = train, family = "binomial")
p_hat <- predict(m2, test, type = "response")

Predicting into test set

m2 <- glm(sex ~ body_mass_g + bill_length_mm,
          data = train, family = "binomial")
p_hat <- predict(m2, test, type = "response")

test |>
  select(sex)

# A tibble: 70 × 1
   sex   
   <fct> 
 1 female
 2 male  
 3 female
 4 male  
 5 male  
 6 female
 7 male  
 8 female
 9 male  
10 male  
# ℹ 60 more rows

Predicting into test set

m2 <- glm(sex ~ body_mass_g + bill_length_mm,
          data = train, family = "binomial")
p_hat <- predict(m2, test, type = "response")

test |>
  select(sex) |>
  mutate(p_hat = p_hat)

# A tibble: 70 × 2
   sex    p_hat
   <fct>  <dbl>
 1 female 0.345
 2 male   0.566
 3 female 0.259
 4 male   0.280
 5 male   0.365
 6 female 0.196
 7 male   0.428
 8 female 0.220
 9 male   0.559
10 male   0.279
# ℹ 60 more rows

Predicting into test set

m2 <- glm(sex ~ body_mass_g + bill_length_mm,
          data = train, family = "binomial")

test |>
  select(sex) |>
  mutate(p_hat = predict(m2, test, type = "response"),
         y_hat = ifelse(p_hat > .5, "male", "female"))

# A tibble: 70 × 3
   sex    p_hat y_hat 
   <fct>  <dbl> <chr> 
 1 female 0.345 female
 2 male   0.566 male  
 3 female 0.259 female
 4 male   0.280 female
 5 male   0.365 female
 6 female 0.196 female
 7 male   0.428 female
 8 female 0.220 female
 9 male   0.559 male  
10 male   0.279 female
# ℹ 60 more rows

Classification errors

test |>
  select(sex) |>
  mutate(p_hat = p_hat,
         y_hat = ifelse(p_hat > .5, "male", "female"),
         FP = sex == "female" & y_hat == "male",
         FN = sex == "male" & y_hat == "female")

# A tibble: 70 × 5
   sex    p_hat y_hat  FP    FN   
   <fct>  <dbl> <chr>  <lgl> <lgl>
 1 female 0.345 female FALSE FALSE
 2 male   0.566 male   FALSE FALSE
 3 female 0.259 female FALSE FALSE
 4 male   0.280 female FALSE TRUE 
 5 male   0.365 female FALSE TRUE 
 6 female 0.196 female FALSE FALSE
 7 male   0.428 female FALSE TRUE 
 8 female 0.220 female FALSE FALSE
 9 male   0.559 male   FALSE FALSE
10 male   0.279 female FALSE TRUE 
# ℹ 60 more rows

Misclassification Rate

test |>
  select(sex) |>
  mutate(p_hat = p_hat,
         y_hat = ifelse(p_hat > .5, "male", "female")) |>
  summarize(MCR = mean(sex != y_hat))

# A tibble: 1 × 1
    MCR
  <dbl>
1 0.371