Guide to manage missing data

## if necessary, install the latest version of seaborn
# %pip install -U seaborn 

# import necessary packages
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

from scipy.stats import mode
from sklearn.impute import KNNImputer
from sklearn.impute import SimpleImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

from sklearn import metrics
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_score
from sklearn.metrics import roc_curve

# display full output in Jupyter and DataBricks and not only last result
from IPython.core.interactiveshell import InteractiveShell  
InteractiveShell.ast_node_interactivity = "all"              

import warnings 
warnings.filterwarnings(action= 'ignore')

### import the dataset ###
df = pd.read_csv('https://bhciaaablob.blob.core.windows.net/featurenegineeringfiles/df_prepared.csv')
df_prep = df.drop(columns=['Unnamed: 0','city', 'experience', 'enrollee_id']) 

### recode some of the categorical variables to numerical for the purpose of this article ###

# making 'company_size'  numerical
print(df_prep['company_size'].value_counts())

cat2val_compsize_dict = { 
  '<10'       :     5, 
  '10/49'     :    30,
  '50-99'     :    75, 
  '100-500'   :   300,
  '500-999'   :   750,
  '1000-4999' :  3000,
  '5000-9999' :  7500,
  '10000+'    : 10000,
  np.nan      : np.nan }

def cat2val_compsize(value): return cat2val_compsize_dict[value]
df_prep['company_size_num'] = df_prep['company_size'].apply(cat2val_compsize)

print(df_prep['company_size_num'].value_counts())

df_prep = df_prep.drop(columns=['company_size']) 


#  making 'last_new_job' numerical
print(df_prep['last_new_job'].value_counts())

cat2val_job_dict = { 
  'never' : 0,
  '1'     : 1, 
  '2'     : 2,
  '3'     : 3, 
  '4'     : 4,
  '>4'    : 5,
  np.nan  : np.nan }

def cat2val_job(value): return cat2val_job_dict[value]
df_prep['last_new_job_num'] = df_prep['last_new_job'].apply(cat2val_job)

print(df_prep['last_new_job_num'].value_counts())

df_prep = df_prep.drop(columns=['last_new_job']) 

# making 'education_level' numerical
print(df_prep['education_level'].value_counts())

cat2val_edu_dict = { 
  'Primary School' : 1,
  'High School'    : 2, 
  'Graduate'       : 3,
  'Masters'        : 4, 
  'Phd'            : 5,
  np.nan           : np.nan }

def cat2val_edu(value): return cat2val_edu_dict[value]
df_prep['education_level_num'] = df_prep['education_level'].apply(cat2val_edu)

print(df_prep['education_level_num'].value_counts())

df_prep = df_prep.drop(columns=['education_level']) 

### make lists of variables types to use later in the analysis ###
categorical_features = df_prep.select_dtypes(exclude="number").columns.tolist()  
categorical_features.remove('city name')                                          # drop city, this variable has too much levels

numeric_features = df_prep.select_dtypes(include="number").columns.tolist()
numeric_features.remove('target')                                                 # target is the dependent variable and not one of the numerical independent/predictor variables

### split the data in a train and test dataset ###

# Define the target vector y
y = df_prep['target']  
X = df_prep.drop(columns=['target']) 

# split X and y into training and testing sets
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.20,random_state=42)

# basic exploration: what does the dataset look like? Do you see missing values?
X_train.head(20)

# which variables have the highest number of missing values?
x = X_train.isna().sum().sort_values(ascending=False).reset_index()
x = x.rename(columns={'index':'name_var',0:'# missings'})

n_records = df_prep.shape[0]
x['% missings']=x['# missings']/n_records
x['% missings'] = x['% missings'].astype(float).map("{:.2%}".format)  

# print the columns with missing values
print(x[x['# missings'] > 0])
print('\n')    

nums_with_miss = [i for i in X_train.columns if X_train[i].isnull().sum()>0]
print('Number of variables with missing values: ',len(nums_with_miss))

# join X and y back together
Xy_train = X_train.join(y_train)

# drop all records with missing values
Xy_train_drop_all_miss = Xy_train.dropna()

# original dataset has 19.158 records
print('Original train dataset has', len(Xy_train), 'records.')
# after dropping records with missing values, the dataset only has 8.955 records left  
print('Dataset after dropping all rows with missing values has', len(Xy_train_drop_all_miss), 'records.')
print('This is only', round(100*len(Xy_train_drop_all_miss)/len(Xy_train),1), '% of the original dataset.')

print('\n')
print('The original train dataset has ',round(Xy_train['target'].mean()*100,1),'% of data scientists looking for a new job.')
print('The dataset after dropping all rows with missing values has ',round(Xy_train_drop_all_miss['target'].mean()*100,1),'% of data scientists looking for a new job.')

# basic exploration: describe distributions of variables in the dataset
def describe(df):
    return round(pd.concat([df.describe().T, 
                            df.skew().rename('skewness'), 
                            df.isnull().sum().rename('# missings'),
                            df.isnull().mean().rename('% missings').map("{:.2%}".format),
                            df[categorical_features].mode().T.rename(columns={0:'mode'})]
                           ,  axis=1).T,1)

describe(X_train)

# plot numerical values with missing values 
# how symmetric or skewed are these variables? Can missings be replaced by the mean or median?
for variable in numeric_features:
   if df_prep[variable].isna().sum() > 0:    # only variables with missing values
      x1 = sns.histplot(X_train,x=variable)
      plt.show(x1)      
      x2 = sns.boxplot(X_train[variable])
      plt.show(x2)

# one-hot encoding for the categorical variables - only keep dummies for categories with more than 100 observations 

# fit and transform onehot-encoding based on the trainingsdataset
one_hot_encode = OneHotEncoder(sparse=False)
one_hot_encode.fit(X_train[categorical_features])
X_train_dummies = pd.DataFrame(one_hot_encode.transform(X_train[categorical_features]))
X_train_dummies.columns = one_hot_encode.get_feature_names(categorical_features)

# only keep the dummy variable for categories that occur more than 100 times (including the category 'nan') 
category_list_keep = [dummy for dummy in X_train_dummies.columns if X_train_dummies[dummy].sum() >= 100]
X_train_dummies = X_train_dummies[category_list_keep]

# apply the onehot-encoding by using transform on the testdataset 
X_test_dummies = pd.DataFrame(one_hot_encode.transform(X_test[categorical_features]))
X_test_dummies.columns = one_hot_encode.get_feature_names(categorical_features)
X_test_dummies = X_test_dummies[category_list_keep]

# One-hot encoding removed the index - with this it is put back! 
X_train_dummies.index = X_train.index 
X_test_dummies.index = X_test.index 

# Printing the resulting dataframe
X_train_dummies

# create a list which numerical variables that have more than x missings 
def list_missing_vars(df, features, x_missings = 100):
    list_vars = []
    for variable in features:
        if df[variable].isna().sum() >= x_missings:
           list_vars.append(variable)
    return list_vars

# for every variable in the list, create indicator variables: 1 if value is missing and 0 if value is not missing
def create_missing_indicators(df, variable_list):
    df_new = pd.DataFrame()  
    for variable in variable_list:
        df_new['miss_'+variable] = df[variable].isnull().astype(int)   
    return df_new
    
# create a list which numerical variables have more than 100 missings values - based on the trainingsdataset
list_indicators_ge100 = list_missing_vars(X_train, numeric_features, 100) 

# create indicator variables for the train 
X_train_miss = create_missing_indicators(X_train, list_indicators_ge100)
# create the same list of indicator variables on the test dataset 
X_test_miss = create_missing_indicators(X_test, list_indicators_ge100)

# Printing the resulting dataframe
X_train_miss

# list of variables with missings in the trainingsdataset
nums_with_miss = [i for i in numeric_features if X_train[i].isnull().sum()>0]
#print(nums_with_miss)

# calculate a dictionary of pairs of variable-names and the replacement values - based on the trainingsdataset
dict_replace = {'experience_num'         : X_train['experience_num'].mean(),         
                'company_size_num'       : X_train['company_size_num'].median(),     
                'last_new_job_num'       : 0,                                        
                'education_level_num'    : X_train['education_level_num'].mode()[0],  
                'city_development_index' : X_train['city_development_index'].mean(), 
                'training_hours'         : X_train['training_hours'].mean(),
                'ind_relevent_experience': X_train['ind_relevent_experience'].mean()}

# define a function to replace the missing values with the replacement values in the dictionary 
def impute_values(df, dict_impute):
    df_new = pd.DataFrame()  
    for variablename, value in dict_impute.items():
        df_new[variablename] = df[variablename].fillna(value)
    return df_new 

# replace the missings of the train and test dataset with replacement values of the dictionary 
X_train_miss_imp = impute_values(X_train[numeric_features], dict_replace)  
X_test_miss_imp = impute_values(X_test[numeric_features], dict_replace)  

# Printing the resulting dataframe
X_train_miss_imp

# confirm that there are no missings left in the data
nums_with_miss = [i for i in X_train_miss_imp.columns if X_train_miss_imp[i].isnull().sum()>0]
print('number of missing values: ',len(nums_with_miss))

# the prepared data consists of three parts: the numerical variables with missings imputed, missing indicator variables and dummies for categorical variables.  
X_train_hand = X_train_miss_imp.join(X_train_miss).join(X_train_dummies)
X_test_hand = X_test_miss_imp.join(X_test_miss).join(X_test_dummies)

# Use default logistic regression to fit the model on the trainingsdata
logreg = LogisticRegression()
logreg.fit(X_train_hand,y_train)

# Predict on the test dataset:
y_pred_dummy = logreg.predict(X_test_hand)
y_pred_proba_dummy = logreg.predict_proba(X_test_hand)[::,1]

# Evaluate the model with the use of cv on the test dataset:
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) 
scores = cross_val_score(logreg, X_test_hand, y_test, cv=cv, scoring = 'roc_auc')
print("roc_auc = %f (%f)" % (scores.mean(), scores.std()))

### replace missings with the MEAN value ###

# first initialize the dataframes
X_train_mean = pd.DataFrame() 
X_train_mean = X_train[numeric_features]
X_test_mean = pd.DataFrame() 
X_test_mean = X_test[numeric_features]

# Use the default SimpleImputer to fit the model on the trainingsdata
imputer = SimpleImputer(strategy='mean')
mean_imputer_fit = imputer.fit(X_train_mean)

# use the fitted imputer to replace missings by the mean and make indicator variables (=transform) for the trainingsdata and test dataset.
meanimputed_transform = mean_imputer_fit.transform(X_train_mean)
X_train_mean_imp = pd.DataFrame(meanimputed_transform, columns=X_train_mean.columns) 
meanimputed_transform = mean_imputer_fit.transform(X_test_mean)
X_test_mean_imp = pd.DataFrame(meanimputed_transform, columns=X_test_mean.columns) 

# One-hot encoding removed the index - with this it is put back! 
X_train_mean_imp.index = X_train.index 
X_test_mean_imp.index = X_test.index 

# Printing the resulting dataframe
X_train_mean_imp

# confirm that there are no missings left in the data
nums_with_miss = [i for i in X_train_mean_imp.columns if X_train_mean_imp[i].isnull().sum()>0]
print('number of missing values: ',len(nums_with_miss))

# join the dummy-categories, missing indicator and imputed numerical variables
X_train_mean2 = X_train_mean_imp.join(X_train_miss).join(X_train_dummies)
X_test_mean2 = X_test_mean_imp.join(X_test_miss).join(X_test_dummies)

# Use default logistic regression to fit the model on the trainingsdata
logreg = LogisticRegression()
logreg.fit(X_train_mean2,y_train)

# Predict on the test dataset:
y_pred_dummy=logreg.predict(X_test_mean2)
y_pred_proba_dummy = logreg.predict_proba(X_test_mean2)[::,1]

# Evaluate the model with the use of cv on the test dataset:
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) 
scores = cross_val_score(logreg, X_test_mean2, y_test, cv=cv, scoring = 'roc_auc')
print("roc_auc = %f (%f)" % (scores.mean(), scores.std()))

### replace missings with the MEDIAN value ###

# first initialize the dataframes
X_train_median = pd.DataFrame() 
X_train_median = X_train[numeric_features]
X_test_median = pd.DataFrame() 
X_test_median = X_test[numeric_features]

# Use the default SimpleImputer to fit the model on the trainingsdata
imputer = SimpleImputer(strategy='median')
median_imputer_fit = imputer.fit(X_train_median)

# use the fitted imputer to replace missings by the median and make indicator variables (=transform) for the trainingsdata and test dataset.
medianimputed_transform = median_imputer_fit.transform(X_train_median)
X_train_median_imp = pd.DataFrame(medianimputed_transform, columns=X_train_median.columns) 
medianimputed_transform = median_imputer_fit.transform(X_test_median)
X_test_median_imp = pd.DataFrame(medianimputed_transform, columns=X_test_median.columns) 

# One-hot encoding removed the index - with this it is put back! 
X_train_median_imp.index = X_train.index 
X_test_median_imp.index = X_test.index 

# Printing the resulting dataframe
X_train_median_imp

# confirm that there are no missings left in the data
nums_with_miss = [i for i in X_train_median_imp.columns if X_train_median_imp[i].isnull().sum()>0]
print('number of missing values: ',len(nums_with_miss))

# join the dummy-categories, missing indicator and imputed numerical variables
X_train_median2 = X_train_median_imp.join(X_train_miss).join(X_train_dummies)
X_test_median2 = X_test_median_imp.join(X_test_miss).join(X_test_dummies)

# Use default logistic regression to fit the model on the trainingsdata
logreg = LogisticRegression()
logreg.fit(X_train_median2,y_train)

# Predict on the test dataset:
y_pred_dummy=logreg.predict(X_test_median2)
y_pred_proba_dummy = logreg.predict_proba(X_test_median2)[::,1]

# Evaluate the model with the use of cv on the test dataset:
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) 
scores = cross_val_score(logreg, X_test_median2, y_test, cv=cv, scoring = 'roc_auc')
print("roc_auc = %f (%f)" % (scores.mean(), scores.std()))

### replace missings with the MOST_FREQUENT value ###

# first initialize the dataframes
X_train_most_frequent = pd.DataFrame() 
X_train_most_frequent = X_train[numeric_features]
X_test_most_frequent = pd.DataFrame() 
X_test_most_frequent = X_test[numeric_features]

# Use the default SimpleImputer to fit the model on the trainingsdata
imputer = SimpleImputer(strategy='most_frequent')
most_frequent_imputer_fit = imputer.fit(X_train_most_frequent)

# use the fitted imputer to replace missings by the most_frequent and make indicator variables (=transform) for the trainingsdata and test dataset.
most_frequentimputed_transform = most_frequent_imputer_fit.transform(X_train_most_frequent)
X_train_most_frequent_imp = pd.DataFrame(most_frequentimputed_transform, columns=X_train_most_frequent.columns) 
most_frequentimputed_transform = most_frequent_imputer_fit.transform(X_test_most_frequent)
X_test_most_frequent_imp = pd.DataFrame(most_frequentimputed_transform, columns=X_test_most_frequent.columns) 

# One-hot encoding removed the index - with this it is put back! 
X_train_most_frequent_imp.index = X_train.index 
X_test_most_frequent_imp.index = X_test.index 

# Printing the resulting dataframe
X_train_most_frequent_imp

# confirm that there are no missings left in the data
nums_with_miss = [i for i in X_train_most_frequent_imp.columns if X_train_most_frequent_imp[i].isnull().sum()>0]
print('number of missing values: ',len(nums_with_miss))

# join the dummy-categories, missing indicator and imputed numerical variables
X_train_mode = X_train_most_frequent_imp.join(X_train_miss).join(X_train_dummies)
X_test_mode = X_test_most_frequent_imp.join(X_test_miss).join(X_test_dummies)

# Use default logistic regression to fit the model on the trainingsdata
logreg = LogisticRegression()
logreg.fit(X_train_mode,y_train)

# Predict on the test dataset:
y_pred_dummy=logreg.predict(X_test_mode)
y_pred_proba_dummy = logreg.predict_proba(X_test_mode)[::,1]

# Evaluate the model with the use of cv on the test dataset:
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) 
scores = cross_val_score(logreg, X_test_mode, y_test, cv=cv, scoring = 'roc_auc')
print("roc_auc = %f (%f)" % (scores.mean(), scores.std()))

### replace missings with the K-means Nearest Neighbour ###

# first replace the categoricals by dummies, but leave the numerical variables with missings --> these will be imputed per group with KNN-means
X_train_KNN = pd.DataFrame() 
X_train_KNN = X_train[numeric_features].join(X_train_miss).join(X_train_dummies)
X_test_KNN = pd.DataFrame() 
X_test_KNN = X_test[numeric_features].join(X_test_miss).join(X_test_dummies)

# Use the default KNNImputer to fit the model on the trainingsdata
imputer = KNNImputer(n_neighbors=5)
KNNimputer_fit = imputer.fit(X_train_KNN)

# use the fitted KNNImputer to transform (=replace missings) the trainingsdata and test dataset.
KNNimputed_transform = KNNimputer_fit.transform(X_train_KNN)
X_train_KNN_imp = pd.DataFrame(KNNimputed_transform, columns=X_train_KNN.columns) 
KNNimputed_transform = KNNimputer_fit.transform(X_test_KNN)
X_test_KNN_imp = pd.DataFrame(KNNimputed_transform, columns=X_test_KNN.columns) 

# Printing the resulting dataframe
X_train_KNN_imp

# confirm that there are no missings left in the data
nums_with_miss = [i for i in X_train_KNN_imp.columns if X_train_KNN_imp[i].isnull().sum()>0]
print('number of missing values: ',len(nums_with_miss))

# Use default logistic regression to fit the model on the trainingsdata
logreg = LogisticRegression()
logreg.fit(X_train_KNN_imp,y_train)

# Predict on the test dataset:
y_pred_dummy=logreg.predict(X_test_KNN_imp)
y_pred_proba_dummy = logreg.predict_proba(X_test_KNN_imp)[::,1]

# Evaluate the model with the use of cv on the test dataset:
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) 
scores = cross_val_score(logreg, X_test_KNN_imp, y_test, cv=cv, scoring = 'roc_auc')
print("roc_auc = %f (%f)" % (scores.mean(), scores.std()))

# first replace the categoricals by dummies, but leave the numerical variables with missings --> these will be imputed per group with KNN-means
X_train_iterative = pd.DataFrame() 
X_test_iterative = pd.DataFrame() 
X_train_iterative = X_train[numeric_features].join(X_train_miss).join(X_train_dummies)
X_test_iterative = X_test[numeric_features].join(X_test_miss).join(X_test_dummies)

# Use the default IterativeImputer to fit the model on the trainingsdata
IterativeImp = IterativeImputer(max_iter=5, verbose=0)
IterativeImp.fit(X_train_iterative)

# use the fitted IterativeImputer to transform (=replace missings) the trainingsdata and test dataset.
df_iterative_imputed = IterativeImp.transform(X_train_iterative)
X_train_iterative_imp = pd.DataFrame(df_iterative_imputed, columns=X_train_iterative.columns)

df_iterative_imputed = IterativeImp.transform(X_test_iterative)
X_test_iterative_imp = pd.DataFrame(df_iterative_imputed, columns=X_test_iterative.columns)

# Printing the resulting dataframe
X_train_iterative_imp

# confirm that there are no missings left in the data
nums_with_miss = [i for i in X_train_iterative_imp.columns if X_train_iterative_imp[i].isnull().sum()>0]
print('number of missing values: ',len(nums_with_miss))

# Instantiate the model (using the default parameters)
logreg = LogisticRegression()
 
# Fit the model with data
logreg.fit(X_train_iterative_imp,y_train)

#Predict:
y_pred_dummy=logreg.predict(X_test_iterative_imp)
y_pred_proba_dummy = logreg.predict_proba(X_test_iterative_imp)[::,1]

# Evaluate the model with the use of cv on the test dataset:
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42) 
scores = cross_val_score(logreg, X_test_iterative_imp, y_test, cv=cv, scoring = 'roc_auc')
print("roc_auc = %f (%f)" % (scores.mean(), scores.std()))