Logistic Regression (LR) with sklearn: Breast Cancer Dataset Case

This is a categorical dataset, meaning that we must adapt the target vectors from qualitative labels to binary where malignant=0 and benign=1 could represent one’s health state. 

Examples of categorical problems can be found on Medical diagnosis, where one can predict the future condition of a patient given previous symptoms and health state; Banking, where it is possible to forecast fraud using location, historical transactions, frequency of categorical transactions, etc.; Advertisement: using the clients’ history, site, and web search tags, we can predict the likelihood of increased engagement and clicks to boost advertisement and revenue; Entertainment industry, platforms such as Spotify and Netflix make extract and process engineered features from songs and movies, respectively, taking into account the customers’ wishlist to provide better recommendations.

First, we will discuss the Linear and Multi-linear Regression models, which are straightforward methods and very useful in some study cases. Often, two quantities x and y, both directly measured, are related by a theoretical expression y_i:

y_i = \beta_0+\beta_1x_i + \epsilon_i

such that \epsilon_i \approx N(0,\sigma^2) is the noise. The latter involves parameters that can be evaluated from the observed dataset. The most common and simplest case is when the linear equation can fit the quantities, where \beta_1 is the slope and \beta_0 represents the intersection with the Y-axis. The simplest problem is to adjust the best set of parameters to the experimental/ground truth by considering a linear equation. If the pairs (x_i, y_i) are considered “true” values, we would be able to fit the line perfectly. Nevertheless, since x_i and y_i are subject to errors, each point’s position is not precisely determined. So, instead of the ideal point, there is an ellipse with axes s_x and s_y (standard deviations), whose centers are not expected to lie on a straight line but distributed on each side. The least-squares method suggests that the best line minimizes the sum of the squares of the distances from the ellipses’ centers to the line. This distance is measured in some appropriate direction, which depends on the relative deviations of x and y‘s, and whether these quantities have the same physical dimension. The simplest case relies on supposing that one of the quantities, say, the variable x_i, is measured precisely, while all errors are focused on the y_i. This situation can be represented graphically by vertical lines centered on the (x_i, y_i). The relationship defines the deviation of y_i as y_i-\overline{y}.

Anyway, our primary goal is to estimate the parameters of the model \beta_0 and \beta_1 that minimizes the residual sum of squares using the training dataset:

E =\sum_{i=1}^{n} (y_i – \overline{y})^2 =\sum_{i=1}^{n} (y_i – \beta_0 -\beta_1x_i)^2 = \sum_{i=1}^{n} y_i^2 +n\beta_0^2 + \beta_1 \sum_{i=1}^{n}x_i^2 – 2\beta_0\sum_{i=1}^{n}(y_i – \beta_1 x_i)- 2\beta_1\sum_{i=1}^{n}(y_ix_i)

The error represents the differences between the model’s prediction and the ground truth values obtained from an experiment or observation, and we need to choose the parameters that minimizes E. To find the parameters of the model that minimizes the sum of the errors, we can use the partial derivates \frac{\partial E}{\partial \beta_0} and \frac{\partial E}{\partial \beta_1} (\frac{\partial E}{\partial \beta_i}=0):

\frac{\partial E}{\partial\beta_1} = 2\beta_1 \sum_{i=1}^{n} x_i^2 + 2\beta_0 \sum_{i=1}^{n} x_i – 2\sum_{i=1}^{n}(x_i y_i) = 0 ~(1)\newline \frac{\partial E}{\partial\beta_0} = 2n\beta_0 – \sum_{i=1}^{n} y_i + 2\beta_1 \sum_{i=1}^{n} x_i = 0 \\ \rightarrow \beta_0 = \sum_{i=1}^{n}y_i -\beta_1 (n^{-1} \sum_{i=1}^{n}x_i) = \overline{y}-\beta_1 \overline{x}~(2)

where~\overline{x}=n^{-1}\sum_{i=1}^{n} x_i~,~\overline{y}=n^{-1}\sum_{i=1}^{n} y_i \newline

and applying equation (2) into (1), we have:

\beta_1 (\sum_{i=1}^{n} x_i^2 – \overline{x}\sum_{i=1}^{n} x_i) = \sum_{i=1}^{n} x_i y_i – \overline{y}\sum_{i=1}^{n}(x_i)

\beta_1 = \frac{n\sum_{i=1}^{n}x_i y_i – \sum_{i=1}^{n}y_i \sum_{i=1}^{n}x_i}{n\sum_{i=1}^{n} x_i^2 – (\sum_{i=1}^{n} x_i)^2 } = \frac{\sum_{i=1}^{n}(x_i -\overline{x})(y_i – \overline{y})}{\sum_{i=1}^{n}(x_i – \overline{x})^2}

and if we substitute \beta_1, \overline{x} and \overline{y} into (2), we obtain:

\beta_0=\frac{\sum_{i=1}^{n}x_i^2 \sum_{i=1}^{n}y_i – \sum_{i=1}^{n}x_i \sum_{i=1}^{n}x_i y_i}{n\sum_{i=1}^{n} x_i^2 – (\sum_{i=1}^{n} x_i)^2 } = \frac{1}{n}(\sum_{i=1}^{n}y_i – \beta_1\sum_{i=1}^{n}x_i)

We can also estimate the intervals with a probability of 95% to reproduce the measurements by using the computed parameters, where \beta_0[95\%] = [\beta_0 \pm 2*\sigma_{\beta_0} ] and \beta_1[95\%] = [\beta_1 \pm 2*\sigma_{\beta_1} ]:

\sigma_{\beta_0}^2 = \left[ \sigma^2 n^{-1} + \frac{\overline{x}}{\sum_{i=1}^{n}(x_i-\overline{x})^2}\right] \\ \sigma_{\beta_1}^2 = \left[ \frac{\sigma^2}{\sum_{i=1}^{n}(x_i-\overline{x})^2} \right]

where \sigma^2 is the variance.

If you want to go deeper, please try to follow the references [1] and [2]. In our course, we will also discuss the hypothesis test and interpretation of linear models. Since we are moving towards the logistic regression, it worth mentioning that if we have more than one feature, which is the case of our dataset, it would be more reasonable to work with multi-linear regression. Therefore, this new model implies:

y_i = \beta_0 + \beta_1 x_{i,1} + \beta_2 x_{i,2}, … + \beta_p x_{i,p} + \epsilon_{i}

where p is the number of features and i represents the number of data for training. The latter can also be represented through a matrix formulation, where we can have a feature vector {X_i}=[x_{i,1},x_{i,2},…,x_{i,p}] for each datapoint.

Furthermore, we can express the vector Y as:

Y = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2}, … + \beta_p X_{p} + \epsilon_{i}

The next step is to use the least-square minimum approach, where:

E(\beta) = \sum_{i=1}^{m}(y_i – \overline{y}(\beta))^2 = \sum_{i=1}^{m}(y_i -\beta_0 – \beta_1 x_{i,1}-\beta_2 x_{i,2}, … – \beta_p x_{i,p})^2  \\ = ||Y-X\beta||^2_2, ~where~ \beta = (\beta_0, \beta_1,…,\beta_p)

However, we do not know if the variables X_i are useful for predicting the outcome of Y, and actually, how many variables would be reasonable to use? Furthermore, how accurate are these predictions? These nuances will be discussed in our machine learning course.

Let us consider the matrix formulation for our multi-linear regression model X\beta=Y, therefore, we can rewrite the latter as (X^TX)\beta = X^TY. By using matrix operations:

(X^TX)^{-1}(X^TX)\beta = (X^TX)^{-1}X^TY \rightarrow (X^TX)^{-1}(X^TX) = \mathbb{1} \\ \beta = (X^TX)^{-1}X^TY

We can find the parameters of the model that minimizes the sum of the square of the errors. However, each patient has a medical statement (Benign or Malignant breast cancer), and by taking the Bayesian approach, we can assign a probability of having a benign or malignant condition for each patient given a set of features (medical conditions). We know P(Y|X) from the training dataset. But our goal is to predict if a patient has cancer, which means that, given an input feature vector for patient i X_i=[x_{i,1},x_{i,2},…,x_{i,p}], we need to predict the response Y_i:

\hat {y} = argmax_{y} ~\hat {P_i}(Y=y | X=X_i)

and since Y is a binary variable (benign or malignant), the probabilities can be expressed as:

\hat {P_i}(Y=0 | X=X_i) = \beta_0 + \beta_1 X_{i,1} + + \beta_2 X_{i,2} + … + + \beta_1 X_{i,p}

The main problem relies on the fact that this model allows the probability to be lower than 0 and higher than 1, and therefore \hat {P}(Y=0 | X=X_i) + \hat {P}(Y=1 | X=X_i) \neq 1 . Furthermore, it would be extremely hard to extend this model to multi-classes beyond binary. Instead, we can model the joint probability as:

\hat {P_i}(Y=1 | X=X_i) = \frac{e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}}{1+e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}} = \frac{e^{\sum_{j=1}^p \beta_j X_{i,j}}}{1+e^{\sum_{j=1}^p \beta_j X_{i,j}}}

and since P(Y=1|X) + P(Y=0|X) = 1, therefore:

\hat {P_i}(Y=1 | X=X_i) = 1-\frac{e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}}{1+e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}} =\frac{1}{1+e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}}= \frac{1}{1+e^{\sum_{j=1}^p \beta_j X_{i,j}}}

which is the same as if we use the logaritmic function:

\log \left[\frac{P_i(Y=1|X)}{P_i(Y=1|X)} \right] = \beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}

and a list of pairs (x_1,y_1),(x_2,y_2),…,(x_n,y_n) as the training dataset. The probability associated with the whole dataset is represented by the likelihood, and our goal is to find the set of parameters \beta_0,beta_1,…,\beta_p that maximize the likelihood. The likelihood L(\beta) can be expressed as:

L = \prod_{i=1}^{n} P(Y=y_i|X=X_i) \\ L=\prod_{i=1}^{n} \left[ \frac{e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}}{1+e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}} \right]^{y_i} \\* \left[ \frac{1}{1+e^{\beta_0 + \beta_1 X_{i,1}+\beta_2 X_{i,2}+..+\beta_p X_{i,p}}} \right]^{1-y_i}

where the probability is represented by \prod_{i=1}^{n} P(Y=1|X=X_i)^1 P(Y=0|X=X_i)^{(1-1)}=\prod_{i=1}^{n} P(Y=1|X=X_i) for y_i = 1, and filtered for the second term inside the product when y_i=0, therefore \prod_{i=1}^{n} P(Y=1|X=X_i)^0 P(Y=0|X=X_i)^{(1-0)}=\prod_{i=1}^{n} P(Y=0|X=X_i). In our machine learning course, you will learn how to solve this problem from scratch using Newton’s algorithm. But for now, you will use the sklearn library to find the parameters of the model that best fits our training data. Furthermore, we will use these parameters to predict the outcome classes of a testing data set of patients. Then, we will compare our results with the ground truth to estimate the accuracy of our model. Let’s start!

We have the index, ID, 9 features from briopsy[:,2:11] and the classification as benign or malign.