9 Logistic Regression*
This chapter introduces the popular classification method logistic regression.
Packages Used in this Chapter
This chapter only uses R’s base functionality and does not need extra packages.
9.1 Introduction
Logistic regression contains the word regression, but it is actually a probabilistic, statistical classification model to predict the probability \(p\) of a binary outcome given a set of features. It is a very powerful model that can be fit very quickly. It is one of the first classification models you should try on new data.
Logistic regression is a generalized linear model with the logit as the link function and a binomial error model. It can be thought of as a linear regression with the log odds ratio (logit) of the binary outcome as the dependent variable:
\[logit(p) = ln\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ...\]
The logit function links the probability p to the linear regression by converting the logit into a number in the probability range \([0,1]\).
logit <- function(p) log(p/(1-p))
p <- seq(0, 1, length.out = 100)
plot(logit(p), p, type = "l")
abline(h = 0.5, lty = 2)
abline(v = 0, lty = 2)
This is equivalent to modeling the probability of the outcome \(p\) by:
\[ p = \frac{e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ...}}{1 + e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ...}} = \frac{1}{1+e^{-(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ...)}}\]
9.2 Data Preparation
We load and shuffle the data. We also add a useless variable to see if the logistic regression removes it.
data(iris)
set.seed(100) # for reproducability
x <- iris[sample(1:nrow(iris)),]
x <- cbind(x, useless = rnorm(nrow(x)))We create a binary classification problem by classify if a flower is of species Virginica.
x$virginica <- x$Species == "virginica"
x$Species <- NULL
plot(x, col=x$virginica + 1)
9.3 A first Logistic Regression Model
Logistic regression is a generalized linear model (GLM) with logit as the link function and a binomial error model.
model <- glm(virginica ~ .,
family = binomial(logit), data = x)
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurredAbout the warning: glm.fit: fitted probabilities numerically 0 or 1 occurred means that the data is possibly linearly separable.
model
##
## Call: glm(formula = virginica ~ ., family = binomial(logit), data = x)
##
## Coefficients:
## (Intercept) Sepal.Length Sepal.Width Petal.Length
## -41.649 -2.531 -6.448 9.376
## Petal.Width useless
## 17.696 0.098
##
## Degrees of Freedom: 149 Total (i.e. Null); 144 Residual
## Null Deviance: 191
## Residual Deviance: 11.9 AIC: 23.9Check which features are significant?
summary(model)
##
## Call:
## glm(formula = virginica ~ ., family = binomial(logit), data = x)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -41.649 26.556 -1.57 0.117
## Sepal.Length -2.531 2.458 -1.03 0.303
## Sepal.Width -6.448 4.794 -1.34 0.179
## Petal.Length 9.376 4.763 1.97 0.049 *
## Petal.Width 17.696 10.632 1.66 0.096 .
## useless 0.098 0.807 0.12 0.903
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 190.954 on 149 degrees of freedom
## Residual deviance: 11.884 on 144 degrees of freedom
## AIC: 23.88
##
## Number of Fisher Scoring iterations: 12AIC (Akaike information criterion) is a measure of how good the model is. Smaller is better. It can be used for model selection.
9.4 Stepwise Variable Selection
model2 <- step(model, data = x)
## Start: AIC=23.88
## virginica ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width +
## useless
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Df Deviance AIC
## - useless 1 11.9 21.9
## - Sepal.Length 1 13.2 23.2
## <none> 11.9 23.9
## - Sepal.Width 1 14.8 24.8
## - Petal.Width 1 22.4 32.4
## - Petal.Length 1 25.9 35.9
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
##
## Step: AIC=21.9
## virginica ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Df Deviance AIC
## - Sepal.Length 1 13.3 21.3
## <none> 11.9 21.9
## - Sepal.Width 1 15.5 23.5
## - Petal.Width 1 23.8 31.8
## - Petal.Length 1 25.9 33.9
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
##
## Step: AIC=21.27
## virginica ~ Sepal.Width + Petal.Length + Petal.Width
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1
## occurred
## Df Deviance AIC
## <none> 13.3 21.3
## - Sepal.Width 1 20.6 26.6
## - Petal.Length 1 27.4 33.4
## - Petal.Width 1 31.5 37.5
summary(model2)
##
## Call:
## glm(formula = virginica ~ Sepal.Width + Petal.Length + Petal.Width,
## family = binomial(logit), data = x)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -50.53 23.99 -2.11 0.035 *
## Sepal.Width -8.38 4.76 -1.76 0.079 .
## Petal.Length 7.87 3.84 2.05 0.040 *
## Petal.Width 21.43 10.71 2.00 0.045 *
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 190.954 on 149 degrees of freedom
## Residual deviance: 13.266 on 146 degrees of freedom
## AIC: 21.27
##
## Number of Fisher Scoring iterations: 12The estimates (\(\beta_0, \beta_1,...\) ) are log-odds and can be converted into odds using \(exp(\beta)\). A negative log-odds ratio means that the odds go down with an increase in the value of the predictor. A predictor with a positive log-odds ratio increases the odds. In this case, the odds of looking at a Virginica iris goes down with Sepal.Width and increases with the other two predictors.
9.5 Calculate the Response
Note: we do here in-sample testing on the data we learned the data from. To get a generalization error estimate you should use a test set or cross-validation!
pr <- predict(model2, x, type = "response")
round(pr[1:10], 2)
## 102 112 4 55 70 98 135 7 43 140
## 1.00 1.00 0.00 0.00 0.00 0.00 0.86 0.00 0.00 1.00The response is the predicted probability of the flower being of species virginica. The probabilities of the first 10 flowers are shown. Below is a histogram of predicted probabilities. The color is used to show the examples that have the true class virginica.
hist(pr, breaks = 20, main = "Predicted Probability vs. True Class")
hist(pr[x$virginica == TRUE], col = "red", breaks = 20, add = TRUE)
9.6 Check Classification Performance
We calculate the predicted class by checking if the probability is larger than .5.
pred <- pr > .5Now we can create a confusion table and calculate the accuracy.
tbl <- table(predicted = pred, actual = x$virginica)
tbl
## actual
## predicted FALSE TRUE
## FALSE 98 1
## TRUE 2 49
sum(diag(tbl))/sum(tbl)
## [1] 0.98We can also use caret’s more advanced function caret::confusionMatrix(). Our code
above uses logical vectors. For caret, we need to make sure that both,
the reference and the predictions are coded as factor.
caret::confusionMatrix(
reference = factor(x$virginica, labels = c("Yes", "No"), levels = c(TRUE, FALSE)),
data = factor(pred, labels = c("Yes", "No"), levels = c(TRUE, FALSE)))
## Confusion Matrix and Statistics
##
## Reference
## Prediction Yes No
## Yes 49 2
## No 1 98
##
## Accuracy : 0.98
## 95% CI : (0.943, 0.996)
## No Information Rate : 0.667
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.955
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.980
## Specificity : 0.980
## Pos Pred Value : 0.961
## Neg Pred Value : 0.990
## Prevalence : 0.333
## Detection Rate : 0.327
## Detection Prevalence : 0.340
## Balanced Accuracy : 0.980
##
## 'Positive' Class : Yes
## We see that the model performs well with a very high accuracy and kappa value.
9.7 Regularized Logistic Regression
Glmnet fits generalized linear models (including logistic regression)
using regularization via penalized maximum likelihood.
The regularization parameter \(\lambda\) is a hyperparameter and
glmnet can use cross-validation to find an appropriate
value. glmnet does not have a function interface, so we have
to supply a matrix for X and a vector of responses for y.
library(glmnet)
## Loaded glmnet 4.1-8
X <- as.matrix(x[, 1:5])
y <- x$virginica
fit <- cv.glmnet(X, y, family = "binomial")
fit
##
## Call: cv.glmnet(x = X, y = y, family = "binomial")
##
## Measure: Binomial Deviance
##
## Lambda Index Measure SE Nonzero
## min 0.00164 59 0.126 0.0456 5
## 1se 0.00664 44 0.167 0.0422 3There are several selection rules for lambda, we look at the coefficients of the logistic regression using the lambda that gives the most regularized model such that the cross-validated error is within one standard error of the minimum cross-validated error.
coef(fit, s = fit$lambda.1se)
## 6 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -16.961
## Sepal.Length .
## Sepal.Width -1.766
## Petal.Length 2.197
## Petal.Width 6.820
## useless .A dot means 0. We see that the predictors Sepal.Length and useless are not used in the prediction giving a models similar to stepwise variable selection above.
A predict function is provided. We need to specify what regularization to use and that we want to predict a class label.
predict(fit, newx = X[1:5,], s = fit$lambda.1se, type = "class")
## s1
## 102 "TRUE"
## 112 "TRUE"
## 4 "FALSE"
## 55 "FALSE"
## 70 "FALSE"Glmnet provides supports many types of generalized linear models. Examples can be found in the article An Introduction to glmnet.
9.8 Exercises
We will again use the Palmer penguin data for the exercises.
library(palmerpenguins)
head(penguins)
## # A tibble: 6 × 8
## species island bill_length_mm bill_depth_mm
## <chr> <chr> <dbl> <dbl>
## 1 Adelie Torgersen 39.1 18.7
## 2 Adelie Torgersen 39.5 17.4
## 3 Adelie Torgersen 40.3 18
## 4 Adelie Torgersen NA NA
## 5 Adelie Torgersen 36.7 19.3
## 6 Adelie Torgersen 39.3 20.6
## # ℹ 4 more variables: flipper_length_mm <dbl>,
## # body_mass_g <dbl>, sex <chr>, year <dbl>Create an R markdown document that performs the following:
- Create a test and a training data set (see section Holdout Method in Chapter 3).
- Create a logistic regression using the training set to predict the variable sex.
- Use stepwise variable selection. What variables are selected?
- What do the parameters for for each of the selected features tell you?
- Predict the sex of the penguins in the test set. Create a confusion table and calculate the accuracy and discuss how well the model works.