Logistic Regression
Contents
- What is Logistic Regression
- Math
- Implementation
- Dataset
- Modeling
- Visualization
- Basic Evaluation
- Optimization
- Evaluation
- ROC Curve
- Area under ROC Curve
- References
What is Logistic Regression
Logistic regression is a type of linear regression. However, it is used for classification only. Huh.. that’s confusing, right ? Let’s dive in.
Let’s take the simple iris dataset. The target variable as you know by now ( from day 9 – Introduction to Classification in Python, where we discussed classification using K Nearest neighbours ) is categorical in nature. Let’s load the data first.
head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species <dbl> <dbl> <dbl> <dbl> <fct> 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa
Let’s simplify the dataset to just 2 species –
- 0 – setosa
- 1 – versi-color
Let’s just take data for 2 of the species ( say setosa and versi-color ) with just the sepal data ( sepal length and sepal width ) and plot it.
iris_data = iris[0:100,]
iris_data = iris[0:100,]
plot(iris_data$Sepal.Length, iris_data$Sepal.Width,
col = iris_data$Species, pch = 19,
xlab = "Sepal Length",
ylab = "Sepal Width")
Let’s simplify this further – say, we wanted to predict the species based on a single parameter – Sepal Length. Let’s first plot it.
plot ( iris_data[,1], iris_data[,5],
pch = 19, col = "blue",
xlab = "Sepal Length",
ylab = "Setosa or Versicolor")
We know that regression is used to predict a continous variable. What about a categorical variable like this ? (species). If we can draw a curve like this,
and for all target values predicted with value > 0.5 put it in one category, and for all target values less than 0.5, put it in the other category – like this.
A linear regression (multilinear in this case) equation looks like this.
Logistic regression is almost similar to linear regression. The difference lies in how the predictor is calculated. Let’s see it in the next section.
Math
The name logistic regression is derived from the logit function. This function is based on odds.
logit function
Let’s take an example. A standard dice roll has 6 outcomes. So, what is the probability of landing a 4 ?
Now, what about odds ? The odds of landing a 4 is
So, when we substitute p into the odds equation, it becomes
OK. Now that we understand Probability and Odds, let’s get to the log of odds.
How exactly is the logistic regression similar to linear regression ? Like so.
Where the predictor ( log of odds ) varies between ( -∞ to +∞ ).
To understand this better, let’s plot the log of odds between a probabilty value of 0 and 1.
x = runif ( 100, min = 0.1, max = 0.9) x = sort ( x ) y = log ( x / (1-x))
# plt.plot(x,y) # plt.grid() plot ( x, y , type="l") grid(5, 5, lwd = 2)
This is the logistic regression curve. It maps a probability value ( 0 to 1 ) to a number ( -∞ to +∞ ). However, we are not looking for a continous variable, right ? The predictor we are looking for is a categorical variable – in our case, we said we would be able to predict this based on probability.
- p >= 0.5 – Category 1
- p < 0.5 – Category 2
In order to calculate those probabilities, we would have to calculate the inverse function of the logit function.
sigmoid function
The inverse of the logit curve is the inverse-logit or sigmoid function ( or expit function). The sigmoid function transforms the numbers ( -∞ to +∞ ) back to values between 0 and 1. Here is the formula for the sigmoid function.
x_new = y y_new = exp(x_new) / (1 + exp(x_new)) plot (x_new, y_new, type="l") grid(5, 5, lwd = 2)
Essentially, if we flip the logit function 900, you get the sigmoid function.
Here is the trick – As long as we are able to find a curve like the one below, although the target (predictor) is a value between 0 and 1 ( probabilities), we can say that all values below 0.5 ( half way mark ) belongs to one category and the remaining ( values above 0.5 ) belong to the next category. This is the essence of logistic regression.
Implementation
Let’s try to implement the logistic regression function in R step by step.
Data & Modeling
Just to keep the same example going, let’s try to fit the sepal length data to try and predict the species as either Setosa or Versicolor.
model = glm(Species ~ Sepal.Length + Sepal.Width,
data=iris_data,
family=binomial(link="logit"))
Warning message: "glm.fit: algorithm did not converge" Warning message: "glm.fit: fitted probabilities numerically 0 or 1 occurred"
y_pred = predict(model, iris_data, type="response")
y_pred
1 2.22044604925031e-16 2 2.22044604925031e-16 3 2.22044604925031e-16 4 2.22044604925031e-16 5 2.22044604925031e-16 6 2.22044604925031e-16 7 2.22044604925031e-16 8 2.22044604925031e-16 9 2.22044604925031e-16 10 2.22044604925031e-16 11 2.22044604925031e-16 12 2.22044604925031e-16 13 2.22044604925031e-16 14 2.22044604925031e-16 15 2.22044604925031e-16 16 2.22044604925031e-16 17 2.22044604925031e-16 18 2.22044604925031e-16 19 7.94264732123352e-13 20 2.22044604925031e-16 21 7.2365355113674e-11 22 2.22044604925031e-16 23 2.22044604925031e-16 24 2.22044604925031e-16 25 2.22044604925031e-16 26 2.22044604925031e-16 27 2.22044604925031e-16 28 2.22044604925031e-16 29 2.22044604925031e-16 30 2.22044604925031e-16 31 2.22044604925031e-16 32 7.2365355113674e-11 33 2.22044604925031e-16 34 2.22044604925031e-16 35 2.22044604925031e-16 36 2.22044604925031e-16 37 6.31794371395771e-10 38 2.22044604925031e-16 39 2.22044604925031e-16 40 2.22044604925031e-16 41 2.22044604925031e-16 42 9.0266451369787e-10 43 2.22044604925031e-16 44 2.22044604925031e-16 45 2.22044604925031e-16 46 2.22044604925031e-16 47 2.22044604925031e-16 48 2.22044604925031e-16 49 2.22044604925031e-16 50 2.22044604925031e-16 51 1 52 1 53 1 54 1 55 1 56 1 57 1 58 0.999999999144556 59 1 60 0.999999999998715 61 1 62 1 63 1 64 1 65 1 66 1 67 1 68 1 69 1 70 1 71 1 72 1 73 1 74 1 75 1 76 1 77 1 78 1 79 1 80 1 81 1 82 1 83 1 84 1 85 0.999999998977491 86 1 87 1 88 1 89 1 90 1 91 1 92 1 93 1 94 1 95 1 96 1 97 1 98 1 99 1 100 1
These are probability values. We need to convert them to the actual factors (setosa & versicolor), because, we are dealing with just 2 classes. We can just use a simple ifelse ( ) syntax to convert all values > 0.5 to a versicolor and < 0.5 to a setosa.
y_pred_levels = as.factor(ifelse(y_pred>0.5,"versicolor","setosa"))
y_pred_levels
1 setosa 2 setosa 3 setosa 4 setosa 5 setosa 6 setosa 7 setosa 8 setosa 9 setosa 10 setosa 11 setosa 12 setosa 13 setosa 14 setosa 15 setosa 16 setosa 17 setosa 18 setosa 19 setosa 20 setosa 21 setosa 22 setosa 23 setosa 24 setosa 25 setosa 26 setosa 27 setosa 28 setosa 29 setosa 30 setosa 31 setosa 32 setosa 33 setosa 34 setosa 35 setosa 36 setosa 37 setosa 38 setosa 39 setosa 40 setosa 41 setosa 42 setosa 43 setosa 44 setosa 45 setosa 46 setosa 47 setosa 48 setosa 49 setosa 50 setosa 51 versicolor 52 versicolor 53 versicolor 54 versicolor 55 versicolor 56 versicolor 57 versicolor 58 versicolor 59 versicolor 60 versicolor 61 versicolor 62 versicolor 63 versicolor 64 versicolor 65 versicolor 66 versicolor 67 versicolor 68 versicolor 69 versicolor 70 versicolor 71 versicolor 72 versicolor 73 versicolor 74 versicolor 75 versicolor 76 versicolor 77 versicolor 78 versicolor 79 versicolor 80 versicolor 81 versicolor 82 versicolor 83 versicolor 84 versicolor 85 versicolor 86 versicolor 87 versicolor 88 versicolor 89 versicolor 90 versicolor 91 versicolor 92 versicolor 93 versicolor 94 versicolor 95 versicolor 96 versicolor 97 versicolor 98 versicolor 99 versicolor 100 versicolor
library(caret) cm = confusionMatrix(y_pred_levels,iris_data[,5]) cm
Warning message in levels(reference) != levels(data): "longer object length is not a multiple of shorter object length" Warning message in confusionMatrix.default(y_pred_levels, iris_data[, 5]): "Levels are not in the same order for reference and data. Refactoring data to match."
Confusion Matrix and Statistics
Reference
Prediction setosa versicolor virginica
setosa 50 0 0
versicolor 0 50 0
virginica 0 0 0
Overall Statistics
Accuracy : 1
95% CI : (0.9638, 1)
No Information Rate : 0.5
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 1
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: setosa Class: versicolor Class: virginica
Sensitivity 1.0 1.0 NA
Specificity 1.0 1.0 1
Pos Pred Value 1.0 1.0 NA
Neg Pred Value 1.0 1.0 NA
Prevalence 0.5 0.5 0
Detection Rate 0.5 0.5 0
Detection Prevalence 0.5 0.5 0
Balanced Accuracy 1.0 1.0 NA
Basic Evaluation
Let’s split up the data into training and test data and model it. This time let’s do the full iris dataset. Since there are 3 Species to be predicted, we cannot use glm with a “binomial” family algorithm. Let’s use another library called nnet. As usual, to evaluate categorical target data, we use a confusion matrix.
library(nnet)
index = sample(1:nrow(iris),nrow(iris)*.8) train = iris[index,] test = iris[-index,]
model = multinom(Species~.,data = train)
# weights: 18 (10 variable) initial value 131.833475 iter 10 value 11.516467 iter 20 value 4.881298 iter 30 value 4.469920 iter 40 value 4.263054 iter 50 value 3.911756 iter 60 value 3.823284 iter 70 value 3.598069 iter 80 value 3.591202 iter 90 value 3.570975 iter 100 value 3.570835 final value 3.570835 stopped after 100 iterations
pred = predict(model,test)
As usual, to evaluate categorical target data, we use a confusion matrix.
library(caret) cm = confusionMatrix(pred, as.factor(test$Species)) cm
Confusion Matrix and Statistics
Reference
Prediction setosa versicolor virginica
setosa 16 0 0
versicolor 0 8 1
virginica 0 2 3
Overall Statistics
Accuracy : 0.9
95% CI : (0.7347, 0.9789)
No Information Rate : 0.5333
P-Value [Acc > NIR] : 1.989e-05
Kappa : 0.8315
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: setosa Class: versicolor Class: virginica
Sensitivity 1.0000 0.8000 0.7500
Specificity 1.0000 0.9500 0.9231
Pos Pred Value 1.0000 0.8889 0.6000
Neg Pred Value 1.0000 0.9048 0.9600
Prevalence 0.5333 0.3333 0.1333
Detection Rate 0.5333 0.2667 0.1000
Detection Prevalence 0.5333 0.3000 0.1667
Balanced Accuracy 1.0000 0.8750 0.8365
That’s a 84% score.
Optimization
Let’s plot the logistic regression curve for the test data set.
Step 1 – Get the data
iris_data = iris[51:150,] iris_data = iris_data[order(iris_data$Sepal.Length),] head(iris_data)
| Sepal.Length | Sepal.Width | Petal.Length | Petal.Width | Species | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <fct> | |
| 58 | 4.9 | 2.4 | 3.3 | 1.0 | versicolor |
| 107 | 4.9 | 2.5 | 4.5 | 1.7 | virginica |
| 61 | 5.0 | 2.0 | 3.5 | 1.0 | versicolor |
| 94 | 5.0 | 2.3 | 3.3 | 1.0 | versicolor |
| 99 | 5.1 | 2.5 | 3.0 | 1.1 | versicolor |
| 60 | 5.2 | 2.7 | 3.9 | 1.4 | versicolor |
Step 2 – Model the data using a classifier
model = glm( Species ~ Sepal.Length, data = iris_data , family = binomial)
Step 3 – Plot the Logit curve.
library(pROC) plot( iris_data$Sepal.Length, iris_data$Species, pch = 19, col = iris_data$Species) points ( iris_data$Sepal.Length, model$fitted.values + 2, col = "orange", type="l", lty=2, lwd=3)
As you can see, still there are quite a bit of mis-classifications. All the false negatives and false positives in the plot above are examples of mis-classification. Irrespective of the algorithm used to calculate the fit, there is only so much that can be done in increasing the classification accuracy given the data as-is. Other terms for True Positive and True Negative are
- Sensitivity ( True Positive )
- Specificity ( True Negative )
There is a specific optimization that can be done – and that is to specifically increase accuracy of one segment of the confusion matrix at the expense of the other segments. For example, if you look at a visual of the confusion matrix for our dataset.
For this dataset, classifying the species as “setosa” is positive and “versi-color” as negative.
- setosa – positive
- versi-color – negative
Let’s actuall calculate the accuracy values. Say the confusion matrix looks like this.[[11 1] [ 1 12]]
tp = (11) / (11 + 1)
fp = (1) / (11 + 1)
fn = (1) / (1 + 12)
tn = (12) / (12 + 1)
cat ("True Positive = ", tp, "\n")
cat ("False Positive = ", fp, "\n")
cat ("True Negative = ", tn, "\n")
cat ("False Negative = ", fn, "\n")
True Positive = 0.9166667 False Positive = 0.08333333 True Negative = 0.9230769 False Negative = 0.07692308
What if we want to predict 100% of setosa ( or a much more accurate classification than 0.9 ). Of course, like we discussed earlier, it will come at a cost. However, there is a usecase for this scenario. For example, if getting a particular classification right is extremely important, then we focus more on that particular classification than the others. Have you seen the Brad Pitt’s movie “World War Z” ? A plague emerges all around the world and an asylum is set up in Israel with a high wall. However, when you enter the wall, they make absolutely sure that you do not have the plague. Say, if you have the plague and if you call that as positive, then essentially you maximize the green box in the picture above.
Or another example would be, if you were to diagonize cancer patients, you would rather want to increase the odds of predicting a cancer patient if he/she really has it (True positive). Even it it comes at a cost of wrongly classifying a non-cancer patient as positive ( false positive ). The former can save a life while the later will just cost the company a patient.
Evaluation
ROC Curve
Receiver Operating Characteristics – also called ROC Curve is a measure of how good the classification is. Scikit Learn has an easy way to create ROC curve and calculate the area under the ROC curve. First off, let’s start with a classifier like Logistic Regression
Step 1 – Get the data
iris_data = iris[51:150,] iris_data = iris_data[order(iris_data$Sepal.Length),] model = glm( Species ~ Sepal.Length, data = iris_data , family = binomial) library(pROC) # iris$Species has 3 classes and hence 3 factors. So, we are converting them to # 0/1 factor using factor (c(iris_data$Species) - 2). # -2 so that 3,2 become 1,0 roc = roc ( factor (c(iris_data$Species) - 2), model$fitted.values, plot = TRUE, legacy.axes = TRUE)
Setting levels: control = 0, case = 1 Setting direction: controls < cases
Area under the curve is an indicator of how accurate our classifier is. You can get it as follows.
roc$auc
0.7896
Reference
- Scikit Learn Implementation details of Logistic Regression
- Receiver Operating Characteristics (ROC) Curve from scratch