<jupyter_start><jupyter_text>Logistic regression----This kind of model tries to deal with binary outcomes. It's a generalized linear model: the basic ideia is that we know well how to model continuous outcome from linear regression, so we'll try to fit our binary outcome with a continuous function that, with some transformations, can be linearized given our preditor variables.Each outcome $y_i$ is believed to follow a Bernoulli distribuiton given by:$$ \Pr(y_i = 1) = p$$Our model will try to estimate a probability $p_i$ given a called $logit$ - which is modeled as log function of $p$ and it's value is obtained by a linear function of the predictors:$$ logit(\hat{p}) = log(\frac{\hat{p}}{1-\hat{p}}) = \beta X $$From the equation above, one can see that $logit^{-1}$ gives us the value of $\hat{p}$ estimated by the model. We can calculate doing some algebraic transformations: $$\frac{\hat{p}}{1-p} = \exp{(\beta X)} \\\hat{p}(1 + \exp{(\beta X)}) = \exp{(\beta X)} \\\hat{p} = \frac{1}{1 + \exp{(-\beta X)}}$$Note that this is a sigmoid function, that is steppest at $logit^{-1}(0.5)$ $-$ one can interpret it thinking about the probability of $y_i = 1$: values above $0.5$ are believed to have a positive outcome since we're try to model the expected value [1] as a "turning point", because that's no such thing as a $0.5$ outcome.[1] $E(p) = \frac{1+0}{2} = 0.5$ References- [<NAME>., & <NAME>. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge university press.](https://www.cambridge.org/core/books/data-analysis-using-regression-and-multilevelhierarchical-models/32A29531C7FD730C3A68951A17C9D983) Table of Contents1&nbsp;&nbsp;References2&nbsp;&nbsp;Example: One predictor<jupyter_code># %load first_cell.py
%reload_ext autoreload
%autoreload 2

from paths import RAW_PATH, TREAT_PATH, OUTPUT_PATH, FIGURES_PATH

from copy import deepcopy
import numpy as np
import pandas as pd
pd.options.display.max_columns = 999
import pandas_profiling

import warnings
warnings.filterwarnings('ignore')

# Plotting
import plotly
import plotly.graph_objs as go
import cufflinks as cf
plotly.offline.init_notebook_mode(connected=True)

cf.go_offline()
cf.set_config_file(offline=False, world_readable=True)<jupyter_output><empty_output><jupyter_text>Example: One predictorLet's simulate some trivial data - in this case, we have a categorical variable $x \in [1, 5]$ and a outcome $y$ that is simply $y_i = 1, x \geq 3; y = 0, otherwise.$<jupyter_code>x = np.random.random_integers(1, 5, 100)
y = [1 if i > 2 else 0 for i in x]

df = pd.DataFrame(data=np.array([x,y]).T, columns=['x', 'y'])
df.head()<jupyter_output><empty_output><jupyter_text>We can visualize this outcomes for each value of the predictor in the graph below (note that the point are _jittered_ - we added random noise to $x$ and $y$ values so that they wont all fall into the same point).<jupyter_code>df['x_jitter'] = df['x'].apply(lambda x: x + np.random.random()/5 - np.random.random()/5)
df['y_jitter'] = df['y'].apply(lambda x: x + np.random.random()/10 - np.random.random()/10)

df[['x_jitter', 'y_jitter']].set_index('x_jitter').iplot(kind='scatter', mode='markers', 
                                                         color='blue', theme='white',
                                                         title='Simulated outcomes for 5-level variable',
                                                         yTitle='Outcome (y = 1 / y = 0)',
                                                         xTitle='Predictor variable')<jupyter_output><empty_output><jupyter_text>Now, we can run a logistic regression model to get the probabilities of the outcome given our predictor variable:<jupyter_code>from sklearn.linear_model import LogisticRegression

model = LogisticRegression(random_state=0)
model.fit(df[['x']], df['y'])<jupyter_output><empty_output><jupyter_text>The model generated the following formula for our data: $logit(p) = -3.66 + 1.66 x$. The independent term, $-3.66$, is our _intercept_ and $1.66$ is the _slope_ for the variable x (you can check the values bellow).<jupyter_code>model.intercept_, model.coef_<jupyter_output><empty_output><jupyter_text>We can then visualize our logistic function below:<jupyter_code>def logit(x):
    return model.intercept_[0] + x*model.coef_[0][0]
    
df['logit'] = df['x'].apply(lambda x: logit(x))

df[['x', 'logit']].set_index('x').iplot(kind='scatter', mode='markers', 
                                        color='orange', theme='white',
                                        title='Logistic regression',
                                        xTitle='Predictor variable', 
                                        yTitle='logit(p)')<jupyter_output><empty_output><jupyter_text>Note that the logistic regression is a linear function that has no specifc range - it can go from $(-\infty, \infty)$.We can check the probability function estimated by the model doing the inverse transformation: $$\hat{p} = \frac{1}{1 + \exp{(-\beta X)}}$$<jupyter_code>def inv_logit(x): 
    beta_x = logit(x)
    return 1 / (1 + np.e**(-beta_x))

df['inv_logit'] = df['x'].apply(lambda x: inv_logit(x))
df['inv_logit'].head()
df[['x', 'inv_logit']].set_index('x').iplot(kind='scatter', mode='markers', 
                                            color='green', theme='white',
                                            title='Estimated probability', 
                                            xTitle='Predictor variable', 
                                            yTitle='p estimated')<jupyter_output><empty_output><jupyter_text>Although $x$ is a categorical variable, we've put in more values in between $[1,5]$ for better visualize the shape of the function:<jupyter_code>xvalues = np.linspace(1, 5, 100)
inv_logit_xvalues = [inv_logit(x) for x in xvalues]

df_ex = pd.DataFrame(data=np.array([xvalues, logit_xvalues, inv_logit_xvalues]).T, 
                     columns=['x', 'logit', 'inv_logit'])

df_ex[['x', 'inv_logit']].set_index('x').iplot(kind='scatter', mode='lines', 
                                               color='green', theme='white',
                                               title='Estimated probability', 
                                               xTitle='Predictor variable', 
                                               yTitle='p estimated')<jupyter_output><empty_output><jupyter_text>This is the usual plot we see as the relationship of the estimated probability and our predictor variable, and we can see that the relationship of those two is not linear, but exponential. * Which is the point that gives us the separation of positive and negative outcomes?$$ logit(\hat{p}) = log(\frac{p}{1-p}) = -\beta X \\ p = 0.5 \rightarrow log(\frac{0.5}{1-0.5}) = log(1) = 0$$So we need to find the point at which the value of $logit(\hat{p}) = \beta X = 0$. Going back to our coeficientts, we can see that:$$-\beta X = 0 \rightarrow -3.66 + 1.66 x = 0 \rightarrow x \approx 2.205$$<jupyter_code>logit(2.205)<jupyter_output><empty_output><jupyter_text>Great! We've found the "turning point", and we can check also that $\hat{p} = logit^{-1}(\beta X) = 0.5$:<jupyter_code>inv_logit(2.205)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup