DTSA 5505 - Ask A Manager Survey Reported Salaries from 2022 and 2023

# Packages
## Various
import pickle
import random

## Dates and Times
import pytz
from dateutil import parser
from datetime import timezone,datetime,timedelta
from dateutil.relativedelta import relativedelta

## String manipulation
import string
import re #regex
import json #json parsing
## Data Generation
import random

# Data Manipulation
import pandas as pd
import numpy as np
## Vector Math
from numpy.linalg import norm
## Math
import math
## Statistical Tests
from scipy.stats import ttest_ind
from scipy.stats import ttest_rel
from scipy.stats import median_test

# Data Visualization
import matplotlib as mpl
import matplotlib.pyplot as plt
import plotly as pl
import plotly.express as px
import plotly.graph_objects as go

# Language Processing
import spacy
nlp = spacy.load("en_core_web_md")

# SciKit Learn
from sklearn.cluster import KMeans
from sklearn.cluster import Birch
from sklearn.cluster import OPTICS
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.feature_extraction import FeatureHasher
from sklearn.model_selection import train_test_split

# Helper Functions
def clean_cat(ind, area, flag):
    if flag == 1:
        if ind == 'NAN':
            if area != 'NAN':
                return area
        return ind
    if flag == 2:
        if area == 'NAN':
            if ind != 'NAN':
                return ind
        return area

# 2022 Form Responses
aam_2022 = pd.read_csv('aam_survey_2022.csv')
aam_2022 = aam_2022.rename(columns = {'Timestamp':'form_timestamp', 'How old are you?':'age_group', 
                                      'What industry is your employer in?':'industry',
       'What is the functional area of your job (this might be different from your companys industry)?':'functional_area',
       'Job title':'job_title',
       'If your job title needs additional context, please clarify here:':'add_job_context',
       'What is your annual salary? This should be your GROSS (pre-tax) income. (Youll indicate the currency in a later question.) If you are part-time or hourly, please enter an annualized equivalent -- what you would earn if you worked the job 40 hours a week, 52 weeks a year.':'gross_salary',
       'How much additional monetary compensation do you get, if any (for example, bonuses or overtime in an average year)? Only include monetary compensation here, not the value of benefits, tuition reimbursement, etc. If your bonus or overtime varies from year to year, use the most recent figures.':'bonus',
       'Please indicate the currency':'currency',
       'If Other, please indicate the currency here: ':'currency_text',
       'If your income needs additional context, please provide it here:':'add_salary_context',
       'What country do you work in? (Countries listed had by far the largest representation last year. Please write in your country if its not listed.)':'country',
       'If youre in the U.S., what state do you work in?':'us_state',
       'What city do you work in?':'city', 'Are you remote or on-site?':'work_type',
       'How many years of professional work experience do you have overall?':'exp_group',
       'How many years of professional work experience do you have in your field?':'exp_field_group',
       'What is your highest level of education completed?':'education',
       'What is your gender?':'gender', 'What is your race? (Choose all that apply.)':'race'})

# 2023 Form Responses
aam_2023 = pd.read_csv('aam_survey_2023.csv')
aam_2023 = aam_2023.rename(columns = {'Timestamp':'form_timestamp', 'How old are you?':'age_group', 
                                      'Industry':'industry', 'Functional area of job':'functional_area',
       'Job title':'job_title', 'Job title - additional context':'add_job_context', 'Annual salary (gross)':'gross_salary',
       'Additional monetary compensation':'bonus', 'Currency':'currency', 'Currency - other':'currency_text',
       'Income - additional context':'add_salary_context', 'Country':'country', 'State':'us_state', 'City':'city',
       'Remote or on-site?':'work_type', 'Years of experience, overall':'exp_group',
       'Years of experience in field':'exp_field_group', 'Highest level of education completed':'education',
       'Gender':'gender', 'Race':'race'})

aam_2022['year'] = 2022
aam_2023['year'] = 2023

# Fill NAs for Categorical Variables
aam_2022[['country','us_state','city','work_type','gender','race']] = aam_2022[['country','us_state','city','work_type','gender','race']].fillna(value='Unknown')
aam_2023[['country','us_state','city','work_type','gender','race']] = aam_2023[['country','us_state','city','work_type','gender','race']].fillna(value='Unknown')

# Clean Categories
aam_2022['race'] = aam_2022['race'].apply(lambda x: str(x).replace('Hispanic, Latino, or Spanish origin', 'Hispanic or Latino or Spanish origin'))
aam_2023['race'] = aam_2023['race'].apply(lambda x: str(x).replace('Hispanic, Latino, or Spanish origin', 'Hispanic or Latino or Spanish origin'))

aam_2022['race'] = aam_2022['race'].apply(lambda x: str(x).replace('Another option not listed here or prefer not to answer', 'Unknown'))
aam_2023['race'] = aam_2023['race'].apply(lambda x: str(x).replace('Another option not listed here or prefer not to answer', 'Unknown'))

aam_2022['gender'] = aam_2022['gender'].apply(lambda x: str(x).replace('Other or prefer not to answer', 'Unknown'))
aam_2023['gender'] = aam_2023['gender'].apply(lambda x: str(x).replace('Other or prefer not to answer', 'Unknown'))

aam_2022['gender'] = aam_2022['gender'].apply(lambda x: str(x).split(',')[0])
aam_2023['gender'] = aam_2023['gender'].apply(lambda x: str(x).split(',')[0])

aam_2022['race'] = aam_2022['race'].apply(lambda x: str(x).split(',')[0])
aam_2023['race'] = aam_2023['race'].apply(lambda x: str(x).split(',')[0])

aam_2022['work_type'] = aam_2022['work_type'].apply(lambda x: str(x).replace("Other/it's complicated", 'Other'))
aam_2023['work_type'] = aam_2023['work_type'].apply(lambda x: str(x).replace("Other/it's complicated", 'Other'))
aam_2022['work_type'] = aam_2022['work_type'].apply(lambda x: str(x).replace("Fully remote, Hybrid", 'Other'))

aam_2022['country'] = aam_2022['country'].apply(lambda x: str(x).split(',')[0])
aam_2022['us_state'] = aam_2022['us_state'].apply(lambda x: str(x).split(',')[0])
aam_2022['industry'] = aam_2022['industry'].apply(lambda x: str(x).replace('Accounting, Banking & Finance', 'Accounting & Banking & Finance'))
aam_2022['industry'] = aam_2022['industry'].apply(lambda x: str(x).replace('Leisure, Sport & Tourism', 'Leisure & Sport & Tourism'))
aam_2022['industry'] = aam_2022['industry'].apply(lambda x: str(x).replace('Marketing, Advertising & PR', 'Marketing & Advertising & PR'))
aam_2022['industry'] = aam_2022['industry'].apply(lambda x: str(x).split(',')[0])
aam_2022['functional_area'] = aam_2022['functional_area'].apply(lambda x: str(x).replace('Accounting, Banking & Finance', 'Accounting & Banking & Finance'))
aam_2022['functional_area'] = aam_2022['functional_area'].apply(lambda x: str(x).replace('Leisure, Sport & Tourism', 'Leisure & Sport & Tourism'))
aam_2022['functional_area'] = aam_2022['functional_area'].apply(lambda x: str(x).replace('Marketing, Advertising & PR', 'Marketing & Advertising & PR'))
aam_2022['functional_area'] = aam_2022['functional_area'].apply(lambda x: str(x).split(',')[0])

aam_2023['industry'] = aam_2023['industry'].apply(lambda x: str(x).replace('Accounting, Banking & Finance', 'Accounting & Banking & Finance'))
aam_2023['industry'] = aam_2023['industry'].apply(lambda x: str(x).replace('Leisure, Sport & Tourism', 'Leisure & Sport & Tourism'))
aam_2023['industry'] = aam_2023['industry'].apply(lambda x: str(x).replace('Marketing, Advertising & PR', 'Marketing & Advertising & PR'))
aam_2023['functional_area'] = aam_2023['functional_area'].apply(lambda x: str(x).replace('Accounting, Banking & Finance', 'Accounting & Banking & Finance'))
aam_2023['functional_area'] = aam_2023['functional_area'].apply(lambda x: str(x).replace('Leisure, Sport & Tourism', 'Leisure & Sport & Tourism'))
aam_2023['functional_area'] = aam_2023['functional_area'].apply(lambda x: str(x).replace('Marketing, Advertising & PR', 'Marketing & Advertising & PR'))

# Fill NAs for Numberic Variables
aam_2022['gross_salary'] = aam_2022['gross_salary'].apply(lambda x: str(x).replace(',',''))
aam_2022['bonus'] = aam_2022['bonus'].apply(lambda x: str(x).replace(',',''))
aam_2022 = aam_2022.astype({'gross_salary':'float','bonus':'float'})
aam_2022['bonus'].fillna(value=0, inplace=True)

aam_2023['gross_salary'] = aam_2023['gross_salary'].apply(lambda x: str(x).replace(',',''))
aam_2023['bonus'] = aam_2023['bonus'].apply(lambda x: str(x).replace(',',''))
aam_2023 = aam_2023.astype({'gross_salary':'float','bonus':'float'})
aam_2023['bonus'].fillna(value=0, inplace=True)

# Drop Potentially Bad or Incomplete Data
print('2022 Total Rows:',aam_2022.shape[0])
print('2023 Total Rows:',aam_2023.shape[0])
aam_2022_rows = aam_2022.shape[0]
aam_2023_rows = aam_2023.shape[0]

aam_2022.dropna(subset=['gross_salary','age_group','industry','functional_area','job_title','education'], inplace=True)
aam_2023.dropna(subset=['gross_salary','age_group','industry','functional_area','job_title','education'], inplace=True)

# Additional Cleaning and Filtering for Valid Education Category
edu_lvls = ['High School', 'Some college', 'College degree', "Master's degree", 'PhD', 'Professional degree (MD, JD, etc.)']

aam_2023.loc[aam_2023['education'].str.contains('Associate', case=False), 'education'] = 'Some college'
aam_2023.loc[aam_2023['education'].str.contains('Vocational', case=False), 'education'] = 'Some college'
aam_2023.loc[aam_2023['education'].str.contains('Technical', case=False), 'education'] = 'Some college'
aam_2023.loc[aam_2023['education'].str.contains('Trade', case=False), 'education'] = 'Some college'
aam_2023.loc[aam_2023['education'].str.contains('Bachelor', case=False), 'education'] = 'College degree'
aam_2023.loc[aam_2023['education'].str.contains('B.A ', case=False), 'education'] = 'College degree'
aam_2023.loc[aam_2023['education'].str.contains('B.S ', case=False), 'education'] = 'College degree'
aam_2023.loc[aam_2023['education'].str.contains('BA ', case=False), 'education'] = 'College degree'
aam_2023.loc[aam_2023['education'].str.contains('BS ', case=False), 'education'] = 'College degree'
aam_2023.loc[aam_2023['education'].str.contains('college degree', case=False), 'education'] = 'College degree'
aam_2023.loc[aam_2023['education'].str.contains("Master's Degree", case=False), 'education'] = "Master's degree"
aam_2023.loc[aam_2023['education'].str.contains("MS ", case=False), 'education'] = "Master's degree"
aam_2023.loc[aam_2023['education'].str.contains("EdD ", case=False), 'education'] = 'Professional degree'

aam_2023 = aam_2023[aam_2023['education'].isin(edu_lvls)]

# Formatting string Fields
aam_2022['industry'] = aam_2022['industry'].apply(lambda x: x.upper().rstrip().lstrip())
aam_2022['functional_area'] = aam_2022['functional_area'].apply(lambda x: x.upper().rstrip().lstrip())
aam_2022['job_title'] = aam_2022['job_title'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())
aam_2022['country'] = aam_2022['country'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())
aam_2022['currency'] = aam_2022['currency'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())
aam_2022['city'] = aam_2022['city'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())

aam_2023['industry'] = aam_2023['industry'].apply(lambda x: x.upper().rstrip().lstrip())
aam_2023['functional_area'] = aam_2023['functional_area'].apply(lambda x: x.upper().rstrip().lstrip())
aam_2023['job_title'] = aam_2023['job_title'].apply(lambda x: x.upper().rstrip().lstrip())
aam_2023['country'] = aam_2023['country'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())
aam_2023['currency'] = aam_2023['currency'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())
aam_2022['city'] = aam_2022['city'].apply(lambda x: x.upper().rstrip().rstrip().lstrip())

# Fill in NAN's for industry and functional_area
aam_2022['industry'] = aam_2022.apply(lambda x: clean_cat(x['industry'], x['functional_area'], 1), axis=1)
aam_2022['functional_area'] = aam_2022.apply(lambda x: clean_cat(x['industry'], x['functional_area'], 2), axis=1)
aam_2023['industry'] = aam_2023.apply(lambda x: clean_cat(x['industry'], x['functional_area'], 1), axis=1)
aam_2023['functional_area'] = aam_2023.apply(lambda x: clean_cat(x['industry'], x['functional_area'], 2), axis=1)

aam_2022 = aam_2022[aam_2022['functional_area'] != 'NAN']
aam_2023 = aam_2023[aam_2023['functional_area'] != 'NAN']
aam_2022 = aam_2022[aam_2022['job_title'].str.count('[A-z]')>0]
aam_2023 = aam_2023[aam_2023['job_title'].str.count('[A-z]')>0]

print('2022 Complete Rows:',aam_2022.shape[0])
print('2023 Complete Rows:',aam_2023.shape[0])
print('2022 Row Loss:',aam_2022_rows - aam_2022.shape[0])
print('2023 Row Loss:',aam_2023_rows - aam_2023.shape[0])

# Ordinal-ization of Categories
exp_groups = pd.DataFrame(np.unique(aam_2022['exp_group']), columns = {'group'})
exp_groups['levels'] = [1,5,2,6,7,8,3,4]

aam_2022.loc[aam_2022['age_group'].str.contains('under 18', case=False), 'age_group'] = '18 or under'
aam_2022['age_group_o'] = pd.factorize(aam_2022['age_group'], sort=True)[0]+1
aam_2022 = aam_2022.merge(exp_groups, how='left', left_on='exp_group', right_on='group')
aam_2022 = aam_2022.drop(columns={'group'})
aam_2022 = aam_2022.rename(columns={'levels':'exp_group_o'})
aam_2022 = aam_2022.merge(exp_groups, how='left', left_on='exp_field_group', right_on='group')
aam_2022 = aam_2022.drop(columns={'group'})
aam_2022 = aam_2022.rename(columns={'levels':'exp_field_group_o'})

aam_2023.loc[aam_2023['age_group'].str.contains('under 18', case=False), 'age_group'] = '18 or under'
aam_2023['age_group_o'] = pd.factorize(aam_2023['age_group'], sort=True)[0]+1
aam_2023 = aam_2023.merge(exp_groups, how='left', left_on='exp_group', right_on='group')
aam_2023 = aam_2023.drop(columns={'group'})
aam_2023 = aam_2023.rename(columns={'levels':'exp_group_o'})
aam_2023 = aam_2023.merge(exp_groups, how='left', left_on='exp_field_group', right_on='group')
aam_2023 = aam_2023.drop(columns={'group'})
aam_2023 = aam_2023.rename(columns={'levels':'exp_field_group_o'})

edu_groups = pd.DataFrame(np.unique(aam_2022['education']), columns = {'group'})
edu_groups['levels'] = [3,1,4,6,5,2]

aam_2022 = aam_2022.merge(edu_groups, how='left', left_on='education', right_on='group')
aam_2022 = aam_2022.drop(columns={'group'})
aam_2022 = aam_2022.rename(columns={'levels':'education_o'})

aam_2023 = aam_2023.merge(edu_groups, how='left', left_on='education', right_on='group')
aam_2023 = aam_2023.drop(columns={'group'})
aam_2023 = aam_2023.rename(columns={'levels':'education_o'})

# Adding Additional Columns
aam_2022['total_gross_salary'] = aam_2022['gross_salary'] + aam_2022['bonus']
aam_2023['total_gross_salary'] = aam_2023['gross_salary'] + aam_2023['bonus']

# Combined Data Frame
aam = pd.concat([aam_2022, aam_2023], ignore_index=True)

print('Columns:', len(aam.columns))

# Category Counts
print('Age Groups:',len(np.unique(aam['age_group'])))
print('Industries:',len(np.unique(aam['industry'])))
print('Functional Areas:',len(np.unique(aam['functional_area'])))
print('Job Titles:',len(np.unique(aam['job_title'])))
print('Max Gross Salary:',max(aam['gross_salary']))
print('Minimum Gross Salary:',min(aam['gross_salary']))
print('Max Bonus Compensation:',max(aam['bonus']))
print('Minimum Bonus Compensation:',min(aam['bonus']))
print('Countries:',len(np.unique(aam['country'])))
print('US States:',len(np.unique(aam['us_state'])))
print('Cities:',len(np.unique(aam['city'])))
print('Work Type:',len(np.unique(aam['work_type'])))
print('Experience Groups:',len(np.unique(aam['exp_group'])))
print('Field Experience Groups:',len(np.unique(aam['exp_field_group'])))
print('Education:',len(np.unique(aam['education'])))
print('Genders:',len(np.unique(aam['gender'])))
print('Race:',len(np.unique(aam['race'])))

# Most Responses Fall Into a Small Subset of All Reported Industries/Functional Areas
unique_industries = aam[['industry','gross_salary']].groupby(['industry'], as_index=False).count().sort_values('gross_salary', ascending=False)
print('Total Industries:', len(unique_industries))
print('Total Counts:', sum(unique_industries['gross_salary']))
print('Industries With Fewer than 10 Responses:', len(unique_industries[unique_industries['gross_salary'] < 10]))
print('Total Counts for Industries With Fewer than 10 Responses:', sum(unique_industries[unique_industries['gross_salary'] < 10]['gross_salary']))
print('Percent of Total:', sum(unique_industries[unique_industries['gross_salary'] < 10]['gross_salary'])/sum(unique_industries['gross_salary'])*100,'%')

print('')
unique_areas = aam[['functional_area','gross_salary']].groupby(['functional_area'], as_index=False).count().sort_values('gross_salary', ascending=False)
print('Total Functional Areas:', len(unique_areas))
print('Total Counts:', sum(unique_areas['gross_salary']))
print('Industries With Fewer than 10 Responses:', len(unique_areas[unique_areas['gross_salary'] < 10]))
print('Total Counts for Industries With Fewer than 10 Responses:', sum(unique_areas[unique_areas['gross_salary'] < 10]['gross_salary']))
print('Percent of Total:', sum(unique_areas[unique_areas['gross_salary'] < 10]['gross_salary'])/sum(unique_areas['gross_salary'])*100,'%')

print('')
unique_jobs = aam[['job_title','gross_salary']].groupby(['job_title'], as_index=False).count().sort_values('gross_salary', ascending=False)
print('Total Job Titles Areas:', len(unique_jobs))
print('Total Counts:', sum(unique_jobs['gross_salary']))
print('Industries With Fewer than 10 Responses:', len(unique_jobs[unique_jobs['gross_salary'] < 10]))
print('Total Counts for Job Titles With Fewer than 10 Responses:', sum(unique_jobs[unique_jobs['gross_salary'] < 10]['gross_salary']))
print('Percent of Total:', sum(unique_jobs[unique_jobs['gross_salary'] < 10]['gross_salary'])/sum(unique_jobs['gross_salary'])*100,'%')

# Removing All Industries and Functional Areas With Less Than 10 Responses
small_industries = np.unique(unique_industries[unique_industries['gross_salary'] < 10]['industry'])
small_areas = np.unique(unique_areas[unique_areas['gross_salary'] < 10]['functional_area'])

aam_2022 = aam_2022[aam_2022['industry'].isin(small_industries) == False]
aam_2023 = aam_2023[aam_2023['industry'].isin(small_industries) == False]
aam = aam[aam['industry'].isin(small_industries) == False]

aam_2022 = aam_2022[aam_2022['functional_area'].isin(small_areas) == False]
aam_2023 = aam_2023[aam_2023['functional_area'].isin(small_areas) == False]
aam = aam[aam['functional_area'].isin(small_areas) == False]

aam_2022 = aam_2022.reset_index(drop=True)
aam_2023 = aam_2023.reset_index(drop=True)
aam = aam.reset_index(drop=True)

print('2022 After Removal of Small Industries/Areas:',aam_2022.shape[0])
print('2023 After Removal of Small Industries/Areas:',aam_2023.shape[0])
print('All Data After Removal of Small Industries/Areas:',aam.shape[0])
print('Industries:',len(np.unique(aam['industry'])))
print('Functional Areas:',len(np.unique(aam['functional_area'])))

# NLP Example
print(nlp(aam.loc[0,'job_title']).text,":",nlp(aam.loc[0,'job_title']).vector_norm)
print('DATA SCIENTIST:',nlp('DATA SCIENTIST').vector_norm)
print('data scientist:',nlp('data scientist').vector_norm)
print('scientist data:',nlp('scientist data').vector_norm)
print('scientist:',nlp('scientist').vector_norm)
print('DEBT ADVISOR:',nlp('DEBT ADVISOR').vector_norm)
print('Senior Data Scientist to Data Scientist Similarity:', nlp(aam.loc[0,'job_title']).similarity(nlp('DATA SCIENTIST')))

# print('DATA SCIENTIST:',nlp('DATA SCIENTIST').vector)

# Add vector to the data_set (SLOW)
aam['job_title_vector'] = aam['job_title'].apply(lambda x: nlp(str(x)).vector_norm)

search_vector_norm = aam.loc[0,'job_title_vector']
similarity_band = 0.001
print("0.001 Similarity:", np.unique(aam[(aam['job_title_vector'] <= search_vector_norm + similarity_band) \
    & (aam['job_title_vector'] >= search_vector_norm - similarity_band)]['job_title']))

print("")
similarity_band = 0.005
print("0.005 Similarity:", np.unique(aam[(aam['job_title_vector'] <= search_vector_norm + similarity_band) \
    & (aam['job_title_vector'] >= search_vector_norm - similarity_band)]['job_title']))

print("")
similarity_band = 0.01
print("0.01 Similarity:", np.unique(aam[(aam['job_title_vector'] <= search_vector_norm + similarity_band) \
    & (aam['job_title_vector'] >= search_vector_norm - similarity_band)]['job_title']))

# Add vector to the data_set (SLOW)
aam['combined_title'] = aam['industry'] + ' ' + aam['functional_area'] + ' ' + aam['job_title']
aam['combined_title_vector'] = aam['combined_title'].apply(lambda x: nlp(str(x)).vector_norm)

search_vector_norm = aam.loc[0,'combined_title_vector']
similarity_band = 0.0005
print("0.0005 Similarity:", np.unique(aam[(aam['combined_title_vector'] <= search_vector_norm + similarity_band) \
    & (aam['combined_title_vector'] >= search_vector_norm - similarity_band)]['job_title']))

print("")
similarity_band = 0.001
print("0.001 Similarity:", np.unique(aam[(aam['combined_title_vector'] <= search_vector_norm + similarity_band) \
    & (aam['combined_title_vector'] >= search_vector_norm - similarity_band)]['job_title']))

# print("")
# similarity_band = 0.005
# print("0.005 Similarity:", np.unique(aam[(aam['combined_title_vector'] <= search_vector_norm + similarity_band) \
#     & (aam['combined_title_vector'] >= search_vector_norm - similarity_band)]['job_title']))

# Preprocessing Data Pipeline
numeric_features = ['job_title_vector']
numeric_transformer = Pipeline( #no missing values, so no need for an imputer
    steps=[('scaler', StandardScaler())]
)

preprocessor = ColumnTransformer(
    transformers = [
        ('num', numeric_transformer, numeric_features)
    ]
)

first_pipeline = Pipeline(
    steps=[('preprocessor', preprocessor)
    ], verbose=True
)

# hmm.Subset Data, Transform
aam_subset = aam.loc[0:999]
aam_cluster = aam_subset
aam_cluster = aam_cluster[['job_title','age_group_o','exp_group_o','exp_field_group_o',
                           'education_o','total_gross_salary','job_title_vector',
                          'combined_title_vector']]
aam_cluster = aam_cluster.groupby(['job_title'], as_index=False).mean()

# Data Pipeline and Fit K-Means Model
test_df1 = first_pipeline.fit_transform(aam_cluster)
first_model = KMeans(n_clusters=20).fit(test_df1)
pickle.dump(first_model, open('kmeans_first_test_model.sav', 'wb')) #save model?
aam_cluster = aam_cluster.assign(cluster = first_model.labels_)

# Creating Cluster Lists
cluster_list = []
for i in range(20):
    exec("cluster_%d = []" % (i))
    exec("cluster_%d = np.unique(aam_cluster[aam_cluster['cluster'] == %d]['job_title'])" % (i,i))
    exec("cluster_list.append(cluster_%d)" % (i))
    
# Checking for Duplicates
overlap = []
for i in range(len(cluster_list)-1):
    for j in range(i+1,20-i):
        if len(set(cluster_list[i]) & set(cluster_list[j])) > 1:
            overlap.append([i,j,set(cluster_list[i]) & set(cluster_list[j])])
overlap

cluster_list

# Preprocessing Data Pipeline
preprocessor = ColumnTransformer(
    transformers = [
#         ('num', numeric_transformer, numeric_features),
        ('str1', FeatureHasher(input_type='string', n_features=1000), 'job_title')
#         ('str2', FeatureHasher(input_type='string', n_features=1000), 'combined_title')
    ]
)

second_pipeline = Pipeline(
    steps=[('preprocessor', preprocessor)
    ], verbose=True
)

# K-Means Cluster
aam_subset = aam.loc[0:999]
aam_cluster2 = aam_subset
aam_cluster2 = aam_cluster2[['job_title','age_group_o','exp_group_o','exp_field_group_o',
                           'education_o','total_gross_salary']]
aam_cluster2 = aam_cluster2.groupby(['job_title'], as_index=False).median()

test_df2 = second_pipeline.fit_transform(aam_cluster2)
cluster_model = KMeans(n_clusters=20).fit(test_df2)
pickle.dump(cluster_model, open('kmeans_second_test_model.sav', 'wb')) #save model?
aam_cluster2 = aam_cluster2.assign(cluster = cluster_model.labels_)

# Creating Cluster Lists
cluster_list = []
for i in range(20):
    exec("cluster_%d = []" % (i))
    exec("cluster_%d = np.unique(aam_cluster2[aam_cluster2['cluster'] == %d]['job_title'])" % (i,i))
    exec("cluster_list.append(cluster_%d)" % (i))
    
# Checking for Duplicates
overlap = []
for i in range(len(cluster_list)-1):
    for j in range(i+1,20-i):
        if len(set(cluster_list[i]) & set(cluster_list[j])) > 1:
            overlap.append([i,j,set(cluster_list[i]) & set(cluster_list[j])])
overlap

cluster_list

# birch
from sklearn.cluster import Birch
op_model = Birch(n_clusters = 20).fit(test_df2)
aam_cluster2 = aam_cluster2.assign(cluster3 = op_model.labels_)

cluster_list = []
for i in range(max(aam_cluster2['cluster3'])):
    exec("cluster_%d = []" % (i))
    exec("cluster_%d = np.unique(aam_cluster2[aam_cluster2['cluster3'] == %d]['job_title'])" % (i,i))
    exec("cluster_list.append(cluster_%d)" % (i))
print(cluster_list)

from sklearn.cluster import OPTICS

op_model = OPTICS(min_samples = 5).fit(test_df2.toarray())
aam_cluster2 = aam_cluster2.assign(cluster3 = op_model.labels_)

cluster_list = []
for i in range(max(aam_cluster2['cluster3'])):
    exec("cluster_%d = []" % (i))
    exec("cluster_%d = np.unique(aam_cluster2[aam_cluster2['cluster3'] == %d]['job_title'])" % (i,i))
    exec("cluster_list.append(cluster_%d)" % (i))
print(cluster_list)
print('Unlabeled:', sum(op_model.labels_ == -1) / len(op_model.labels_)*100,'%')

# this takes a very long time
aam_cluster_final = aam
aam_cluster_final = aam_cluster_final[['job_title','age_group_o','exp_group_o','exp_field_group_o',
                           'education_o','total_gross_salary']]
aam_cluster_final = aam_cluster_final.groupby(['job_title'], as_index=False).median()

aam_transform = second_pipeline.fit_transform(aam_cluster_final)
cluster_model = KMeans(n_clusters=2000).fit(aam_transform)
pickle.dump(cluster_model, open('kmeans_final_model.sav', 'wb')) #save model?
aam_cluster_final = aam_cluster_final.assign(cluster = cluster_model.labels_)
aam_cluster_final.to_csv('aam_usd_cluster.csv', index=False, encoding='utf-8')

cluster_list = []
for i in range(max(aam_cluster_final['cluster'])):
    exec("cluster_%d = []" % (i))
    exec("cluster_%d = np.unique(aam_cluster_final[aam_cluster_final['cluster'] == %d]['job_title'])" % (i,i))
    exec("cluster_list.append(cluster_%d)" % (i))
    
# Checking for Duplicates
overlap = []
for i in range(len(cluster_list)-1):
    for j in range(i+1,len(cluster_list)-i):
        if len(set(cluster_list[i]) & set(cluster_list[j])) > 1:
            overlap.append([i,j,set(cluster_list[i]) & set(cluster_list[j])])
overlap
# cluster_list

# takes even longer
aam_cluster_final = aam
aam_cluster_final = aam_cluster_final[['job_title','age_group_o','exp_group_o','exp_field_group_o',
                           'education_o','total_gross_salary']]
aam_cluster_final = aam_cluster_final.groupby(['job_title'], as_index=False).median()

aam_transform = second_pipeline.fit_transform(aam_cluster_final)
cluster_model = OPTICS(min_samples = 3).fit(aam_transform.toarray())
pickle.dump(cluster_model, open('optics_final_model.sav', 'wb')) #save model?
aam_cluster_final = aam_cluster_final.assign(cluster2 = cluster_model.labels_)
aam_cluster_final.to_csv('aam_usd_cluster_optics.csv', index=False, encoding='utf-8')

cluster_list = []
for i in range(max(aam_cluster_final['cluster2'])):
    exec("cluster_%d = []" % (i))
    exec("cluster_%d = np.unique(aam_cluster_final[aam_cluster_final['cluster2'] == %d]['job_title'])" % (i,i))
    exec("cluster_list.append(cluster_%d)" % (i))
    
# Checking for Duplicates
overlap = []
for i in range(len(cluster_list)-1):
    for j in range(i+1,len(cluster_list)-i):
        if len(set(cluster_list[i]) & set(cluster_list[j])) > 1:
            overlap.append([i,j,set(cluster_list[i]) & set(cluster_list[j])])
overlap
max(aam_cluster_final['cluster2'])

# Preprocessing Skip
aam = pd.read_csv('aam_clean.csv')
aam_cluster = pd.read_csv('aam_usd_cluster.csv')
aam_cluster = aam_cluster.astype({'cluster':'string'})

# Join with Cluster Association
aam = aam.merge(aam_cluster[['job_title','cluster']], how='left', on=['job_title'])
print('Missing Job Title Clusters:',len(aam[aam['cluster'].isna()]))

# K-Means Cluster
cluster_counts = aam_cluster[['job_title','cluster']].groupby(['cluster'], as_index=False).count()
print("Clusters With Only 1 Member:", len(cluster_counts[cluster_counts['job_title']==1]))

# OPTICS
aam_cluster2 = pd.read_csv('aam_usd_cluster_optics.csv')
# aam_cluster2 = aam_cluster2.astype({'cluster2':'string'})

cluster_counts = aam_cluster2[['job_title','cluster2']].groupby(['cluster2'], as_index=False).count()
print("Total Clusters:", max(cluster_counts['cluster2']))
print("Unlabeled Rows:", len(aam_cluster2[aam_cluster2['cluster2'] == -1]) / len(aam_cluster2)*100,'%')

random.seed(30)
k_means_rand = random.choices(np.unique(aam_cluster['cluster']), k=10)
optics_list = []

print('K-Means Clusters:')
for i in k_means_rand:
    print('Cluster:',i, np.unique(aam_cluster[aam_cluster['cluster'] == i]['job_title']))
    optics_list.append(np.unique(aam_cluster[aam_cluster['cluster'] == i]['job_title'])[0])

print('')
print('OPTICS Clusters:')
for j in optics_list:
    o_c = int(aam_cluster2[aam_cluster2['job_title'] == j]['cluster2'])
    if o_c != -1:
        print('Cluster:',o_c, np.unique(aam_cluster2[aam_cluster2['cluster2'] == o_c]['job_title']))

df = aam[['form_timestamp','country']].groupby(['country'], as_index=False).count()
df.loc[df['form_timestamp'] <= 50, 'country'] = 'OTHER'
df = df.groupby(['country'], as_index=False).sum()
fig = px.pie(df, names='country', values='form_timestamp', 
             title='Breakdown by Country', labels={'country':'Country', 'form_timestamp':'Responses'})
fig.update_traces(textposition='inside', textinfo='percent+label')
# fig.show()
fig.write_image("fig1.png")

aam['currency'] = aam['currency'].apply(lambda x: x.upper())
df = aam[['form_timestamp','currency']].groupby(['currency'], as_index=False).count()
df.loc[df['form_timestamp'] <= 50, 'currency'] = 'OTHER'
df = df.groupby(['currency'], as_index=False).sum()
fig = px.pie(df, names='currency', values='form_timestamp', 
             title='Breakdown by Currency', labels={'currency':'Currency', 'form_timestamp':'Responses'})
fig.update_traces(textposition='inside', textinfo='percent+label')
# fig.show()
fig.write_image("fig2.png")

aam_usd = aam[aam['currency'] == 'USD']

df = aam_usd[['form_timestamp','industry']].groupby(['industry'], as_index=False).count()
df.loc[df['form_timestamp'] <= 500, 'industry'] = 'OTHER'
df = df.groupby(['industry'], as_index=False).sum()
fig = px.pie(df, names='industry', values='form_timestamp', 
             title='Breakdown by Industry', labels={'industry':'Industry', 'form_timestamp':'Responses'})
fig.update_traces(textposition='inside', textinfo='percent+label')
# fig.show()
fig.write_image("fig3.png")

df = aam_usd[['form_timestamp','functional_area']].groupby(['functional_area'], as_index=False).count()
df.loc[df['form_timestamp'] <= 500, 'functional_area'] = 'OTHER'
df = df.groupby(['functional_area'], as_index=False).sum()
fig = px.pie(df, names='functional_area', values='form_timestamp', 
             title='Breakdown by Functional Area', labels={'functional_area':'Functional Area', 'form_timestamp':'Responses'})
fig.update_traces(textposition='inside', textinfo='percent+label')
# fig.show()
fig.write_image("fig4.png")

df = aam_usd[(aam_usd['total_gross_salary'] <= 600000)]
fig = px.violin(df, x='year',y='total_gross_salary', points='outliers', title='Box Plot of Total Salary by Year', 
            color='year', box=True, labels={'year':'Year','total_gross_salary':'Total Salary'})
# fig.show()
fig.write_image("fig5.png")

df = aam_usd[(aam_usd['total_gross_salary'] <= 600000)]
print('2022 Quantiles:')
print('Q3: $%.2f' %(df[df['year'] == 2022]['total_gross_salary'].quantile([0.75])))
print('Median: $%.2f' %(df[df['year'] == 2022]['total_gross_salary'].quantile([0.5])))
print('Q1: $%.2f' %(df[df['year'] == 2022]['total_gross_salary'].quantile([0.25])))
print('')
print('2023 Quantiles:')
print('Q3: $%.2f' %(df[df['year'] == 2023]['total_gross_salary'].quantile([0.75])))
print('Median: $%.2f' %(df[df['year'] == 2023]['total_gross_salary'].quantile([0.5])))
print('Q1: $%.2f' %(df[df['year'] == 2023]['total_gross_salary'].quantile([0.25])))

df = aam_usd[(aam_usd['total_gross_salary'] >= 0)]
df = df.astype({'year':'string'})
fig = px.histogram(df, x='year',y='total_gross_salary', title='Industry Average Total Salary in 2022 vs. 2023',
            color='year', labels={'industry':'Industry','total_gross_salary':'Total Salary', 'year':'Year'}, 
            histfunc='avg', animation_frame='industry', animation_group='industry', range_y = [0,200000])
# fig.show()
fig.write_html('fig6.html')

df = aam_usd[(aam_usd['total_gross_salary'] >= 0)]
df = df.astype({'year':'string'})
fig = px.histogram(df, x='year',y='total_gross_salary', title='Functional Area Average Gross Salary in 2022 vs. 2023',
            color='year', labels={'functional_area':'Functional Area','total_gross_salary':'Total Salary', 'year':'Year'}, 
            histfunc='avg', animation_frame='functional_area', animation_group='functional_area', range_y = [0,200000])
# fig.show()
fig.write_html('fig7.html')

df = aam_usd[(aam_usd['total_gross_salary'] <= 600000)]
fig = px.histogram(df, x='total_gross_salary', title='Total Salary Distribution', color='year',
                   labels={'total_gross_salary':'Total Salary', 'year':'Year'})
# fig.show()
fig.write_image("fig8.png")

# Log Normal Fit
df = aam_usd[(aam_usd['total_gross_salary'] <= 600000)]
s, loc, scale = stats.lognorm.fit(df['total_gross_salary'], floc=0)
z_x = list(range(0,600000,100))

z = stats.lognorm.pdf(z_x, s, scale=scale)

plt.hist(df['total_gross_salary'], bins = 100, density=True, alpha=0.75)
plt.plot(z_x, z)
plt.savefig('fig16.png')

print('Scale:',scale)
print('Shape',s)

# Student's t-test
df = aam_usd[(aam_usd['total_gross_salary'] <= 600000)]
salaries_23 = df[df['year'] == 2023]['total_gross_salary'].to_list()
salaries_22 = df[df['year'] == 2022]['total_gross_salary'].to_list()

stat, p = ttest_ind(salaries_23, salaries_22)
print('Test Statistic: %.3f, p-value: %.10f' % (stat, p))
print('Mean of Salaries from 2023: $%.2f' %(df[df['year'] == 2023]['total_gross_salary'].mean()))
print('Mean of Salaries from 2022: $%.2f' %(df[df['year'] == 2022]['total_gross_salary'].mean()))

# Mood's Median Test
stat, p, med, tabl = median_test(salaries_23, salaries_22)
print('Test Statistic: %.3f, p-value: %.17f' % (stat, p))
print('Median of Salaries from 2023: $%.2f' %(df[df['year'] == 2023]['total_gross_salary'].median()))
print('Median of Salaries from 2022: $%.2f' %(df[df['year'] == 2022]['total_gross_salary'].median()))

aam_usd23 = aam_usd[aam_usd['year'] == 2023]

df = aam_usd23[['form_timestamp','us_state']].groupby(['us_state'], as_index=False).count()
df.loc[df['form_timestamp'] <= 350, 'us_state'] = 'Other'
df = df.groupby(['us_state'], as_index=False).sum()
fig = px.pie(df, names='us_state', values='form_timestamp', 
             title='Breakdown by State', labels={'us_state':'State', 'form_timestamp':'Responses'})
fig.update_traces(textposition='inside', textinfo='percent+label')
# fig.show()
fig.write_image("fig9.png")

df = aam_usd23
df.loc[df['city'].str.contains('NEW YORK', case=False), 'city'] = 'NEW YORK'
df = df[['form_timestamp','city']].groupby(['city'], as_index=False).count()
df.loc[df['form_timestamp'] <= 150, 'city'] = 'OTHER'
df = df.groupby(['city'], as_index=False).sum()
fig = px.pie(df, names='city', values='form_timestamp', 
             title='Breakdown by City', labels={'city':'State', 'form_timestamp':'Responses'})
fig.update_traces(textposition='inside', textinfo='percent+label')
# fig.show()
fig.write_image("fig10.png")

df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
fig = px.box(df, x='age_group',y='total_gross_salary', points='outliers', title='Box Plot of Total Salary by Age', 
            color='age_group', labels={'age_group':'Age Group','total_gross_salary':'Total Salary'}, 
            category_orders={'age_group':{y: x for x, y in aam_usd23[['age_group','age_group_o']].groupby(['age_group','age_group_o'], as_index=False).count().to_dict()['age_group'].items()}})
# fig.show()
fig.write_image("fig11.png")

df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
fig = px.box(df, x='exp_group',y='total_gross_salary', points='outliers', title='Box Plot of Total Salary by Work Experience', 
            color='exp_group', labels={'exp_group':'Experience Group','total_gross_salary':'Total Salary'}, 
            category_orders={'exp_group': list(aam_usd23[['exp_group','exp_group_o']].groupby(['exp_group','exp_group_o'], as_index=False).count().sort_values('exp_group_o')['exp_group'])})
# fig.show()
fig.write_image("fig12.png")

df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
fig = px.box(df, x='exp_field_group',y='total_gross_salary', points='outliers', title='Box Plot of Total Gross Salary by Field Experience', 
            color='exp_field_group', labels={'exp_field_group':'Field Experience Group','total_gross_salary':'Total Salary'}, 
            category_orders={'exp_field_group': list(aam_usd23[['exp_field_group','exp_field_group_o']].groupby(['exp_field_group','exp_field_group_o'], as_index=False).count().sort_values('exp_field_group_o')['exp_field_group'])})
# fig.show()
fig.write_image("fig13.png")

df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
fig = px.box(df, x='work_type',y='total_gross_salary', points='outliers', title='Box Plot of Total Salary by Work Type', 
            color='work_type', labels={'work_type':'Work Type','total_gross_salary':'Total Salary'})
# fig.show()
fig.write_image("fig14.png")

# Mood's Median Test
df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
onsite = df[df['work_type'] == 'On-site']['total_gross_salary'].to_list()
remote = df[df['work_type'] == 'Fully remote']['total_gross_salary'].to_list()
hybrid = df[df['work_type'] == 'Hybrid']['total_gross_salary'].to_list()
other = df[df['work_type'] == 'Other']['total_gross_salary'].to_list()

stat, p, med, tabl = median_test(onsite, remote, hybrid, other)
print('Test Statistic: %.3f, p-value: %.10f' % (stat, p))
print('Median of Onsite Salaries: $%.2f' %(df[df['work_type'] == 'On-site']['total_gross_salary'].median()))
print('Median of Remote Salaries: $%.2f' %(df[df['work_type'] == 'Fully remote']['total_gross_salary'].median()))
print('Median of Hybrid Salaries: $%.2f' %(df[df['work_type'] == 'Hybrid']['total_gross_salary'].median()))
print('Median of Other Salaries: $%.2f' %(df[df['work_type'] == 'Other']['total_gross_salary'].median()))

df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
fig = px.box(df, x='gender',y='total_gross_salary', points='outliers', title='Box Plot of Total Salary by Gender', 
            color='gender', labels={'gender':'Gender','total_gross_salary':'Total Salary'})
# fig.show()
fig.write_image("fig15.png")

# Mood's Median Test
df = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
man = df[df['gender'] == 'Man']['total_gross_salary'].to_list()
woman = df[df['gender'] == 'Woman']['total_gross_salary'].to_list()

stat, p, med, tabl = median_test(man, woman)
print('Test Statistic: %.3f, p-value: %.10f' % (stat, p))
print('Median of Male Salaries: $%.2f' %(df[df['gender'] == 'Man']['total_gross_salary'].median()))
print('Median of Female Salaries: $%.2f' %(df[df['gender'] == 'Woman']['total_gross_salary'].median()))

# Preprocessing Skip
aam = pd.read_csv('aam_clean.csv')
aam_cluster = pd.read_csv('aam_usd_cluster.csv')
aam_cluster = aam_cluster.astype({'cluster':'string'})

# Join with Cluster Association
aam = aam.merge(aam_cluster[['job_title','cluster']], how='left', on=['job_title'])

# Only 2023 Data in USD under $600k
aam_usd = aam[aam['currency'] == 'USD']
aam_usd23 = aam_usd[aam_usd['year'] == 2023]
aam_usd23_f = aam_usd23[(aam_usd23['total_gross_salary'] <= 600000)]
aam_usd23_f = aam_usd23_f.drop(columns={'form_timestamp','add_job_context','currency',
                                       'currency_text','add_job_context','add_salary_context',
                                       'city'})

# Load K-Means Model
k_means_final_model = pickle.load(open('kmeans_final_model.sav', 'rb'))

# Pipeline
numeric_features = ['age_group_o','exp_group_o','exp_field_group_o', 'education_o']
numeric_transformer = Pipeline( #no missing values, so no need for an imputer
    steps=[('scaler', StandardScaler())]
)

categorical_features = ['industry','functional_area','work_type','us_state','gender','race','cluster']
categorical_transformer = Pipeline(
    steps=[('encoder',OneHotEncoder(handle_unknown='error'))]
)

preprocessor = ColumnTransformer(
    transformers = [
        ('num', numeric_transformer, numeric_features),
        ('cat', categorical_transformer, categorical_features)
    ]
)

salary_pipeline = Pipeline(
    steps=[('preprocessor', preprocessor)
    ], verbose=True
)

# 2nd Pipeline (no Job Cluster)
numeric_features = ['age_group_o','exp_group_o','exp_field_group_o', 'education_o']
numeric_transformer = Pipeline( #no missing values, so no need for an imputer
  steps=[('scaler', StandardScaler())]
)

categorical_features = ['industry','functional_area','work_type','us_state','gender','race']
categorical_transformer = Pipeline(
    steps=[('encoder',OneHotEncoder(handle_unknown='error'))]
)

preprocessor = ColumnTransformer(
    transformers = [
        ('num', numeric_transformer, numeric_features),
        ('cat', categorical_transformer, categorical_features)
    ]
)

salary_pipeline2 = Pipeline(
    steps=[('preprocessor', preprocessor)
    ], verbose=True
)