<jupyter_start><jupyter_text>Let's look at some iron ore data<jupyter_code>import pandas
import seaborn
import matplotlib.pyplot as plt
import os
import numpy as np
seaborn.set()

df = pandas.read_csv(os.path.abspath('iron_ore_study.csv'))

# Splits from oscar Fe>60%, SiO2<9, Al2O3<2, P<0.08
split_points = [
    ('FE', 60, [False, True]),
    ('SIO2', 9, [True, False]),
    ('AL2O3', 2, [True, False]),
    ('P', 0.08, [True, False]),  
]

# It's ore if everything is True
df['is_ore'] = np.vstack([
    pandas.cut(df[elem], bins=[0, split, 100], labels=is_ore)
    for elem, split, is_ore in split_points
]).sum(axis=0) == 4

# Take a look
seaborn.pairplot(df.iloc[::5], hue='is_ore', plot_kws={'alpha': 0.5})
seaborn.countplot('is_ore', data=df)<jupyter_output><empty_output><jupyter_text>Plot our x data<jupyter_code>seaborn.jointplot('SIO2', 'FE', df, joint_kws={'alpha': 0.4, 'marker': '.'})
seaborn.jointplot('P', 'FE', df, joint_kws={'alpha': 0.4, 'marker': '.'})<jupyter_output>C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\axes\_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\axes\_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "<jupyter_text>Logistic regressionYou can think of the logistic function as a function that takes a real number (as comes out of the linear regression) and 'squashes' it into a 0, 1 label. It's defined as$$g(z) = \frac{1}{1 + e^{-z}}$$<jupyter_code>def logistic(z):    
    return 1 / (1 + np.exp(-z))

zs = np.linspace(-10, 10)
logistic_data = pandas.DataFrame( 
    {'z': zs, 'logistic': logistic(zs)}
)
logistic_data.plot('z', 'logistic')<jupyter_output><empty_output><jupyter_text>If we can generate a linear model fit with one parameter - call it $f$:$$z = f(x) = a + b x$$where $a$ is the intercept, $b$ the coefficient and $x$ is the input features. Then we get label predictions$$\mathrm{label} = g(f(x)) > threshold$$  ProblemLet's pretend we have a crappy sensor which only measures Al. Can we still make good predictions of ore/not ore using just this feature?We should look at transforming our aluminium data so that we go from (0, inf) -> (-inf, inf). We'll do this in a hacky sense by using a log function but we should really use a log-ratio transform here!<jupyter_code>seaborn.distplot(df['SIO2'])
seaborn.distplot(np.log(df['SIO2']))<jupyter_output>C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\axes\_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "<jupyter_text>We'll do this using a scikit-learn pipeline - this lets us chain transformations and predictions into one object which makes life a lot easier.<jupyter_code>from sklearn import preprocessing, pipeline, linear_model

# Make up our pipeline where we transform the aluminium first to make it more gaussian!
regressor = pipeline.Pipeline([
    ('transform', preprocessing.FunctionTransformer(np.log, validate=True)),
    ('model', linear_model.LogisticRegression())
])<jupyter_output><empty_output><jupyter_text>Next we map the data into the `y ~ f(X)` format that scikit-learn wants<jupyter_code>X = df[['AL2O3']]
y = df['is_ore']<jupyter_output><empty_output><jupyter_text>Now fitting the model is as simple as<jupyter_code>regressor.fit(X, y)<jupyter_output><empty_output><jupyter_text>Once we've fitted the model we can make predictions straight away<jupyter_code>predict_df = pandas.DataFrame(
    {'test_al2o3_values': [0.5, 1, 2, 3, 4]} # Are these values ore?
)  
regressor.predict(predict_df)<jupyter_output><empty_output><jupyter_text>To see what's going on in a bit more depth, we can pull the coefficients out of the scikit-learn pipeline<jupyter_code>model = regressor.named_steps.model
model.intercept_, model.coef_<jupyter_output><empty_output><jupyter_text>and rewrite our logistic function to include the linear model<jupyter_code>def logistic(regressor, X):
    """
    Plot our logistic model given input values x
    
    We're doing this so that we can see the output of the logistic function - normally
    you'd just do `regressor.predict(x)` to get actual 1, 0 labels for your data.
    
    Parameters:
        regressor - a fitted logistic regression pipeline
        x - the values to evaulate the function at
    """
    # We can pull the model and transforms from our pipeline
    model = regressor.named_steps.model
    tf = regressor.named_steps.transform
    
    # Next we replay the steps in the pipeline to make a prediction
    z = model.intercept_ + model.coef_[0][0] * tf.transform(X)
    return 1 / (1 + np.exp(-z)).ravel()<jupyter_output><empty_output><jupyter_text>Now we can evaulate our logistic function for our test values<jupyter_code>logistic(regressor, predict_df)<jupyter_output><empty_output><jupyter_text>With these in hand lets generate some plots<jupyter_code>f, ax = plt.subplots(1, 1)

# Some aluminium values to predict from
al_compositions = pandas.DataFrame(
    {'test_al2o3_values': np.linspace(0.1, 3)}
)

# An offset to stop everything plotting on top of everything else
offset = 0.02

# shows predictions given contents
predictions = regressor.predict(al_compositions)
ax.plot(al_compositions, predictions + offset, '.', alpha=0.7, label='predicted (+ offset)')  

# shows measured values plus jitter
jitter = np.random.normal(scale=0.01, size=len(df))
ax.plot(df['AL2O3'], df['is_ore'] + jitter - offset, '.', alpha=0.1, label='measured (+ jitter - offset)')

# shows logistic function fitted from regressor
ax.plot(al_compositions, logistic(regressor, al_compositions), '--', label='fitted logistic function')

# Generate the logistic curve showing the location of 
ax.set_ylim(-0.1, 1.1)
ax.legend()
ax.set_title('Logistic regression with scikit-learn')
f.tight_layout()<jupyter_output><empty_output><jupyter_text>Try this using one of the other variables!  Measuring model performanceWe don't get everything right! How can we get a feeling for the model performance? What are some of the issues that we might need to take into account?<jupyter_code>from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = \
    train_test_split(X, y, test_size=0.33)<jupyter_output><empty_output><jupyter_text>Now we can train on just the training set, predict on the test set and see how we do!<jupyter_code>regressor.fit(X_train, y_train)

y_predict = regressor.predict(X_test)
results = pandas.DataFrame({
    'AL2O3': X_test['AL2O3'],
    'is_ore_actual': y_test,
    'is_ore_predicted': y_predict,
    'count': 1
})
results<jupyter_output><empty_output><jupyter_text>Now we can see the number of false positives and false negatives using pivot_table from last week<jupyter_code>confusion = results.pivot_table(values='count', index='is_ore_actual', columns='is_ore_predicted', aggfunc='sum')
confusion<jupyter_output><empty_output><jupyter_text>and plot with seaborn<jupyter_code>seaborn.heatmap(confusion, cmap='Reds')<jupyter_output><empty_output><jupyter_text>Alternatively we can unstack the array into true and false negatives and positives<jupyter_code>confusion.unstack()<jupyter_output><empty_output><jupyter_text>And we can convert these values to fractions of the total<jupyter_code>true_neg, false_neg, true_pos, false_pos = confusion.unstack() / confusion.unstack().sum()
true_neg, false_neg, true_pos, false_pos<jupyter_output><empty_output><jupyter_text>When are true positives and false positives important? ThresholdWe haven't done anything with the threshold yet - how should we pick the value for this?We've already got a confusion matrix - we can take the ratio of the true_positive vs the false_positive rates and compare the two<jupyter_code>confusion
correct = np.diag(confusion)
incorrect = np.diag(np.roll(confusion, 1, axis=1))
correct, incorrect
total_correct = correct.sum()
total_incorrect = incorrect.sum()
total_correct, total_incorrect<jupyter_output><empty_output><jupyter_text>And we can plot this for each threshold. We need a way of adjusting the class weights in the model. Scikit-learn doesn't let you specify a threshold directly but<jupyter_code>def fit_with_class_threshold(threshold):
    "Fit a logistic regression to get an ROC value for a given threshold"
    # Transform our threshold into class weights
    class_weights = {True: threshold, False: 1 - threshold}
    
    # Make a regressor
    regressor = pipeline.Pipeline([
        ('transform', preprocessing.FunctionTransformer(np.log, validate=True)),
        ('model', linear_model.LogisticRegression(class_weight=class_weights))
    ])
    
    # Fit it
    regressor.fit(X_train, y_train)
    
    # Make some predictions, see how we did
    results = pandas.DataFrame({
        'AL2O3': X_test['AL2O3'],
        'is_ore_actual': y_test,
        'is_ore_predicted': regressor.predict(X_test),
        'count': 1
    })
    confusion = results.pivot_table(
        values='count', 
        index='is_ore_actual', 
        columns='is_ore_predicted', 
        aggfunc='sum')
    true_neg, false_neg, true_pos, false_pos = confusion.unstack()
    correct = np.diag(confusion)
    incorrect = np.diag(np.roll(confusion, 1, axis=1))
    
    # Return results as a dictionary
    return {
        'threshold': threshold,
        'total_correct': correct.sum(), 
        'total_incorrect': incorrect.sum(), 
        'true_negative': true_neg,
        'false_negative': false_neg,
        'true_positive': true_pos,
        'false_positive': false_pos
    }<jupyter_output><empty_output><jupyter_text>Now we can generate results using our metric<jupyter_code>fit_with_class_threshold(0.1)
fit_with_class_threshold(0.9)<jupyter_output><empty_output><jupyter_text>Now we can iterate over all our thresholds and see what does the best<jupyter_code>results = pandas.DataFrame.from_records(
    [fit_with_class_threshold(t) for t in np.linspace(0.1, 0.9)],
    index='threshold'
)
results.head()<jupyter_output><empty_output><jupyter_text>We'll generate a few plots<jupyter_code>results.plot(y=['total_correct', 'total_incorrect'])
results.plot(y=['false_negative', 'true_negative'])
results.plot(y=['false_positive', 'true_positive'])<jupyter_output><empty_output><jupyter_text>Extension - modelling with statsmodelsFor what it's worth we can also generate these using statsmodels<jupyter_code>import statsmodels.api as sm
from numpy import log<jupyter_output><empty_output><jupyter_text>We have to handle preprocessing ourselves<jupyter_code>from statsmodels.tools import add_constant

def preprocess(x):
    "Our preprocessing pipeline for Al2O3"
    return add_constant(np.log(x)) # add_constant adds an intercept to the fit<jupyter_output><empty_output><jupyter_text>Statsmodels uses stats jargon - endog -> endogenous variable -> y- exog -> exogenous variable -> X<jupyter_code>endog = df.is_ore
exog = preprocess(df['AL2O3'])<jupyter_output><empty_output><jupyter_text>Fitting the model is pretty similar though<jupyter_code>model = sm.Logit(endog, exog)
results = model.fit()
results.summary()<jupyter_output><empty_output><jupyter_text>Statsmodels does a bit more statistical testing/automated confidence intervals for us at the cost of having to manage crossvalidation etc ourselves. Depending on what you're trying to achieve this could be a viable way to go.As before we can immediately make some predictions - statsmodels gives us the value of the logistic function<jupyter_code>preds = model.predict(results.params, preprocess(test_al_values))
preds
preds > 0.5<jupyter_output><empty_output><jupyter_text>and compare to scikit-learn<jupyter_code>regressor.predict(test_al_values.reshape(-1, 1))<jupyter_output><empty_output><jupyter_text>We can also compare the logistic values get slightly different answers here - probably down to the solver used under the hood<jupyter_code>print(' statsmodels:', model.predict(results.params, preprocess(test_al_values)))
print('scikit-learn:', logistic(regressor, test_al_values))<jupyter_output><empty_output><jupyter_text>We can generate the same plot again though<jupyter_code>f, ax = plt.subplots(1, 1)

# Some aluminium values to predict from
al_compositions = np.linspace(0.1, 3)

# An offset to stop everything plotting on top of everything else
offset = 0.02

# shows predictions given contents
predictions = model.predict(results.params, preprocess(al_compositions)) > 0.5
ax.plot(al_compositions, predictions + offset, '.', alpha=0.7, label='predicted (+ offset)')  

# shows measured values plus jitter
jitter = np.random.normal(scale=0.01, size=len(df))
ax.plot(df['AL2O3'], df['is_ore'] + jitter - offset, '.', alpha=0.1, label='measured (+ jitter - offset)')

# shows logistic function fitted from regressor
ax.plot(al_compositions, model.predict(results.params, preprocess(al_compositions)), '--', label='fitted logistic function')

# Generate the logistic curve showing the location of 
ax.set_ylim(-0.1, 1.1)
ax.legend()
ax.set_title('Logistic regression with statsmodels')
f.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11041932338159177, 'fraction_numerical': 0.01185207411296977, 'mean_word_length': 5.506588579795022, 'pattern_counts': {'":': 0, '<': 126, '<?xml version=': 0, '>': 132, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '356122', 'license': 'MIT', 'n_tokens_mistral': 4532, 'n_tokens_neox': 4319, 'n_words': 1226, 'path': 'notebooks/logistic_regression_example.ipynb', 'repo_name': 'bussomick/03-simple-predictions', 'size': '13331'}