
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.gridspec as gridspec
import colorcet as cc
import seaborn as sns
from datetime                   import datetime
from scipy.stats                import pearsonr
from scipy.stats                import pointbiserialr
from scipy.stats                import chi2_contingency
from scipy.stats                import kruskal

# Frequent pattern mining  
from mlxtend.preprocessing      import TransactionEncoder
from mlxtend.frequent_patterns  import apriori, association_rules

# Data processing
from sklearn.preprocessing      import LabelEncoder, OneHotEncoder

# Correlation Matrix
from scipy                      import stats

# Machine learning
import sklearn
from sklearn.cluster            import KMeans, AgglomerativeClustering
from sklearn.neighbors          import KNeighborsClassifier
from sklearn.linear_model       import LogisticRegression
import xgboost

# Imbalanced learning
from imblearn.over_sampling     import SMOTE
from imblearn.under_sampling    import RandomUnderSampler
from imblearn.pipeline          import make_pipeline

# Model selection and evaluation
from sklearn                    import metrics
from sklearn.metrics            import classification_report, accuracy_score, cohen_kappa_score, f1_score, recall_score, precision_score, confusion_matrix, roc_curve
from sklearn.metrics            import mean_squared_error, mean_absolute_error, r2_score
from sklearn.model_selection    import train_test_split

# Suppress all warnings
import warnings
warnings.filterwarnings('ignore')


# constant

RANDOM_STATE = 42


# calculate months between 2 dates

def calculate_month_interval(start_time, end_time):
    return (end_time.year - start_time.year)*12 + (end_time.month - start_time.month)


def cf_matrix_labels(cf_matrix):
    
    group_names = ['True Negatives','False Positives','False Negatives','True Positives']
    group_counts = ['{0:0.0f}'.format(value) for value in cf_matrix.flatten()]
    group_percentages = ['{0:.2%}'.format(value) for value in cf_matrix.flatten()/np.sum(cf_matrix)]
    labels = [f'{v1}\n{v2}\n{v3}' for v1, v2, v3 in zip(group_names,group_counts,group_percentages)]
    labels = np.asarray(labels).reshape(2,2)
    
    return labels


# customize classification report

def get_classification_report(y_true, proba_pred):
    precision, recall, threshold = metrics.precision_recall_curve(y_true, proba_pred)
    report_df = pd.DataFrame({
        'threshold': threshold,
        'precision': precision[:-1],
        'recall': recall[:-1]
    }).assign(
        f1 = lambda x: 2 * x['recall'] * x['precision'] / (x['recall'] + x['precision']),
        f2 = lambda x: (5 * x['precision'] * x['recall']) / (4 * x['precision'] + x['recall'])
    )
    fig = plt.figure(figsize=(7,7))
    plt.plot(report_df['threshold'], report_df['precision'], label = 'precision')
    plt.plot(report_df['threshold'], report_df['recall'], label = 'recall')
    # plt.plot(report_df['threshold'], report_df['f1'], '--', color='brown', label = 'f1')
    plt.plot(report_df['threshold'], report_df['f2'], '--', color='green', label = 'f2')
    plt.xlabel('threshold')
    plt.ylabel('score')
    plt.legend(loc='upper center')
    plt.show()
    return report_df


# put all data files in folder "dataset"

DATA_LOC = 'dataset'


loyalty_file = 'Customer Loyalty History.csv'
loyalty_df = pd.read_csv(os.path.join(DATA_LOC, loyalty_file))
loyalty_df.columns

Index(['Loyalty Number', 'Country', 'Province', 'City', 'Postal Code',
       'Gender', 'Education', 'Salary', 'Marital Status', 'Loyalty Card',
       'CLV', 'Enrollment Type', 'Enrollment Year', 'Enrollment Month',
       'Cancellation Year', 'Cancellation Month'],
      dtype='object')


# remove all space character in column name

loyalty_new_cols = [col.replace(' ', '') for col in loyalty_df.columns]
loyalty_df = loyalty_df.rename(columns=dict(zip(loyalty_df.columns, loyalty_new_cols)))

display(loyalty_df.info(), loyalty_df.head())
print('No. of unique customers:', loyalty_df['LoyaltyNumber'].nunique())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 16737 entries, 0 to 16736
Data columns (total 16 columns):
 #   Column             Non-Null Count  Dtype  
---  ------             --------------  -----  
 0   LoyaltyNumber      16737 non-null  int64  
 1   Country            16737 non-null  object 
 2   Province           16737 non-null  object 
 3   City               16737 non-null  object 
 4   PostalCode         16737 non-null  object 
 5   Gender             16737 non-null  object 
 6   Education          16737 non-null  object 
 7   Salary             12499 non-null  float64
 8   MaritalStatus      16737 non-null  object 
 9   LoyaltyCard        16737 non-null  object 
 10  CLV                16737 non-null  float64
 11  EnrollmentType     16737 non-null  object 
 12  EnrollmentYear     16737 non-null  int64  
 13  EnrollmentMonth    16737 non-null  int64  
 14  CancellationYear   2067 non-null   float64
 15  CancellationMonth  2067 non-null   float64
dtypes: float64(4), int64(3), object(9)
memory usage: 2.0+ MB

None

No. of unique customers: 16737


flight_file = 'Customer Flight Activity.csv'
flight_df = pd.read_csv(os.path.join(DATA_LOC, flight_file))
flight_df.columns

Index(['Loyalty Number', 'Year', 'Month', 'Flights Booked',
       'Flights with Companions', 'Total Flights', 'Distance',
       'Points Accumulated', 'Points Redeemed', 'Dollar Cost Points Redeemed'],
      dtype='object')


# remove all space character in column name

flight_new_cols = [col.replace(' ', '') for col in flight_df.columns]
flight_df = flight_df.rename(columns=dict(zip(flight_df.columns, flight_new_cols)))\
    .rename(columns={'FlightswithCompanions': 'FlightsWithCompanions'})
    
display(flight_df.info(), flight_df.head())
print('No. of unique customers:', flight_df['LoyaltyNumber'].nunique())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 405624 entries, 0 to 405623
Data columns (total 10 columns):
 #   Column                    Non-Null Count   Dtype  
---  ------                    --------------   -----  
 0   LoyaltyNumber             405624 non-null  int64  
 1   Year                      405624 non-null  int64  
 2   Month                     405624 non-null  int64  
 3   FlightsBooked             405624 non-null  int64  
 4   FlightsWithCompanions     405624 non-null  int64  
 5   TotalFlights              405624 non-null  int64  
 6   Distance                  405624 non-null  int64  
 7   PointsAccumulated         405624 non-null  float64
 8   PointsRedeemed            405624 non-null  int64  
 9   DollarCostPointsRedeemed  405624 non-null  int64  
dtypes: float64(1), int64(9)
memory usage: 30.9 MB

None

No. of unique customers: 16737


# convert data type

loyalty_df = loyalty_df.assign(
    LoyaltyNumber = lambda x: x['LoyaltyNumber'].astype(str),
    CancellationYear = lambda x: x['CancellationYear'].astype('Int64'),
    CancellationMonth = lambda x: x['CancellationMonth'].astype('Int64')
)


loyalty_df.describe()


# records with negative salary

print('Number of records with negative salary:', loyalty_df[loyalty_df['Salary'].le(0)].shape[0])
loyalty_df[loyalty_df['Salary'].le(0)]

Number of records with negative salary: 20


# records with missing salary

print(
    'Number of records with missing salary:', loyalty_df['Salary'].isna().sum(),
    '\nPercentage of records with missing salary:', round(loyalty_df['Salary'].isna().sum()*100/loyalty_df.shape[0],2), '%'
)

Number of records with missing salary: 4238 
Percentage of records with missing salary: 25.32 %


# treat negative salary as missing values, fill all missing salary with 0

loyalty_df['Salary'] = loyalty_df['Salary'].fillna(0).where(loyalty_df['Salary'].gt(0), 0)
print('Number of records with salary = 0:', loyalty_df['Salary'].eq(0).sum())

Number of records with salary = 0: 4258


# EnrollmentDate (type: datetime), EnrollmentYearMonth (type: string, format: YYYY-MM)
# CancellationDate (type: datetime), CancellationYearMonth (type: string, format: YYYY-MM)

time_membership = ['Enrollment', 'Cancellation']

for i in time_membership: 
    loyalty_df[f'{i}Date'] = pd.to_datetime(
        loyalty_df[[f'{i}Year', f'{i}Month']].astype('Float64')\
        .rename(columns={f'{i}Year': 'YEAR', f'{i}Month': 'MONTH'}).assign(DAY=1)
    )
    loyalty_df[f'{i}YearMonth'] = loyalty_df[f'{i}Date'].dt.strftime('%Y-%m')
    
loyalty_df.head()


loyalty_df['ChurnIndicator'] = np.where(loyalty_df['CancellationMonth'].isna().eq(False), 1, 0)


loyalty_df.head()


# convert data type

flight_df['LoyaltyNumber'] = flight_df['LoyaltyNumber'].astype(str)


flight_df.describe()


flight_df['Date'] = pd.to_datetime(flight_df[['Year', 'Month']].assign(DAY=1))
flight_df['YearMonth'] = flight_df['Date'].dt.strftime('%Y-%m')


# loyalty history
display(loyalty_df.info(), loyalty_df.head(), loyalty_df.describe())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 16737 entries, 0 to 16736
Data columns (total 21 columns):
 #   Column                 Non-Null Count  Dtype         
---  ------                 --------------  -----         
 0   LoyaltyNumber          16737 non-null  object        
 1   Country                16737 non-null  object        
 2   Province               16737 non-null  object        
 3   City                   16737 non-null  object        
 4   PostalCode             16737 non-null  object        
 5   Gender                 16737 non-null  object        
 6   Education              16737 non-null  object        
 7   Salary                 16737 non-null  float64       
 8   MaritalStatus          16737 non-null  object        
 9   LoyaltyCard            16737 non-null  object        
 10  CLV                    16737 non-null  float64       
 11  EnrollmentType         16737 non-null  object        
 12  EnrollmentYear         16737 non-null  int64         
 13  EnrollmentMonth        16737 non-null  int64         
 14  CancellationYear       2067 non-null   Int64         
 15  CancellationMonth      2067 non-null   Int64         
 16  EnrollmentDate         16737 non-null  datetime64[ns]
 17  EnrollmentYearMonth    16737 non-null  object        
 18  CancellationDate       2067 non-null   datetime64[ns]
 19  CancellationYearMonth  2067 non-null   object        
 20  ChurnIndicator         16737 non-null  int32         
dtypes: Int64(2), datetime64[ns](2), float64(2), int32(1), int64(2), object(12)
memory usage: 2.6+ MB

None


# flight activity
display(flight_df.info(), flight_df.head(), flight_df.describe())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 405624 entries, 0 to 405623
Data columns (total 12 columns):
 #   Column                    Non-Null Count   Dtype         
---  ------                    --------------   -----         
 0   LoyaltyNumber             405624 non-null  object        
 1   Year                      405624 non-null  int64         
 2   Month                     405624 non-null  int64         
 3   FlightsBooked             405624 non-null  int64         
 4   FlightsWithCompanions     405624 non-null  int64         
 5   TotalFlights              405624 non-null  int64         
 6   Distance                  405624 non-null  int64         
 7   PointsAccumulated         405624 non-null  float64       
 8   PointsRedeemed            405624 non-null  int64         
 9   DollarCostPointsRedeemed  405624 non-null  int64         
 10  Date                      405624 non-null  datetime64[ns]
 11  YearMonth                 405624 non-null  object        
dtypes: datetime64[ns](1), float64(1), int64(8), object(2)
memory usage: 37.1+ MB

None


# Number of months of membership

loyalty_df = loyalty_df.assign(
    ObservationDate = lambda x: x['CancellationDate'].where(x['CancellationDate'].isna().eq(False), flight_df['Date'].max()),
    MonthsOfMembership = lambda x: x.apply(lambda z: calculate_month_interval(z['EnrollmentDate'], z['ObservationDate']), axis=1)
).drop(columns='ObservationDate')


# Aggregation on flight activity (sum)

flight_cols = [col for col in flight_df.columns if flight_df[col].dtype in ['int64', 'float64']][2:]
temp_total_flight = flight_df.groupby(by='LoyaltyNumber', sort=False)[flight_cols].sum()
temp_total_flight.columns = [f'Total{col}' for col in flight_cols]
temp_total_flight


# Percentage of Redeemed Points / Accumulated Points

temp_total_flight['PercentagePointsRedeemed'] = temp_total_flight['TotalPointsRedeemed'].where(
    temp_total_flight['TotalPointsRedeemed'].eq(0),
    temp_total_flight['TotalPointsRedeemed'].div(temp_total_flight['TotalPointsAccumulated']).round(2)
)
temp_total_flight['PercentagePointsRedeemed'].describe()

count    16737.000000
mean         0.231936
std          0.277969
min          0.000000
25%          0.000000
50%          0.190000
75%          0.350000
max          5.200000
Name: PercentagePointsRedeemed, dtype: float64


# Aggregation on flight activity (Monthly average)

temp_monthly_flight = temp_total_flight[['TotalTotalFlights', 'TotalPointsAccumulated']].merge(
    loyalty_df[['LoyaltyNumber', 'MonthsOfMembership']].set_index('LoyaltyNumber'),
    how='left', left_index=True, right_index=True
).assign(
    MonthlyTotalFlights = lambda x: x['MonthsOfMembership'].where(x['MonthsOfMembership'].eq(0), 
                                                                  x['TotalTotalFlights'].div(x['MonthsOfMembership']).round(2)), 
    MonthlyPointsAccumulated = lambda x: x['MonthsOfMembership'].where(x['MonthsOfMembership'].eq(0), 
                                                                       x['TotalPointsAccumulated'].div(x['MonthsOfMembership']).round(2))
).iloc[:, 3:]

temp_monthly_flight.describe()


loyalty_df = loyalty_df.merge(temp_total_flight.reset_index(), how='left', on='LoyaltyNumber')\
    .merge(temp_monthly_flight.reset_index(), how='left', on='LoyaltyNumber')


loyalty_df.head()


flight_df = flight_df.sort_values(by=['Date', 'LoyaltyNumber'])#.set_index('LoyaltyNumber')
TIME_INDEX = pd.DatetimeIndex(flight_df['Date'])


# Aggregation of flight stats for 1 year, 6 months, 3 months and 1 month lookback

lookback_period = {
    '1Y': '366D',
    '6M': '183D',
    '3M': '93D',
    '1M': '31D',
}

temp_lookback_flight_dfs = []

for last, day in lookback_period.items():
    df = flight_df.groupby(by='LoyaltyNumber', sort=False)[flight_cols]\
        .rolling(window=day, on=TIME_INDEX, closed='left').sum()\
        .reset_index().sort_values(by=['Date', 'LoyaltyNumber'])\
        .drop(columns='level_1').set_index(['LoyaltyNumber', 'Date'])
    df.columns = [f'Last{last}_{col}' for col in flight_cols]
    temp_lookback_flight_dfs.append(df)

lookback_flight_df = flight_df[['LoyaltyNumber', 'Date']]

for df in temp_lookback_flight_dfs:
    lookback_flight_df = lookback_flight_df.merge(df.reset_index(), how='left', on=['LoyaltyNumber', 'Date'])

lookback_flight_df


# Months of membership

temp_mom_lookback = flight_df[['LoyaltyNumber', 'Date']].merge(loyalty_df[['LoyaltyNumber', 'EnrollmentDate', 'CancellationDate']], how='left', on='LoyaltyNumber')\
    .assign(MonthsOfMembership = lambda x: x.apply(lambda z: calculate_month_interval(z['EnrollmentDate'], z['Date']), axis=1))
    
temp_mom_lookback


cust_eda_loyalty_df = loyalty_df.copy()


# Confirm that dataset only contains Canadian customers
cust_eda_loyalty_df["Country"].value_counts()

Country
Canada    16737
Name: count, dtype: int64


cust_eda_loyalty_df['ChurnIndicator'] = cust_eda_loyalty_df['ChurnIndicator'].astype("category")


flights_summary = cust_eda_loyalty_df.groupby('MaritalStatus').agg(
    {'TotalTotalFlights':['sum'],
     'TotalFlightsWithCompanions': ['sum']}
).reset_index()

# Calculate the percentage of flights with companions
flights_summary['PercentageWithCompanions'] = (flights_summary['TotalFlightsWithCompanions'] / flights_summary['TotalTotalFlights']) * 100

# Display the summary DataFrame
print(flights_summary)

  MaritalStatus TotalTotalFlights TotalFlightsWithCompanions  \
                              sum                        sum   
0      Divorced            311069                      62289   
1       Married           1213745                     243295   
2        Single            562875                     112941   

  PercentageWithCompanions  
                            
0                20.024175  
1                20.044985  
2                20.065023


# Distribution of salary 

plt.figure(figsize=(6, 9))
sns.boxplot(data=cust_eda_loyalty_df, y='Salary', palette='viridis')
plt.title('Salary Distribution')
plt.ylabel('Salary')
plt.xticks(rotation=45)
plt.show()


# Check distribution by gender

x = 'Gender'
label = 'Gender'
order = ["Female","Male"]
rows, cols = 5, 3

# Calculate total flights booked by gender
total_flights = cust_eda_loyalty_df.groupby(x)['TotalTotalFlights'].sum().reset_index()
total_flights.Gender = total_flights.Gender.astype("category")
total_flights.Gender = total_flights.Gender.cat.set_categories(order)
total_flights.sort_values([x])

# Calculate proportion of churned customers by education level
churn_proportions = cust_eda_loyalty_df.groupby(x)['ChurnIndicator'].value_counts(normalize=True).unstack().fillna(0)
churn_proportions.columns = ['Not Churn', 'Churn']
churn_proportions['Not Churn'] *= 100
churn_proportions['Churn'] *= 100
churn_proportions.sort_values('Churn', ascending=False, inplace=True)

# Focus on months of membership between 0 to 12 months
mom12 = cust_eda_loyalty_df[(loyalty_df['MonthsOfMembership'] >= 0) & (loyalty_df['MonthsOfMembership'] <= 12)]
mom12[x] = pd.Categorical(mom12[x], categories=order, ordered=True)
mom12 = mom12.sort_values(x)

# Create custom colour mapping to each gender
colour =  {'Female': 'salmon', 'Male': 'skyblue'}

# Define plot area
fig, axs = plt.subplots(nrows = rows, ncols = cols, figsize=(20,25))

# Create plots
plt.subplot(rows, cols, 1)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue=x, palette=colour)
plt.title(f'Frequency of Customers by {label}')
plt.xlabel(label)
plt.ylabel('Count')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 2)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='Salary', hue=x, palette=colour)
plt.title(f'Distribution of Salary by {label}')
plt.xlabel(label)
plt.ylabel('Salary')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 3)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='CLV', hue=x, palette=colour)
plt.title(f'Distribution of CLV by {label}')
plt.xlabel(label)
plt.ylabel('CLV')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 4)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='MonthsOfMembership', hue=x, palette=colour)
plt.title(f'Distribution of Membership Length (Months) by {label}')
plt.xlabel(label)
plt.ylabel('Membership Length (Months)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 5)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsBooked', hue=x, palette=colour)
plt.title(f'Distribution of Individual Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Individual Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 6)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsWithCompanions', hue=x, palette=colour)
plt.title(f'Distribution of Flights with Companions by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights with Companions Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 7)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalTotalFlights', hue=x, palette=colour)
plt.title(f'Distribution of Total Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 8)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalDistance', hue=x, palette=colour)
plt.title(f'Distribution of Total Distance Travelled by {label}')
plt.xlabel(label)
plt.ylabel('Total Distance')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 9)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='PercentagePointsRedeemed', hue=x, order=order, palette=colour)
plt.title(f'Distribution of Points Redeemed by {label}')
plt.xlabel(label)
plt.ylabel('Points Redeemed (%)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 10)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue='ChurnIndicator', palette="viridis")
plt.title(f'Frequency of Churn by {label}')
plt.xlabel(label)
plt.ylabel('Churn Frequency')
plt.legend(title='Churn Indicator', loc='upper right')

churn_proportions.plot(kind='bar', stacked=True, ax = axs[3,1], color=['steelblue', 'mediumseagreen'])
axs[3,1].set(title = f'Proportion of Churn by {label}',
             xlabel = label,
             ylabel = 'Churn Proportion (%)')
axs[3,1].set_xticklabels(axs[3,1].get_xticklabels(), rotation=0)

gs = gridspec.GridSpec(rows, cols, figure=fig)
ax_big = fig.add_subplot(gs[4, 0:])
sns.countplot(data=mom12, x='MonthsOfMembership', hue=x, ax=ax_big, palette=colour)
ax_big.set(title=f'Count of Months of Membership by {label}',
           xlabel="Months of Membership",
           ylabel='Count')

axs[3,2].set_axis_off()
axs[4,0].set_axis_off()
axs[4,1].set_axis_off()
axs[4,2].set_axis_off()
plt.tight_layout()

plt.show()


# Check distribution by education

x = 'Education'
label = 'Education Level'
order = ["Doctor", "Master", "Bachelor", "College", "High School or Below"]
rows, cols = 5, 3

# Calculate total flights booked by education level
total_flights = cust_eda_loyalty_df.groupby(x)['TotalTotalFlights'].sum().reset_index()
total_flights.Education = total_flights.Education.astype("category")
total_flights.Education = total_flights.Education.cat.set_categories(order)
total_flights.sort_values([x])

# Sort education level in cust_eda_loyalty_df in descending order
cust_eda_loyalty_df.Education = cust_eda_loyalty_df.Education.astype("category")
cust_eda_loyalty_df.Education = cust_eda_loyalty_df.Education.cat.set_categories(order)
cust_eda_loyalty_df.sort_values([x])

# Calculate proportion of churned customers by education level
churn_proportions = cust_eda_loyalty_df.groupby(x)['ChurnIndicator'].value_counts(normalize=True).unstack().fillna(0)
churn_proportions.columns = ['Not Churn', 'Churn']
churn_proportions['Not Churn'] *= 100
churn_proportions['Churn'] *= 100
churn_proportions.sort_values('Churn', ascending=False, inplace=True)

# Focus on months of membership between 0 to 12 months
mom12 = cust_eda_loyalty_df[(cust_eda_loyalty_df['MonthsOfMembership'] >= 0) & (cust_eda_loyalty_df['MonthsOfMembership'] <= 12)]
mom12[x] = pd.Categorical(mom12[x], categories=order, ordered=True)
mom12 = mom12.sort_values(x)

# Define plot area
fig, axs = plt.subplots(nrows = rows, ncols = cols, figsize=(20,25))

# Create plots
plt.subplot(rows, cols, 1)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue=x)
plt.title(f'Frequency of Customers by {label}')
plt.xlabel(label)
plt.ylabel('Count')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 2)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='Salary', hue=x)
plt.title(f'Distribution of Salary by {label}')
plt.xlabel(label)
plt.ylabel('Salary')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 3)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='CLV', hue=x)
plt.title(f'Distribution of CLV by {label}')
plt.xlabel(label)
plt.ylabel('CLV')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 4)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='MonthsOfMembership', hue=x)
plt.title(f'Distribution of Membership Length (Months) by {label}')
plt.xlabel(label)
plt.ylabel('Membership Length (Months)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 5)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsBooked', hue=x)
plt.title(f'Distribution of Individual Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Individual Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 6)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsWithCompanions', hue=x)
plt.title(f'Distribution of Flights with Companions by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights with Companions Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 7)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalTotalFlights', hue=x)
plt.title(f'Distribution of Total Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 8)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalDistance', hue=x)
plt.title(f'Distribution of Total Distance Travelled by {label}')
plt.xlabel(label)
plt.ylabel('Total Distance')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 9)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='PercentagePointsRedeemed', hue=x, order=order)
plt.title(f'Distribution of Points Redeemed by {label}')
plt.xlabel(label)
plt.ylabel('Points Redeemed (%)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 10)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue='ChurnIndicator', palette="viridis")
plt.title(f'Frequency of Churn by {label}')
plt.xlabel(label)
plt.ylabel('Churn Frequency')
plt.legend(title='Churn Indicator', loc='upper right')

churn_proportions.plot(kind='bar', stacked=True, ax = axs[3,1], color=['steelblue', 'mediumseagreen'])
axs[3,1].set(title = f'Proportion of Churn by {label}',
             xlabel = label,
             ylabel = 'Churn Proportion (%)')
axs[3,1].set_xticklabels(axs[3,1].get_xticklabels(), rotation=0)

gs = gridspec.GridSpec(rows, cols, figure=fig)
ax_big = fig.add_subplot(gs[4, 0:])
sns.countplot(data=mom12, x='MonthsOfMembership', hue=x, ax=ax_big)
ax_big.set(title=f'Count of Months of Membership by {label}',
           xlabel="Months of Membership",
           ylabel='Count')

axs[3,2].set_axis_off()
axs[4,0].set_axis_off()
axs[4,1].set_axis_off()
axs[4,2].set_axis_off()
plt.tight_layout()

plt.show()


# Check distribution by marital status

x = 'MaritalStatus'
label = "Marital Status"
order = ["Married","Divorced","Single"]
rows, cols = 5, 3

# Calculate total flights booked by marital status
total_flights = cust_eda_loyalty_df.groupby(x)['TotalTotalFlights'].sum().reset_index()
total_flights.MaritalStatus = total_flights.MaritalStatus.astype("category")
total_flights.MaritalStatus = total_flights.MaritalStatus.cat.set_categories(order)
total_flights.sort_values([x])

# Calculate proportion of churned customers by marital status
churn_proportions = cust_eda_loyalty_df.groupby(x)['ChurnIndicator'].value_counts(normalize=True).unstack().fillna(0)
churn_proportions.columns = ['Not Churn', 'Churn']
churn_proportions['Not Churn'] *= 100
churn_proportions['Churn'] *= 100
churn_proportions.sort_values('Churn', ascending=False, inplace=True)

# Focus on months of membership between 0 to 12 months
mom12 = cust_eda_loyalty_df[(cust_eda_loyalty_df['MonthsOfMembership'] >= 0) & (cust_eda_loyalty_df['MonthsOfMembership'] <= 12)]
mom12[x] = pd.Categorical(mom12[x], categories=order, ordered=True)
mom12 = mom12.sort_values(x)

# Deine plot area
fig, axs = plt.subplots(nrows = rows, ncols = cols, figsize=(20,25))

# Create plots
plt.subplot(rows, cols, 1)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue=x)
plt.title(f'Frequency of Customers by {label}')
plt.xlabel(label)
plt.ylabel('Count')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 2)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='Salary', hue=x)
plt.title(f'Distribution of Salary by {label}')
plt.xlabel(label)
plt.ylabel('Salary')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 3)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='CLV', hue=x)
plt.title(f'Distribution of CLV by {label}')
plt.xlabel(label)
plt.ylabel('CLV')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 4)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='MonthsOfMembership', hue=x)
plt.title(f'Distribution of Membership Length (Months) by {label}')
plt.xlabel(label)
plt.ylabel('Membership Length (Months)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 5)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsBooked', hue=x)
plt.title(f'Distribution of Individual Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Individual Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 6)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsWithCompanions', hue=x)
plt.title(f'Distribution of Flights with Companions by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights with Companions Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 7)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalTotalFlights', hue=x)
plt.title(f'Distribution of Total Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 8)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalDistance', hue=x)
plt.title(f'Distribution of Total Distance Travelled by {label}')
plt.xlabel(label)
plt.ylabel('Total Distance')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 9)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='PercentagePointsRedeemed', hue=x, order=order)
plt.title(f'Distribution of Points Redeemed by {label}')
plt.xlabel(label)
plt.ylabel('Points Redeemed (%)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 10)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue='ChurnIndicator', palette="viridis")
plt.title(f'Frequency of Churn by {label}')
plt.xlabel(label)
plt.ylabel('Churn Frequency')
plt.legend(title='Churn Indicator', loc='upper right')

churn_proportions.plot(kind='bar', stacked=True, ax = axs[3,1], color=['steelblue', 'mediumseagreen'])
axs[3,1].set(title = f'Proportion of Churn by {label}',
             xlabel = label,
             ylabel = 'Churn Proportion (%)')
axs[3,1].set_xticklabels(axs[3,1].get_xticklabels(), rotation=0)

gs = gridspec.GridSpec(rows, cols, figure=fig)
ax_big = fig.add_subplot(gs[4, 0:])
sns.countplot(data=mom12, x='MonthsOfMembership', hue=x, ax=ax_big)
ax_big.set(title=f'Count of Months of Membership by {label}',
           xlabel="Months of Membership",
           ylabel='Count')

axs[3,2].set_axis_off()
axs[4,0].set_axis_off()
axs[4,1].set_axis_off()
axs[4,2].set_axis_off()
plt.tight_layout()

plt.show()


# Check distribution by loyalty card type

x = 'LoyaltyCard'
label = "Loyalty Card"
order = ["Star", "Nova", "Aurora"]
rows, cols = 5, 3

# Calculate total flights booked by loyalty card type
total_flights = cust_eda_loyalty_df.groupby(x)['TotalTotalFlights'].sum().reset_index()
total_flights.LoyaltyCard = total_flights.LoyaltyCard.astype("category")
total_flights.LoyaltyCard = total_flights.LoyaltyCard.cat.set_categories(order)
total_flights.sort_values([x])

# Calculate proportion of churned customers by loyalty card type
churn_proportions = cust_eda_loyalty_df.groupby(x)['ChurnIndicator'].value_counts(normalize=True).unstack().fillna(0)
churn_proportions.columns = ['Not Churn', 'Churn']
churn_proportions['Not Churn'] *= 100
churn_proportions['Churn'] *= 100
churn_proportions.sort_values('Churn', ascending=False, inplace=True)

# Focus on months of membership between 0 to 12 months
mom12 = cust_eda_loyalty_df[(cust_eda_loyalty_df['MonthsOfMembership'] >= 0) & (cust_eda_loyalty_df['MonthsOfMembership'] <= 12)]
mom12[x] = pd.Categorical(mom12[x], categories=order, ordered=True)
mom12 = mom12.sort_values(x)

# Define plot area
fig, axs = plt.subplots(nrows = rows, ncols = cols, figsize=(20,25))

# Create plots
plt.subplot(rows, cols, 1)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue=x)
plt.title(f'Frequency of Customers by {label}')
plt.xlabel(label)
plt.ylabel('Count')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 2)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='Salary', hue=x)
plt.title(f'Distribution of Salary by {label}')
plt.xlabel(label)
plt.ylabel('Salary')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 3)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='CLV', hue=x)
plt.title(f'Distribution of CLV by {label}')
plt.xlabel(label)
plt.ylabel('CLV')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 4)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='MonthsOfMembership', hue=x)
plt.title(f'Distribution of Membership Length (Months) by {label}')
plt.xlabel(label)
plt.ylabel('Membership Length (Months)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 5)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsBooked', hue=x)
plt.title(f'Distribution of Individual Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Individual Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 6)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsWithCompanions', hue=x)
plt.title(f'Distribution of Flights with Companions by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights with Companions Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 7)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalTotalFlights', hue=x)
plt.title(f'Distribution of Total Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 8)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalDistance', hue=x)
plt.title(f'Distribution of Total Distance Travelled by {label}')
plt.xlabel(label)
plt.ylabel('Total Distance')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 9)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='PercentagePointsRedeemed', hue=x, order=order)
plt.title(f'Distribution of Points Redeemed by {label}')
plt.xlabel(label)
plt.ylabel('Points Redeemed (%)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 10)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue='ChurnIndicator', palette="viridis")
plt.title(f'Frequency of Churn by {label}')
plt.xlabel(label)
plt.ylabel('Churn Frequency')
plt.legend(title='Churn Indicator', loc='upper right')

churn_proportions.plot(kind='bar', stacked=True, ax = axs[3,1], color=['steelblue', 'mediumseagreen'])
axs[3,1].set(title = f'Proportion of Churn by {label}',
             xlabel = label,
             ylabel = 'Churn Proportion (%)')
axs[3,1].set_xticklabels(axs[3,1].get_xticklabels(), rotation=0)

gs = gridspec.GridSpec(rows, cols, figure=fig)
ax_big = fig.add_subplot(gs[4, 0:])
sns.countplot(data=mom12, x='MonthsOfMembership', hue=x, ax=ax_big)
ax_big.set(title=f'Count of Months of Membership by {label}',
           xlabel="Months of Membership",
           ylabel='Count')

axs[3,2].set_axis_off()
axs[4,0].set_axis_off()
axs[4,1].set_axis_off()
axs[4,2].set_axis_off()
plt.tight_layout()

plt.show()


# Check distribution by city

x = 'City'
label = "City"
rows, cols = 5, 3

# Create new df to sort in descending order of median of CLV
median_clv_df = loyalty_df.copy()
median_clv_df["City"] = pd.Categorical(median_clv_df["City"], categories=median_clv_df.groupby("City")['CLV'].median().sort_values(ascending=False).index, ordered=True)

# Calculate frequency of each city and sort in descending order
city_counts = cust_eda_loyalty_df[x].value_counts().reset_index()
city_counts.columns = [x, 'Count']
city_counts = city_counts.sort_values('Count', ascending=False)
order = city_counts[x].to_list()

# Sort city in cust_eda_loyalty_df in descending order of frequency
cust_eda_loyalty_df['City'] = pd.Categorical(cust_eda_loyalty_df['City'], categories=order, ordered=True)
cust_eda_loyalty_df = cust_eda_loyalty_df.sort_values('City')

# Calculate total flights booked by each city
total_flights = cust_eda_loyalty_df.groupby(x)['TotalTotalFlights'].sum().reset_index()
total_flights.City = total_flights.City.astype("category")
total_flights.City = total_flights.City.cat.set_categories(order)
total_flights.sort_values([x])

# Calculate proportion of churned customers by city
churn_proportions = cust_eda_loyalty_df.groupby(x)['ChurnIndicator'].value_counts(normalize=True).unstack().fillna(0)
churn_proportions.columns = ['Not Churn', 'Churn']
churn_proportions['Not Churn'] *= 100
churn_proportions['Churn'] *= 100
churn_proportions.sort_values('Churn', ascending=False, inplace=True)

# Focus on months of membership between 0 to 12 months
mom12 = cust_eda_loyalty_df[(cust_eda_loyalty_df['MonthsOfMembership'] >= 0) & (cust_eda_loyalty_df['MonthsOfMembership'] <= 12)]
mom12[x] = pd.Categorical(mom12[x], categories=order, ordered=True)
mom12 = mom12.sort_values(x)

# Create custom colour mapping to each city
palette = cc.cm.glasbey_dark
colour = [cc.cm.glasbey_dark(i / len(cust_eda_loyalty_df["City"])) for i in range(len(cust_eda_loyalty_df["City"]))]
col_map = {city: col for city, col in zip(cust_eda_loyalty_df['City'], colour)}

# Define plot area
fig, axs = plt.subplots(nrows = rows, ncols = cols, figsize=(30,40))

# Create plots
plt.subplot(rows, cols, 1)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue=x, palette=col_map)
plt.title(f'Frequency of Customers by {label}')
plt.xlabel(label)
plt.ylabel('Count')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 2)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='Salary', hue=x, palette=col_map)
plt.title(f'Distribution of Salary by {label}')
plt.xlabel(label)
plt.ylabel('Salary')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 3)
sns.boxplot(data=median_clv_df, x=x, y='CLV', hue=x, palette=col_map)
plt.title(f'Distribution of CLV by {label}')
plt.xlabel(label)
plt.ylabel('CLV')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 4)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='MonthsOfMembership', hue=x, palette=col_map)
plt.title(f'Distribution of Membership Length (Months) by {label}')
plt.xlabel(label)
plt.ylabel('Membership Length (Months)')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 5)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsBooked', hue=x, palette=col_map)
plt.title(f'Distribution of Individual Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Individual Flights Booked')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 6)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsWithCompanions', hue=x, palette=col_map)
plt.title(f'Distribution of Flights with Companions by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights with Companions Flights Booked')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 7)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalTotalFlights', hue=x, palette=col_map)
plt.title(f'Distribution of Total Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights Booked')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 8)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalDistance', hue=x, palette=col_map)
plt.title(f'Distribution of Total Distance Travelled by {label}')
plt.xlabel(label)
plt.ylabel('Total Distance')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 9)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='PercentagePointsRedeemed', hue=x, order=order, palette=col_map)
plt.title(f'Distribution of Points Redeemed by {label}')
plt.xlabel(label)
plt.ylabel('Points Redeemed (%)')
plt.legend([], [], frameon=False)
plt.xticks(rotation=80)

plt.subplot(rows, cols, 10)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue='ChurnIndicator', palette="viridis")
plt.title(f'Frequency of Churn by {label}')
plt.xlabel(label)
plt.ylabel('Churn Frequency')
plt.legend(title='Churn Indicator', loc='upper right')
plt.xticks(rotation=80)

churn_proportions.plot(kind='bar', stacked=True, ax = axs[3,1], color=['steelblue', 'mediumseagreen'])
axs[3,1].set(title = f'Proportion of Churn by {label}',
             xlabel = label,
             ylabel = 'Churn Proportion (%)')
axs[3,1].set_xticklabels(axs[3,1].get_xticklabels(), rotation=80)

gs = gridspec.GridSpec(rows, cols, figure=fig)
ax_big = fig.add_subplot(gs[4, 0:])
sns.countplot(data=mom12, x='MonthsOfMembership', hue=x, ax=ax_big, palette=col_map)
ax_big.set(title=f'Count of Months of Membership by {label}',
           xlabel="Months of Membership",
           ylabel='Count')

axs[3,2].set_axis_off()
axs[4,0].set_axis_off()
axs[4,1].set_axis_off()
axs[4,2].set_axis_off()
plt.tight_layout()

plt.show()


# Distribution by Enrollment Type

x = 'EnrollmentType'
label = "Enrollment Type"
order = ["Standard","2018 Promotion"]
rows, cols = 5, 3

# Calculate total flights booked by enrollment type
total_flights = cust_eda_loyalty_df.groupby(x)['TotalTotalFlights'].sum().reset_index()
total_flights.EnrollmentType = total_flights.EnrollmentType.astype("category")
total_flights.EnrollmentType = total_flights.EnrollmentType.cat.set_categories(order)
total_flights.sort_values([x])

# Calculate proportion of churned customers by enrollment type
churn_proportions = cust_eda_loyalty_df.groupby(x)['ChurnIndicator'].value_counts(normalize=True).unstack().fillna(0)
churn_proportions.columns = ['Not Churn', 'Churn']
churn_proportions['Not Churn'] *= 100
churn_proportions['Churn'] *= 100
churn_proportions.sort_values('Churn', ascending=False, inplace=True)

# Focus on months of membership between 0 to 12 months
mom12 = cust_eda_loyalty_df[(cust_eda_loyalty_df['MonthsOfMembership'] >= 0) & (cust_eda_loyalty_df['MonthsOfMembership'] <= 12)]
mom12[x] = pd.Categorical(mom12[x], categories=order, ordered=True)
mom12 = mom12.sort_values(x)

# Define plot area
fig, axs = plt.subplots(nrows = rows, ncols = cols, figsize=(20,25))

# Create plots
plt.subplot(rows, cols, 1)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue=x)
plt.title(f'Frequency of Customers by {label}')
plt.xlabel(label)
plt.ylabel('Count')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 2)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='Salary', hue=x)
plt.title(f'Distribution of Salary by {label}')
plt.xlabel(label)
plt.ylabel('Salary')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 3)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='CLV', hue=x)
plt.title(f'Distribution of CLV by {label}')
plt.xlabel(label)
plt.ylabel('CLV')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 4)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='MonthsOfMembership', hue=x)
plt.title(f'Distribution of Membership Length (Months) by {label}')
plt.xlabel(label)
plt.ylabel('Membership Length (Months)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 5)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsBooked', hue=x)
plt.title(f'Distribution of Individual Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Individual Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 6)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalFlightsWithCompanions', hue=x)
plt.title(f'Distribution of Flights with Companions by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights with Companions Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 7)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalTotalFlights', hue=x)
plt.title(f'Distribution of Total Flights by {label}')
plt.xlabel(label)
plt.ylabel('Total Flights Booked')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 8)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='TotalDistance', hue=x)
plt.title(f'Distribution of Total Distance Travelled by {label}')
plt.xlabel(label)
plt.ylabel('Total Distance')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 9)
sns.boxplot(data=cust_eda_loyalty_df, x=x, y='PercentagePointsRedeemed', hue=x, order=order)
plt.title(f'Distribution of Points Redeemed by {label}')
plt.xlabel(label)
plt.ylabel('Points Redeemed (%)')
plt.legend([], [], frameon=False)

plt.subplot(rows, cols, 10)
sns.countplot(data=cust_eda_loyalty_df, x=x, hue='ChurnIndicator', palette="viridis")
plt.title(f'Frequency of Churn by {label}')
plt.xlabel(label)
plt.ylabel('Churn Frequency')
plt.legend(title='Churn Indicator', loc='upper right')

churn_proportions.plot(kind='bar', stacked=True, ax = axs[3,1], color=['steelblue', 'mediumseagreen'])
axs[3,1].set(title = f'Proportion of Churn by {label}',
             xlabel = label,
             ylabel = 'Churn Proportion (%)')
axs[3,1].set_xticklabels(axs[3,1].get_xticklabels(), rotation=0)

gs = gridspec.GridSpec(rows, cols, figure=fig)
ax_big = fig.add_subplot(gs[4, 0:])
sns.countplot(data=mom12, x='MonthsOfMembership', hue='EnrollmentType', ax=ax_big)
ax_big.set(title=f'Count of Months of Membership by {label}',
           xlabel="Months of Membership",
           ylabel='Count')

axs[3,2].set_axis_off()
axs[4,0].set_axis_off()
axs[4,1].set_axis_off()
axs[4,2].set_axis_off()
plt.tight_layout()

plt.show()


# Check distribution by CLV

plt.figure(figsize=(20,25))

plt.subplot(3, 3, 1)
sns.lineplot(data=cust_eda_loyalty_df, x='EnrollmentYear', y='CLV', ci=None, marker='o', sort=True)
plt.title('Customer Lifetime Value by Enrollment Year')
plt.xlabel('Enrollment Year')
plt.ylabel('Customer Lifetime Value')

plt.subplot(3, 3, 2)
sns.lineplot(data=cust_eda_loyalty_df, x='City', y='CLV', ci=None, marker='o', sort=True)
plt.title('Customer Lifetime Value by City')
plt.xlabel('City')
plt.ylabel('Customer Lifetime Value')
plt.xticks(rotation=80)

plt.subplot(3, 3, 4)
corr_coeff, _ = pearsonr(cust_eda_loyalty_df['Salary'], cust_eda_loyalty_df['CLV'])
sns.regplot(data=cust_eda_loyalty_df, x='Salary', y='CLV', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Customer Lifetime Value by Salary')
plt.xlabel('Salary')
plt.ylabel('Customer Lifetime Value')
plt.xticks(rotation=80)

plt.subplot(3, 3, 5)
corr_coeff, _ = pearsonr(cust_eda_loyalty_df['TotalTotalFlights'], cust_eda_loyalty_df['CLV'])
sns.regplot(data=cust_eda_loyalty_df, x='TotalTotalFlights', y='CLV', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Customer Lifetime Value by Total Flights Booked')
plt.xlabel('Total Flights Booked')
plt.ylabel('Customer Lifetime Value')
plt.xticks(rotation=80)

plt.subplot(3, 3, 6)
corr_coeff, _ = pearsonr(cust_eda_loyalty_df['TotalDistance'], cust_eda_loyalty_df['CLV'])
sns.regplot(data=cust_eda_loyalty_df, x='TotalDistance', y='CLV', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Customer Lifetime Value by Total Distance Travelled')
plt.xlabel('Total Distance Travelled')
plt.ylabel('Customer Lifetime Value')
plt.xticks(rotation=80)

plt.subplot(3, 3, 7)
corr_coeff, _ = pearsonr(cust_eda_loyalty_df['MonthsOfMembership'], cust_eda_loyalty_df['CLV'])
sns.regplot(data=cust_eda_loyalty_df, x='MonthsOfMembership', y='CLV', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Customer Lifetime Value by Months of Membership')
plt.xlabel('Months of Membership')
plt.ylabel('Customer Lifetime Value')
plt.xticks(rotation=80)

plt.subplot(3, 3, 8)
corr_coeff, _ = pearsonr(cust_eda_loyalty_df['TotalPointsAccumulated'], cust_eda_loyalty_df['CLV'])
sns.regplot(data=cust_eda_loyalty_df, x='TotalPointsAccumulated', y='CLV', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Customer Lifetime Value by Total Points Accumulated')
plt.xlabel('Total Points Accumulated')
plt.ylabel('Customer Lifetime Value')
plt.xticks(rotation=80)

plt.show()


# Assuming 'continuous_var' is your continuous variable and 'nominal_var' is your nominal variable
def kruskal_wallis(df, cont_var, nom_var):
    groups = [df[df[f'{nom_var}'] == category][f'{cont_var}'] for category in df[f'{nom_var}'].unique()]
    kruskal_result = stats.kruskal(*groups)
    print(f"Kruskal-Wallis H-statistic for {nom_var}: {kruskal_result.statistic}, p-value: {kruskal_result.pvalue}")


for i in ['MaritalStatus', 'Education', 'LoyaltyCard']:
    kruskal_wallis(loyalty_df, 'CLV', i)

Kruskal-Wallis H-statistic for MaritalStatus: 32.67197728755935, p-value: 8.042114934653465e-08
Kruskal-Wallis H-statistic for Education: 82.73214822370325, p-value: 4.5915624634674436e-17
Kruskal-Wallis H-statistic for LoyaltyCard: 2204.896056033002, p-value: 0.0


### 3.2.4 Relationship between Points Accumulated and Points Redeemed
corr_coeff, _ = pearsonr(cust_eda_loyalty_df['TotalPointsAccumulated'], cust_eda_loyalty_df['TotalPointsRedeemed'])
sns.regplot(data=cust_eda_loyalty_df, x='TotalPointsAccumulated', y='TotalPointsRedeemed', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Points Accumulated by Points Redeemed')
plt.xlabel('Total Points Accumulated')
plt.ylabel('Total Points Redeemed')
plt.xticks(rotation=80)

(array([-2000.,     0.,  2000.,  4000.,  6000.,  8000., 10000., 12000.]),
 [Text(-2000.0, 0, '−2000'),
  Text(0.0, 0, '0'),
  Text(2000.0, 0, '2000'),
  Text(4000.0, 0, '4000'),
  Text(6000.0, 0, '6000'),
  Text(8000.0, 0, '8000'),
  Text(10000.0, 0, '10000'),
  Text(12000.0, 0, '12000')])


corr_coeff, _ = pearsonr(cust_eda_loyalty_df['TotalPointsAccumulated'], cust_eda_loyalty_df['TotalPointsRedeemed'])
sns.regplot(data=cust_eda_loyalty_df, x='TotalPointsAccumulated', y='TotalPointsRedeemed', 
            scatter_kws={'s': 50, 'alpha': 0.9, 'edgecolor': 'white'}, 
            line_kws={'color': 'red'}, marker="o")
plt.text(0.1, 0.9, f'Correlation Coefficient: {corr_coeff:.2f}', 
         ha='left', va='baseline', transform=plt.gca().transAxes)
plt.title('Points Accumulated by Points Redeemed')
plt.xlabel('Total Points Accumulated')
plt.ylabel('Total Points Redeemed')
plt.xticks(rotation=80)

(array([-2000.,     0.,  2000.,  4000.,  6000.,  8000., 10000., 12000.]),
 [Text(-2000.0, 0, '−2000'),
  Text(0.0, 0, '0'),
  Text(2000.0, 0, '2000'),
  Text(4000.0, 0, '4000'),
  Text(6000.0, 0, '6000'),
  Text(8000.0, 0, '8000'),
  Text(10000.0, 0, '10000'),
  Text(12000.0, 0, '12000')])


year_enrollment = cust_eda_loyalty_df.groupby('EnrollmentYear')['EnrollmentYear'].value_counts().reset_index()
year_enrollment['growth'] = year_enrollment['count'].pct_change()*100
year_enrollment = year_enrollment.fillna(0)
year_enrollment['EnrollmentYear'] = year_enrollment['EnrollmentYear'].astype('str')
year_enrollment


fig, axs = plt.subplots(nrows = 1, ncols = 3, figsize = (15,5))


sns.lineplot(year_enrollment[~year_enrollment['EnrollmentYear'].isin(['2012','2013'])], x = 'EnrollmentYear', y='growth',label='growth', ax=axs[0])
axs[0].set_title('Percentage change in # of customers')
for x, y in zip(year_enrollment[~year_enrollment['EnrollmentYear'].isin(['2012','2013'])]['EnrollmentYear'], year_enrollment[~year_enrollment['EnrollmentYear'].isin(['2012','2013'])]['growth']):
    axs[0].text(x, y, f'{y:.2f}%', ha='center', va='bottom')

sns.lineplot(year_enrollment[~year_enrollment['EnrollmentYear'].isin(['2012','2013'])], x = 'EnrollmentYear', y='count',label='count', ax=axs[1])
axs[1].set_title('Change in # of customers')
for x, y in zip(year_enrollment[~year_enrollment['EnrollmentYear'].isin(['2012','2013'])]['EnrollmentYear'], year_enrollment[~year_enrollment['EnrollmentYear'].isin(['2012','2013'])]['count']):
    axs[1].text(x, y, f'{y}', ha='center', va='bottom')

pivot_plot = loyalty_df.groupby(['EnrollmentYear','EnrollmentType']).size().reset_index().pivot(columns='EnrollmentType',index = 'EnrollmentYear', values = 0)
# Reorder the columns
desired_order = ['Standard', '2018 Promotion']
pivot_plot = pivot_plot[desired_order]
# sns.countplot(loyalty_df[~loyalty_df['EnrollmentYear'].isin(['2012','2013'])], x = 'EnrollmentYear', hue = 'EnrollmentType', multiple='stack', hue_order=['2018 Promotion', 'Standard'])
pivot_plot.plot(kind='bar', stacked=True, ax = axs[2])
axs[2].set_title('Change in # of customers')

plt.tight_layout()


# Cancellation by year and month

cancellations_heatmap = cust_eda_loyalty_df.pivot_table(values='LoyaltyNumber', index='CancellationYear', columns='CancellationMonth', aggfunc='count', fill_value=0)

plt.figure(figsize=(12, 8))
sns.heatmap(cancellations_heatmap, annot=True, fmt='d', cmap='viridis')
plt.title('Heatmap of Cancellations by Year and Month')
plt.xlabel('Cancellation Month')
plt.ylabel('Cancellation Year')
plt.show()


enrollments = cust_eda_loyalty_df['EnrollmentDate'].dt.to_period('M').value_counts().sort_index()
cancellations = cust_eda_loyalty_df['CancellationDate'].dropna().dt.to_period('M').value_counts().sort_index()

plt.figure(figsize=(14, 8))
plt.plot(enrollments.index.astype(str), enrollments.values, label='Enrollments', marker='o')
plt.plot(cancellations.index.astype(str), cancellations.values, label='Cancellations', marker='o')
plt.title('Number of Enrollments and Cancellations Over Time')
plt.xlabel('Date (Year-Month)')
plt.ylabel('Count')
plt.xticks(rotation=90)
plt.legend()
plt.show()


enrollments_2018 = cust_eda_loyalty_df.groupby([loyalty_df['EnrollmentDate'].dt.to_period('M'), 'EnrollmentType']).size().unstack(fill_value=0)
enrollment_order = ['Standard', '2018 Promotion']
enrollments_2018 = enrollments_2018[enrollment_order]

plt.figure(figsize=(14, 8))
for column in enrollments_2018.columns:
    plt.plot(enrollments_2018.index.astype(str), enrollments_2018[column], label=column, marker='o')
plt.xlabel('Date')
plt.ylabel('Enrollments')
plt.title('Enrollments by Enrollment Type')
plt.legend(title='Enrollment Type')
plt.xticks(rotation=90)
plt.show()


loyalty_df2018 = cust_eda_loyalty_df[cust_eda_loyalty_df['EnrollmentYear'] == 2018]


loyalty_df2018.info()

<class 'pandas.core.frame.DataFrame'>
Index: 3010 entries, 10966 to 13406
Data columns (total 32 columns):
 #   Column                         Non-Null Count  Dtype         
---  ------                         --------------  -----         
 0   LoyaltyNumber                  3010 non-null   object        
 1   Country                        3010 non-null   object        
 2   Province                       3010 non-null   object        
 3   City                           3010 non-null   category      
 4   PostalCode                     3010 non-null   object        
 5   Gender                         3010 non-null   object        
 6   Education                      3010 non-null   category      
 7   Salary                         3010 non-null   float64       
 8   MaritalStatus                  3010 non-null   object        
 9   LoyaltyCard                    3010 non-null   object        
 10  CLV                            3010 non-null   float64       
 11  EnrollmentType                 3010 non-null   object        
 12  EnrollmentYear                 3010 non-null   int64         
 13  EnrollmentMonth                3010 non-null   int64         
 14  CancellationYear               159 non-null    Int64         
 15  CancellationMonth              159 non-null    Int64         
 16  EnrollmentDate                 3010 non-null   datetime64[ns]
 17  EnrollmentYearMonth            3010 non-null   object        
 18  CancellationDate               159 non-null    datetime64[ns]
 19  CancellationYearMonth          159 non-null    object        
 20  ChurnIndicator                 3010 non-null   category      
 21  MonthsOfMembership             3010 non-null   int64         
 22  TotalFlightsBooked             3010 non-null   int64         
 23  TotalFlightsWithCompanions     3010 non-null   int64         
 24  TotalTotalFlights              3010 non-null   int64         
 25  TotalDistance                  3010 non-null   int64         
 26  TotalPointsAccumulated         3010 non-null   float64       
 27  TotalPointsRedeemed            3010 non-null   int64         
 28  TotalDollarCostPointsRedeemed  3010 non-null   int64         
 29  PercentagePointsRedeemed       3010 non-null   float64       
 30  MonthlyTotalFlights            3010 non-null   float64       
 31  MonthlyPointsAccumulated       3010 non-null   float64       
dtypes: Int64(2), category(3), datetime64[ns](2), float64(6), int64(9), object(10)
memory usage: 721.7+ KB


# Calculate count and percentage for each enrollment type
enrollment_counts = loyalty_df2018['EnrollmentType'].value_counts()
enrollment_percentage = (enrollment_counts / len(loyalty_df2018)) * 100

# Combine count and percentage into a DataFrame
enrollment_summary = pd.DataFrame({'Count': enrollment_counts, 'Percentage': enrollment_percentage})

print(enrollment_summary)

                Count  Percentage
EnrollmentType                   
Standard         2039   67.740864
2018 Promotion    971   32.259136


# Plot churn indicator by enrollment type
plt.figure(figsize=(10, 6))
ax = sns.countplot(data=loyalty_df2018, x='EnrollmentType', hue='ChurnIndicator')

# Calculate total count of each enrollment type
total_counts = loyalty_df2018['EnrollmentType'].value_counts()

# Add percentage labels
for container in ax.containers:
    total_height = sum([patch.get_height() for patch in container.patches])
    for patch, en_type in zip(container.patches, total_counts.index):
        height = patch.get_height()
        percentage = height / total_counts[en_type]
        ax.annotate(f'{percentage:.1%}', 
                    xy=(patch.get_x() + patch.get_width() / 2, height), 
                    xytext=(0, 3),  # 3 points vertical offset
                    textcoords="offset points",
                    ha='center', va='bottom', fontsize=8)

plt.title('Churn Indicator by Enrollment Type')
plt.xlabel('Enrollment Type')
plt.ylabel('Count')
plt.legend(title='Churn Indicator', loc='upper right')
plt.show()


from scipy.stats import chi2_contingency

# Create a contingency table of churn indicator by enrollment type
contingency_table = pd.crosstab(loyalty_df2018['EnrollmentType'], loyalty_df2018['ChurnIndicator'])

# Perform the chi-square test
chi2, p, dof, expected = chi2_contingency(contingency_table)

print("Chi-square statistic:", chi2)
print("P-value:", p)
print("Degrees of freedom:", dof)
print("Expected frequencies table:")
print(expected)

Chi-square statistic: 121.39824875372153
P-value: 3.126313190708848e-28
Degrees of freedom: 1
Expected frequencies table:
[[ 919.70797342   51.29202658]
 [1931.29202658  107.70797342]]


plt.figure(figsize=(10, 10))
sns.boxplot(data=loyalty_df2018, x='EnrollmentType', y='CLV')
plt.title('CLV by Enrollment Type in 2018')
plt.xlabel('Enrollment Type')
plt.ylabel('CLV')
plt.show()


from scipy.stats import ttest_ind

# Extract CLV values for each EnrollmentType group
group1 = loyalty_df2018[loyalty_df2018['EnrollmentType'] == 'Standard']['CLV']
group2 = loyalty_df2018[loyalty_df2018['EnrollmentType'] == '2018 Promotion']['CLV']

# Perform independent t-test
t_statistic, p_value = ttest_ind(group1, group2, equal_var=False)

print(f'T-statistic: {t_statistic}')
print(f'p-value: {p_value}')

T-statistic: -0.13740317545689854
p-value: 0.8907276449119362


# get all churn customers

churn_monthsofmembership = loyalty_df[loyalty_df['ChurnIndicator'].eq(1)]['MonthsOfMembership'].value_counts().reset_index()\
    .sort_values(by='MonthsOfMembership')
    
# distribution of months of membership

plt.figure(figsize=(7,5))
sns.lineplot(data=churn_monthsofmembership, x='MonthsOfMembership', y='count')
plt.xlabel('Months of membership')
plt.ylabel('Number of customers')
plt.yticks(np.arange(0, 1400, 200))
plt.title('Distribution of all customers by months of membership')

Text(0.5, 1.0, 'Distribution of all customers by months of membership')


# get all churn customers with Standard Enrollment Type

churn_monthsofmembership_standard = loyalty_df[loyalty_df['ChurnIndicator'].eq(1) & loyalty_df['EnrollmentType'].ne('2018 Promotion')]\
    ['MonthsOfMembership'].value_counts().reset_index()\
    .sort_values(by='MonthsOfMembership')

# distribution of months of membership

plt.figure(figsize=(7,5))
sns.lineplot(data=churn_monthsofmembership_standard, x='MonthsOfMembership', y='count')
plt.xlabel('Months of membership')
plt.ylabel('Number of customers')
plt.yticks(np.arange(0, 1400, 200))
plt.title('Distribution of standard customers by months of membership')

Text(0.5, 1.0, 'Distribution of standard customers by months of membership')


continuous_vars = [
    'Salary', 'CLV', 'MonthsOfMembership', 'TotalFlightsBooked',
    'TotalFlightsWithCompanions', 'TotalTotalFlights', 'TotalDistance',
    'TotalPointsAccumulated', 'TotalPointsRedeemed', 'TotalDollarCostPointsRedeemed',
    'PercentagePointsRedeemed'
]

categorical_vars = ['Gender', 'Education', 'MaritalStatus', 'LoyaltyCard', 'EnrollmentType', 'Province', 'City']

# Convert `ChurnIndicator` to boolean for point-biserial correlation
loyalty_df['ChurnIndicator_bool'] = loyalty_df['ChurnIndicator'].astype(bool)


# Point-biserial correlation between ChurnIndicator and continuous variables
for var in continuous_vars:
    correlation, p_value = pointbiserialr(loyalty_df['ChurnIndicator'], loyalty_df[var])
    print(f'Point-biserial correlation between ChurnIndicator and {var}: {correlation:.2f}, p-value: {p_value}')

Point-biserial correlation between ChurnIndicator and Salary: -0.00, p-value: 0.8618576547393644
Point-biserial correlation between ChurnIndicator and CLV: 0.01, p-value: 0.3118873451918476
Point-biserial correlation between ChurnIndicator and MonthsOfMembership: -0.30, p-value: 0.0
Point-biserial correlation between ChurnIndicator and TotalFlightsBooked: -0.48, p-value: 0.0
Point-biserial correlation between ChurnIndicator and TotalFlightsWithCompanions: -0.42, p-value: 0.0
Point-biserial correlation between ChurnIndicator and TotalTotalFlights: -0.48, p-value: 0.0
Point-biserial correlation between ChurnIndicator and TotalDistance: -0.49, p-value: 0.0
Point-biserial correlation between ChurnIndicator and TotalPointsAccumulated: -0.49, p-value: 0.0
Point-biserial correlation between ChurnIndicator and TotalPointsRedeemed: -0.27, p-value: 1.0250023372563784e-286
Point-biserial correlation between ChurnIndicator and TotalDollarCostPointsRedeemed: -0.27, p-value: 5.016401724085526e-287
Point-biserial correlation between ChurnIndicator and PercentagePointsRedeemed: -0.11, p-value: 7.426373817580738e-48


# Calculate point-biserial correlations and p-values
correlation_results = []

for var in continuous_vars:
    correlation, p_value = pointbiserialr(loyalty_df['ChurnIndicator'], loyalty_df[var])
    correlation_results.append({'Variable': var, 'Correlation': correlation, 'P-value': p_value})

# Convert the results into a DataFrame
correlation_df = pd.DataFrame(correlation_results)
correlation_df.sort_values(by='Correlation', inplace=True, ascending=False)

# Plot the correlations
plt.figure(figsize=(12, 8))
sns.barplot(x='Correlation', y='Variable', data=correlation_df, palette='coolwarm', orient='h')
plt.axvline(x=0, color='gray', linestyle='--')
plt.title('Point-Biserial Correlations between ChurnIndicator and Continuous Variables')
plt.xlabel('Correlation Coefficient')
plt.ylabel('Continuous Variable')

# Add p-value annotations 
for i, row in correlation_df.iterrows():
    # Format p-values in scientific notation with three decimal places
    p_value_str = "{:.3e}".format(row["P-value"])
    color = 'black' 
    x_coord = -0.005
    plt.text(x_coord, i, f'p={p_value_str}', color=color, va='center', ha='right', fontsize=8)

plt.show()


def cramers_v(confusion_matrix):
    chi2 = chi2_contingency(confusion_matrix)[0]
    n = confusion_matrix.sum()
    r, k = confusion_matrix.shape
    return np.sqrt(chi2 / (n * (min(r, k) - 1)))

for var in categorical_vars:
    contingency_table = pd.crosstab(loyalty_df['ChurnIndicator'], loyalty_df[var])
    cramers_v_value = cramers_v(contingency_table.to_numpy())
    print(f"Cramér's V between ChurnIndicator and {var}: {cramers_v_value:.2f}")

Cramér's V between ChurnIndicator and Gender: 0.01
Cramér's V between ChurnIndicator and Education: 0.01
Cramér's V between ChurnIndicator and MaritalStatus: 0.01
Cramér's V between ChurnIndicator and LoyaltyCard: 0.02
Cramér's V between ChurnIndicator and EnrollmentType: 0.00
Cramér's V between ChurnIndicator and Province: 0.03
Cramér's V between ChurnIndicator and City: 0.03


def cramers_v(confusion_matrix):
    chi2 = chi2_contingency(confusion_matrix)[0]
    n = confusion_matrix.sum()
    r, k = confusion_matrix.shape
    return np.sqrt(chi2 / (n * (min(r, k) - 1)))

# Assuming categorical_vars contains the list of categorical variables
cramers_v_values = []

for var in categorical_vars:
    contingency_table = pd.crosstab(loyalty_df['ChurnIndicator'], loyalty_df[var])
    cramers_v_value = cramers_v(contingency_table.to_numpy())
    print(f"Cramér's V between ChurnIndicator and {var}: {cramers_v_value:.2f}")
    cramers_v_values.append(cramers_v_value)

# Plot all Cramér's V values in one plot
plt.figure(figsize=(10, 6))
bars = plt.barh(categorical_vars, cramers_v_values, color='skyblue')
plt.xlabel("Cramér's V Value")
plt.title("Cramér's V Values between ChurnIndicator and Categorical Variables")
plt.xlim(0, 1)  # Adjust the limit if needed
plt.grid(axis='x')

# Add value labels to the bars
for bar, value in zip(bars, cramers_v_values):
    plt.text(value, bar.get_y() + bar.get_height()/2, f'{value:.3f}', ha='left', va='center', fontsize=10)

plt.show()

Cramér's V between ChurnIndicator and Gender: 0.01
Cramér's V between ChurnIndicator and Education: 0.01
Cramér's V between ChurnIndicator and MaritalStatus: 0.01
Cramér's V between ChurnIndicator and LoyaltyCard: 0.02
Cramér's V between ChurnIndicator and EnrollmentType: 0.00
Cramér's V between ChurnIndicator and Province: 0.03
Cramér's V between ChurnIndicator and City: 0.03


# total number of customers in dataset

loyalty_df['LoyaltyNumber'].nunique()

16737


# calculate total flights in 2017 and 2018

total_flight_1718 = flight_df[['LoyaltyNumber', 'TotalFlights']].groupby(by='LoyaltyNumber', sort=False).sum()['TotalFlights']\
    .reset_index()
total_flight_1718


# get active customers in 2017 and 2018

active_cx_1718 = total_flight_1718[total_flight_1718['TotalFlights'].gt(0)]['LoyaltyNumber']
len(active_cx_1718)

15236


# get active customers in 2017 and 2018 enrolling through Standard Enrollment

active_standard_cx_1718 = loyalty_df[
    loyalty_df['LoyaltyNumber'].isin(active_cx_1718) 
    & loyalty_df['EnrollmentType'].ne('2018 Promotion')
]['LoyaltyNumber']
len(active_standard_cx_1718)

14303


active_loyalty_df = loyalty_df[loyalty_df['LoyaltyNumber'].isin(active_standard_cx_1718)]
active_loyalty_df.shape

(14303, 33)


# Filter active standard customers with: 
# - Year 2018 
# - Customers with more than 6 months of membership 
# - Remove records after customers cancel memberships

predict_retention_df = lookback_flight_df.merge(temp_mom_lookback, how='left', on=['LoyaltyNumber', 'Date'])\
    .pipe(lambda x: x[x['Date'].ge('2018-01-01') & x['Date'].ge(x['EnrollmentDate'])])\
    .pipe(lambda x: x[~x['Date'].gt(x['CancellationDate'])])\
    .assign(ChurnIndicator = lambda x: np.where(x['Date'].eq(x['CancellationDate']), 1, 0))\
    .pipe(lambda x: x[x['LoyaltyNumber'].isin(active_standard_cx_1718)])
    
predict_retention_df


predict_retention_df['LoyaltyNumber'].nunique()

13814


# Join demographic profiles

predict_retention_df = predict_retention_df.merge(loyalty_df.iloc[:, :9], how='left', on='LoyaltyNumber')\
    .drop(columns=['EnrollmentDate', 'CancellationDate', 'Country'])
    
predict_retention_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 250425 entries, 0 to 250424
Data columns (total 39 columns):
 #   Column                           Non-Null Count   Dtype         
---  ------                           --------------   -----         
 0   LoyaltyNumber                    250425 non-null  object        
 1   Date                             250425 non-null  datetime64[ns]
 2   Last1Y_FlightsBooked             250425 non-null  float64       
 3   Last1Y_FlightsWithCompanions     250425 non-null  float64       
 4   Last1Y_TotalFlights              250425 non-null  float64       
 5   Last1Y_Distance                  250425 non-null  float64       
 6   Last1Y_PointsAccumulated         250425 non-null  float64       
 7   Last1Y_PointsRedeemed            250425 non-null  float64       
 8   Last1Y_DollarCostPointsRedeemed  250425 non-null  float64       
 9   Last6M_FlightsBooked             250425 non-null  float64       
 10  Last6M_FlightsWithCompanions     250425 non-null  float64       
 11  Last6M_TotalFlights              250425 non-null  float64       
 12  Last6M_Distance                  250425 non-null  float64       
 13  Last6M_PointsAccumulated         250425 non-null  float64       
 14  Last6M_PointsRedeemed            250425 non-null  float64       
 15  Last6M_DollarCostPointsRedeemed  250425 non-null  float64       
 16  Last3M_FlightsBooked             250425 non-null  float64       
 17  Last3M_FlightsWithCompanions     250425 non-null  float64       
 18  Last3M_TotalFlights              250425 non-null  float64       
 19  Last3M_Distance                  250425 non-null  float64       
 20  Last3M_PointsAccumulated         250425 non-null  float64       
 21  Last3M_PointsRedeemed            250425 non-null  float64       
 22  Last3M_DollarCostPointsRedeemed  250425 non-null  float64       
 23  Last1M_FlightsBooked             250425 non-null  float64       
 24  Last1M_FlightsWithCompanions     250425 non-null  float64       
 25  Last1M_TotalFlights              250425 non-null  float64       
 26  Last1M_Distance                  250425 non-null  float64       
 27  Last1M_PointsAccumulated         250425 non-null  float64       
 28  Last1M_PointsRedeemed            250425 non-null  float64       
 29  Last1M_DollarCostPointsRedeemed  250425 non-null  float64       
 30  MonthsOfMembership               250425 non-null  int64         
 31  ChurnIndicator                   250425 non-null  int32         
 32  Province                         250425 non-null  object        
 33  City                             250425 non-null  object        
 34  PostalCode                       250425 non-null  object        
 35  Gender                           250425 non-null  object        
 36  Education                        250425 non-null  object        
 37  Salary                           250425 non-null  float64       
 38  MaritalStatus                    250425 non-null  object        
dtypes: datetime64[ns](1), float64(29), int32(1), int64(1), object(7)
memory usage: 73.6+ MB


# Label encoding for categorical features

label_encoder = LabelEncoder()

for col in predict_retention_df.columns[1:]:
    if predict_retention_df[col].dtype == 'object':
        predict_retention_df[col] = label_encoder.fit_transform(predict_retention_df[col])


predict_retention_df = predict_retention_df.reset_index()
predict_retention_df


predict_retention_df = predict_retention_df.set_index(['LoyaltyNumber', 'Date'])
predict_retention_df['ChurnIndicator'].value_counts(normalize=True) # churn rate

ChurnIndicator
0    0.996642
1    0.003358
Name: proportion, dtype: float64


predict_retention_df.index.get_level_values('LoyaltyNumber').nunique() # number of customers

13814


# get X and y

X_retention_data = predict_retention_df.drop(columns='ChurnIndicator')
y_retention_data = predict_retention_df[['ChurnIndicator']]


# get training and test sets

X_retention_train, X_retention_test, y_retention_train, y_retention_test = train_test_split(
    X_retention_data,
    y_retention_data,
    test_size=0.3,
    random_state=RANDOM_STATE,
    stratify=y_retention_data
)

print(
    'Train set:',
    '\n x_train:', X_retention_train.shape, 
    '\n y_train:', y_retention_train.shape, 
    '\n train size:', y_retention_train.shape[0]/y_retention_data.shape[0], 
    '\n churn rate:', y_retention_train.sum()/y_retention_train.count(),
    '\nTest set:',
    '\n x_test:', X_retention_test.shape, 
    '\n y_test:', y_retention_test.shape, 
    '\n test size:', y_retention_test.shape[0]/y_retention_data.shape[0], 
    '\n churn rate:', y_retention_test.sum()/y_retention_test.count()
)

Train set: 
 x_train: (175297, 37) 
 y_train: (175297, 1) 
 train size: 0.6999980033942298 
 churn rate: ChurnIndicator    0.00336
dtype: float64 
Test set: 
 x_test: (75128, 37) 
 y_test: (75128, 1) 
 test size: 0.3000019966057702 
 churn rate: ChurnIndicator    0.003354
dtype: float64


# define random under-sampling

rus = RandomUnderSampler(
    sampling_strategy=0.05,
    random_state=RANDOM_STATE
)

# define XGBoost algorithm

algo_xgb = xgboost.XGBClassifier(
    max_depth=6,
    learning_rate=0.3,
    objective='binary:logistic',
    tree_method='hist',
    sub_sample=1.0,
    random_state=RANDOM_STATE
)

# create pipeline with random sampling

pipe_rus_xgb = make_pipeline(rus, algo_xgb)


# fit model on training set

pipe_rus_xgb.fit(X_retention_train, y_retention_train)

Pipeline(steps=[('randomundersampler',
                 RandomUnderSampler(random_state=42, sampling_strategy=0.05)),
                ('xgbclassifier',
                 XGBClassifier(base_score=None, booster=None, callbacks=None,
                               colsample_bylevel=None, colsample_bynode=None,
                               colsample_bytree=None, device=None,
                               early_stopping_rounds=None,
                               enable_categorical=False, eval_metric=None,
                               feature_types=None, gamma=None, grow_policy=None,
                               importance_type=None,
                               interaction_constraints=None, learning_rate=0.3,
                               max_bin=None, max_cat_threshold=None,
                               max_cat_to_onehot=None, max_delta_step=None,
                               max_depth=6, max_leaves=None,
                               min_child_weight=None, missing=nan,
                               monotone_constraints=None, multi_strategy=None,
                               n_estimators=None, n_jobs=None,
                               num_parallel_tree=None, random_state=42, ...))])

Pipeline(steps=[('randomundersampler',
                 RandomUnderSampler(random_state=42, sampling_strategy=0.05)),
                ('xgbclassifier',
                 XGBClassifier(base_score=None, booster=None, callbacks=None,
                               colsample_bylevel=None, colsample_bynode=None,
                               colsample_bytree=None, device=None,
                               early_stopping_rounds=None,
                               enable_categorical=False, eval_metric=None,
                               feature_types=None, gamma=None, grow_policy=None,
                               importance_type=None,
                               interaction_constraints=None, learning_rate=0.3,
                               max_bin=None, max_cat_threshold=None,
                               max_cat_to_onehot=None, max_delta_step=None,
                               max_depth=6, max_leaves=None,
                               min_child_weight=None, missing=nan,
                               monotone_constraints=None, multi_strategy=None,
                               n_estimators=None, n_jobs=None,
                               num_parallel_tree=None, random_state=42, ...))])

RandomUnderSampler(random_state=42, sampling_strategy=0.05)

XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric=None, feature_types=None,
              gamma=None, grow_policy=None, importance_type=None,
              interaction_constraints=None, learning_rate=0.3, max_bin=None,
              max_cat_threshold=None, max_cat_to_onehot=None,
              max_delta_step=None, max_depth=6, max_leaves=None,
              min_child_weight=None, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=None, n_jobs=None,
              num_parallel_tree=None, random_state=42, ...)


# get classification report on test set at maximum F2

y_retention_prob = pipe_rus_xgb.predict_proba(X_retention_test)[:, 1]
xgb_test = get_classification_report(y_retention_test, y_retention_prob)
xgb_test[xgb_test['f1'].eq(xgb_test['f1'].max())]


# get churn score threshold at maximum F1 score

lifelong_score_threshold = xgb_test[xgb_test['f1'].eq(xgb_test['f1'].max())]['threshold'].iloc[0]
lifelong_score_threshold

0.98810136


# get probability prediction for each customer

lifelong_churn_score_df = X_retention_data.copy()
lifelong_churn_score_df['predict_proba'] = pipe_rus_xgb.predict_proba(lifelong_churn_score_df)[:, 1]


# get customers predicted as churn - customers who have score (probability) > churn threshold

lifelong_high_churn_score_df = lifelong_churn_score_df[lifelong_churn_score_df['predict_proba'].ge(lifelong_score_threshold)]
lifelong_churn_cx = lifelong_high_churn_score_df.index.get_level_values('LoyaltyNumber').unique().to_list()
len(lifelong_churn_cx)

36


# get data of active standard customers

active_loyalty_df = loyalty_df[loyalty_df['LoyaltyNumber'].isin(active_standard_cx_1718)]
active_loyalty_df.shape

(14303, 33)


# get all churn active customers

all_active = loyalty_df[loyalty_df['LoyaltyNumber'].isin(active_cx_1718)]

churn_monthsofmembership = all_active[all_active['ChurnIndicator'].eq(1)]['MonthsOfMembership'].value_counts().reset_index()\
    .sort_values(by='MonthsOfMembership')
    
# distribution of months of membership

plt.figure(figsize=(7,5))
sns.lineplot(data=churn_monthsofmembership, x='MonthsOfMembership', y='count')
plt.xlabel('Months of membership')
plt.ylabel('Number of customers')
plt.yticks(np.arange(0, 700, 100))
plt.title('Distribution of all active customers')

Text(0.5, 1.0, 'Distribution of all active customers')


# get all churn active standard customers

churn_monthsofmembership_standard = active_loyalty_df[active_loyalty_df['ChurnIndicator'].eq(1)]\
    ['MonthsOfMembership'].value_counts().reset_index()\
    .sort_values(by='MonthsOfMembership')

# distribution of months of membership

plt.figure(figsize=(7,5))
sns.lineplot(data=churn_monthsofmembership_standard, x='MonthsOfMembership', y='count')
plt.xlabel('Months of membership')
plt.ylabel('Number of customers')
plt.yticks(np.arange(0, 700, 100))
plt.title('Distribution of active standard customers')

Text(0.5, 1.0, 'Distribution of active standard customers')


# choose features for frequent pattern mining

fpm_8m = active_loyalty_df.pipe(lambda x: x[x['MonthsOfMembership'].le(8)]).pipe(lambda x: x[[
    'ChurnIndicator', 
    'Education', 'Salary', 'MaritalStatus', 'LoyaltyCard',
    'MonthlyTotalFlights', 'PercentagePointsRedeemed'
]])

fpm_8m.describe()


fpm_8m[fpm_8m['ChurnIndicator'].eq(1)].describe()


# bin numeric features

salary_bins = [-np.inf, 0, 15000, 60000, 80000, np.inf] # min, mean-std, mean, 75%, max
monthlyflights_bins = [-np.inf, 0, 4, 6, 8, np.inf] # min, 25%, mean, max 
pctpointsredeemed_bins = [-np.inf, 0, 0.4, np.inf] # min, 25%, 50%, >75%, pct=1, max

fpm_8m = fpm_8m.assign(
    SalaryBins = lambda x: pd.cut(x=x['Salary'], bins=salary_bins, right=True),
    # MonthsOfMembershipBins = lambda x: pd.cut(x=x['MonthsOfMembership'], bins=membership_bins, right=True),
    MonthlyFlightsBins = lambda x: pd.cut(x=x['MonthlyTotalFlights'], bins=monthlyflights_bins, right=True),
    PctPointsRedeemedBins = lambda x: pd.cut(x=x['PercentagePointsRedeemed'], bins=pctpointsredeemed_bins, right=True)
).astype(str).replace({'(-inf, 0.0]': '0'})\
.drop(columns=['Salary', 'MonthlyTotalFlights', 'PercentagePointsRedeemed'])\

fpm_8m.head()


# one-hot encoding

fpm_churn_6m = pd.get_dummies(fpm_8m[fpm_8m['ChurnIndicator'].eq('1')].drop(columns='ChurnIndicator'))
fpm_nonchurn_6m = pd.get_dummies(fpm_8m[fpm_8m['ChurnIndicator'].eq('0')].drop(columns='ChurnIndicator'))
display(fpm_churn_6m.shape, fpm_nonchurn_6m.shape)

(528, 22)

(1547, 23)


# get churn customers' unique frequent patterns for various scenarios with minimum support from 0.1 to 0.9

for i in np.arange(0.1, 1, 0.1):
    min_support = round(i, ndigits=2)
    patterns_churn = apriori(fpm_churn_6m, min_support=min_support, use_colnames=True).assign(length = lambda x: x['itemsets'].str.len())
    patterns_nonchurn = apriori(fpm_nonchurn_6m, min_support=min_support, use_colnames=True).assign(length = lambda x: x['itemsets'].str.len())
    unique_patterns_churn = patterns_churn[~patterns_churn['itemsets'].isin(patterns_nonchurn['itemsets'])]
    print('min support =', min_support)
    display(unique_patterns_churn)

min support = 0.1

min support = 0.2

min support = 0.3

min support = 0.4

min support = 0.5

min support = 0.6

min support = 0.7

min support = 0.8

min support = 0.9


# Filter active standard customers with: 
# - Year 2018 
# - Customers with more than 6 months of membership 
# - Remove records after customers cancel memberships

predict_retention_df = lookback_flight_df.merge(temp_mom_lookback, how='left', on=['LoyaltyNumber', 'Date'])\
    .pipe(lambda x: x[x['Date'].ge('2018-01-01') & x['Date'].ge(x['EnrollmentDate'])])\
    .pipe(lambda x: x[~x['Date'].gt(x['CancellationDate'])])\
    .assign(ChurnIndicator = lambda x: np.where(x['Date'].eq(x['CancellationDate']), 1, 0))\
    .pipe(lambda x: x[x['MonthsOfMembership'].gt(6)])\
    .pipe(lambda x: x[x['LoyaltyNumber'].isin(active_standard_cx_1718)])
    
predict_retention_df


predict_retention_df['LoyaltyNumber'].nunique()

12480


# Join demographic profiles

predict_retention_df = predict_retention_df.merge(loyalty_df.iloc[:, :9], how='left', on='LoyaltyNumber')\
    .drop(columns=['EnrollmentDate', 'CancellationDate', 'Country'])
    
predict_retention_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 231031 entries, 0 to 231030
Data columns (total 39 columns):
 #   Column                           Non-Null Count   Dtype         
---  ------                           --------------   -----         
 0   LoyaltyNumber                    231031 non-null  object        
 1   Date                             231031 non-null  datetime64[ns]
 2   Last1Y_FlightsBooked             231031 non-null  float64       
 3   Last1Y_FlightsWithCompanions     231031 non-null  float64       
 4   Last1Y_TotalFlights              231031 non-null  float64       
 5   Last1Y_Distance                  231031 non-null  float64       
 6   Last1Y_PointsAccumulated         231031 non-null  float64       
 7   Last1Y_PointsRedeemed            231031 non-null  float64       
 8   Last1Y_DollarCostPointsRedeemed  231031 non-null  float64       
 9   Last6M_FlightsBooked             231031 non-null  float64       
 10  Last6M_FlightsWithCompanions     231031 non-null  float64       
 11  Last6M_TotalFlights              231031 non-null  float64       
 12  Last6M_Distance                  231031 non-null  float64       
 13  Last6M_PointsAccumulated         231031 non-null  float64       
 14  Last6M_PointsRedeemed            231031 non-null  float64       
 15  Last6M_DollarCostPointsRedeemed  231031 non-null  float64       
 16  Last3M_FlightsBooked             231031 non-null  float64       
 17  Last3M_FlightsWithCompanions     231031 non-null  float64       
 18  Last3M_TotalFlights              231031 non-null  float64       
 19  Last3M_Distance                  231031 non-null  float64       
 20  Last3M_PointsAccumulated         231031 non-null  float64       
 21  Last3M_PointsRedeemed            231031 non-null  float64       
 22  Last3M_DollarCostPointsRedeemed  231031 non-null  float64       
 23  Last1M_FlightsBooked             231031 non-null  float64       
 24  Last1M_FlightsWithCompanions     231031 non-null  float64       
 25  Last1M_TotalFlights              231031 non-null  float64       
 26  Last1M_Distance                  231031 non-null  float64       
 27  Last1M_PointsAccumulated         231031 non-null  float64       
 28  Last1M_PointsRedeemed            231031 non-null  float64       
 29  Last1M_DollarCostPointsRedeemed  231031 non-null  float64       
 30  MonthsOfMembership               231031 non-null  int64         
 31  ChurnIndicator                   231031 non-null  int32         
 32  Province                         231031 non-null  object        
 33  City                             231031 non-null  object        
 34  PostalCode                       231031 non-null  object        
 35  Gender                           231031 non-null  object        
 36  Education                        231031 non-null  object        
 37  Salary                           231031 non-null  float64       
 38  MaritalStatus                    231031 non-null  object        
dtypes: datetime64[ns](1), float64(29), int32(1), int64(1), object(7)
memory usage: 67.9+ MB


# Label encoding for categorical features

label_encoder = LabelEncoder()

for col in predict_retention_df.columns[1:]:
    if predict_retention_df[col].dtype == 'object':
        predict_retention_df[col] = label_encoder.fit_transform(predict_retention_df[col])


predict_retention_df = predict_retention_df.reset_index()
predict_retention_df


predict_retention_df = predict_retention_df.set_index(['LoyaltyNumber', 'Date'])
predict_retention_df['ChurnIndicator'].value_counts(normalize=True)

ChurnIndicator
0    0.996438
1    0.003562
Name: proportion, dtype: float64


predict_retention_df.index.get_level_values('LoyaltyNumber').nunique()

12480


# get X and y

X_retention_data = predict_retention_df.drop(columns='ChurnIndicator')
y_retention_data = predict_retention_df[['ChurnIndicator']]


# get training and test sets

X_retention_train, X_retention_test, y_retention_train, y_retention_test = train_test_split(
    X_retention_data,
    y_retention_data,
    test_size=0.3,
    random_state=RANDOM_STATE,
    stratify=y_retention_data
)

print(
    'Train set:',
    '\n x_train:', X_retention_train.shape, 
    '\n y_train:', y_retention_train.shape, 
    '\n train size:', y_retention_train.shape[0]/y_retention_data.shape[0], 
    '\n churn rate:', y_retention_train.sum()/y_retention_train.count(),
    '\nTest set:',
    '\n x_test:', X_retention_test.shape, 
    '\n y_test:', y_retention_test.shape, 
    '\n test size:', y_retention_test.shape[0]/y_retention_data.shape[0], 
    '\n churn rate:', y_retention_test.sum()/y_retention_test.count()
)

Train set: 
 x_train: (161721, 37) 
 y_train: (161721, 1) 
 train size: 0.6999969701035792 
 churn rate: ChurnIndicator    0.003562
dtype: float64 
Test set: 
 x_test: (69310, 37) 
 y_test: (69310, 1) 
 test size: 0.3000030298964208 
 churn rate: ChurnIndicator    0.003564
dtype: float64


# define random under-sampling

rus = RandomUnderSampler(
    sampling_strategy=0.05,
    random_state=RANDOM_STATE
)

# define XGBoost algorithm

algo_xgb = xgboost.XGBClassifier(
    max_depth=6,
    learning_rate=0.3,
    objective='binary:logistic',
    tree_method='hist',
    sub_sample=1.0,
    random_state=RANDOM_STATE
)

# create pipeline with random sampling

pipe_rus_xgb = make_pipeline(rus, algo_xgb)


# fit model on training set

pipe_rus_xgb.fit(X_retention_train, y_retention_train)

Pipeline(steps=[('randomundersampler',
                 RandomUnderSampler(random_state=42, sampling_strategy=0.05)),
                ('xgbclassifier',
                 XGBClassifier(base_score=None, booster=None, callbacks=None,
                               colsample_bylevel=None, colsample_bynode=None,
                               colsample_bytree=None, device=None,
                               early_stopping_rounds=None,
                               enable_categorical=False, eval_metric=None,
                               feature_types=None, gamma=None, grow_policy=None,
                               importance_type=None,
                               interaction_constraints=None, learning_rate=0.3,
                               max_bin=None, max_cat_threshold=None,
                               max_cat_to_onehot=None, max_delta_step=None,
                               max_depth=6, max_leaves=None,
                               min_child_weight=None, missing=nan,
                               monotone_constraints=None, multi_strategy=None,
                               n_estimators=None, n_jobs=None,
                               num_parallel_tree=None, random_state=42, ...))])

Pipeline(steps=[('randomundersampler',
                 RandomUnderSampler(random_state=42, sampling_strategy=0.05)),
                ('xgbclassifier',
                 XGBClassifier(base_score=None, booster=None, callbacks=None,
                               colsample_bylevel=None, colsample_bynode=None,
                               colsample_bytree=None, device=None,
                               early_stopping_rounds=None,
                               enable_categorical=False, eval_metric=None,
                               feature_types=None, gamma=None, grow_policy=None,
                               importance_type=None,
                               interaction_constraints=None, learning_rate=0.3,
                               max_bin=None, max_cat_threshold=None,
                               max_cat_to_onehot=None, max_delta_step=None,
                               max_depth=6, max_leaves=None,
                               min_child_weight=None, missing=nan,
                               monotone_constraints=None, multi_strategy=None,
                               n_estimators=None, n_jobs=None,
                               num_parallel_tree=None, random_state=42, ...))])

RandomUnderSampler(random_state=42, sampling_strategy=0.05)

XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=None, device=None, early_stopping_rounds=None,
              enable_categorical=False, eval_metric=None, feature_types=None,
              gamma=None, grow_policy=None, importance_type=None,
              interaction_constraints=None, learning_rate=0.3, max_bin=None,
              max_cat_threshold=None, max_cat_to_onehot=None,
              max_delta_step=None, max_depth=6, max_leaves=None,
              min_child_weight=None, missing=nan, monotone_constraints=None,
              multi_strategy=None, n_estimators=None, n_jobs=None,
              num_parallel_tree=None, random_state=42, ...)


# get classification report on test set at maximum F2

y_retention_prob = pipe_rus_xgb.predict_proba(X_retention_test)[:, 1]
xgb_test = get_classification_report(y_retention_test, y_retention_prob)
xgb_test[xgb_test['f1'].eq(xgb_test['f1'].max())]


xgb_test[xgb_test['f1'].eq(xgb_test['f1'].max())].to_clipboard()


# get customers predicted as churn in stage 1 (with minimum_support = 0.3)
stage1_churn_cx = active_loyalty_df.pipe(lambda x: x[x['MonthsOfMembership'].le(8)]).pipe(lambda x: x[
    x['MaritalStatus'].eq('Single')
])['LoyaltyNumber'].to_list()
len(stage1_churn_cx)

569


# get churn threshold at maximum F1 score

stage2_threshold = xgb_test[xgb_test['f1'].eq(xgb_test['f1'].max())]['threshold'].iloc[0]
stage2_threshold

0.9878924


# get probability prediction for each customer

stage2_churn_score_df = X_retention_data.copy()
stage2_churn_score_df['predict_proba'] = pipe_rus_xgb.predict_proba(stage2_churn_score_df)[:, 1]


# get customers predicted as churn in stage 2 - customers who have score (probability) > churn threshold

stage2_high_churn_score_df = stage2_churn_score_df[stage2_churn_score_df['predict_proba'].ge(stage2_threshold)]
stage2_churn_cx = stage2_high_churn_score_df.index.get_level_values('LoyaltyNumber').unique().to_list()
len(stage2_churn_cx)

80


# all customers predicted as churn

dual_phase_churn_cx = stage1_churn_cx + stage2_churn_cx


# Original churn rate

display(
    active_loyalty_df[['ChurnIndicator']].value_counts(normalize=False),
    active_loyalty_df[['ChurnIndicator']].value_counts(normalize=True)
)

ChurnIndicator
0                 13288
1                  1015
Name: count, dtype: int64

ChurnIndicator
0                 0.929036
1                 0.070964
Name: proportion, dtype: float64


# add columns of churn prediction from 2 approaches

active_loyalty_df = active_loyalty_df.assign(
    lifelong_ChurnPrediction = lambda x: np.where(x['LoyaltyNumber'].isin(lifelong_churn_cx), 1, 0),
    dualphase_ChurnPrediction = lambda x: np.where(x['LoyaltyNumber'].isin(dual_phase_churn_cx), 1, 0)
)


# churn_ind_list = ['ChurnIndicator', 'lifelong_ChurnPrediction', 'dualphase_ChurnPrediction']
total_number_of_customers = active_loyalty_df['LoyaltyNumber'].nunique()
retention_strategies = {
    'all': 'ChurnIndicator',
    'lifelong': 'lifelong_ChurnPrediction',
    'dualphase': 'dualphase_ChurnPrediction'
}

retention_stats_dict = {}

for strategy, col in retention_strategies.items():
    # number of unique customers predicted as churn
    if strategy == 'all':
        predicted_churn = total_number_of_customers
    else:
        predicted_churn = active_loyalty_df[col].sum()
    # true positive - number of churn correctly predicted as churn (to receive promotion -> convert to nonchurn)
    true_positive = active_loyalty_df[active_loyalty_df['ChurnIndicator'].eq(1) & active_loyalty_df[col].eq(1)].shape[0]
    # false negative - number of predicted nonchurn but are actual churn (missing alert)
    false_negative = active_loyalty_df[active_loyalty_df['ChurnIndicator'].eq(1) & active_loyalty_df[col].eq(0)].shape[0]
    # churn rate after promotion = false negative / total number of customers
    churn_rate_after_promo = round(false_negative * 100 / total_number_of_customers, ndigits=1)
    # cost per successful retention
    cost_efficiency = round(predicted_churn / true_positive, ndigits=1)
    
    retention_stats_dict[strategy] = [predicted_churn, true_positive, false_negative, churn_rate_after_promo, cost_efficiency]
    

retention_stats_dict_df = pd.DataFrame(
    retention_stats_dict, 
    index=['predicted_churn', 'true_positive', 'false_negative', 'churn_rate_after_promo', 'cost_efficiency']    
)    
retention_stats_dict_df


loyalty_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 16737 entries, 0 to 16736
Data columns (total 33 columns):
 #   Column                         Non-Null Count  Dtype         
---  ------                         --------------  -----         
 0   LoyaltyNumber                  16737 non-null  object        
 1   Country                        16737 non-null  object        
 2   Province                       16737 non-null  object        
 3   City                           16737 non-null  object        
 4   PostalCode                     16737 non-null  object        
 5   Gender                         16737 non-null  object        
 6   Education                      16737 non-null  object        
 7   Salary                         16737 non-null  float64       
 8   MaritalStatus                  16737 non-null  object        
 9   LoyaltyCard                    16737 non-null  object        
 10  CLV                            16737 non-null  float64       
 11  EnrollmentType                 16737 non-null  object        
 12  EnrollmentYear                 16737 non-null  int64         
 13  EnrollmentMonth                16737 non-null  int64         
 14  CancellationYear               2067 non-null   Int64         
 15  CancellationMonth              2067 non-null   Int64         
 16  EnrollmentDate                 16737 non-null  datetime64[ns]
 17  EnrollmentYearMonth            16737 non-null  object        
 18  CancellationDate               2067 non-null   datetime64[ns]
 19  CancellationYearMonth          2067 non-null   object        
 20  ChurnIndicator                 16737 non-null  int32         
 21  MonthsOfMembership             16737 non-null  int64         
 22  TotalFlightsBooked             16737 non-null  int64         
 23  TotalFlightsWithCompanions     16737 non-null  int64         
 24  TotalTotalFlights              16737 non-null  int64         
 25  TotalDistance                  16737 non-null  int64         
 26  TotalPointsAccumulated         16737 non-null  float64       
 27  TotalPointsRedeemed            16737 non-null  int64         
 28  TotalDollarCostPointsRedeemed  16737 non-null  int64         
 29  PercentagePointsRedeemed       16737 non-null  float64       
 30  MonthlyTotalFlights            16737 non-null  float64       
 31  MonthlyPointsAccumulated       16737 non-null  float64       
 32  ChurnIndicator_bool            16737 non-null  bool          
dtypes: Int64(2), bool(1), datetime64[ns](2), float64(6), int32(1), int64(9), object(12)
memory usage: 4.1+ MB


loyalty_df.describe()


loyalty_df.head()


c_df = loyalty_df.drop(loyalty_df.columns[12:20], axis=1)
c_df = c_df.sort_index()


c_df_k1 = c_df.copy()


kmeans = KMeans(n_clusters=4)
labels =['low','mid-low','mid-high', 'high']

X = np.array(c_df_k1['CLV']).reshape(-1,1)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)
print(y_kmeans)

clist=np.array([y_kmeans[i] for i in c_df_k1['CLV'].sort_index().sort_values().index])
print(clist)

values, counts = np.unique(y_kmeans, return_counts=True)
print(values,counts)

mapping = np.append(clist[0], clist[1:][clist[1:] != clist[:-1]])
c_df_k1['CLV_cluster']=[labels[int(np.where(mapping == i)[0])] for i in y_kmeans]
c_df_k1.head()

[0 0 0 ... 3 2 3]
[0 0 0 ... 3 3 3]
[0 1 2 3] [9570 1437 5386  344]


from sklearn import datasets

X = np.array(c_df_k1['CLV']).reshape(-1,1)

# Calculate within-cluster sum of squares (WCSS) for different values of k
wcss = []
for i in range(1, 11):
    kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=0)
    kmeans.fit(X)
    wcss.append(kmeans.inertia_)

# Plot the elbow method curve
plt.plot(range(1, 11), wcss)
plt.title('Elbow Method')
plt.xlabel('Number of Clusters')
plt.ylabel('WCSS')
plt.show()


c_df_k2 = c_df.copy()


c_df_k2['AvgDistPerFlight'] = c_df_k2['TotalDistance'] / c_df_k2['TotalTotalFlights']
c_df_k2['AvgDistPerFlight'].fillna(0, inplace=True)


kmeans = KMeans(n_clusters=5)
labels =['short','mid-short', 'mid', 'mid-high', 'long']

X = np.array(c_df_k2['AvgDistPerFlight']).reshape(-1,1)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)
print(y_kmeans)

clist=np.array([y_kmeans[i] for i in c_df_k2['AvgDistPerFlight'].sort_values().index])
print(clist)

values, counts = np.unique(y_kmeans, return_counts=True)
print(values,counts)

mapping = np.append(clist[0], clist[1:][clist[1:] != clist[:-1]])
c_df_k2['AvgDist_cluster']=[labels[int(np.where(mapping == i)[0])] for i in y_kmeans]
c_df_k2.head()

[4 0 0 ... 0 4 2]
[2 2 2 ... 1 1 1]
[0 1 2 3 4] [9736   13 1578   92 5318]


X = np.array(c_df_k2['AvgDistPerFlight']).reshape(-1,1)

# Calculate within-cluster sum of squares (WCSS) for different values of k
wcss = []
for i in range(1, 11):
    kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=0)
    kmeans.fit(X)
    wcss.append(kmeans.inertia_)

# Plot the elbow method curve
plt.plot(range(1, 11), wcss)
plt.title('Elbow Method')
plt.xlabel('Number of Clusters')
plt.ylabel('WCSS')
plt.show()


c_df_k3 = c_df.copy()


c_df_k3['PointUsage'] = c_df_k3['TotalPointsRedeemed'] / c_df_k3['TotalPointsAccumulated']
c_df_k3['PointUsage'].fillna(0, inplace=True)


X = np.array(c_df_k3['PointUsage']).reshape(-1,1)

# Calculate within-cluster sum of squares (WCSS) for different values of k
wcss = []
for i in range(1, 11):
    kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=0)
    kmeans.fit(X)
    wcss.append(kmeans.inertia_)

# Plot the elbow method curve
plt.plot(range(1, 11), wcss)
plt.title('Elbow Method')
plt.xlabel('Number of Clusters')
plt.ylabel('WCSS')
plt.show()


kmeans = KMeans(n_clusters=3)
labels =['low', 'mid', 'high']

X = np.array(c_df_k3['PointUsage']).reshape(-1,1)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)
print(y_kmeans)

clist=np.array([y_kmeans[i] for i in c_df_k3['PointUsage'].sort_values().index])
print(clist)

values, counts = np.unique(y_kmeans, return_counts=True)
print(values,counts)

mapping = np.append(clist[0], clist[1:][clist[1:] != clist[:-1]])
c_df_k3['PU_cluster']=[labels[int(np.where(mapping == i)[0])] for i in y_kmeans]
c_df_k3.head()

[1 2 1 ... 2 1 1]
[1 1 1 ... 0 0 0]
[0 1 2] [  106 10740  5891]


# Group by CLV_cluster and Gender
grouped = c_df_k1.groupby(["CLV_cluster", "Gender"])

# Count occurrences of each gender within each cluster
cluster_gender_counts = grouped.size().unstack()

# Plot bar chart
cluster_gender_counts.plot(kind='bar', color=['b', 'r'], alpha=0.7)

# Set labels and title
plt.xlabel('CLV Cluster')
plt.ylabel('Count')
plt.title('Distribution of CLV Cluster by Gender')

# Show plot
plt.legend(title='Gender')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


c = np.array(['b', 'g', 'r', 'c'])
labels =['low','mid-low','mid-high', 'high']

grouped = c_df_k1.groupby("CLV_cluster")

for name, group in grouped:
    plt.plot(group["Salary"],color = c[np.array(labels) == name][0], marker='o', 
             linestyle ='', data = c_df_k1)


grouped = c_df_k1.groupby(["CLV_cluster", "Education"])

cluster_edu_counts = grouped.size().unstack()

cluster_edu_counts.plot(kind='bar', color=['b', 'r', 'g', 'c', 'm'], alpha=0.7)

plt.xlabel('CLV Cluster')
plt.ylabel('Count')
plt.title('Distribution of CLV Cluster by Education')

plt.legend(title='Education')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k1.groupby(["CLV_cluster", "LoyaltyCard"])

cluster_loyalty_counts = grouped.size().unstack()

cluster_loyalty_counts.plot(kind='bar', color=['b', 'r', 'g'], alpha=0.7)

plt.xlabel('CLV Cluster')
plt.ylabel('Count')
plt.title('Distribution of CLV Cluster by LoyaltyCard')

plt.legend(title='LoyaltyCard')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k1.groupby(["CLV_cluster", "MaritalStatus"])

cluster_loyalty_counts = grouped.size().unstack()

cluster_loyalty_counts.plot(kind='bar', color=['b', 'r', 'g'], alpha=0.7)

plt.xlabel('CLV Cluster')
plt.ylabel('Count')
plt.title('Distribution of CLV Cluster by MaritalStatus')

plt.legend(title='LoyaltyCard')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


# Group by AvgDist_cluster and Gender
grouped = c_df_k2.groupby(["AvgDist_cluster", "Gender"])

cluster_gender_counts = grouped.size().unstack()
cluster_gender_counts.plot(kind='bar', color=['b', 'r'], alpha=0.7)
plt.xlabel('AvgDist Cluster')
plt.ylabel('Count')
plt.title('Distribution of AvgDist Cluster by Gender')
plt.legend(title='Gender')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


c = np.array(['b', 'g', 'r', 'c', 'm'])
labels =['short','mid-short', 'mid', 'mid-high', 'long']

grouped = c_df_k2.groupby("AvgDist_cluster")

for name, group in grouped:
    plt.plot(group["Salary"],color = c[np.array(labels) == name][0], marker='o', 
             linestyle ='', data = c_df_k1)


grouped = c_df_k2.groupby(["AvgDist_cluster", "Education"])

cluster_edu_counts = grouped.size().unstack()

cluster_edu_counts.plot(kind='bar', color=['b', 'r', 'g', 'c', 'm'], alpha=0.7)

plt.xlabel('AvgDist Cluster')
plt.ylabel('Count')
plt.title('Distribution of AvgDist Cluster by Education')

plt.legend(title='Education')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k2.groupby(["AvgDist_cluster", "LoyaltyCard"])

cluster_loyalty_counts = grouped.size().unstack()

cluster_loyalty_counts.plot(kind='bar', color=['b', 'r', 'g'], alpha=0.7)

plt.xlabel('AvgDist Cluster')
plt.ylabel('Count')
plt.title('Distribution of AvgDist Cluster by LoyaltyCard')

plt.legend(title='LoyaltyCard')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k2.groupby(["AvgDist_cluster", "EnrollmentType"])

cluster_enrol_counts = grouped.size().unstack()

cluster_enrol_counts.plot(kind='bar', color=['b', 'r'], alpha=0.7)

plt.xlabel('AvgDist Cluster')
plt.ylabel('Count')
plt.title('Distribution of AvgDist Cluster by EnrollmentType')

plt.legend(title='EnrollmentType')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k3.groupby(["PU_cluster", "LoyaltyCard"])

cluster_loyalty_counts = grouped.size().unstack()

cluster_loyalty_counts.plot(kind='bar', color=['b', 'r', 'g'], alpha=0.7)

plt.xlabel('PU Cluster')
plt.ylabel('Count')
plt.title('Distribution of PU Cluster by LoyaltyCard')

plt.legend(title='LoyaltyCard')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k3.groupby(["PU_cluster", "EnrollmentType"])

cluster_enrol_counts = grouped.size().unstack()

cluster_enrol_counts.plot(kind='bar', color=['b', 'r'], alpha=0.7)

plt.xlabel('PU Cluster')
plt.ylabel('Count')
plt.title('Distribution of PU Cluster by EnrollmentType')

plt.legend(title='EnrollmentType')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k3.groupby(["PU_cluster", "Education"])

cluster_edu_counts = grouped.size().unstack()

cluster_edu_counts.plot(kind='bar', color=['b', 'r', 'g', 'c', 'm'], alpha=0.7)

plt.xlabel('PU Cluster')
plt.ylabel('Count')
plt.title('Distribution of PU Cluster by Education')

plt.legend(title='Education')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


grouped = c_df_k3.groupby(["PU_cluster", "MaritalStatus"])

cluster_edu_counts = grouped.size().unstack()

cluster_edu_counts.plot(kind='bar', color=['b', 'r', 'g', 'c', 'm'], alpha=0.7)

plt.xlabel('PU Cluster')
plt.ylabel('Count')
plt.title('Distribution of PU Cluster by Marital Status')

plt.legend(title='Education')
plt.xticks(rotation=0)  # Rotate x-axis labels if needed
plt.tight_layout()  # Adjust layout to prevent clipping of labels
plt.show()


from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler


loyalty_df_cluster = loyalty_df.copy()


loyalty_df_cluster.columns

Index(['LoyaltyNumber', 'Country', 'Province', 'City', 'PostalCode', 'Gender',
       'Education', 'Salary', 'MaritalStatus', 'LoyaltyCard', 'CLV',
       'EnrollmentType', 'EnrollmentYear', 'EnrollmentMonth',
       'CancellationYear', 'CancellationMonth', 'EnrollmentDate',
       'EnrollmentYearMonth', 'CancellationDate', 'CancellationYearMonth',
       'ChurnIndicator', 'MonthsOfMembership', 'TotalFlightsBooked',
       'TotalFlightsWithCompanions', 'TotalTotalFlights', 'TotalDistance',
       'TotalPointsAccumulated', 'TotalPointsRedeemed',
       'TotalDollarCostPointsRedeemed', 'PercentagePointsRedeemed',
       'MonthlyTotalFlights', 'MonthlyPointsAccumulated',
       'ChurnIndicator_bool'],
      dtype='object')


# Extracting only the numerical features
cluster_features = ['Salary','CLV','MonthsOfMembership','TotalFlightsBooked',
       'TotalFlightsWithCompanions', 'TotalTotalFlights', 'TotalDistance',
       'TotalPointsAccumulated', 'TotalPointsRedeemed',
       'TotalDollarCostPointsRedeemed']


# Select features for clustering
X = loyalty_df_cluster.loc[:,cluster_features]


# Fit standard scaler
scaler = StandardScaler()
scaled_features = scaler.fit_transform(X)


# Apply PCA to reduce dimensionality to 2D
pca = PCA(n_components=3)
X_pca = pca.fit_transform(scaled_features)

# Get the components
components = pca.components_

# Create a DataFrame to view component values in order to identify more important features for clustering
component_df = pd.DataFrame(components, columns=cluster_features, index=['PC1', 'PC2', 'PC3'])
component_df


# Create a graph to determine the number of Principal Components and the features which contribute to its explained variance
explained_variance = pca.explained_variance_ratio_
plt.plot(range(1, len(explained_variance) + 1), explained_variance.cumsum(), marker='o')
plt.xlabel('Number of Principal Components')
plt.ylabel('Cumulative Explained Variance')
plt.show()


# Keeping only certain features based on their importance in the PCA
cluster_features = ['TotalFlightsBooked', 'TotalTotalFlights', 'TotalDistance',
       'TotalPointsAccumulated', 'TotalPointsRedeemed',
       'TotalDollarCostPointsRedeemed']


# Select features for clustering
X = loyalty_df_cluster.loc[:,cluster_features]


# Fit standard scaler
scaler = StandardScaler()
scaled_features = scaler.fit_transform(X)


# Calculate within-cluster sum of squares (WCSS) for different values of k
wcss = []
for i in range(1, 11):
    kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=2024)
    kmeans.fit(scaled_features)
    wcss.append(kmeans.inertia_)

# Plot the elbow method curve
plt.plot(range(1, 11), wcss)
plt.title('Elbow Method')
plt.xlabel('Number of Clusters')
plt.ylabel('WCSS')
plt.show()


# Calculate  silhouette coefficient for different values of k
silhoutte_coeff = []
for i in range(2, 11):
    kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=2024)
    kmeans.fit(scaled_features)
    silhoutte_coeff.append(silhouette_score(X,kmeans.labels_))

# Plot the elbow method curve
plt.plot(range(2, 11), silhoutte_coeff)
plt.title('Elbow Method')
plt.xlabel('Number of Clusters')
plt.ylabel('silhouette')
plt.show()


# Number of clusters
k = 3

# Initialize model
kmean = KMeans(n_clusters=k, init='k-means++', max_iter=300, n_init=10, random_state=2024)

# fit model
kmean.fit(scaled_features)

# get labels
cluster_labels = kmean.labels_

# Add cluster labels to the DF and adjust the cluster number so that they are more intuitive (cluster 0 will become cluster 1, for example)
loyalty_df_cluster['cluster'] = cluster_labels+1


# Print the count of points in each cluster
pd.concat([loyalty_df_cluster['cluster'].value_counts().sort_index().reset_index(),
loyalty_df_cluster['cluster'].value_counts(normalize=True).sort_index().reset_index()],axis = 1)


# Graph by loyalty card
fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(20, 20));

c = 0
for i in range(0,3):
    for j in range(0, 2):
        # This determines the x_lim and y_lim by considering the values of the features
        cluster_0 = loyalty_df_cluster[loyalty_df_cluster['cluster'] == 1][cluster_features[c]]
        cluster_1 = loyalty_df_cluster[loyalty_df_cluster['cluster'] == 2][cluster_features[c]]
        cluster_2 = loyalty_df_cluster[loyalty_df_cluster['cluster'] == 3][cluster_features[c]]
        cluster_3 = loyalty_df_cluster[loyalty_df_cluster['cluster'] == 4][cluster_features[c]]
        cluster_4 = loyalty_df_cluster[loyalty_df_cluster['cluster'] == 5][cluster_features[c]]
        cluster_5 = loyalty_df_cluster[loyalty_df_cluster['cluster'] == 6][cluster_features[c]]
        min_c = min(cluster_0.min(),cluster_1.min(),cluster_2.min(),cluster_3.min(),cluster_4.min(),cluster_5.min())
        max_c = max(cluster_0.max(),cluster_1.max(),cluster_2.max(),cluster_3.max(),cluster_4.max(),cluster_5.max())

        # This automatically divides the feature into bins based on the entire dataset
        num_bin = list(np.histogram(loyalty_df_cluster[cluster_features[c]],bins=10)[1])
        
        sns.histplot(loyalty_df_cluster[loyalty_df_cluster['cluster']==1], x = cluster_features[c], ax = axs[i,j], stat = 'probability', alpha = 0.3, bins = num_bin, label = 'cluster_1')
        sns.histplot(loyalty_df_cluster[loyalty_df_cluster['cluster']==2], x = cluster_features[c], ax = axs[i,j], stat = 'probability', alpha = 0.3, bins = num_bin, label = 'cluster_2')
        sns.histplot(loyalty_df_cluster[loyalty_df_cluster['cluster']==3], x = cluster_features[c], ax = axs[i,j], stat = 'probability', alpha = 0.3, bins = num_bin, label = 'cluster_3')
        sns.histplot(loyalty_df_cluster[loyalty_df_cluster['cluster']==4], x = cluster_features[c], ax = axs[i,j], stat = 'probability', alpha = 0.3, bins = num_bin, label = 'cluster_4')
        sns.histplot(loyalty_df_cluster[loyalty_df_cluster['cluster']==5], x = cluster_features[c], ax = axs[i,j], stat = 'probability', alpha = 0.3, bins = num_bin, label = 'cluster_5')
        sns.histplot(loyalty_df_cluster[loyalty_df_cluster['cluster']==6], x = cluster_features[c], ax = axs[i,j], stat = 'probability', alpha = 0.3, bins = num_bin, label = 'cluster_6')
        axs[i, j].set_xlim(min_c,max_c)
        axs[i, j].legend(title='Cluster')
        c += 1
                

plt.tight_layout()


# numerical features for aggregation
num_features = ['Salary','CLV','MonthsOfMembership','TotalFlightsBooked',
       'TotalFlightsWithCompanions', 'TotalTotalFlights', 'TotalDistance',
       'TotalPointsAccumulated', 'TotalPointsRedeemed',
       'TotalDollarCostPointsRedeemed']


median_clusters = loyalty_df_cluster.groupby('cluster')[num_features].median()


# Function to apply colors based on value positions
def highlight_positions(s):
    is_min = s.isin(s.nsmallest(1))
    is_max = s.isin(s.nlargest(1))
    is_mid = ~is_min & ~is_max

    result = ['background-color: #dddcdc' if v else '' for v in is_min]
    result = ['background-color: #ffaf93' if v else res for v, res in zip(is_max, result)]
    result = ['background-color: #93fff5' if v else res for v, res in zip(is_mid, result)]
    
    return result


median_clusters = median_clusters.style.apply(highlight_positions)


median_clusters


fig, axs = plt.subplots(nrows = 3, ncols= 2, figsize = (20,20))

c = 0

for i in range(0,3):
    for j in range(0,2):
        if i*j == 2:
            sns.scatterplot(loyalty_df_cluster, x = cluster_features[0], y = cluster_features[-1], hue='cluster', ax = axs[i,j])
        else:    
            sns.scatterplot(loyalty_df_cluster, x = cluster_features[c], y = cluster_features[c+1], hue='cluster', ax = axs[i,j])
        c += 1

plt.show()
plt.tight_layout()

<Figure size 640x480 with 0 Axes>


loyalty_df_metric = loyalty_df_cluster.copy()


# Remove all customers with 0 months of membership and 0 flight flown
loyalty_df_metric = loyalty_df_metric[loyalty_df_metric['TotalDistance']!=0]


# Calculate Customers' Average Monthly CLV
loyalty_df_metric['MonthsOfMembership'] = loyalty_df_metric['MonthsOfMembership'].apply(lambda x: 1 if x== 0 else x)
loyalty_df_metric['AvgCLV'] = loyalty_df_metric['CLV']/loyalty_df_metric['MonthsOfMembership']


# Calculate Customers' Unit Revenue based on the Total Flights over the last 2 years
loyalty_df_metric['UnitRevenue'] = (loyalty_df_metric['AvgCLV']*24)/loyalty_df_metric['TotalTotalFlights']

loyalty_df_metric['UnitRevenue'] = np.where(
    (loyalty_df_metric['MonthsOfMembership'] == 1) & (loyalty_df_metric['TotalTotalFlights']==1),
    loyalty_df_metric['AvgCLV'],
    loyalty_df_metric['UnitRevenue']
)


# This function will determine the optimal customer index for the optimal coupon value
def promo_ob3_v2(customers):  
    marketing_cost = 0.2
    l = np.array(customers['UnitRevenue'].sort_values())
    return l[np.argmax(np.add.accumulate(l) - l*marketing_cost*len(customers))]*marketing_cost


# Calculate profit and revenue based on no segmentation
m_cost = 0.2

revenue_unsegmented = sum(loyalty_df_metric['UnitRevenue'][loyalty_df_metric['UnitRevenue']*m_cost <= promo_ob3_v2(loyalty_df_metric)])
profit_unsegmented = revenue_unsegmented - promo_ob3_v2(loyalty_df_metric)*len(loyalty_df_metric)

print(f"Without customer segmentation:\nProfit is {profit_unsegmented}") 
print(f"Profit margin is {profit_unsegmented/revenue_unsegmented}")
print(f'Total revenue is {revenue_unsegmented}')
print(f'Total cost is:', promo_ob3_v2(loyalty_df_metric)*len(loyalty_df_metric))

Without customer segmentation:
Profit is 38505.985186105885
Profit margin is 0.17774345978361278
Total revenue is 216637.9861907919
Total cost is: 178132.00100468603


loyalty_df_metric['target_customer'] = loyalty_df_metric['Education'] + loyalty_df_metric['LoyaltyCard'] + loyalty_df_metric['MaritalStatus']
loyalty_df_metric['target_customer_groups'] = loyalty_df_metric['target_customer'].apply(lambda x: 'target' if x == 'BachelorAuroraMarried' else 'other')


labels = list(set(loyalty_df_metric['target_customer']))


# Calculate profit and revenue based on grouped segmentations

def promo_ob3_v3(customers, group_label): 
    marketing_cost = 0.2
    l = np.array(customers[customers['target_customer']==group_label]['UnitRevenue'].sort_values())
    return l[np.argmax(np.add.accumulate(l) - l*marketing_cost*len(customers[customers['target_customer']==group_label]))]*marketing_cost


profits = []
revenue = []
cost = []
for i in labels:
    m_cost = 0.2
    revenue_segmented = sum(loyalty_df_metric['UnitRevenue'][(loyalty_df_metric['UnitRevenue']*m_cost <= promo_ob3_v3(loyalty_df_metric,i)) & (loyalty_df_metric['target_customer']==i)])
    profit_segmented = revenue_segmented - promo_ob3_v3(loyalty_df_metric,i)*len(loyalty_df_metric[loyalty_df_metric['target_customer']==i])
    cost.append(promo_ob3_v3(loyalty_df_metric,i)*len(loyalty_df_metric[loyalty_df_metric['target_customer']==i]))
    revenue.append(revenue_segmented)
    profits.append(profit_segmented)

print(f"With customer segmentation:\nProfit is {sum(profits)}") 
print(f"Profit margin is {sum(profits)/sum(revenue)}")
print(f'Total revenue is {sum(revenue)}')
print(f'Total cost is {sum(cost)}')

With customer segmentation:
Profit is 50930.078002717484
Profit margin is 0.2185888835412032
Total revenue is 232994.82195816856
Total cost is 182064.74395545103


labels = [1,2,3]


# Calculate profit and revenue based on grouped segmentations

def promo_ob3_v3(customers, cluster_num): 
    marketing_cost = 0.2
    l = np.array(customers[customers['cluster']==cluster_num]['UnitRevenue'].sort_values())
    return l[np.argmax(np.add.accumulate(l) - l*marketing_cost*len(customers[customers['cluster']==cluster_num]))]*marketing_cost


profits = []
revenue = []
cost = []
for i in labels:
    m_cost = 0.2
    revenue_segmented = sum(loyalty_df_metric['UnitRevenue'][(loyalty_df_metric['UnitRevenue']*m_cost <= promo_ob3_v3(loyalty_df_metric,i)) & (loyalty_df_metric['cluster']==i)])
    profit_segmented = revenue_segmented - promo_ob3_v3(loyalty_df_metric,i)*len(loyalty_df_metric[loyalty_df_metric['cluster']==i])
    cost.append(promo_ob3_v3(loyalty_df_metric,i)*len(loyalty_df_metric[loyalty_df_metric['cluster']==i]))
    revenue.append(revenue_segmented)
    profits.append(profit_segmented)

print(f"With customer segmentation:\nProfit is {sum(profits)}") 
print(f"Profit margin is {sum(profits)/sum(revenue)}")
print(f'Total revenue is {sum(revenue)}')
print(f'Total cost is {sum(cost)}')

With customer segmentation:
Profit is 185374.4748100365
Profit margin is 0.18146255851987536
Total revenue is 1021557.7049175833
Total cost is 836183.2301075468


test_df = pd.DataFrame({'method':['m1','m2','m3'],
                        'profit':[38505,50930,185147],
                        'revenue':[216637,232994,1021297]})


test_df['log'] = np.log(50*test_df['profit'])
test_df['margin'] = (test_df['profit']/test_df['revenue'])
weight = 2
test_df['margin_profit_score'] = (np.log10(test_df['profit'])**weight)*((test_df['profit'])/test_df['revenue'])
test_df


satisfaction_file = 'satisfaction.csv'
satisfaction_df = pd.read_csv(os.path.join(DATA_LOC, satisfaction_file))
satisfaction_df.columns

Index(['Unnamed: 0', 'id', 'Gender', 'Customer Type', 'Age', 'Type of Travel',
       'Class', 'Flight Distance', 'Inflight wifi service',
       'Departure/Arrival time convenient', 'Ease of Online booking',
       'Gate location', 'Food and drink', 'Online boarding', 'Seat comfort',
       'Inflight entertainment', 'On-board service', 'Leg room service',
       'Baggage handling', 'Checkin service', 'Inflight service',
       'Cleanliness', 'Departure Delay in Minutes', 'Arrival Delay in Minutes',
       'satisfaction'],
      dtype='object')


satisfaction_df['Customer Type']

0            Loyal Customer
1         disloyal Customer
2            Loyal Customer
3            Loyal Customer
4            Loyal Customer
                ...        
129875    disloyal Customer
129876       Loyal Customer
129877       Loyal Customer
129878       Loyal Customer
129879       Loyal Customer
Name: Customer Type, Length: 129880, dtype: object


# we only need the columns with demo characteristics and survey results in satisfaction_df, we can discard the churn
keep_columns = ['id', 'Gender', 'Customer Type', 'Age', 'Type of Travel',
       'Class', 'Flight Distance', 'Inflight wifi service',
       'Departure/Arrival time convenient', 'Ease of Online booking',
       'Gate location', 'Food and drink', 'Online boarding', 'Seat comfort',
       'Inflight entertainment', 'On-board service', 'Leg room service',
       'Baggage handling', 'Checkin service', 'Inflight service',
       'Cleanliness', 'Departure Delay in Minutes', 'Arrival Delay in Minutes','satisfaction']

satisfaction_df = satisfaction_df.loc[:,keep_columns]


# Create copy of loyalty_df_cluster (with cluster labels) for manipulation
loyalty_df_copy = loyalty_df_cluster.copy()


# Identify the Average distance per flight in loyalty df
loyalty_df_copy['AverageDistancePerFlight'] = loyalty_df_copy['TotalDistance']/loyalty_df_copy['TotalFlightsBooked']


# keep only average distances in satisfaction with a distance that is the same or lesser than the max average distance in loyalty_df
satisfaction_df_match = satisfaction_df[satisfaction_df['Flight Distance']<=max(loyalty_df_copy['AverageDistancePerFlight'])]


# Fill all missing values with 0
loyalty_df_copy['AverageDistancePerFlight'].fillna(0, inplace=True)


# Get the bins of average flight distance in the loyalty df
np.histogram(loyalty_df_copy['AverageDistancePerFlight'])

(array([14422,  2227,    45,    15,    12,     3,     2,     8,     1,
            2], dtype=int64),
 array([   0. ,  352.7,  705.4, 1058.1, 1410.8, 1763.5, 2116.2, 2468.9,
        2821.6, 3174.3, 3527. ]))


# Create a new column to put the average flight distance into different bins based on the bins identified in loyalty_df
bins = [ 0 ,  352.7,  705.4, 1058.1, 1410.8, 1763.5, 2116.2, 2468.9,
        2821.6, 3174.3, 3527. ]
loyalty_df_copy['AverageDistancePerFlight_bin'] = pd.cut(loyalty_df_copy['AverageDistancePerFlight'],bins=bins, include_lowest=True)
satisfaction_df_match['AverageDistancePerFlight_bin'] = pd.cut(satisfaction_df_match['Flight Distance'],bins=bins, include_lowest=True)


# Define the age group for the people in the loyalty program dataframe
# If someone is Married or Divorced, they must be at least 16 years of age
# If someone has finished their bachelor degree, it is reasonable to assume they are at least 16 years of age?

# Define a function to apply the conditions
def age_group(row):
    if row['MaritalStatus'] in ['Married', 'Divorced'] or row['Education'] in ['Bachelor', 'Master', 'Doctor']:
        return '>=16'
    else:
        return 'Unknown'

# Apply the function to create the AgeGroup column
loyalty_df_copy['AgeGroup'] = loyalty_df_copy.apply(age_group, axis=1)


# We can considered customers who have churned as 'neutral or dissatisfied' in the satisfaction file
def satisfaction_loyalty(row):
    if row['ChurnIndicator'] == 0:
        return 'churned'
    else:
        return 'not-churned'

def satisfaction_satisfaction(row):
    if row['Customer Type'] == 'Loyal Customer':
        return 'churned'
    else:
        return 'not-churned'    
    
# Apply the function to create the satisfaction column in loyalty_df
loyalty_df_copy['churn_type'] = loyalty_df_copy.apply(satisfaction_loyalty, axis=1)
satisfaction_df_match['churn_type'] = satisfaction_df_match.apply(satisfaction_satisfaction, axis=1)


# Define the age group for the people in the satisfaction survey
satisfaction_df_match['AgeGroup'] = satisfaction_df_match['Age'].apply(lambda x: '>=16' if x >= 16 else 'Unknown')


# This function will combine the gender, avg flight distance, satisfaction, and age group into one column, just for sampling later on
def age_flight_gender_func(row):
    return str(row['Gender'] + str(row['AverageDistancePerFlight_bin']) + row['churn_type'] + row['AgeGroup'])

loyalty_df_copy['gender_age_flight'] = loyalty_df_copy.apply(age_flight_gender_func,axis = 1)
satisfaction_df_match['gender_age_flight'] = satisfaction_df_match.apply(age_flight_gender_func,axis = 1)


# Extract the proportions of each type of Gender + Flight bin + satisfaction + Age from the loyalty card dataframe
age_flight_gender = loyalty_df_copy.groupby('gender_age_flight')['gender_age_flight'].value_counts().reset_index()
age_flight_gender = dict(zip(age_flight_gender['gender_age_flight'],age_flight_gender['count']))

# Sample from the satisfaction df
result_df = []
for i,u in age_flight_gender.items():
    subset_df = satisfaction_df_match[satisfaction_df_match['gender_age_flight']==i]
    sampled_df = subset_df.sample(n = u, replace=True, random_state=2024)
    result_df.append(sampled_df)

sampled_df = pd.concat(result_df)
sampled_df = sampled_df.drop(['id', 'Gender', 'Customer Type', 'Age', 'Type of Travel',
       'Class', 'Flight Distance', 'AverageDistancePerFlight_bin', 'churn_type',
       'AgeGroup'], axis=1)

# Merge the survey results from the sampled df with the loyalty df
sampled_df = sampled_df.sort_values('gender_age_flight').reset_index()
loyalty_df_copy = loyalty_df_copy.sort_values('gender_age_flight').reset_index()

loyalty_satisfaction_df = pd.concat([loyalty_df_copy,sampled_df], axis=1)


loyalty_satisfaction_df.columns

Index(['index', 'LoyaltyNumber', 'Country', 'Province', 'City', 'PostalCode',
       'Gender', 'Education', 'Salary', 'MaritalStatus', 'LoyaltyCard', 'CLV',
       'EnrollmentType', 'EnrollmentYear', 'EnrollmentMonth',
       'CancellationYear', 'CancellationMonth', 'EnrollmentDate',
       'EnrollmentYearMonth', 'CancellationDate', 'CancellationYearMonth',
       'ChurnIndicator', 'MonthsOfMembership', 'TotalFlightsBooked',
       'TotalFlightsWithCompanions', 'TotalTotalFlights', 'TotalDistance',
       'TotalPointsAccumulated', 'TotalPointsRedeemed',
       'TotalDollarCostPointsRedeemed', 'PercentagePointsRedeemed',
       'MonthlyTotalFlights', 'MonthlyPointsAccumulated',
       'ChurnIndicator_bool', 'cluster', 'AverageDistancePerFlight',
       'AverageDistancePerFlight_bin', 'AgeGroup', 'churn_type',
       'gender_age_flight', 'index', 'Inflight wifi service',
       'Departure/Arrival time convenient', 'Ease of Online booking',
       'Gate location', 'Food and drink', 'Online boarding', 'Seat comfort',
       'Inflight entertainment', 'On-board service', 'Leg room service',
       'Baggage handling', 'Checkin service', 'Inflight service',
       'Cleanliness', 'Departure Delay in Minutes', 'Arrival Delay in Minutes',
       'satisfaction', 'gender_age_flight'],
      dtype='object')


loyalty_satisfaction_lean = loyalty_satisfaction_df.drop(['index', 'gender_age_flight'], axis=1)
# loyalty_satisfaction_lean = loyalty_satisfaction_lean.drop_duplicates()
# loyalty_satisfaction_lean = loyalty_satisfaction_lean.iloc[:,:-1]


loyalty_satisfaction_lean


loyalty_satisfaction_lean.columns

Index(['LoyaltyNumber', 'Country', 'Province', 'City', 'PostalCode', 'Gender',
       'Education', 'Salary', 'MaritalStatus', 'LoyaltyCard', 'CLV',
       'EnrollmentType', 'EnrollmentYear', 'EnrollmentMonth',
       'CancellationYear', 'CancellationMonth', 'EnrollmentDate',
       'EnrollmentYearMonth', 'CancellationDate', 'CancellationYearMonth',
       'ChurnIndicator', 'MonthsOfMembership', 'TotalFlightsBooked',
       'TotalFlightsWithCompanions', 'TotalTotalFlights', 'TotalDistance',
       'TotalPointsAccumulated', 'TotalPointsRedeemed',
       'TotalDollarCostPointsRedeemed', 'PercentagePointsRedeemed',
       'MonthlyTotalFlights', 'MonthlyPointsAccumulated',
       'ChurnIndicator_bool', 'cluster', 'AverageDistancePerFlight',
       'AverageDistancePerFlight_bin', 'AgeGroup', 'churn_type',
       'Inflight wifi service', 'Departure/Arrival time convenient',
       'Ease of Online booking', 'Gate location', 'Food and drink',
       'Online boarding', 'Seat comfort', 'Inflight entertainment',
       'On-board service', 'Leg room service', 'Baggage handling',
       'Checkin service', 'Inflight service', 'Cleanliness',
       'Departure Delay in Minutes', 'Arrival Delay in Minutes',
       'satisfaction'],
      dtype='object')


survey_columns = ['Inflight wifi service', 'Departure/Arrival time convenient', 'Ease of Online booking',
       'Gate location', 'Food and drink', 'Online boarding', 
       'Seat comfort', 'Inflight entertainment', 'On-board service', 
       'Leg room service', 'Baggage handling', 'Checkin service', 
       'Inflight service', 'Cleanliness', 'satisfaction']


sns.set_style('whitegrid', {'axes.grid' : False})
sns.set_palette('colorblind')


# Graph by of distribution of answers for satisfaction
fig, axs = plt.subplots(nrows=3, ncols=5, figsize=(20, 10))

c = 0
for i in range(3):
    for j in range(5):
        if c < len(survey_columns):  # Ensure we don't go out of bounds

            tally_survey = loyalty_satisfaction_lean[survey_columns[c]].value_counts(normalize=False).sort_index()
            axs[i,j].bar(tally_survey.index, tally_survey.values, alpha = 0.5)
            axs[i,j].set_title(survey_columns[c])
            
            if survey_columns[c] == 'satisfaction':
                unique_values = sorted(loyalty_satisfaction_lean[survey_columns[c]].unique())
                axs[i, j].set_title('Overall Satisfaction')
                axs[i, j].set_xticks(unique_values)
                axs[i, j].set_xticklabels(unique_values)
            else:
                axs[i, j].set_xticks([0, 1, 2, 3, 4, 5])
                axs[i, j].set_xticklabels([0,1,2,3,4,5])
       
            c += 1
        else:
            axs[i, j].axis('off')

plt.tight_layout()
plt.show()


# Retain relevant columns for analysis
col = ['cluster', 'CLV', 'Inflight wifi service', 'Departure/Arrival time convenient', 'Ease of Online booking', 'Gate location', 'Food and drink', 'Online boarding',
      'Seat comfort', 'Inflight entertainment', 'On-board service', 'Leg room service', 'Baggage handling', 'Checkin service', 'Inflight service', 'Cleanliness']

loyalty_satisfaction_rf = loyalty_satisfaction_lean[col]
loyalty_satisfaction_rf


def corr_sig(df=None):
    p_matrix = np.zeros(shape=(df.shape[1],df.shape[1]))
    for col in df.columns:
        for col2 in df.drop(col,axis=1).columns:
            _ , p = stats.pearsonr(df[col],df[col2])
            p_matrix[df.columns.to_list().index(col),df.columns.to_list().index(col2)] = p
    return p_matrix

def plot_cor_matrix(corr, mask=None, annot_fontsize=9, label_fontsize=9):
    f, ax = plt.subplots(figsize=(11, 9))
    sns.heatmap(corr, ax=ax,
                mask=mask,
                # cosmetics
                annot=True, vmin=-1, vmax=1, center=0,
                cmap='PuBu_r', #linewidths=0.1, linecolor='grey', 
                cbar_kws={'orientation': 'vertical', 'shrink': 0.5, 'aspect': 30, 'pad': 0.02},
                annot_kws={"size": annot_fontsize})
    
    # Adjust axis labels font size
    ax.set_xticklabels(ax.get_xticklabels(), fontsize=label_fontsize)
    ax.set_yticklabels(ax.get_yticklabels(), fontsize=label_fontsize)


def corr_spear(df=None):
    p_matrix = np.zeros(shape=(df.shape[1],df.shape[1]))
    for col in df.columns:
        for col2 in df.drop(col,axis=1).columns:
            _ , p = stats.spearmanr(df[col],df[col2])
            p_matrix[df.columns.to_list().index(col),df.columns.to_list().index(col2)] = p
    return p_matrix


# Plotting without significance filtering
corr = loyalty_satisfaction_rf.corr()
mask = np.triu(corr)
plot_cor_matrix(corr,mask)
plt.show()


# Plotting with significance filter
corr = loyalty_satisfaction_rf.corr()                            # get correlation
p_values = corr_sig(loyalty_satisfaction_rf)                     # get p-Value
mask = np.invert(np.tril(p_values<0.05))    # mask - only get significant corr
plot_cor_matrix(corr,mask)


loyalty_satisfaction_rf.describe()


loyalty_satisfaction_clean = loyalty_satisfaction_lean.copy()

# Replace 0 with NaN for Survey Questions
survey_columns = ['Inflight wifi service',
       'Departure/Arrival time convenient', 'Ease of Online booking',
       'Gate location', 'Food and drink', 'Online boarding', 'Seat comfort',
       'Inflight entertainment', 'On-board service', 'Leg room service',
       'Baggage handling', 'Checkin service', 'Inflight service',
       'Cleanliness']
loyalty_satisfaction_clean[survey_columns] = loyalty_satisfaction_clean[survey_columns].replace(0, np.nan)

# Retain relevant columns
keep_columns = [ 
       'Inflight wifi service','Departure/Arrival time convenient', 'Ease of Online booking',
       'Gate location', 'Food and drink', 'Online boarding', 'Seat comfort',
       'Inflight entertainment', 'On-board service', 'Leg room service',
       'Baggage handling', 'Checkin service', 'Inflight service',
       'Cleanliness', 'satisfaction']
loyalty_satisfaction_clean = loyalty_satisfaction_clean[keep_columns]

# Drop rows where response = 0
loyalty_satisfaction_clean = loyalty_satisfaction_clean.dropna()

# Create Target Variable 'Overall Satisfaction'
loyalty_satisfaction_clean.loc[loyalty_satisfaction_clean['satisfaction'] == 'satisfied', 'Overall Satisfaction'] = 1
loyalty_satisfaction_clean.loc[loyalty_satisfaction_clean['satisfaction'] == 'neutral or dissatisfied', 'Overall Satisfaction'] = 0
loyalty_satisfaction_clean['Overall Satisfaction'] = loyalty_satisfaction_clean['Overall Satisfaction'].astype(int)
loyalty_satisfaction_clean = loyalty_satisfaction_clean.drop(columns=['satisfaction'])

loyalty_satisfaction_clean


# Checking for distribution of target variable
targetcount = loyalty_satisfaction_clean["Overall Satisfaction"].value_counts()
targetcount.plot(kind="bar", title="Class distribution of the target variable")

<Axes: title={'center': 'Class distribution of the target variable'}, xlabel='Overall Satisfaction'>


#Split Dataset for Training and Testing
trainDF_lr, testDF_lr = train_test_split(loyalty_satisfaction_clean, 
                                   test_size=0.2, 
                                   random_state=1234, 
                                   stratify=loyalty_satisfaction_clean[["Overall Satisfaction"]])

X_train_lr = trainDF_lr.iloc[:, trainDF_lr.columns != "Overall Satisfaction"]
y_train_lr = trainDF_lr.iloc[:, trainDF_lr.columns == "Overall Satisfaction"]
X_test_lr = testDF_lr.iloc[:, testDF_lr.columns != "Overall Satisfaction"]
y_test_lr = testDF_lr.iloc[:, testDF_lr.columns == "Overall Satisfaction"]


# Train a Logistic Regression model
lr_model = LogisticRegression(max_iter=1000)
lr_model.fit(X_train_lr, y_train_lr)

y_pred_lr = lr_model.predict(X_test_lr)
    
#  Model Evaluation
print(classification_report(y_test_lr, y_pred_lr))

              precision    recall  f1-score   support

           0       0.88      0.90      0.89      1942
           1       0.79      0.74      0.77       937

    accuracy                           0.85      2879
   macro avg       0.84      0.82      0.83      2879
weighted avg       0.85      0.85      0.85      2879


# Predict probabilities on the test set
baseline_probabilities = lr_model.predict_proba(X_test_lr)[:, 1]  # Predicted probabilities of satisfaction
baseline_avg_satisfaction = np.mean(baseline_probabilities)
print(f"Baseline Average Satisfaction: {baseline_avg_satisfaction:.4f}")

Baseline Average Satisfaction: 0.3268


# Get coefficients
lr_coefficients = lr_model.coef_[0]

# Create a DataFrame for better visualization
lr_coefficients_df = pd.DataFrame({
    'Feature': X_train_lr.columns,
    'Coefficient': lr_coefficients
})

# Sort the DataFrame by absolute value of coefficients
lr_coefficients_df = lr_coefficients_df.reindex(lr_coefficients_df.Coefficient.abs().sort_values(ascending=False).index)

# Plot the coefficients
colors = ['lightblue' if coef > 0 else 'pink' for coef in lr_coefficients_df['Coefficient']]

plt.figure(figsize=(12, 8))
bars = plt.barh(lr_coefficients_df['Feature'], lr_coefficients_df['Coefficient'], color=colors)
plt.xlabel('Coefficient')
plt.title('Feature Coefficients from Logistic Regression Model')
plt.gca().invert_yaxis()  # Invert y-axis to have the most significant feature on top

# Add values at the end of the bars
for bar in bars:
    width = bar.get_width()
    plt.text(width, bar.get_y() + bar.get_height() / 2, f'{width:.2f}', 
             va='center', ha='left', color='black')

plt.show()

# Display the coefficients DataFrame
lr_coefficients_df


# 1. Allocate Budget
# Define Equal Budget Allocation
num_features = 14
total_budget = 1000000

# Budget per feature
budget_per_feature_equal = total_budget / num_features


# 2. Apply improvements

# Ensure that the features used in the model are correctly referenced
features = X_train_lr.columns

# Define an assumed improvement value per unit budget (based on domain knowledge or assumptions)
unit_budget_to_improvement = 0.1  # Assume each $100,000 investment increases the feature by 0.1 unit

# Calculate the improvement factor based on this assumption
improvement_factor = unit_budget_to_improvement / 100_000

print(f"Derived Improvement Factor: {improvement_factor}")

# Create budget allocation dictionary for equal allocation
budget_allocation_equal = {feature: budget_per_feature_equal for feature in X_train_lr.columns}

# Function to apply improvements based on budget allocation
def apply_improvements(X, budget_allocation, improvement_factor):
    X_improved = X.copy()
    
    for feature, budget in budget_allocation.items():
        improvement_value = budget * improvement_factor
        if feature in X_improved.columns:
            X_improved[feature] = np.clip(X_improved[feature] + improvement_value, 1, 5)  # Assuming scale of 1 to 5
    
    return X_improved

# Apply improvements to the test set
X_test_improved_equal = apply_improvements(X_test_lr, budget_allocation_equal, improvement_factor)

Derived Improvement Factor: 1e-06


# 3. Measure Impact
# Calculate new predicted probabilities after feature improvements
predicted_probabilities_improved_equal = lr_model.predict_proba(X_test_improved_equal)[:, 1]

# Calculate new average satisfaction
new_avg_satisfaction_equal = np.mean(predicted_probabilities_improved_equal)
print(f"New Average Satisfaction (Equal Allocation): {new_avg_satisfaction_equal:.4f}")

# Measure the impact
impact_equal = new_avg_satisfaction_equal - baseline_avg_satisfaction
print(f"Impact of Equal Allocation on Satisfaction: {impact_equal:.4f}")

New Average Satisfaction (Equal Allocation): 0.3443
Impact of Equal Allocation on Satisfaction: 0.0175


# Remove 'Overall Satisfaction' column
loyalty_satisfaction_clean_filtered = loyalty_satisfaction_clean.drop(columns=['Overall Satisfaction'])

# Calculate the mean score for each feature excluding 'Overall Satisfaction'
mean_scores = loyalty_satisfaction_clean_filtered.mean()

# Create a DataFrame for better visualization
mean_scores_df = pd.DataFrame(mean_scores, columns=['Mean Score'])

# Sort the DataFrame by mean scores in ascending order
mean_scores_df = mean_scores_df.sort_values(by='Mean Score')

mean_scores_df.style.set_caption("Features and their Mean Scores").background_gradient(cmap='PiYG')


# 1. Allocate Budget
# Rank features based on mean scores (ascending order)
ranked_features = mean_scores.sort_values()

# Select the lowest 3 performing features (for example)
lowest_features = ranked_features.index[:3]

# Budget per lowest performing feature
budget_per_feature_lowest = total_budget / len(lowest_features)


# 2. Apply Improvements
# Create budget allocation dictionary for lowest performing features
budget_allocation_lowest = {feature: budget_per_feature_lowest for feature in lowest_features}

# Apply improvements to the test set
X_test_improved_lowest = apply_improvements(X_test_lr, budget_allocation_lowest, improvement_factor)


# 3. Measure Impact
# Calculate new predicted probabilities after feature improvements
predicted_probabilities_improved_lowest = lr_model.predict_proba(X_test_improved_lowest)[:, 1]

# Calculate new average satisfaction
new_avg_satisfaction_lowest = np.mean(predicted_probabilities_improved_lowest)
print(f"New Average Satisfaction (Lowest Performing Features): {new_avg_satisfaction_lowest:.4f}")

# Measure the impact
impact_lowest = new_avg_satisfaction_lowest - baseline_avg_satisfaction
print(f"Impact of Targeting Lowest Performing Features on Satisfaction: {impact_lowest:.4f}")

New Average Satisfaction (Lowest Performing Features): 0.3625
Impact of Targeting Lowest Performing Features on Satisfaction: 0.0357


# 1. Calculate Feature Importance

# Calculate the absolute values of the coefficients
lr_coefficients_df['Abs_Coefficient'] = lr_coefficients_df['Coefficient'].abs()

# Sort the DataFrame by absolute value of coefficients in descending order
# Separate positive and negative coefficients
positive_coef_df = lr_coefficients_df[lr_coefficients_df['Coefficient'] > 0].sort_values(by='Coefficient', ascending=False)
negative_coef_df = lr_coefficients_df[lr_coefficients_df['Coefficient'] < 0].sort_values(by='Coefficient')
print("Feature importance based on logistic regression coefficients:")
print(positive_coef_df)

Feature importance based on logistic regression coefficients:
                   Feature  Coefficient  Abs_Coefficient
0    Inflight wifi service     0.877880         0.877880
5          Online boarding     0.786087         0.786087
7   Inflight entertainment     0.741962         0.741962
9         Leg room service     0.417182         0.417182
8         On-board service     0.334594         0.334594
11         Checkin service     0.220951         0.220951
2   Ease of Online booking     0.137866         0.137866
3            Gate location     0.058092         0.058092
6             Seat comfort     0.054074         0.054074


# 2. Allocate Budget

# Calculate total sum of absolute positive coefficients
total_positive_abs_coef = positive_coef_df['Coefficient'].sum()

# Calculate budget allocation for each positive feature based on its importance
positive_coef_df['Budget_Allocation'] = (positive_coef_df['Coefficient'] / total_positive_abs_coef) * total_budget

# Display the budget allocation
print("Budget Allocation based on Positive Feature Importance:")
print(positive_coef_df[['Feature', 'Budget_Allocation']])

Budget Allocation based on Positive Feature Importance:
                   Feature  Budget_Allocation
0    Inflight wifi service      241927.682828
5          Online boarding      216631.319705
7   Inflight entertainment      204471.181712
9         Leg room service      114967.691932
8         On-board service       92207.988369
11         Checkin service       60889.933058
2   Ease of Online booking       37993.469594
3            Gate location       16008.961634
6             Seat comfort       14901.771168


# 3. Apply Improvements
# Create budget allocation dictionary for positive features
budget_allocation_positive = positive_coef_df.set_index('Feature')['Budget_Allocation'].to_dict()

# Apply improvements to the test set
X_test_improved_positive = apply_improvements(X_test_lr, budget_allocation_positive, improvement_factor)


# 4. Measure Impact
# Calculate new predicted probabilities after feature improvements
predicted_probabilities_improved_positive = lr_model.predict_proba(X_test_improved_positive)[:, 1]

# Calculate new average satisfaction
new_avg_satisfaction_positive = np.mean(predicted_probabilities_improved_positive)
print(f"New Average Satisfaction (Positive Features): {new_avg_satisfaction_positive:.4f}")

# Measure the impact
impact_positive = new_avg_satisfaction_positive - baseline_avg_satisfaction
print(f"Impact of Targeting Positive Features on Satisfaction: {impact_positive:.4f}")

New Average Satisfaction (Positive Features): 0.3868
Impact of Targeting Positive Features on Satisfaction: 0.0600

	Salary	CLV	EnrollmentYear	EnrollmentMonth	CancellationYear	CancellationMonth
count	12499.000000	16737.000000	16737.000000	16737.000000	2067.0	2067.0
mean	79245.609409	7988.896536	2015.253211	6.669116	2016.503145	6.962748
std	35008.297285	6860.982280	1.979111	3.398958	1.380743	3.455297
min	-58486.000000	1898.010000	2012.000000	1.000000	2013.0	1.0
25%	59246.500000	3980.840000	2014.000000	4.000000	2016.0	4.0
50%	73455.000000	5780.180000	2015.000000	7.000000	2017.0	7.0
75%	88517.500000	8940.580000	2017.000000	10.000000	2018.0	10.0
max	407228.000000	83325.380000	2018.000000	12.000000	2018.0	12.0

	LoyaltyNumber	Country	Province	City	PostalCode	Gender	Education	Salary	MaritalStatus	LoyaltyCard	CLV	EnrollmentType	EnrollmentYear	EnrollmentMonth	CancellationYear	CancellationMonth
1082	542976	Canada	Quebec	Montreal	H2Y 4R4	Male	High School or Below	-49830.0	Divorced	Star	24127.50	2018 Promotion	2018	3	<NA>	<NA>
1894	959977	Canada	British Columbia	Vancouver	V5R 1W3	Female	Bachelor	-12497.0	Married	Aurora	9453.00	2018 Promotion	2018	3	<NA>	<NA>
2471	232755	Canada	British Columbia	Vancouver	V1E 4R6	Female	Bachelor	-46683.0	Single	Nova	4787.81	2018 Promotion	2018	3	<NA>	<NA>
3575	525245	Canada	British Columbia	Victoria	V10 6T5	Male	Bachelor	-45962.0	Married	Star	2402.33	2018 Promotion	2018	3	<NA>	<NA>
3932	603070	Canada	British Columbia	West Vancouver	V6V 8Z3	Female	Bachelor	-19325.0	Single	Star	2893.74	2018 Promotion	2018	3	<NA>	<NA>
4712	491242	Canada	British Columbia	Dawson Creek	U5I 4F1	Male	Bachelor	-43234.0	Married	Star	7597.91	2018 Promotion	2018	3	<NA>	<NA>
6560	115505	Canada	Newfoundland	St. John's	A1C 6H9	Male	Bachelor	-10605.0	Married	Nova	5860.17	2018 Promotion	2018	4	<NA>	<NA>
6570	430398	Canada	Newfoundland	St. John's	A1C 6H9	Male	Bachelor	-17534.0	Married	Nova	49423.80	2018 Promotion	2018	3	<NA>	<NA>
7373	152016	Canada	Ontario	Toronto	P1J 8T7	Female	Bachelor	-58486.0	Married	Aurora	5067.21	2018 Promotion	2018	2	<NA>	<NA>
8576	194065	Canada	Ontario	Sudbury	M5V 1G5	Female	Bachelor	-31911.0	Married	Nova	2888.85	2018 Promotion	2018	2	<NA>	<NA>
8767	212128	Canada	Ontario	Toronto	P2T 6G3	Male	Bachelor	-49001.0	Married	Nova	3130.68	2018 Promotion	2018	2	<NA>	<NA>
10232	790475	Canada	Ontario	Trenton	K8V 4B2	Female	Bachelor	-34079.0	Married	Nova	12913.50	2018 Promotion	2018	2	<NA>	<NA>
11635	366599	Canada	Ontario	Toronto	M1R 4K3	Female	Bachelor	-9081.0	Married	Star	6915.73	2018 Promotion	2018	4	<NA>	<NA>
12596	436087	Canada	Quebec	Montreal	H2T 9K8	Male	Bachelor	-46470.0	Married	Star	4786.89	2018 Promotion	2018	4	2018	8
13564	364596	Canada	Quebec	Tremblant	H5Y 2S9	Female	Bachelor	-26322.0	Single	Aurora	16710.84	2018 Promotion	2018	4	2018	12
14327	239955	Canada	Quebec	Hull	J8Y 3Z5	Female	Bachelor	-47310.0	Married	Nova	6366.23	2018 Promotion	2018	3	<NA>	<NA>
14355	347013	Canada	Quebec	Quebec City	G1B 3L5	Male	Bachelor	-39503.0	Married	Nova	6446.71	2018 Promotion	2018	3	<NA>	<NA>
15416	729561	Canada	Quebec	Quebec City	G1B 3L5	Female	Bachelor	-19332.0	Divorced	Star	5308.29	2018 Promotion	2018	2	<NA>	<NA>
16431	734647	Canada	Saskatchewan	Regina	S1J 3C5	Male	Bachelor	-46303.0	Married	Nova	11280.73	2018 Promotion	2018	4	<NA>	<NA>
16735	906428	Canada	Yukon	Whitehorse	Y2K 6R0	Male	Bachelor	-57297.0	Married	Star	10018.66	2018 Promotion	2018	4	<NA>	<NA>

	Year	Month	FlightsBooked	FlightsWithCompanions	TotalFlights	Distance	PointsAccumulated	PointsRedeemed	DollarCostPointsRedeemed
count	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000
mean	2017.500000	6.500000	4.115052	1.031805	5.146858	1208.880059	123.692721	30.696872	2.484503
std	0.500001	3.452057	5.225518	2.076869	6.521227	1433.155320	146.599831	125.486049	10.150038
min	2017.000000	1.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	2017.000000	3.750000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
50%	2017.500000	6.500000	1.000000	0.000000	1.000000	488.000000	50.000000	0.000000	0.000000
75%	2018.000000	9.250000	8.000000	1.000000	10.000000	2336.000000	239.000000	0.000000	0.000000
max	2018.000000	12.000000	21.000000	11.000000	32.000000	6293.000000	676.500000	876.000000	71.000000

	Salary	CLV	EnrollmentYear	EnrollmentMonth	CancellationYear	CancellationMonth	EnrollmentDate	CancellationDate	ChurnIndicator
count	16737.000000	16737.000000	16737.000000	16737.000000	2067.0	2067.0	16737	2067	16737.000000
mean	59222.180618	7988.896536	2015.253211	6.669116	2016.503145	6.962748	2015-09-22 11:45:55.977773824	2016-12-30 23:07:03.222060800	0.123499
min	0.000000	1898.010000	2012.000000	1.000000	2013.0	1.0	2012-04-01 00:00:00	2013-01-01 00:00:00	0.000000
25%	0.000000	3980.840000	2014.000000	4.000000	2016.0	4.0	2014-01-01 00:00:00	2016-01-01 00:00:00	0.000000
50%	63654.000000	5780.180000	2015.000000	7.000000	2017.0	7.0	2015-11-01 00:00:00	2017-04-01 00:00:00	0.000000
75%	82940.000000	8940.580000	2017.000000	10.000000	2018.0	10.0	2017-07-01 00:00:00	2018-03-01 00:00:00	0.000000
max	407228.000000	83325.380000	2018.000000	12.000000	2018.0	12.0	2018-12-01 00:00:00	2018-12-01 00:00:00	1.000000
std	45781.737001	6860.982280	1.979111	3.398958	1.380743	3.455297	NaN	NaN	0.329019

	Year	Month	FlightsBooked	FlightsWithCompanions	TotalFlights	Distance	PointsAccumulated	PointsRedeemed	DollarCostPointsRedeemed	Date
count	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624.000000	405624
mean	2017.500000	6.500000	4.115052	1.031805	5.146858	1208.880059	123.692721	30.696872	2.484503	2017-12-15 23:59:59.999999744
min	2017.000000	1.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	2017-01-01 00:00:00
25%	2017.000000	3.750000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	2017-06-23 12:00:00
50%	2017.500000	6.500000	1.000000	0.000000	1.000000	488.000000	50.000000	0.000000	0.000000	2017-12-16 12:00:00
75%	2018.000000	9.250000	8.000000	1.000000	10.000000	2336.000000	239.000000	0.000000	0.000000	2018-06-08 12:00:00
max	2018.000000	12.000000	21.000000	11.000000	32.000000	6293.000000	676.500000	876.000000	71.000000	2018-12-01 00:00:00
std	0.500001	3.452057	5.225518	2.076869	6.521227	1433.155320	146.599831	125.486049	10.150038	NaN

	LoyaltyNumber	Country	Province	City	PostalCode	Gender	Education	Salary	MaritalStatus	LoyaltyCard	CLV	EnrollmentType	EnrollmentYear	EnrollmentMonth	CancellationYear	CancellationMonth
0	480934	Canada	Ontario	Toronto	M2Z 4K1	Female	Bachelor	83236.0	Married	Star	3839.14	Standard	2016	2	NaN	NaN
1	549612	Canada	Alberta	Edmonton	T3G 6Y6	Male	College	NaN	Divorced	Star	3839.61	Standard	2016	3	NaN	NaN
2	429460	Canada	British Columbia	Vancouver	V6E 3D9	Male	College	NaN	Single	Star	3839.75	Standard	2014	7	2018.0	1.0
3	608370	Canada	Ontario	Toronto	P1W 1K4	Male	College	NaN	Single	Star	3839.75	Standard	2013	2	NaN	NaN
4	530508	Canada	Quebec	Hull	J8Y 3Z5	Male	Bachelor	103495.0	Married	Star	3842.79	Standard	2014	10	NaN	NaN

	LoyaltyNumber	Year	Month	FlightsBooked	FlightsWithCompanions	TotalFlights	Distance	PointsAccumulated
0	100018	2017	1	3	0	3	1521	152.0
1	100102	2017	1	10	4	14	2030	203.0
2	100140	2017	1	6	0	6	1200	120.0
3	100214	2017	1	0	0	0	0	0.0
4	100272	2017	1	0	0	0	0	0.0

	TotalFlightsBooked	TotalFlightsWithCompanions	TotalTotalFlights	TotalDistance	TotalPointsAccumulated	TotalPointsRedeemed	TotalDollarCostPointsRedeemed
LoyaltyNumber
100018	157	35	192	50682	5376.00	1513	123
100102	173	42	215	40222	4115.25	1195	96
100140	152	38	190	41252	4184.25	593	48
100214	79	17	96	33982	3426.00	861	70
100272	127	36	163	40872	4108.04	1007	82
...	...	...	...	...	...	...	...
998972	80	24	104	31676	3164.00	2201	178
999304	46	5	51	10662	1065.00	809	65
999731	57	16	73	15598	1555.00	407	33
999788	42	16	58	15784	1575.00	361	29
999891	1	1	2	2526	252.00	414	34

	LoyaltyNumber	Date	Last1Y_FlightsBooked	Last1Y_FlightsWithCompanions	Last1Y_TotalFlights	Last1Y_Distance	Last1Y_PointsAccumulated	Last1Y_PointsRedeemed	Last1Y_DollarCostPointsRedeemed	Last6M_FlightsBooked	...	Last3M_PointsAccumulated	Last3M_PointsRedeemed	Last3M_DollarCostPointsRedeemed	Last1M_FlightsBooked	Last1M_FlightsWithCompanions	Last1M_TotalFlights	Last1M_Distance	Last1M_PointsAccumulated	Last1M_PointsRedeemed	Last1M_DollarCostPointsRedeemed
0	100018	2017-01-01	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
1	100102	2017-01-01	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	100140	2017-01-01	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
3	100214	2017-01-01	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
4	100272	2017-01-01	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
528019	999902	2018-12-01	87.0	15.0	102.0	28122.0	3149.5	384.0	31.0	37.0	...	120.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
528020	999911	2018-12-01	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
528021	999940	2018-12-01	51.0	18.0	69.0	15929.0	1668.0	672.0	54.0	33.0	...	615.0	0.0	0.0	12.0	5.0	17.0	3859.0	385.0	0.0	0.0
528022	999982	2018-12-01	22.0	2.0	24.0	5948.0	594.0	0.0	0.0	22.0	...	535.0	0.0	0.0	3.0	0.0	3.0	1560.0	156.0	0.0	0.0
528023	999986	2018-12-01	49.0	10.0	59.0	13665.0	1495.0	1199.0	97.0	40.0	...	274.0	0.0	0.0	11.0	4.0	15.0	2040.0	204.0	0.0	0.0

	LoyaltyNumber	Date	EnrollmentDate	CancellationDate	MonthsOfMembership
0	100018	2017-01-01	2016-08-01	NaT	5
1	100102	2017-01-01	2013-03-01	NaT	46
2	100140	2017-01-01	2016-07-01	NaT	6
3	100214	2017-01-01	2015-08-01	NaT	17
4	100272	2017-01-01	2014-01-01	NaT	36
...	...	...	...	...	...
405619	999902	2018-12-01	2014-05-01	NaT	55
405620	999911	2018-12-01	2012-08-01	NaT	76
405621	999940	2018-12-01	2017-07-01	NaT	17
405622	999982	2018-12-01	2018-07-01	NaT	5
405623	999986	2018-12-01	2018-02-01	NaT	10

	EnrollmentYear	count	growth
0	2012	1686	0.000000
1	2013	2397	42.170819
2	2014	2370	-1.126408
3	2015	2331	-1.645570
4	2016	2456	5.362505
5	2017	2487	1.262215
6	2018	3010	21.029353

	LoyaltyNumber	Date	Last1Y_FlightsBooked	Last1Y_FlightsWithCompanions	Last1Y_TotalFlights	Last1Y_Distance	Last1Y_PointsAccumulated	Last1Y_PointsRedeemed	Last1Y_DollarCostPointsRedeemed	Last6M_FlightsBooked	...	Last1M_FlightsWithCompanions	Last1M_TotalFlights	Last1M_Distance	Last1M_PointsAccumulated	Last1M_PointsRedeemed	Last1M_DollarCostPointsRedeemed	EnrollmentDate	CancellationDate	MonthsOfMembership	ChurnIndicator
332052	100018	2018-01-01	81.0	16.0	97.0	26663.0	2664.00	1128.0	92.0	44.0	...	0.0	6.0	1908.0	190.0	0.0	0.0	2016-08-01	NaT	17	0
332053	100102	2018-01-01	83.0	18.0	101.0	16933.0	1691.00	1195.0	96.0	49.0	...	3.0	14.0	3542.0	354.0	510.0	41.0	2013-03-01	NaT	58	0
332054	100140	2018-01-01	72.0	18.0	90.0	24681.0	2466.00	0.0	0.0	41.0	...	0.0	15.0	4500.0	450.0	0.0	0.0	2016-07-01	NaT	18	0
332055	100214	2018-01-01	37.0	3.0	40.0	15259.0	1523.00	861.0	70.0	30.0	...	0.0	11.0	4334.0	433.0	0.0	0.0	2015-08-01	NaT	29	0
332056	100272	2018-01-01	66.0	19.0	85.0	22711.0	2268.00	393.0	32.0	15.0	...	6.0	17.0	4539.0	453.0	0.0	0.0	2014-01-01	NaT	48	0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
664094	999631	2018-12-01	44.0	19.0	63.0	17336.0	1756.92	0.0	0.0	29.0	...	5.0	10.0	2160.0	216.0	0.0	0.0	2014-01-01	NaT	59	0
664096	999758	2018-12-01	15.0	4.0	19.0	7822.0	782.00	573.0	46.0	15.0	...	0.0	10.0	3570.0	357.0	0.0	0.0	2018-08-01	NaT	4	0
664099	999902	2018-12-01	87.0	15.0	102.0	28122.0	3149.50	384.0	31.0	37.0	...	0.0	0.0	0.0	0.0	0.0	0.0	2014-05-01	NaT	55	0
664101	999940	2018-12-01	51.0	18.0	69.0	15929.0	1668.00	672.0	54.0	33.0	...	5.0	17.0	3859.0	385.0	0.0	0.0	2017-07-01	NaT	17	0
664102	999982	2018-12-01	22.0	2.0	24.0	5948.0	594.00	0.0	0.0	22.0	...	0.0	3.0	1560.0	156.0	0.0	0.0	2018-07-01	NaT	5	0

	ChurnIndicator	Salary	MonthlyTotalFlights	PercentagePointsRedeemed
count	2075.000000	2075.000000	2075.000000	2075.000000
mean	0.254458	58952.838072	11.597942	0.276236
std	0.435661	44711.063915	23.614560	0.544584
min	0.000000	0.000000	0.000000	0.000000
25%	0.000000	0.000000	4.670000	0.000000
50%	0.000000	64381.000000	7.200000	0.000000
75%	1.000000	82338.000000	10.330000	0.390000
max	1.000000	298851.000000	223.000000	5.200000

	ChurnIndicator	Salary	MonthlyTotalFlights	PercentagePointsRedeemed
count	528.0	528.000000	528.000000	528.000000
mean	1.0	58272.257576	6.142140	0.281761
std	0.0	46077.435863	2.988565	0.501163
min	1.0	0.000000	0.120000	0.000000
25%	1.0	0.000000	4.120000	0.000000
50%	1.0	64162.500000	6.120000	0.000000
75%	1.0	81685.500000	8.152500	0.422500
max	1.0	270797.000000	17.000000	5.200000

	support	itemsets	length
27	0.143939	(MonthlyFlightsBins_(0.0, 4.0], Education_Bach...	2
43	0.136364	(MonthlyFlightsBins_(0.0, 4.0], MaritalStatus_...	2
44	0.145833	(MonthlyFlightsBins_(4.0, 6.0], MaritalStatus_...	2
47	0.104167	(PctPointsRedeemedBins_(0.0, 0.4], MaritalStat...	2
50	0.100379	(LoyaltyCard_Nova, MaritalStatus_Single)	2
53	0.102273	(MaritalStatus_Single, PctPointsRedeemedBins_(...	2
56	0.111742	(LoyaltyCard_Nova, SalaryBins_(80000.0, inf])	2
61	0.106061	(MonthlyFlightsBins_(0.0, 4.0], LoyaltyCard_Star)	2
62	0.123106	(LoyaltyCard_Star, MonthlyFlightsBins_(6.0, 8.0])	2
74	0.106061	(LoyaltyCard_Aurora, Education_Bachelor, Marit...	3
80	0.102273	(MonthlyFlightsBins_(4.0, 6.0], Education_Bach...	3
87	0.117424	(MonthlyFlightsBins_(0.0, 4.0], Education_Bach...	3
94	0.104167	(MaritalStatus_Married, MonthlyFlightsBins_(0....	3

	support	itemsets	length
4	0.215909	(LoyaltyCard_Aurora)	1
10	0.240530	(MonthlyFlightsBins_(0.0, 4.0])	1
11	0.255682	(MonthlyFlightsBins_(4.0, 6.0])	1
12	0.240530	(MonthlyFlightsBins_(6.0, 8.0])	1

	all	lifelong	dualphase
predicted_churn	14303.0	36.0	643.0
true_positive	1015.0	10.0	180.0
false_negative	0.0	1005.0	835.0
churn_rate_after_promo	0.0	7.0	5.8
cost_efficiency	14.1	3.6	3.6

	Salary	CLV	MonthsOfMembership	TotalFlightsBooked	TotalFlightsWithCompanions	TotalTotalFlights	TotalDistance	TotalPointsAccumulated	TotalPointsRedeemed	TotalDollarCostPointsRedeemed
PC1	0.001779	-0.002746	0.256215	0.391827	0.370512	0.399315	0.395435	0.394964	0.292972	0.293061
PC2	-0.040422	0.018458	-0.217314	-0.201056	-0.061576	-0.175131	-0.183925	-0.183603	0.635855	0.635660
PC3	0.713647	0.699477	0.027271	-0.008601	-0.012483	-0.009772	-0.005293	-0.000071	0.013082	0.013134

	Salary	CLV	MonthsOfMembership	TotalFlightsBooked	TotalFlightsWithCompanions	TotalTotalFlights	TotalDistance	TotalPointsAccumulated	TotalPointsRedeemed	TotalDollarCostPointsRedeemed
cluster
1	63985.000000	5768.460000	43.000000	122.000000	28.000000	151.000000	35968.000000	3674.920000	523.000000	42.000000
2	63863.000000	5780.180000	47.000000	137.000000	36.000000	174.000000	40788.000000	4168.250000	1552.000000	126.000000
3	62866.000000	5796.580000	8.000000	21.000000	4.000000	27.000000	6152.000000	616.500000	0.000000	0.000000

	method	profit	revenue	log	margin	margin_profit_score
0	m1	38505	216637	14.470566	0.177740	3.737327
1	m2	50930	232994	14.750230	0.218589	4.842978
2	m3	185147	1021297	16.040928	0.181286	5.030098

	LoyaltyNumber	Country	Province	City	PostalCode	Gender	Education	Salary	MaritalStatus	LoyaltyCard	...	Inflight entertainment	On-board service	Leg room service	Baggage handling	Checkin service	Inflight service	Cleanliness	Departure Delay in Minutes	Arrival Delay in Minutes	satisfaction
0	811333	Canada	Ontario	Toronto	M1R 4K3	Female	Bachelor	82937.0	Married	Aurora	...	5	3	4	5	5	5	5	0	20.0	neutral or dissatisfied
1	882874	Canada	Ontario	Toronto	M2M 7K8	Female	Bachelor	53559.0	Married	Star	...	4	4	5	1	1	3	4	5	4.0	satisfied
2	504525	Canada	Quebec	Montreal	H2Y 4R4	Female	Bachelor	52808.0	Married	Nova	...	5	5	4	5	5	5	3	0	10.0	satisfied
3	574330	Canada	Quebec	Hull	J8Y 3Z5	Female	Bachelor	64968.0	Married	Star	...	5	5	5	5	1	5	4	30	19.0	satisfied
4	607405	Canada	Quebec	Tremblant	H5Y 2S9	Female	Bachelor	72820.0	Married	Nova	...	2	2	1	3	4	3	2	0	0.0	neutral or dissatisfied
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
16732	430221	Canada	Manitoba	Winnipeg	R3R 3T4	Male	Bachelor	90071.0	Divorced	Star	...	5	3	2	5	1	4	5	0	0.0	satisfied
16733	867284	Canada	Nova Scotia	Halifax	B3C 2M8	Male	College	0.0	Single	Nova	...	1	3	4	5	2	4	1	0	0.0	neutral or dissatisfied
16734	341904	Canada	Manitoba	Winnipeg	R2C 0M5	Male	College	0.0	Single	Star	...	2	1	3	4	1	3	2	0	0.0	neutral or dissatisfied
16735	124377	Canada	British Columbia	Whistler	V6T 1Y8	Male	College	0.0	Single	Nova	...	5	4	2	2	2	2	5	0	1.0	neutral or dissatisfied
16736	971608	Canada	Manitoba	Winnipeg	R2C 0M5	Male	Bachelor	73945.0	Single	Star	...	1	5	5	4	3	4	1	0	0.0	neutral or dissatisfied

	cluster	CLV	Inflight wifi service	Departure/Arrival time convenient	Ease of Online booking	Gate location	Food and drink	Online boarding	Seat comfort	Inflight entertainment	On-board service	Leg room service	Baggage handling	Checkin service	Inflight service	Cleanliness
0	2	20946.19	4	0	4	3	5	4	5	5	3	4	5	5	5	5
1	3	2599.31	1	1	1	1	4	4	4	4	4	5	1	1	3	4
2	1	2999.34	4	4	4	4	3	1	4	5	5	4	5	5	5	3
3	1	2599.97	5	4	4	4	1	5	4	5	5	5	5	1	5	4
4	1	2994.35	2	3	2	3	2	2	3	2	2	1	3	4	3	2
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
16732	3	4626.80	5	1	1	1	5	4	5	5	3	2	5	1	4	5
16733	3	13144.49	2	4	2	1	1	2	1	1	3	4	5	2	4	1
16734	3	2254.40	2	1	2	3	2	2	4	2	1	3	4	1	3	2
16735	3	8175.32	3	2	3	3	5	3	5	5	4	2	2	2	2	5
16736	3	12157.33	4	4	4	3	1	4	1	1	5	5	4	3	4	1

	Feature	Coefficient
0	Inflight wifi service	0.877880
5	Online boarding	0.786087
7	Inflight entertainment	0.741962
1	Departure/Arrival time convenient	-0.542590
9	Leg room service	0.417182
8	On-board service	0.334594
11	Checkin service	0.220951
4	Food and drink	-0.219847
2	Ease of Online booking	0.137866
12	Inflight service	-0.088342
13	Cleanliness	-0.086174
3	Gate location	0.058092
6	Seat comfort	0.054074
10	Baggage handling	-0.035825

	Mean Score
Inflight wifi service	2.828135
Ease of Online booking	2.838694
Gate location	2.977075
Online boarding	3.156999
Checkin service	3.176659
Food and drink	3.181938
Cleanliness	3.209726
Leg room service	3.265787
On-board service	3.272873
Inflight entertainment	3.306982
Departure/Arrival time convenient	3.316568
Seat comfort	3.325113
Baggage handling	3.537478
Inflight service	3.554498

Reaching New Heights: Develop and Retain Valuable Airline Customers

1 Getting started¶

1.1 Import libraries¶

1.2 Define helper functions & constants¶

1.3 Import data¶

1.3.1 Loyalty history (loyalty_df)¶

1.3.2 Flight activity (flight_df)¶

2 Data Preparation¶

2.1 Preliminary Data Cleaning¶

2.1.1 Loyalty history (loyalty_df)¶

2.1.1.1 Correcting data type¶

2.1.1.2 Checking numeric columns¶

2.1.1.3 Time column¶

2.1.1.4 Churn indicator¶

2.1.2 Flight activity (flight_df)¶

2.1.2.1 Converting data type¶

2.1.2.2 Checking numeric column¶

2.1.2.3 Time column¶

2.1.3 Data cleansing result¶

2.2 Feature engineering¶

2.2.1 Customer level¶

Merge all feature sets to loyalty_df¶

2.2.2 Flight level¶

3 Exploratory Data Analysis¶

3.1 Customer Distribution¶

3.1.1 Salary Distribution¶

3.1.2 Gender Distribution¶

3.1.3 Education Distribution¶

3.1.4 Marital Status Distribution¶

3.1.5 Loyalty Card Type Distribution¶

3.1.6 City Distribution¶

3.1.7 Enrollment Type Distribution¶

3.1.8 CLV Distribution¶

3.1.8.1 Kruskal-Wallis Test for Categorical Features and CLV¶

3.1.9 Relationship between Points Accumulated and Points Redeemed¶

3.2 Potential impact of 2018 Promotion¶

3.2.1 Number of Enrollments and Cancellations¶

3.2.2 Number of Enrollments with/without 2018 Promotions¶

3.2.3 Additional Exploration on 2018 Promotion¶

3.3 Churn customers¶

3.3.1 Additional Exploration on Churn¶

4 Maintaining Loyalty with Perks¶

4.1 Filter active customers in 2017-2018 with Standard Enrollment¶

4.2 Lifelong Loyalty Scoring System¶

4.2.1 Pre-process data on flight level¶

4.2.2 Split data¶

4.2.3 Train model¶

4.2.4 Customer churn prediction¶

4.3 Dual-phase retention framework¶

4.3.1 Stage 0 - Exploratory data analysis¶

4.3.2 Stage 1 - Frequent pattern mining¶

4.3.3 Stage 2 - Predictive modelling on customer retention¶

4.3.3.1 Pre-process data on flight level¶

4.3.3.2 Split data¶

4.3.3.3 Train model¶

4.3.4 Customer churn prediction¶

4.3.4.1 Stage 1¶

4.3.4.2 Stage 2¶

4.3.4.3 Both stages¶

4.4 Customer retention statistics¶

5 Encouraging Additional Expenditure with Coupon & Customer Segmentation¶

5.1 Defining the number of clusters¶

5.2 K-means Clustering¶

5.2.1 Distribution of CLV Cluster by Gender¶

5.2.2 Distribution of CLV Cluster by Salary¶

5.2.3 Distribution of CLV Cluster by Education¶

5.2.4 Distribution of CLV Cluster by Loyalty Card¶

5.2.5 Distribution of CLV Cluster by Marital Status¶

5.2.6 Distribution of Average Distance Cluster by Gender¶

5.2.7 Distribution of Average Distance Cluster by Salary¶

5.2.8 Distribution of Average Distance Cluster by Education¶

5.2.9 Distribution of Average Distance Cluster by Loyalty Card¶

5.2.10 Distribution of Average Distance Cluster by Enrollment Type¶

5.2.11 Distribution of Point Usage Cluster by Loyalty Card¶

5.2.12 Distribution of Point Usage Cluster by Enrollment Type¶

5.2.13 Distribution of Point Usage Cluster by Education¶

5.2.14 Distribution of Point Usage Cluster by Marital Status¶

5.3 Clustering based on all features¶

5.4 Determining Optimal Coupon Values by different segmentation methods¶

5.4.1 Promotions based on no differentiation (one large customer group)¶