<jupyter_start><jupyter_text>====================================================== MA477 - Theory and Applications of Data Science   Lesson 11: Logistic Regression    Dr. <NAME> United States Military Academy, West Point AY20-2====================================================== Lecture Outline    Why not Linear Regression?     What is Logistic Regression?    Estimating Regression Coefficients     Logistic Regression with more than 2 response classes     Python Implementation      Why not Linear Regression?We will illustrate below that in the cases when the response variable is a qualitative variable, then linear regression may not be appropriate. So, why is that? Let's illustrate this using the `iris` dataset. In the iris dataset one is tasked to determine whether a plant is of type `setosa, versicolor` or `virginica` based on `sepal length, width` and `petal length, width`. So, the response variable $Y$ is clearly qualitative. Now, before we can fit a Linear Regression model (or most models for that matter) we need to first encode the outcomes using numbers. So, this is exactly where the dilemma lies. In other words, how do we encode a qualitative response variable with more than two possibilities? We will demonstrate shortly, that the model we obtain is highly sensitive to how you encode the qualitative variable. That is, if two different teams encoded the qualitative variables with more than two outcomes differently, they could end up with very different models which in turn would produce very different predictions on the new, previously unseen by the model, data. Obviously, this is a very undesirable outcome. For the purpose of this demonstration we will use the following two encodings for the response variable: Case 1$$Y_1=\begin{cases}0 & \text{ setosa }\\1 & \text{ versicolor }\\2 & \text{ virginica }\end{cases}$$ Case 2 $$Y_2=\begin{cases}2 & \text{ setosa }\\0 & \text{ versicolor }\\1 & \text{ virginica }\end{cases}$$We'll import our standard libraries as well as the `iris` dataset wich may be found in `sklearn.datasets`.<jupyter_code>import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
iris=load_iris()
iris.target_names
iris.feature_names
print(iris.DESCR)
df1=pd.DataFrame(iris.data,columns=iris.feature_names)
df1.head()
df1['target']=iris.target
df1.head()<jupyter_output><empty_output><jupyter_text>Recall that the current encodings of the response variable are as follows: $$Y_1=\begin{cases}0 & \text{ setosa }\\1 & \text{ versicolor }\\2 & \text{ virginica }\end{cases}$$In what follows, we will create a second dataframe `df2` with the response variable encoded as $$Y_1=\begin{cases}2 & \text{ setosa }\\0 & \text{ versicolor }\\1 & \text{ virginica }\end{cases}$$<jupyter_code>df2=df1.copy()
df2['target']=df2['target'].apply(lambda x: 2 if x==0 else 0 if x==1 else 1)
df2.head()<jupyter_output><empty_output><jupyter_text>Next, we will fit two Linear Regression models, and compare the coefficients of the model.<jupyter_code>from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
lg1=LinearRegression()
lg2=LinearRegression()
X1=df1.drop("target",axis=1)
y1=df1['target']

X2=df2.drop("target",axis=1)
y2=df2['target']
X1_train,X1_test,y1_train,y1_test=train_test_split(X1,y1,test_size=0.1,random_state=1)

X2_train,X2_test,y2_train,y2_test=train_test_split(X2,y2,test_size=0.1,random_state=1)
X1_train.head()
X2_train.head()<jupyter_output><empty_output><jupyter_text>One thing to note is that since we used the same random state for both splits, `X1_train` and `X2_train` are identical, as it is desired. The only thing that differs is `y's`, since we used different encodings. Now, we fit two linear models using these two datasets.<jupyter_code>lg1.fit(X1_train,y1_train)

lg2.fit(X2_train,y2_train)<jupyter_output><empty_output><jupyter_text>Now, let's get the coefficients and compare them.<jupyter_code>coef1=lg1.coef_
coef2=lg2.coef_
df_coef=pd.DataFrame()
df_coef['coef1']=coef1
df_coef['coef2']=coef2
df_coef<jupyter_output><empty_output><jupyter_text>Clearly, these coefficients are very different, and it is to be expected that they would produce different outcomes for the same input data. Let's test this.<jupyter_code>pred1=lg1.predict(X1_test)

pred2=lg2.predict(X2_test)<jupyter_output><empty_output><jupyter_text>Next, let's put everything in a dataframe. Namely, the predictions of the two models, along with the actual classes of the plants.<jupyter_code>df_pred=pd.DataFrame(index=y1_test.index)

df_pred['pred1']=pred1
df_pred['pred2']=pred2
df_pred['actual1']=y1_test
df_pred['actual2']=y2_test
df_pred
df_pred['pred class mod1']=df_pred['pred1'].apply(lambda x: 'setosa' if x<0.5 else 'versicolor' if 0.5<x<1.5 else 'virginca')
df_pred['pred class model2']=df_pred['pred2'].apply(lambda x: 'setosa' if x>1.5 else 'versicolor' if x<0.5 else 'virginca')
df_pred<jupyter_output><empty_output><jupyter_text>As we can see, the two models could yield very different predictions. In fact, in this case, they don't match on a single output.The other issue with this approach is the fact that the outcome of the linear regression cannot quite be interpreted as a probability since they may take values above one and below zero.   Logistic Regression ModelFor the remainder of this lesson, we will focus on the case of  binary response variables. We will use the wisconsin breast cancer data set we previously used when we discussed KNN Classifier. Before we begin discussing the innerworkings of Logistic Regression, let's import the dataset.<jupyter_code>from sklearn.datasets import load_breast_cancer
cancer=load_breast_cancer()
print(cancer.DESCR)
df=pd.DataFrame(cancer.data,columns=cancer.feature_names)
df['target']=cancer.target
df.head()<jupyter_output><empty_output><jupyter_text>Recall that the target variable is a qualitative variable with two categories:  Malignant and  Benign. We will use the following encoding for the response variable:$Y=\begin{cases} 0 & \text{ Malignant }\\1 & \text{ Benign }\end{cases}$ How does Logistic Regression Work?So, the response variable $Y$ can either be $0$ (Malignant) or $1$ (Benign). Instead of modeling directly the response variable $Y$, logistic regression, models instead the  probability  taht $Y$ belongs to a particular category. Specifically, if we denote by $X=(X_1,\dots, X_{30})$ all the $30$ cancer measurements, then logistic regression will instead model the probability that $Y=0$ or ($Y=$ Malignant) given the measurements $X$; that is $$p(X):= P\left(Y=0|X\right)$$The values of $p(X)$ will be between $0$ and $1$. So, for example, we may predict  Malignant if $p(X)>0.5$, or if we don't want to take a risk in misclassifying a malignant cell as benign, we can be more conservative and predict Malignant if $p(X)>0.1$ and predict  Benign only if $p(X)\leq 0.1$. This is something that you as a Data Scientist have to decide depending on the goal at hand. Our goal is to find a good way to model the conditional probability $p(X)$. Modeling $p(X)$How can we model $p(X)?$.We could try modeling these probabilities via linear regression $$p(X)=\beta_0+\sum_{i=1}^{30}\beta_iX_i$$ Discussion: Is this a good approach? Explain!What is one of the main conditions that $p(X)$ must satisfy?Well, since $p(X)$ represent probabilities their values must be between $0$ and $1$. So, we must hunt for functions whose output satisfies this condition. While there many functions that satisfy this condition, in the case of Logistic Regression we use the  logistic function  (in Neural Networks this function is typically referred to as the  Sigmoid function):$$S(z)=\frac{1}{1+e^{-z}}=\frac{e^z}{e^z+1}$$ ExercisePlot $S(z)$ in the range $(-10,10)$<jupyter_code>def sigmoid(z):
    return 1/(1+np.exp(-z))
sns.set_style('whitegrid')
plt.figure(figsize=(10,6))

z=np.linspace(-10,10,50)

plt.yticks(np.arange(0,1.01,0.1))

plt.plot(z,sigmoid(z), 'r--', lw=3, label='$S(z)=e^z / (e^z+1)$')
plt.legend(fontsize=14)
plt.show()<jupyter_output><empty_output><jupyter_text>So, we would model $p(X)$ as follows:$$p(X)=P(Y=Malignant|X)=\frac{e^{\beta_0+\sum_{i=1}^{30}\beta_iX_i}}{1+e^{\beta_0+\sum_{i=1}^{30}\beta_iX_i}}$$ Estimating the coefficientsThe goal now is to estimate the regression coefficients $\beta_0,\dots,\beta_n$. We will do so via a method known as the  maximum likelihood method. The maximum likelihood method is very general and is typically used to fit many of the non-linear models. The basic intuition behind the maximum likelihood method is as follows: Let $X^{(i)}$ denote the $30$ measurements for the $i^{th}$ breast cancer tissue. Then, we seek to find coefficients $\hat{\beta_0},\hat{\beta_1},\dots,\hat{\beta}_{30}$ such that $p\left(X^{(i)}\right)$ is very close to $1$ if the breast tissue is Malignant ($Y^{(i)}=0$) and it is very close to $0$ if the breast tissue is Benign ($Y^{(i)}=1)$.Formally, the function that has this property is known as the  likelihood function :$$L(\beta_0,\beta_1,\dots,\beta_{30})=\prod_{i: Y^{(i)}=0}p\left(X^{(i)}\right)\prod_{j:Y^{(j)}=1}\left(1-p\left(X^{(j)}\right)\right)$$The estimates $\hat{\beta_0},\hat{\beta_1},\dots,\hat{\beta}_{30}$ are chosen to maximize $L(\cdot)$  Extensions of Logistic Regression to Multiple-class Response VariablesThere are ways to use Logisitc Regression to perform classification in the setting wehre the response variable can belong to three or more categories. We briefly mention here two such methods:     One-Versus-One     Suppose that the response variable can belong to one of the $K\geq 3$ classes. In the one-versus-one approach, we construct $\frac{K(K-1)}{2}$ Logistic Regression models, each of which compares a pair of classes. The simplest way to classify a new, previously unseen by the model, observation is via a hard majority vote. Specifically, we use each of the $\frac{K(K-1)}{2}$ models to classify the new observation and tally the number of times that the new observation is assigned to each of the $K$ classes. At the end, the new observation is assigned to the class to which it was most frequently assigned overall.     One-Versus-All    One-versus-all is a similar approach. In this case we fit a total of $K$ models, each time comparing one of the $K$ classes to the remaining $K-1$ classes.     There are other methods that accomodate multiple-class responses, one of the most popular one is Linear Discriminant Analysis (LDA). Python ImplementationFor this portion we will use the wisconsin breast cancer dataset.<jupyter_code>from sklearn.linear_model import LogisticRegression<jupyter_output><empty_output><jupyter_text>Logistic Regression has many parameters. The two most important ones that we will experiment with are `penalty` and `C`. Similar to the Lasso and Ridge methods, we can apply a regularization to the regression coefficients, namely, we can apply an $\ell^1$ or $\ell^2$ regularization, or an `elastic` regularization, which is a weighted average of the two. The constant $C$ is the inverse of the strength of the regularization. In other words, the smaller the value of $C$ the stronger the regularization.<jupyter_code>lr=LogisticRegression(solver='liblinear')<jupyter_output><empty_output><jupyter_text>For now, we will chose the validation-set method to asses the prediction power of the model.<jupyter_code>X=df.drop('target',axis=1)
y=df['target']
X.head()
from sklearn.preprocessing import StandardScaler
scaler=StandardScaler()
scaled=scaler.fit_transform(X)
X_sc=pd.DataFrame(scaled, columns=X.columns)
X_sc.head()<jupyter_output><empty_output><jupyter_text>Next we split the data into a training and a test set and then we fit the model.<jupyter_code>X_train,X_test,y_train,y_test=train_test_split(X_sc,y,test_size=0.3,random_state=11)
lr.fit(X_train,y_train)<jupyter_output><empty_output><jupyter_text>Now the model is fit! What does this mean? It means that the coefficients $\beta_0,\dots,\beta_{30}$ have been estimated. We get them below:<jupyter_code>df_coef=pd.DataFrame(lr.coef_,columns=X_sc.columns,index=['Coefficients'])
df_coef<jupyter_output><empty_output><jupyter_text>ExerciseCarry out the following tasks:Part 1-Make predictions-Display the confusion matrix and classification report-Plot the ROC curve and compute the AUCPart 2Tune the model to get the highest accuracy and the best recall for Malignant class (class 0)<jupyter_code>#'pred' will output the actual class...it will be either a 0 or a 1 in our case
pred=lr.predict(X_test)

#pred_prob will output the actual probabilities for each class
pred_prob=lr.predict_proba(X_test)
pred
pred_prob<jupyter_output><empty_output><jupyter_text>Next we display some of the metrics, such as precision, recall, accuracy etc.<jupyter_code>print(confusion_matrix(y_test,pred))
recall_score(y_test,pred,pos_label=1)
print(classification_report(y_test,pred))<jupyter_output>precision    recall  f1-score   support

           0       0.95      0.93      0.94        61
           1       0.96      0.97      0.97       110

    accuracy                           0.96       171
   macro avg       0.96      0.95      0.96       171
weighted avg       0.96      0.96      0.96       171<jupyter_text>We have an accuracy of $96\%$, a recall rate for the Malignant class of $93\%$ and for Benign class of $97\%$. A natural question is if we can tune our model in such a way as to get better results? Before we do so, let's plot the ROC curve.<jupyter_code>fp,tp,_=roc_curve(y_test,pred_prob[:,1])
area=auc(fp,tp)
plt.figure(figsize=(8,6))
plt.plot(fp,tp, 'r--',lw=3,label='AUC={:.3f}'.format(area))
plt.xlabel("False Positive Rate",fontsize=14)
plt.ylabel("True Positive Rate",fontsize=14)
plt.legend(fontsize=14)
plt.show()<jupyter_output><empty_output><jupyter_text>Our model, as is, is doing a great job in differentiating between the two classes...it is better than what we could achieve with KNN Classifier. Next, we will tune the model...our goal will be to improve the recall score for the Malignant class, all the while trying to maintain as high of an accuracy and AUC as possible. To do so we will search for different $C$ values for each of the different regularizations. We will also try different values for the threshold.<jupyter_code>def custom_predict(pred,thresh):
    predictions=[]
    for item in pred:
        if item>thresh:
            predictions.append(0)
        elif item<=thresh:
            predictions.append(1)
    return predictions
from collections import defaultdict
from sklearn.metrics import accuracy_score
reg={'l1','l2'}
param_c=np.arange(0.1,5,0.1)

def tune_log_reg(param_c,reg,thresh=0.5,iterations=10):

    
    recall=defaultdict(list)
    accuracy=defaultdict(list)
    
    for penalty in reg:
        for c in param_c:
            recall_c=[]
            accuracy_c=[]
            auc_c=[]
            for i in range(iterations):
                X_train,X_test,y_train,y_test=train_test_split(X_sc,y,test_size=0.3,random_state=i+1)
                
                
                lg=LogisticRegression(penalty=penalty,C=c,solver='liblinear')
                lg.fit(X_train,y_train)
                pred=lg.predict_proba(X_test)
                pred=custom_predict(pred[:,0],thresh)
                #print(pred)
                #print(recall_score(y_test,pred,pos_label=0))
                recall_c.append(recall_score(y_test,pred,pos_label=0))
                accuracy_c.append(accuracy_score(y_test,pred))
            #print(recall_c)
            recall[penalty].append(np.array(recall_c).mean())
            accuracy[penalty].append(np.array(accuracy_c).mean())
    return recall,accuracy
            
    
    
recall,accuracy=tune_log_reg(param_c,reg,iterations=100,thresh=0.32)
plt.figure(figsize=(10,6))
plt.plot(param_c,recall['l1'],'r--',label='Recall l1')
plt.plot(param_c,recall['l2'],'b-.',label='Recall l2')
plt.plot(param_c,accuracy['l1'],'r-',label='Accuracy l1')
plt.plot(param_c,accuracy['l2'],'b-',label='Accuracy l2')
plt.yticks(np.arange(.93,.985,0.002))
plt.xticks(np.arange(0,5,0.15))
plt.xlim(0,3)
plt.ylim(0.95,0.98)
plt.xlabel('C',fontsize=14)
plt.legend(fontsize=12)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1128129421172418, 'fraction_numerical': 0.02214430706772467, 'mean_word_length': 4.983805668016195, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 77, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '440306', 'license': 'MIT', 'n_tokens_mistral': 5317, 'n_tokens_neox': 5046, 'n_words': 1747, 'path': 'MA477 - Theory and Applications of Data Science/Lessons/Lesson 11 - Logistic Regression/Lesson 11 - LogisticRegression.ipynb', 'repo_name': 'jkstarling/MA477-copy', 'size': '16263'}