<jupyter_start><jupyter_text>Gender Pay Gap in the UK Business UnderstandingIn the UK, employers with 250 or more employees are required to calculate a series of metrics about their gender pay gap and publish them in a yearly report. This [data](https://gender-pay-gap.service.gov.uk/) is publicly available.The reports of 2018/19 will be analyzed since it is the last fully complete dataset available. The steps of CRISP-DM process are followed and indicated with subtitles across this notebook.**Investigation Questions:**1. How balanced are payment amounts by gender?2. What economic activities have largest gender pay gap?3. Which employer characteristics explain the gender pay gap?Apart from exploring answers to the previous questions, further exploratory analysis and alternatives are presented.  Data UnderstandingCompanies must report the following gender pay gap figures, the calculation details can be found [here](https://www.gov.uk/guidance/gender-pay-gap-reporting-make-your-calculations).| Figure                                                    | Column(s)                                                                                  ||:-----------------------------------------------------------|:------------------------------------------------------------------------------------------|| Mean gender pay gap in hourly pay                         | `DiffMeanHourlyPercent`                                                                    || Median gender pay gap in hourly pay                       | `DiffMedianHourlyPercent`                                                                  || Mean bonus gender pay gap                                 | `DiffMeanBonusPercent`                                                                     || Median bonus gender pay gap                               | `DiffMedianBonusPercent`                                                                   || Proportion of males and females receiving a bonus payment |  `(M/F)BonusPercent`                                                                || Proportion of males and females in each pay quartile      | `(M/F)LowQuartile, (M/F)LowerMiddleQuartile, (M/F)UpperMiddleQuartile, (M/FUpperQuartile)` |<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
plt.rcParams["font.family"] = "serif"
raw = pd.read_csv("data/UK-Gender-Pay-Gap-Data-2018-2019.csv")
rows, cols = raw.shape
print(rows, "x", cols)
raw.head()
print(list(raw.columns))<jupyter_output>['EmployerName', 'Address', 'CompanyNumber', 'SicCodes', 'DiffMeanHourlyPercent', 'DiffMedianHourlyPercent', 'DiffMeanBonusPercent', 'DiffMedianBonusPercent', 'MaleBonusPercent', 'FemaleBonusPercent', 'MaleLowerQuartile', 'FemaleLowerQuartile', 'MaleLowerMiddleQuartile', 'FemaleLowerMiddleQuartile', 'MaleUpperMiddleQuartile', 'FemaleUpperMiddleQuartile', 'MaleTopQuartile', 'FemaleTopQuartile', 'CompanyLinkToGPGInfo', 'ResponsiblePerson', 'EmployerSize', 'CurrentName', 'SubmittedAfterTheDeadline', 'DueDate', 'DateSubmitted']<jupyter_text>Apart from gender pay gap figures, the dataset includes company information such as name, size (`EmployerSize`), the industry sectors they belong to (`SicCodes`) and a URL to their written report (`CompanyLinkToGPGInfo`). Additionally, submission metadata is included too.<jupyter_code>columns = raw.columns
non_numeric_cols = ['EmployerName', 'Address', 'CompanyNumber', 'SicCodes',
                    'CompanyLinkToGPGInfo', 'ResponsiblePerson', 'EmployerSize',
                    'CurrentName', 'SubmittedAfterTheDeadline', 'DueDate', 'DateSubmitted']
raw[non_numeric_cols].describe()
numeric_cols = [col for col in columns if col not in non_numeric_cols]
raw[numeric_cols].describe()<jupyter_output><empty_output><jupyter_text>Bonus data have minimum and maximum outliers. This will not be fixed since bonuses are not in the scope of the business questions.<jupyter_code>raw.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 10828 entries, 0 to 10827
Data columns (total 25 columns):
 #   Column                     Non-Null Count  Dtype  
---  ------                     --------------  -----  
 0   EmployerName               10828 non-null  object 
 1   Address                    10828 non-null  object 
 2   CompanyNumber              9486 non-null   object 
 3   SicCodes                   10328 non-null  object 
 4   DiffMeanHourlyPercent      10828 non-null  float64
 5   DiffMedianHourlyPercent    10828 non-null  float64
 6   DiffMeanBonusPercent       9072 non-null   float64
 7   DiffMedianBonusPercent     9072 non-null   float64
 8   MaleBonusPercent           10828 non-null  float64
 9   FemaleBonusPercent         10828 non-null  float64
 10  MaleLowerQuartile          10828 non-null  float64
 11  FemaleLowerQuartile        10828 non-null  float64
 12  MaleLowerMiddleQuartile    10828 non-null  float64
 13  FemaleLowerMiddleQuartile  10828 non-null  flo[...]<jupyter_text>Crucial colummns do not miss values (employer size, hourly pay gaps and payment quartiles). `SicCodes` has 4.6% values missing (500/10.828), this issue will be address in the data preparation stage. **Mean vs Mean Hourly Pay Gap**Both `DiffMedianHourlyPercent` and `DiffMeanHourlyPercent` indicate the gender pay gap of the company.<jupyter_code>pay_gap_indicators = raw[['DiffMeanHourlyPercent', 'DiffMedianHourlyPercent']]
pay_gap_indicators.plot.hist(bins=100, alpha=0.5, figsize=(14,7),
                             grid=True, title="Pay gap indicators histogram")<jupyter_output><empty_output><jupyter_text>Note that negative values indicate the company has a gender pay gap that favours woman. This might be due to gender representation or discrimination.<jupyter_code>pay_gap_indicators.describe()<jupyter_output><empty_output><jupyter_text>The average of the mean curve is 14.92 and the median value of the median curve is 9.5. For further analysis the **Mean** gap will not be used because that metric is very sensitive to statistical outliers (very low of very high salaries) which might be missleading. the **Median** draws a better representation of data like raw salaries.  Data Preparation**Remove Irrelevant columns**<jupyter_code>df = raw.drop(columns=['Address','CompanyNumber', 'CompanyLinkToGPGInfo',
                       'CurrentName', 'SubmittedAfterTheDeadline',
                       'DueDate', 'DateSubmitted', 'ResponsiblePerson'
                      ])
df.shape<jupyter_output><empty_output><jupyter_text>**Company Size** The `EmployerSize` column contain 1 of 7 categories.<jupyter_code>emp_size_dist = df.EmployerSize.value_counts().sort_values()
ax = emp_size_dist.plot.barh(title='EmployerSize distribution')
ax.set_xlabel("Number of Companies")<jupyter_output><empty_output><jupyter_text>`EmployerSize` categorical texts will be replaced with the midpoint of the ranges they represent to have numerical values.<jupyter_code>to_replace     = ['Not Provided', 'Less than 250',  '250 to 499',  '500 to 999',
                  '1000 to 4999', '5000 to 19,999', '20,000 or more']
replace_values = [1, 150, 375, 750, 3000, 12500, 40000]
df.EmployerSize = df.EmployerSize.replace(to_replace, replace_values).astype(int)
df.EmployerSize.value_counts().sort_index()<jupyter_output><empty_output><jupyter_text>The difference between the initial categories (i.e. 1, 150) and the last categories (i.e. 12.500, 40.000) are significantly big. Hence, a logarithmic representation is added in order to obtain a linear representation of the sizes, this might be helfpul for the modelling phase.<jupyter_code>df['EmployerSizeLog'] = pd.Series(np.log(df.EmployerSize), index=df.index)
fig, axs = plt.subplots(1, 2, figsize=(11,4))
axs[0].set_title("EmployerSize")
axs[0].set_ylabel("Regular scale")
axs[0].set_xticks(range(8))
axs[0].set_xticklabels(to_replace, rotation=40, ha="right")
axs[0].grid()

axs[0].plot(np.sort(df.EmployerSize.unique()), marker=".")
axs[1].set_title("EmployerSizeLog")
axs[1].set_ylabel("Logarigthmic scale")
axs[1].set_xticks(range(8))
axs[1].set_xticklabels(to_replace, rotation=40, ha="right")
axs[1].grid()
axs[1].plot(range(0,7), np.sort(df.EmployerSizeLog.unique()), marker=".")<jupyter_output><empty_output><jupyter_text>**Q1: Payment Balance**A [violin plot](https://seaborn.pydata.org/generated/seaborn.violinplot.html) is suitable to represent how companies are balancing salaries based in building four quartiles. Data granularization is required to plot this.<jupyter_code>def extract_quartiles_by_gender(gender_str, df):
    """Returns a data frame only with a gender specific quartile percentages.
    """
    raw_cols = ['LowerQuartile', 'LowerMiddleQuartile', 'UpperMiddleQuartile', 'TopQuartile']
    quartile_cols = [gender_str + quartile for quartile in raw_cols]
    quartiles = df[quartile_cols].copy()
    quartiles.loc[:, 'Sex'] = gender_str
    quartiles.rename(columns={
        quartile_cols[0]: "Lower",
        quartile_cols[1]: "LowerMiddle",
        quartile_cols[2]: "UpperMiddle",
        quartile_cols[3]: "Top",
    }, inplace=True)
    return quartiles
male_quartile   = extract_quartiles_by_gender("Male", df)
female_quartile = extract_quartiles_by_gender("Female", df)
full_quartiles  = male_quartile.append(female_quartile)
print(full_quartiles.shape)
full_quartiles.sample(5)
melt_quartiles = pd.melt(full_quartiles, id_vars=['Sex'], var_name='Quartile',
    value_vars=['Lower', 'LowerMiddle', 'UpperMiddle', 'Top'],
    value_name='Percentage'
)
print(melt_quartiles.shape)
melt_quartiles.sample(5)
fig, ax = plt.subplots(figsize=(12, 5))
ax.set_title("Quartile percentage distribution by gender", fontsize=14, pad=20)
# sns.despine(left=True, bottom=True)
ax = sns.violinplot(data=melt_quartiles, x="Quartile", y="Percentage", ax=ax, hue="Sex", 
                     cut=0, split=True, inner=None, scale="count", palette="Set2")
plt.tick_params(axis='y', labelsize=12)
plt.tick_params(axis='x', labelsize=12)
plt.grid()
plt.legend(loc='lower center', fontsize=11, framealpha=1)
ax.set_xlabel(None)
ax.set_ylabel('% Density distribution' , fontsize=12)
fig.savefig("plots/quartile-percentage-dist.png", dpi=100)<jupyter_output><empty_output><jupyter_text>**Q1 Answer** A balance is found when both curves draw rainbow arcs, meaning that male and female employees are equally present in that payment quartile (smooth mode at 50%). From left to right (also better to worse), the lower quartile is the most balanced one with a slight preference on woman. The greatest imbalance is in the high wage quarile, where men take 80 to 100% of the salaries in the top quartile. **Q2 Economic activities** Sic Codes are five digit encodings used to classify companies into 21 economic sections which are subdivided into 731 economic sectors. `SicCodes` column contains a list of these codes, with some missing, clutter and incorrect values (`1`).<jupyter_code>df.sample(10).SicCodes<jupyter_output><empty_output><jupyter_text>In order to understand the industrial context of the companies, this data should be extended with human interpretable descriptions. For that, [this](https://datahub.io/core/uk-sic-2007-condensed) dataset is available in `data/sic_codes.csv`, it includes the mapping between industrial sections and sectors for company classification.<jupyter_code>sic_data = pd.read_csv('data/sic_codes.csv')
sic_data<jupyter_output><empty_output><jupyter_text>In order to manage SIC sectors and section codes, an independent module was implemented.<jupyter_code>%load_ext autoreload
import sic_codes as sc
%autoreload 1
%aimport sic_codes
df2 = sc.drop_sic_codes_na(df) # drops rows with NaN in the SicCodes column.
print(df2.shape)
df2 = sc.clean_sic_codes(df2) # Transforms cluttered SicCodes into clean lists.
df2 = sc.add_sections(df2)    # Extends the dataset with a SicSection column.
df2[["SicCodes", "SicSections"]].sample(5)<jupyter_output>(10328, 18)<jupyter_text>`SicCodes` correspond to sectors (731 elements), which are grouped in sections (21 elements). For this reason, some `SicSections` lists might be shorter than codes.<jupyter_code>def plot_section_dist(df, title, filename, save=False):
    """Plots Industry sections distribution of companies"""
    exploded_sections = sc.explode_sections(df)
    print(exploded_sections.shape)

    fig, ax = plt.subplots(1,1, figsize=(9,8))
    ax.set_title(title, fontsize='14', loc="left", pad=18)
    ax.spines['top'].set_visible(False); ax.spines['right'].set_visible(False)
    ax.spines['left'].set_visible(False); ax.spines['bottom'].set_visible(False)
    plt.tick_params(axis='both', which='both', bottom=False, top=False,
        labelbottom=False, right=False, left=False, labelleft=True, labelsize=12)

    section_bars = exploded_sections.SectDesc.value_counts().sort_values()

    section_bars.plot(kind='barh', ax=ax)

    for container in ax.containers:
        plt.setp(container, height=0.75)

    for i, v in enumerate(section_bars):
        ax.text(v + 4, i, str(v))
    if save:
        fig.savefig(f"plots/{filename}-industry-dist.png", bbox_inches="tight", dpi=100)

plot_section_dist(df2, title="Company distribution by Industry Section", filename="general")<jupyter_output>(11300, 20)<jupyter_text>The previous plot shows the distribution of all companies. The total sum of the companies is bigger thant the original rows because companies with 2 or more industrial sections have been exploded.The following sub-distributions have been extracted with a sample of 1000 companies:- Companies with largest pay Gap in favour of man- Companies with largest pay Gap in favour of woman- Companies without median pay gap<jupyter_code>sample_size = 1000
sorted_median_paygap = df2.sort_values(by=['DiffMedianHourlyPercent']).copy().reset_index()
male_biased   = sorted_median_paygap.tail(sample_size)
plot_section_dist(male_biased, "Industry section distribution of 1000 companies with largest pay gap", filename="male")<jupyter_output>(1084, 21)<jupyter_text>**Q2 Answer** Education sector conforms the 30% of the "bad" companies. The following four groups conform 10 to 12%. Three involve technical activities such as finances, insurances and construction. Also administrative roles but the reason might be because of the large general number of occurrences. The rest of the companies groups are less than 10%.<jupyter_code>female_biased = sorted_median_paygap.head(sample_size)
plot_section_dist(female_biased, title="Industry sections of 1000 companies with largest pay gap in favour of women", filename="female")
# Find how many reports claim 0 median pay gap.
zero_median_gap = sorted_median_paygap[sorted_median_paygap.DiffMedianHourlyPercent == 0]
zero_median_gap_count = zero_median_gap.shape[0]
print(zero_median_gap_count, "have zero median pay gap")

# Include companies close to zero median gap 
# to build a new index with 1000 elements.
side_pad = int((sample_size - zero_median_gap_count) / 2) + 1
balanced_companies_idx = range(zero_median_gap.index[0] - side_pad,
                               zero_median_gap.index[-1]+ side_pad)

print("Sample balanced companies index", balanced_companies_idx)
balanced = sorted_median_paygap.iloc[balanced_companies_idx] 
balanced.DiffMedianHourlyPercent.describe()
plot_section_dist(balanced, title="Industry sections of 1000 companies with less pay gap", filename="balance")<jupyter_output>(1081, 21)<jupyter_text>The services area is the section with more companies that claim having no median gap. Interestingly, Human and social oriented activities are present in high proportions in this sample of "good" companies. **Split Industry Sections**The `sc.split_sections` function distributes the list of industry codes and extends the DataFrame with 21 new columns `SectA`, `SectB` … `SectU` which indicate whether the company belongs to a particular section (1) or not (0).<jupyter_code>df3 = sc.split_sectors(df2) # Add 21 new columns for Industry sections.
print(df3.shape)
df3.head()<jupyter_output>(10328, 40)<jupyter_text>Data Modelling **Q3: Which company characteristics influence the pay gap** `DiffHourlyMedianPercent` will be modelled with Linear Regression and Random Forest algorithms in order to find the most influencial characteristics of company when estimating the pay gap.<jupyter_code>from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
pd.set_option('display.float_format', lambda x: '%.3f' % x)
features = ['MaleLowerQuartile', 'FemaleLowerQuartile',
#             'DiffMeanBonusPercent', 'DiffMedianBonusPercent',
#             'MaleBonusPercent', 'FemaleBonusPercent',
            'MaleLowerMiddleQuartile', 'FemaleLowerMiddleQuartile',
            'MaleUpperMiddleQuartile', 'FemaleUpperMiddleQuartile',
            'MaleTopQuartile', 'FemaleTopQuartile', 
            'EmployerSize', 'EmployerSizeLog',
            'SectA', 'SectB', 'SectC', 'SectD', 'SectE', 'SectF', 'SectG',
            'SectH', 'SectI', 'SectJ', 'SectK', 'SectL', 'SectM', 'SectN',
            'SectO', 'SectP', 'SectQ', 'SectR', 'SectS', 'SectT', 'SectU']

target = 'DiffMedianHourlyPercent'
X = df3[features]
y = df3[target]
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)<jupyter_output><empty_output><jupyter_text>Note: Modelling with Bonuses data does not improve the results. However, further analysis could be carried out. **Linear Regression**<jupyter_code>linear_reg = LinearRegression(normalize=True)
linear_reg.fit(X_train, y_train)

y_preds = linear_reg.predict(X_test)
score = r2_score(y_test, y_preds)
score
coeffs = pd.DataFrame({'Column': list(X_train.columns), 'Coeff': linear_reg.coef_})
coeffs['AbsCoeff'] = abs(coeffs['Coeff'])
coeffs.sort_values(by=["AbsCoeff"], ascending=False)<jupyter_output><empty_output><jupyter_text>**Random Forest**<jupyter_code>rand_for = RandomForestRegressor(random_state=1)
rand_for.fit(X_train, y_train)

y_preds = rand_for.predict(X_test)
score = r2_score(y_test, y_preds)
score
columns = list(X_train.columns)
columns[25] = "EducationSector"
coeffs2 = pd.DataFrame({'Feature': list(columns),
                        "Importance (%)": rand_for.feature_importances_ * 100})
coeffs2
coeffs2.sort_values(by=["Importance (%)"], ascending=False) \
    .reset_index().plot.barh(x="Feature", y="Importance (%)", figsize=(10,9))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1095794139907017, 'fraction_numerical': 0.02427289436695859, 'mean_word_length': 4.294504865483686, 'pattern_counts': {'":': 1, '<': 75, '<?xml version=': 0, '>': 75, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '124077', 'license': 'MIT', 'n_tokens_mistral': 5990, 'n_tokens_neox': 5396, 'n_words': 1570, 'path': 'UK-Gender-Pay-Gap.ipynb', 'repo_name': 'leanguardia/gender-pay-gap-uk', 'size': '18498'}