# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Introduction to Applied Machine Learning
# ## Operations example based on real-world application
# In this notebook you will implement a machine learning algorithm to model a complex system that has been anonymized from a real problem within the Operations space. You will understand how to look beyond single factor analysis and utilize machine learning to model the system and also discover key factors that contribute to failure modes.
#
# The data has been anonymized, randomized, and any identifying information removed. However, it is still heavily derived from types of operations data you may have seen before and is a great way to begin thinking about how you can use data mining methods and machine learning to solve problems in your space.
#
# ## Contents
# 1. Import dataset
# 2. Data exploration
# 3. Feature engineering
# 4. Modeling
# 5. Model Evaluation
#
# ## How to use this notebook
# - To execute any single block of text or markdown, use ctrl+enter, shift+enter or press the run arrow on the left of the box (only in Colaboratory)
# - To reset the notebook select "Factory reset runtime" from the Runtime tab at the top of Colaboratory

# ## 1. Import dataset

# +
# First let's import our data. The data is split into 2 files
# data_url contains manufacturing-related data that describes the conditions in which the product was manufactured
# label_url contains information related to the outcome of each product subset - and we will discuss later what constitutes bad vs. good product
import pandas as pd

data_url = 'https://raw.githubusercontent.com/jzhangab/DS101/master/1_Data/data.csv'
label_url = 'https://raw.githubusercontent.com/jzhangab/DS101/master/1_Data/outcome.csv'
df = pd.read_csv(data_url, sep = ',')
df_label = pd.read_csv(label_url, sep = ',')
# -

# Let's look at the first 5 rows to begin understanding what factors are available
df.head()

# Now let's look at the outcome dataset
df_label.head()

# You can begin thinking/answering some questions you may have about this data.
#
# 1. How is each data point separated? What denotes a unique data point?
# 2. How many data points do we have?
# 3. Do we have any missing data?
# 4. Do we only have numerical data?
# 5. How many features do we potentially have?
# 6. What is the question we are trying to answer with this data, can this data answer this question?
# 7. etc...

# ## 2. Data Exploration
# Let's now begin exploring the data we have to help ourselves understand the problem and see what approach we might want to take.
#
# Before we can begin plotting our data, we have to deal with 2 problems with data cleanliness
#
# 1. There is a value "Low Value" that is a replacement for undetectable signal. We have to replace this with a number.
# 2. There are null values in our parameter columns that we need to fill.

# For cleanliness, missing data is very important, let's check how much missing data there is for each factor
for col in list(df):
    num_na = len(df[col]) - df[col].count()
    print ("Percent null in column " + col + " is:", 100*num_na/len(df[col]))

# And let's check the number of "Low Value" in each column
for col in list(df):
    num_lv = len(df.loc[df[col] == 'Low Value'][col])
    print ("Percent Low Value in column " + col + " is:", 100*num_lv/len(df[col]))

# ### Important observations
# 1. Thanks domain knowledge, we know that "Low Value" is a placeholder for any measurement that is less than 0.5. Therefore we will replace this by the average of values 0-0.5 (0.25)
# 2. A-priori we don't care about the DATE of manufacture because we know that our process should be in control regardless of the date so we will not consider it for this analysis. This does not mean that in your own analysis on a separate problem that DATE will not be important.
# 3. We also do not care for the DESCRIPTOR_1 and SOURCE columns because we know they are unrelated to manufacturing process. (Again, for just this particular problem. For your own work unstructured data may apply.)

# +
# Let's do the following to clean our data
# 1. Replace all null values with 0
df.fillna(0, inplace=True)

# 2. Replace Low Value with 0.025
df = df.replace("Low Value", 0.025)

# 3. Convert parameter columns to float64 from string objects, we have to do this because "Low Value" defaults each column to string objects, but the machine learning model expects numerical data
for c in list(df):
    if "PARAM" in c:
        df[c] = df[c].astype(float)
# -

# We can use dtypes to check the type of each column to make sure our parameters are now numerical
df.dtypes

# +
# Now let's remove the non-parametric columns from our dataframe. These are columns we will not use for modeling.

# This notation is a generator, a one-line for loop in python. It is used to make the code cleaner and use fewer new variable declarations.
df = df[[c for c in list(df) if c not in ['SOURCE', 'DATE', 'DESCRIPTOR_1']]]

# Alternatively it can take the form of a multi-line for loop as shown in the commented block below.
'''
column_list = []
for c in list(df):
    if c not in ['SOURCE', 'DATE', 'DESCRIPTOR_1']:
        column_list.append(c)
df = df[column_list]
'''
# -

# ### Important observation
# Another potential issue for us is determining what to do with duplicate data. We know based on domain knowledge that each data point should be unique by the NAME column, but how do we handle this if each NAME has multiple SUB_NAMES? This is very common when we deal with products that have multiple measurements for the same attribute. We need to decide if we should treat each measurement separately or aggregate them somehow using mean, max, median, or another measure of center or extremes.
#
# For this particular problem we will use the aggregation by mean. We need to do the following.
#
# 1. Figure out if there are duplicate entires by the NAME column
# 2. If there are, we need to average all data by SUB_NAME for each NAME

# Are there duplicate names? We can check for this by doing a value count for the NAME column
# The function value_counts allows us to count the number of each unique value in the column NAME
df['NAME'].value_counts()

# +
# It looks like there are some NAMEs with up to 32 SUB-NAMEs, we need to average the PARAM values for each NAME
# The groupby function aggregates columns of the dataframe using a list of index columns. Here we want to aggregate by NAME only so we get mean of each parameter by NAME
df = df.groupby(['NAME']).mean().reset_index()

# Remove SUB_NAME from columns because we no longer care about the SUB_NAMEs after averaging by NAME
df = df[[c for c in list(df) if c not in ['SUB_NAME']]]
# -

# Check the dataframe
df.head()

# +
# Let's take a look at the histograms of the dataframe to understand each factor.
import matplotlib.pyplot as plt

# %matplotlib inline
fig = plt.figure(figsize = (15,15))
ax = fig.gca()
df.hist(ax = ax)
# -

# That's a lot of histograms! We can make some quick observations.
#
# 1. Some of the data is normally distributed (PARAM_1, PARAM_10, PARAM_2, PARAM_8)
# 2. Some of the data is skewed (PARAM_3, PARAM_4)
# 3. Some of the data is multi-modal (PARAM_5, PARAM_6, PARAM_7, PARAM_9)
# 4. Some of the parameters are between 0 and 1 while others are between 0 and 10, is this going to be a problem?

# The non-normally distributed data is very interesting, perhaps they are related to each other?
# %matplotlib inline
df.plot.scatter(x = 'PARAM_5',
                y = 'PARAM_6')

# ## 3. Feature Engineering
# We've already performed a lot of data cleaning in the exploration phase. Now we need to finish the job by combining the feature and label dataframes as well as creating our classification label
#
# 1. Merge feature data and label data
# 2. Decide which features to use
# 3. Create a label (for binary classification)

# The labels are located in the dataframe df_label
# df_label has the same issue as the parameter dataframe in that there are NAMEs with duplicates by SUB_NAME
df_label.head()

# Check to see if there are duplicate NAME using value_counts
df_label['NAME'].value_counts()

# +
# Now let's aggregate on NAME using mean
df_label = df_label.groupby(['NAME']).mean().reset_index()

# Remove SUB_NAME from columns because we no longer care about the SUB_NAMEs after averaging by NAME
df_label = df_label[[c for c in list(df_label) if c != 'SUB_NAME']]
# -

df_label.head()

# ### Important observation
# How do we determine if a particular data point represents bad product? Well thanks to a-priori knowledge we know that there is an existing process-monitoring and controls action limit at VALUE = 2.0. Let's explore the label using this action limit

# First let's plot a histogram to visualize the label measure distribution
# %matplotlib inline
fig = plt.figure(figsize = (5, 5))
ax = fig.gca()
df_label.hist(ax = ax)

# +
# Hmm, by the way how many standard deviations from the mean is 2.0?

# First calculate the standard deviation of a column
sd = df_label['MEASURE'].std()

# Then calculate the mean
mean = df_label['MEASURE'].mean()

# Finally print # SD between mean and action limit (2.0)
print("Action limit # SD from mean:", (2-mean)/sd)
# -

# ### Important observation
#
# 1. The label measure itself appears to be normally distributed.
# 2. The action limit is 2.85 SD from mean which means our system does not quite meet the 3.0 SD threshold normally found in most preventive capability analysis systems.
# 3. We have an imbalanced dataset, there are far more "good" datapoints than "bad" datapoints, as bad datapoints are only those above the action limit.

# +
# For a classification label, we need to convert the continuous column MEASURE into binary column (0 or 1) LABEL

# First create a new column and set all values to 0 (not a failure)
df_label['LABEL'] = 0

# Set data point values in the new column LABEL to 1 that is greater than or equal to the action limit (2.0) in the continuous column MEASURE
df_label.loc[df_label['MEASURE'] >= 2, 'LABEL'] = 1
# -

# ### Multicollinearity
# The idea of collinearity is that if certain input factors are closely correlated, they will bias the output of the model by amplifying their particular effects. We need to understand if some of our factors are highly collinear and then reduce bias by removing all but 1 of the collinear factors from the dataframe.
#
# 1. Let's check if we have a collinearity problem in our parameter dataframe

# We can check the correlation (R-square) between variables using a correlation matrix
df.corr()

# +
# To quantify multicollinearity, we will use variance inflation factor (VIF)
# Rule of thumb, VIF above 10 indicates a particular variable ought to be removed
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

# Also - VIF for a constant term should be high because the intercept is a proxy for the constant.
# A constant term needs to be added to accurately measure VIF for the other terms
df_c = add_constant(df[[c for c in list(df) if 'PARAM' in c]])

# inline Generator on a pandas series
pd.Series([variance_inflation_factor(df_c.values, i) 
               for i in range(df_c.shape[1])], 
              index=df_c.columns)
# -

# ### Lastly we have to merge the feature and label dataframes for the modeling phase of the analysis

# Left merge label onto features (look up SQL left join if unsure what this does)
df = pd.merge(df, df_label, on = 'NAME', how = 'left')

# ## 4. Modeling
# In the modeling step we will train a supervised machine learning model to understand relationships in the operations data set. We will then evaluate the model to see how well it predicts.
#
# The particular model that we will use is the random forest algorithm, a classical machine learning algorithm very useful for classification tasks.
#
# 1. Split dataset into training and validation datasets
# 2. Train model
# 3. Predict outcomes of validation dataset
# 4. Calculate accuracy of validation dataset

# +
# We will split the data 80%/20% using 80% of the data to train the model and 20% to validate the accuracy of the model
# We can use pre-built functions from the machine learning package sci-kit learn to do this task
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, confusion_matrix

# The NAME, MEASURE, and LABEL columns are not part of input features so we will use a generator to create a new list and call it "features"
features = [col for col in list(df) if col not in ['NAME', 'MEASURE', 'LABEL']]

# The X input is df[features] which is all columns in the dataframe of the list features we created using the generator
# The Y input is df['LABEL'] which is the binary label column we created in data exploration
X = df[features]
y = df['LABEL']

# Generate the 4 datasets we need
# X_train and y_train to train the model
# X_test to generate predictions
# y_test to evaluate the accuracy of the predictions
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
# -

# Declare and fit model
# We have to use the hyperparameter class_weight because of the very imbalanced classes (very few 1's compared to 0's)
# Balanced class weight uses the inverse of frequency to weight each class
model = RandomForestClassifier(random_state=0, class_weight="balanced")
model.fit(X_train, y_train)

# Predict using test set. These are binary predictions of 1's and 0's
y_pred = model.predict(X_test)
y_pred

# ## 5. Model Evaluation
# We will use several techniques to evaluate the strength of the Model
#
# 1. Accuracy
# 2. Confusion Matrix (false positive, true positive, false negative, true negative)
# 3. Receiver operating characteristic
# 4. Feature importance

# Compare y_test (true values) to y_pred (predicted values)
accuracy_score(y_test, y_pred)

# Let's take a look at the confusion matrix, which shows us false positives and false negatives
# According to the confusion matrix we have only a single false negative and selected correctly the other 2 positive classes
confusion_matrix(y_test, y_pred)

# +
# Another method of evaluating a classifier is using the Receiver Operating Characteristic (ROC)
# ROC is a plot of true positive vs. false positive. We calculate the area under the curve (AUC)
# AUC = 1 indicates a perfect classifier, AUC = 0.5 means the classifier is no better than a coin flip
from sklearn.metrics import roc_curve, roc_auc_score

# %matplotlib inline
y_pred_proba = model.predict_proba(X_test)[::,1]
falseposrate, trueposrate, _ = roc_curve(y_test, y_pred_proba)
auc = roc_auc_score(y_test, y_pred_proba)
plt.plot(falseposrate,trueposrate,label="ROC curve, auc="+str(auc))
plt.legend(loc=4)
plt.show()

# +
import numpy as np
# With tree models we can use feature importance to get some insight into root cause

# Mean feature importance across all trees
mean_i = model.feature_importances_

# Standard deviation of feature importances across all trees
st_i = np.std([tree.feature_importances_ for tree in model.estimators_],
             axis=0)

# Features
features_i = np.argsort(mean_i)[::-1]

# Print the feature ranking
print("Feature ranking:")

for f in range(X_train.shape[1]):
    print("%d. PARAM %d (%f)" % (f + 1, features_i[f]+1, mean_i[features_i[f]]))
# -

# ### Preliminary Conclusions
# It appears we can make some preliminary analysis about the root cause of the products that fall above the action limit. Next steps for such an analysis would be to investigate the physical relationship between PARAM_8 and the manufacturing process.



StarcoderdataJupyterScriptsDedupFiltered