k-prototypes for customer segmentation

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import re
!pip install venn &> /dev/null
from venn import venn
!pip install kmodes &> /dev/null
from kmodes.kprototypes import KPrototypes
import plotly.express as px
from sklearn import preprocessing

doctors = pd.read_csv('data/doctors.csv')
orders = pd.read_csv('data/orders.csv')
complaints = pd.read_csv('data/complaints.csv')
instructions = pd.read_csv('data/instructions.csv')

ndoc_doc = set(doctors.DoctorID)
ndoc_inst = set(instructions.DoctorID)
ndoc_ord = set(orders.DoctorID)
ndoc_comp = set(complaints.DoctorID)
doc_data = {
    "doctors": ndoc_doc,
    "instructions": ndoc_inst,
    "orders": ndoc_ord,
    "complaints": ndoc_comp
}

venn(doc_data);
print(f"Total of doctors on register :{len(ndoc_doc.union(ndoc_inst,ndoc_ord,ndoc_comp))}")

# Satisfaction must be float, so we will use nan so that empty values do not interfere with the average.
def isnumb(x):
    '''
    Function  is used to accept a string and convert it into a float. If the input string does not contain a numeral value or If the first character of thestring is not a number 	then it returns NaN.
    '''
    return re.sub('[-.]+','',x).isnumeric()

# Satisfaction to float 
doctors['Satisfaction'] = doctors['Satisfaction'].apply(lambda x : float(x) if isnumb(x) else np.nan)
sns.pairplot(doctors)
plt.show();

###
cor_doc = doctors.corr()
mask = np.triu(np.ones_like(cor_doc, dtype=bool))
sns.heatmap(cor_doc,
            annot = True,
            fmt='.2f',
           mask = mask);

###
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(20,6))

doctors.Category.value_counts().plot(kind='bar',title = "Count by Category", ax=axes[0])
doctors.Rank.value_counts().plot(kind='bar',title = "Count Doctors by Rank", ax=axes[1])
fig.show()

doc_by_reg = doctors.groupby(['Region']).agg({'DoctorID':'count',
                                              'Purchases': 'mean',
                                              'R rate':'mean',
                                              'Satisfaction':'mean'
                                             }).sort_values(by='DoctorID',ascending=False)
doc_by_reg.rename(columns = {'DoctorID':'Doctor'}, inplace = True)
display(doc_by_reg.head())

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(20,9))
sns.set_theme(style="whitegrid")
plt.xticks(rotation=90)
for i,k in enumerate(doc_by_reg.columns):
    axes[i//2,i%2].tick_params(labelrotation=90)
    sns.barplot(x = doc_by_reg.index, y = k,
                data = doc_by_reg, color="c",
                ax = axes[i//2,i%2]).set(title='Distribution of '+k+' for region')
fig.tight_layout()
plt.show()

print('Show doctor order by Purchase')
display(doctors.sort_values(by ='Purchases', ascending=False).head(5))

dfcond = orders.sum()[3:].sort_values(ascending  = False)
plt.figure(figsize=[10,5])
dfcond.plot(kind='bar', title = 'Oders by condition')
plt.show()

dfcond = (dfcond/dfcond.sum()).to_frame()
dfcond['cumsum'] = dfcond.cumsum()
dfcond.rename(columns={0:'Percentual'}, inplace = True)
dfcond

# Set index to DoctorID to join
doctors2 = doctors.set_index('DoctorID')

# "Qty" - number of complaints per complaint *type* per doctor. 'Complaint Type'
complaints['Complaint Type'].value_counts().plot(kind='bar')
complaints2 = complaints.pivot('DoctorID','Complaint Type','Qty').fillna(0)
complaints2 = complaints2[complaints2.columns[1:]]
complaints2.insert(0,'Total_complaint',complaints2.T.sum())

# Grouping the complaints of each type by doctor
select_cols = ['Total_complaint', 'Correct']

complaints3 = pd.concat([doctors2,complaints2], axis =1) # Join doctors and complaints

complaints_purch = complaints3[select_cols+['Purchases']]

# Pairplot for quantitative variables
sns.pairplot(complaints_purch)
plt.show();

# Heatmap of correlations
cor_compl = complaints_purch.corr()

mask = np.triu(np.ones_like(cor_compl, dtype=bool))
sns.heatmap(cor_compl,
            annot = True,
            fmt='.2f',
           mask = mask);
plt.show()

fig = px.scatter(complaints3.dropna(), y="Total_complaint", x="Purchases", color="Rank")
fig.update_layout(title_text="Relationship between the number of Total_complaint and purchases by Rank")
fig.update_traces(marker_size=10)
fig.show()

# Change yes and no values to 1 and 0, and set DoctorId to index
instructions2 = instructions.copy()
instructions2.Instructions = instructions.Instructions.map({'Yes':1,'No':0})
instructions2.set_index('DoctorID',inplace=True)

# Group orders by doctor, I will not consider the conditions, to better distribute the weights of the features in the clustering.
ord_by_doct = orders.groupby('DoctorID').agg({'OrderID':'count'})

# Join orders, doctors, complaints, instructions
# I will not consider type of complaints, to better distribute the weights of the features in the clustering
dfodc = ord_by_doct.join([complaints3.iloc[:,:9].dropna(how = 'all'),instructions2], how='outer') 

dctagg = {'OrderID':['sum','mean'],'Incidence rate':'mean',
          'R rate':'mean','Satisfaction':'mean','Experience':'mean',
          'Purchases':['sum','mean'],'Total_complaint':'mean'}

print('Specialists are the ones who make the most orders on total and average too.')
display(dfodc.groupby('Category').agg(dctagg).sort_values(by=('Purchases','mean'), ascending=False))

print('Identifying the categories that have the highest volume of purchases.')
display(dfodc.groupby('Rank').agg(dctagg).sort_values(by=[('Purchases','sum'),('Purchases','mean')], ascending=False))

fig = px.scatter(dfodc.fillna(0), y="OrderID", x="Purchases", color="Rank")
fig.update_layout(title_text="Relationship between the number of orders and purchases by doctor")
fig.update_traces(marker_size=10)
fig.show()

# As we have many categorical values, I will use k prototypes.
# K-Means is one of the conventional clustering methods commonly used in clustering techniques. However, its method is not good and suitable for data that contains categorical variables. 

# Categorical columns
catcols = ['Region','Category', 'Rank']
iloccatcols = [dfodc.columns.get_loc(k) for k in catcols]

# Numeric columns
numcols = [k for k in dfodc.columns if k not in catcols]
dfodc.fillna(0, inplace=True)

for k in catcols:
    dfodc.loc[:,k] = dfodc.loc[:,k].astype(str)

# To eliminate possible problems of scale difference between the variables, we will apply the standardization.
customers_norm = dfodc.copy()
scaler = preprocessing.StandardScaler()
# Apply standardization to numeric columns
customers_norm[numcols] = scaler.fit_transform(customers_norm[numcols])
    
costs = []
n_clust = []

# Plot for Elbow Method for best choice of number of clusters.

# Random state 0 for reproducing the problem the same every time it is run
# Selecting ‘Huang’ as the init, the model will select the first k distinct objects from the data set as initial k-modes
for i in range(3, 20):
    try:
        kproto = KPrototypes(n_clusters = i, init='Huang',random_state = 0)
        clusters = kproto.fit_predict(customers_norm, categorical = iloccatcols )
        costs.append(kproto.cost_)
        n_clust.append(i)
    except:
        pass
    
plt.figure(figsize=(8,4))
plt.plot(n_clust,costs)
plt.title('Elbow Method')
plt.xlabel('Number of clusters')
plt.ylabel('Inertia')
plt.show()

def optimal_number_of_clusters(sse,n_clust):
    x1, y1 = n_clust[0], sse[0]
    x2, y2 = n_clust[-1], sse[-1]
    distances = []
    for x0,y0 in zip(n_clust,sse):
        numerator = abs((y2-y1)*x0 - (x2-x1)*y0 + x2*y1 - y2*x1)
        denominator = np.sqrt((y2 - y1)**2 + (x2 - x1)**2)
        distances.append(numerator/denominator)
    
    return distances

distances = optimal_number_of_clusters(costs,n_clust)
n_max_dist = n_clust[distances.index(max(distances))]

plt.figure(figsize=(8,4))
plt.plot(n_clust, distances)
plt.plot([n_max_dist,n_max_dist],[0,max(distances)], label = 'Optimal number of clusters')
plt.xticks(n_clust)
plt.title('Elbow Method')
plt.xlabel('Number of clusters')
plt.ylabel('distances')
plt.legend(loc='best')
plt.show()

print("The Optimal number of clusters is :",n_max_dist)

kproto = KPrototypes( n_clusters= n_max_dist , init='Huang',random_state = 0)
clusters = kproto.fit_predict(customers_norm, categorical = iloccatcols )
#join data with labels 
labeledCustomers = dfodc.copy()
labeledCustomers['labels'] = clusters

plt.figure(figsize=(8,4))
labeledCustomers['labels'].value_counts().plot(kind='bar')
plt.title("Size of clusters")
plt.show()


# Groupby labels for mean and mode
mode = lambda x:x.value_counts().index[0]
aggfunc = list(map(lambda x: 'mean' if x not in catcols else mode,labeledCustomers.columns[:-1]))
dctag = {k:j for k,j in zip(labeledCustomers.columns,aggfunc)}
display(labeledCustomers.groupby('labels').agg(dctag).sort_values(by = 'Purchases', ascending = False))

# Groupby labels for count and total of Orders, Purchases and complaints
display(labeledCustomers.groupby('labels').agg({'Rank':'count','OrderID':'sum','Purchases':'sum','Total_complaint':'sum'}).sort_values(by = 'Purchases', ascending = False))

sns.catplot(data=labeledCustomers, x="OrderID", y="Purchases", hue="labels",height = 4 , aspect = 2.5)
plt.show()

for i,j in enumerate(catcols+['OrderID','Purchases']):
	sns.catplot(x = j, y='labels', hue='labels', data=labeledCustomers,height = 4 , aspect = 2.5)
	if j in ['Region','Purchases']:
		plt.xticks(rotation=90)
	plt.show()