Table of Contents
- 1 Customer Segmentation - Introduction
- 2 Imports
- 3 Data Loading
- 4 Preprocessing
- 4.1 Assess Missing Data
- 4.2 Select and Re-Encode Features
- 4.3 Create a Cleaning Function
- 5 Feature Transformation
- 6 Clustering
Customer Segmentation - Introduction¶
In this project, I apply unsupervised learning techniques to identify segments of the population that form the core customer base for a mail-order sales company in Germany. These segments can then be used to direct marketing campaigns towards audiences that will have the highest expected rate of returns.
Imports¶
In [1]:
# import libraries here; add more as necessary
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from IPython.core.display import HTML
from scipy.stats import itemfreq, chisquare
from sklearn.preprocessing import Imputer
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
# Set base plotting style
plt.style.use('seaborn-ticks')
# # Set base plotting size
plt.rcParams['figure.figsize'] = 12, 9
# Magic word for producing visualizations in notebook
%matplotlib inline
# Increase figure resolution for high dpi screens
%config InlineBackend.figure_format = 'retina'
# Autoreload modules
%load_ext autoreload
%autoreload 2
Data Loading¶
There are four files associated with this project:
Demographics.csv: Demographics data for the general population of Germany; 891211 persons (rows) x 85 features (columns).Customers.csv: Demographics data for customers of a mail-order company; 191652 persons (rows) x 85 features (columns).Data_Dictionary.md: Detailed information file about the features in the provided datasets.Feature_Summary.csv: Summary of feature attributes for demographics data; 85 features (rows) x 4 columns
Each row of the demographics files represents a single person, but also includes information outside of individuals, including information about their household, building, and neighborhood. I will use this information to cluster the general population into groups with similar demographic properties. Then, I will see how the people in the customers dataset fit into those created clusters. The hope here is that certain clusters are over-represented in the customers data, as compared to the general population; those over-represented clusters will be assumed to be part of the core userbase. This information can then be used for further applications, such as targeting for a marketing campaign.
In [3]:
# Load in the general demographics data.
azdias = pd.read_csv('Demographics.csv', delimiter=';')
# Load in the feature summary file.
feat_info = pd.read_csv('Feature_Summary.csv', delimiter=';')
In [4]:
feat_info.head()
Out[4]:
In [5]:
azdias.head()
Out[5]:
Preprocessing¶
Assess Missing Data¶
Identify missing or unknown data values and convert them to NaNs.¶
In [6]:
df = azdias.copy()
df.head()
Out[6]:
In [7]:
df.describe()
Out[7]:
In [8]:
df.info()
In [9]:
# Percentage of nans per row before precessing
nan_perc = (df.isnull().mean() * 100).sort_values(ascending=False)
nan_perc[:20]
Out[9]:
KK_KUNDENTYP has the highest percentage of missing values (65%) with the next KBA05_ANTG1 at 14.9% at a different scale
In [10]:
#How many rows don't have nans?
nan_perc[nan_perc == 0].shape[0]
Out[10]:
In [11]:
# as a percentage
nan_perc[nan_perc > 0].shape[0]/ df.shape[0] * 100
Out[11]:
In [12]:
# How many nans before processing?
print('Total percentage of Nans before processing: ',
round(df.isnull().sum().sum() / np.product(df.shape) * 100, 2),
'%')
Convert Missing Value Codes to NaNs¶
The fourth column of the feature attributes summary (loaded in above as feat_info) documents the codes from the data dictionary that indicate missing or unknown data. While the file encodes this as a list (e.g. [-1,0]), this will get read in as a string object.
In [13]:
# find out the different kind of encodings
feat_info.missing_or_unknown.unique()
Out[13]:
In [14]:
# check for the XX type
feat_info[feat_info.missing_or_unknown == '[XX]']
Out[14]:
In [15]:
def convert_to_list_for_X_or_XX(s):
'''Takes the column with the missing value
encodings and converts to list, specially processing
the casses where X or XX is in the encodings
TODO : Generalize for any new string'''
if 'XX' in s:
if len(s) == 4:
return [s[1:-1]]
else:
return [-1, 'XX']
elif 'X' in s:
return [-1, 'X']
else:
return eval(s)
In [16]:
# Process the missing values encodings
feat_info.missing_or_unknown = feat_info.missing_or_unknown.apply(convert_to_list_for_X_or_XX)
In [17]:
# Create a dictionary to map column names to missing data encodings
missing_data_encode = dict(zip(feat_info.attribute, feat_info.missing_or_unknown))
In [18]:
# test
missing_data_encode['AGER_TYP']
Out[18]:
In [19]:
# test
missing_data_encode['CAMEO_DEU_2015']
Out[19]:
In [20]:
# Convert to nans using the missing_data_encode dictionary
for col in missing_data_encode.keys():
nans_before = df[col].isnull().sum()
df.loc[df[col].isin(missing_data_encode[col]), col] = np.nan
nans_after = df[col].isnull().sum()
if nans_after != nans_before:
print(col, 'nans before: ', nans_before)
print(col, 'nans after: ', nans_after, '\n')
Assess Missing Data in Each Column¶
- How much missing data is present in each column?
- For the remaining features, are there any patterns in which columns have, or share, missing data?
In [21]:
# How many nans after processing?
print('Total percentage of Nans after processing: ',
round(df.isnull().sum().sum() / np.product(df.shape) * 100, 2),
'%')
In [22]:
# A Series with the proportion of nans per row after precessing
nan_perc = (df.isnull().mean()).sort_values(ascending=False)
nan_perc[:20]
Out[22]:
In [23]:
#How many coluns don't have nans?
nan_perc[nan_perc == 0].shape[0]
Out[23]:
In [24]:
# as a percentage
nan_perc[nan_perc == 0].shape[0]/ df.shape[1] * 100
Out[24]:
In [25]:
def hist_box_plot(x, x_label, y_label, bin_incr):
'''Take an array as input and draw a histogram with a boxblot above it'''
f, (ax_box, ax_hist) = plt.subplots(2,
sharex=True,
gridspec_kw={
"height_ratios": (.15, .85)},
figsize=(14, 6))
sns.boxplot(x, ax=ax_box)
bins = np.arange(0, x.max() + bin_incr, bin_incr)
x.hist(grid=False, bins=bins)
ax_box.set(yticks=[])
ax_hist.set_ylabel(y_label)
ax_hist.set_xlabel(x_label)
sns.despine(ax=ax_hist)
sns.despine(ax=ax_box, left=True)
Fig1 : Histogram of missing values per Feature¶
In [26]:
hist_box_plot(nan_perc, x_label='Proportion of missing values', y_label='% of features', bin_incr=0.05);
The majority of features have less than 20% missing values. There are six features with a considerably higher number of missing values.
In [27]:
# Which are the columns with the higher NaN values?
high_nan = nan_perc[nan_perc>0.2]
high_nan
Out[27]:
Investigate patterns in the amount of missing data in each column.¶
In [28]:
# Merge feat_info with nan_perc
nan_perc_df = pd.DataFrame(nan_perc, columns=['perc_missing'])
feat_nan = (feat_info
.merge(nan_perc_df, left_on='attribute', right_index=True)
.drop(columns='missing_or_unknown', axis=1)
)
# Split the data into three categories:
# High percentage of nans (>20%)
# Moderate percentage of nans (5-20%)
# Low percentage of nans (0-5%)
# No nans
feat_nan['nan_cat'] = pd.cut(feat_nan.perc_missing,
bins=[-np.inf, 0, 0.1, 0.2, np.inf],
labels=['zero', 'low', 'moderate', 'high'],
)
feat_nan.head()
Out[28]:
Fig2 : CountPlot of missing value category per information level.¶
In [29]:
fig, ax1 = plt.subplots(figsize=(14,8))
order = ['person', 'household', 'building', 'postcode', 'community',
'microcell_rr3', 'microcell_rr4', 'region_rr1', 'macrocell_plz8']
hue_order = ['zero', 'low', 'moderate', 'high']
sns.countplot(data=feat_nan,
x='information_level',
hue='nan_cat',
order=order,
hue_order=hue_order,
ax=ax1)
plt.legend(loc='upper right')
sns.despine(ax=ax1);
Observations
- Only the
personinformational level has features with zeroNaNs. - Low
NaNfeatures (>0% - 10%) are distributed amongperson,householdandregioninformation levels. - Moderate
NaNfeatures (10% - 20%) is the only category distributed among all information levels. - High
NaNfeatures (>20%) are distributed amongperson,householdandmicrocell_rr3information levels.
Fig3 : SwarmPlot of missing value category per information level.¶
In [30]:
fig, ax1 = plt.subplots(figsize=(14,8))
order = ['person', 'household', 'building', 'postcode', 'community',
'microcell_rr3', 'microcell_rr4', 'region_rr1', 'macrocell_plz8']
hue_order = ['zero', 'low', 'moderate', 'high']
sns.swarmplot(data=feat_nan,
x='information_level',
y='perc_missing',
hue='nan_cat',
order=order,
hue_order=hue_order,
ax=ax1)
plt.legend(loc='upper right')
sns.despine(ax=ax1);
Observations
- The moderate - NaN features (10% - 20%) have all the same values for each from the information levels:
postcode,community,microcell_rr3 and 4andmacrocell_plz8. This suggests that they exist some methodical bias in creating theNaNsfor each information level. - The
householdinformation level has the most evenly distributed Nan categories. - All the information levels, except
personandhouseholdconsist mainly of moderate NaN categories. - The high NaN categories are clear outliers.
In [31]:
# Look at the similar values of microcell_rr3
feat_nan[(feat_nan.information_level == 'microcell_rr3') & (feat_nan.nan_cat == 'moderate')]
Out[31]:
Confirms the observation that features of the same attribute type (KBA05 classifies according to the number of family houses and buildings in the microcell) have the exact same percentage of missing values.
Remove the outlier columns from the dataset.¶
The outlying columns with a high percentage of NaNs will be removed.
In [32]:
cols_to_remove = feat_nan[feat_nan.nan_cat == 'high'].attribute
cols_to_remove
Out[32]:
In [33]:
df = df.drop(columns=cols_to_remove, axis=1)
df.shape
Out[33]:
Summary¶
The total percentage of Nans after the converting missing value encoding to NaNs is 11.05 %
Fig1 observations: The majority of features (79 out of 85) have less than 20% missing values. There are six features with a considerably higher number of missing values. These features are:
TITEL_KZ: Academic title flag (1: Dr., 4: Prof. Dr., 5: other)- 99.8% missing, making it unusable.
AGER_TYP: Best-ager typology (1: passive elderly, 3: experience-driven elderly)- 76.9% missing values.This could potentially be an interesting feature for identifying subgroups within older persons groups but the very high level of missing values would make the information unreliable
KK_KUNDENTYP: Consumer pattern over past 12 months (1: regular customer, 6: passive customer)- 65.6% missing values. This could also potentially be an interesting feature for identifying consumer behavior but the fact that that there are already other features that encode shopping behavior like
SHOPPER_TYPand the high level of missing values would suggest that this feature can be dropped without too much loss of information.
- 65.6% missing values. This could also potentially be an interesting feature for identifying consumer behavior but the fact that that there are already other features that encode shopping behavior like
KBA05_BAUMAX: Most common building type within the microcell (1: mainly 1-2 family homes in the microcell, 4: mainly 10+ family homes in the microcell, 5: mainly business buildings in the microcell)- 53.5% missing values. Similar information to the other
KBA05_features and the high level of missing values suggest that this feature can be dropped without too much information loss.
- 53.5% missing values. Similar information to the other
GEBURTSJAHR: Year of birth (Numeric)- 44.0% missing. The high level of missing values for a numeric variable and the fact that that there exist also other features that encode age like
ALTERSKATEGORIE_GROBsuggest that this feature can be dropped without too much information loss.
- 44.0% missing. The high level of missing values for a numeric variable and the fact that that there exist also other features that encode age like
ALTER_HH: Birthdate of head of household (1: 1895-01-01 to 1899-12-31, 21: 1995-01-01 to 1999-12-31)- 34.8% missing values: This feature could provide valuable information about to the age groupings. We will drop it in this phase due to the high level of missing values would make the information unreliable. It could in a later iteration be imputed and tested if it's inclusion leads to better results.
Fig2 observations:
- Only the
personinformational level has features with zeroNaNs. - Low
NaNfeatures (>0% - 10%) are distributed amongperson,householdandregioninformation levels. - Moderate
NaNfeatures (10% - 20%) is the only category distributed among all information levels. - High
NaNfeatures (>20%) are distributed amongperson,householdandmicrocell_rr3information levels.
Fig3 observations:
- The moderate - NaN features (10% - 20%) have all the same values for each from the information levels:
postcode,community,microcell_rr3 and 4andmacrocell_plz8.
For example all theKBA05-type features with a moderate level ofNaNshave 14.6% of missing data. This suggests that there may exist some methodical bias in creating theNaNsfor each information level. - The
householdinformation level has the most evenly distributed Nan categories. - All the information levels, except
personandhouseholdconsist mainly of moderate NaN categories. - The high NaN categories are clear outliers.
All six outlier columns with a high percentage of NaNs (>20%) will be removed.
Assess Missing Data in Each Row¶
How much data is missing in each row of the dataset?¶
Fig4 : Histogram of missing values per row¶
In [34]:
df['row_nan_perc'] = df.isnull().mean(axis=1)
hist_box_plot(df['row_nan_perc'], x_label='Proportion of missing values', y_label='% of rows', bin_incr=0.01)
- The majority of the rows is below the 10% threshhold of missing values.
- The iqr is around the 5 % threshhold.
We will use 10% as the cut of threshhold to split the data into two separate subsets and compare their distributions.
Divide the data into two subsets based on the number of missing values in each row.¶
In [35]:
threshhold = 0.1
df_lowna = df.loc[df['row_nan_perc'] < threshhold, :]
df_highna = df.loc[df['row_nan_perc'] > threshhold, :]
In [36]:
df_lowna.shape
Out[36]:
In [37]:
# percentage of low na rows
df_lowna.shape[0] / df.shape[0]
Out[37]:
In [38]:
df_highna.shape
Out[38]:
In [39]:
# percentage of high na rows
df_highna.shape[0] / df.shape[0]
Out[39]:
Compare the distribution of values for at least five columns where there are no or few missing values, between the two subsets.¶
In [40]:
# find columns with no nans at lowna
df_lowna.dropna(axis=1).shape
Out[40]:
In [41]:
# find columns with no nans at highna
df_highna.dropna(axis=1).shape
Out[41]:
In [42]:
# get the commmon columns
common_cols = list(set(df_lowna.dropna(axis=1).columns).intersection(df_highna.dropna(axis=1).columns))
common_cols
Out[42]:
In [43]:
# Drop the 'row_nan_perc'
common_cols.remove('row_nan_perc')
Fig5 : CountPlots of High and Low NaN features¶
In [44]:
# compare distributions
for col in common_cols[:6]:
fig, axes = plt.subplots(1,2, figsize=(14, 6), sharey=True)
sns.countplot(df_lowna[col], ax=axes[0], color='b')
sns.countplot(df_highna[col], ax=axes[1], color='r')
Observations
- The distributions are discrete and look distinctly different.
In [45]:
# Statistical comparison of distributions
for col in common_cols:
print(col)
print(chisquare(df_highna[col].value_counts().sort_index(),
df_lowna[col].value_counts().sort_index(),
axis=0),
'\n')
The p value is 0 for all the chisquared tests which means that we saw the expected frequencies based on the observed distribution (High NaN features) about 0% of the time which is statistically significant and confirms the assumption that the two distributions are different [5]
Summary¶
Fig4 observations:
- The majority of the rows (83.8%) is below the 10% of missing values.
We used 10% as the cut of threshhold to split the data into two separate subsets and compare their distributions taking 5 sample columns that do no have any NaN values.
Fig5 observations:
- These distributions are discrete and look distinctly different, an assumption that was confirmed by performing a Chisquared test for testing their proportions.
We will continue the analysis with the dataset that contains rows with less than 10% NaNs
Select and Re-Encode Features¶
In [46]:
# Keep only the rows with NaNs perc below 10%
df = df_lowna
# Drop the row nan percentage row
df = df[df.columns[:-1]]
df.head()
Out[46]:
In [47]:
# Drop the missing_or_unknown column from the
# processed feat_info Df
feat_info = feat_info[feat_info.columns[:-1]]
feat_info.head()
Out[47]:
In [48]:
# filter feat_info to the remaining columns
feat_info = feat_info[feat_info.attribute.isin(df.columns)]
feat_info.shape
Out[48]:
How many features are there of each data type?¶
Fig6 : CountPlot of Features per Data Type¶
In [49]:
fig, ax1 = plt.subplots(figsize=(14,8))
sns.countplot(data=feat_info,
x='type',
color='c',
ax=ax1)
ax1.set_ylabel('No of features')
ax1.set_xlabel('Data Type')
for p in ax1.patches:
ax1.annotate('{:.0f}'.format(p.get_height()), (p.get_x()+0.35, p.get_height()+1))
sns.despine(ax=ax1);
- 49 out of 79 features (62%) are ordinal.
- The categorical and mixed features, that need to be processed, represent 18 + 6 (30.38%) of the features.
Re-Encode Categorical Features¶
In [50]:
# Get categorical features
feat_info[feat_info.type == 'categorical']
Out[50]:
The majority of the categorical features is of the person information_level.
In [51]:
# How many of each information_level
feat_info[feat_info.type == 'categorical'].groupby('information_level').count().iloc[:,-1]
Out[51]:
In [52]:
# Get the categorical feature names
df_cat_cols = feat_info[feat_info.type == 'categorical'].attribute
df_cat_cols
Out[52]:
In [53]:
# Find the number of levels
df[df_cat_cols].nunique()
Out[53]:
In [54]:
# Cound the number of levels
df[df_cat_cols].nunique().sort_values()
Out[54]:
In [55]:
# how many binary
sum(df[df_cat_cols].nunique() == 2)
Out[55]:
In [56]:
# how many multilevel
sum(df[df_cat_cols].nunique() != 2)
Out[56]:
In [57]:
df[df_cat_cols].dtypes
Out[57]:
The last three columns are of dtype object. Let's have a look at them
In [58]:
df[df_cat_cols[-3:]].sample(10)
Out[58]:
OST_WEST_KZ: Building location via former East / West Germany (GDR / FRG) (O: East (GDR), W: West (FRG))- is a binary categorical variable and is emcoded with string values and and we will reencode the string catgories in numerical types so that it can be used with the algorithms.
CAMEO_DEUG_2015German CAMEO: Wealth / Life Stage Typology, rough scale (1: upper class, 9: urban working class)- is numeric but encoded as
objectdtype. Since it is a multilevel categorical feature and will be transformed to a dummy variable or dropped it can stay like this.
- is numeric but encoded as
CAMEO_DEU_2015should be afloat64data type but is encoded as object.CAMEO_DEU_2015is a more detailed versionCAMEO_DEUG_2015and would probably be highly correlated. Further investigation should be done.
Re-encode categorical variable(s) to be kept in the analysis.¶
In [59]:
# Reencode the binary OST_WEST_KZ to numeric categories
df = df.copy()
df['OST_WEST_KZ'] = df['OST_WEST_KZ'].replace({'O': 0,
'W': 1})
In [60]:
# test
df[df_cat_cols[-3:]].sample(10)
Out[60]:
In [61]:
# Choose the multilevel categorical features
cat_nunique = df[df_cat_cols].nunique()
cat_to_encode = cat_nunique[cat_nunique>2]
cat_to_encode.sort_values()
Out[61]:
In [62]:
# DROPING THE MULTI-LEVEL CATEGORICAL COLUMNS IN THIS ITERATION
df = df.drop(columns=cat_to_encode.index, axis=0)
df.info()
Summary¶
The dataset contains 18 categorical variables, 5 of which are binary and 13 are multilevel.
OST_WEST_KZ: Building location via former East / West Germany (GDR / FRG) (O: East (GDR), W: West (FRG))- is a binary categorical variable and is emcoded with string values and and we will reencode the string catgories in numerical types so that it can be used with the algorithms.
CAMEO_DEUG_2015: German CAMEO: Wealth / Life Stage Typology, rough scale (1: upper class, 9: urban working class)
The rest of the binary features will be kept as is. The multilevel categorical features need to be one hot encoded in order to be used for feature transformation and training our model. They are:
| Feature | # of Levels |
|---|---|
| NATIONALITAET_KZ | 3 |
| SHOPPER_TYP | 4 |
| LP_FAMILIE_GROB | 5 |
| LP_STATUS_GROB | 5 |
| CJT_GESAMTTYP | 6 |
| FINANZTYP | 6 |
| ZABEOTYP | 6 |
| GEBAEUDETYP | 7 |
| CAMEO_DEUG_2015 | 9 |
| LP_STATUS_FEIN | 10 |
| LP_FAMILIE_FEIN | 11 |
| GFK_URLAUBERTYP | 12 |
| CAMEO_DEU_2015 | 44 |
Since some of the information they contain is covered by other features in the dataset like CAMEO_DEUG_2015 for CAMEO_DEU_2015 or FINANZ_ for FINANZTYP LP_STATUS_FEIN and LP_STATUS_GROB for LP_FAMILIE_FEIN and LP_FAMILIE_GROB, we will choose to drop them in this iteration of the analysis and keep things more straightforward and not increase the number of variables too much at this phase.
Engineer Mixed-Type Features¶
In [63]:
# Get mixed dtype features
feat_info[feat_info.type == 'mixed']
Out[63]:
Investigate "PRAEGENDE_JUGENDJAHRE" and engineer two new variables.¶
Data Dictionary:
PRAEGENDE_JUGENDJAHRE: Dominating movement of person's youth (avantgarde vs. mainstream; east vs. west)
- -1: unknown
- 0: unknown
- 1: 40s - war years (Mainstream, E+W)
- 2: 40s - reconstruction years (Avantgarde, E+W)
- 3: 50s - economic miracle (Mainstream, E+W)
- 4: 50s - milk bar / Individualisation (Avantgarde, E+W)
- 5: 60s - economic miracle (Mainstream, E+W)
- 6: 60s - generation 68 / student protestors (Avantgarde, W)
- 7: 60s - opponents to the building of the Wall (Avantgarde, E)
- 8: 70s - family orientation (Mainstream, E+W)
- 9: 70s - peace movement (Avantgarde, E+W)
- 10: 80s - Generation Golf (Mainstream, W)
- 11: 80s - ecological awareness (Avantgarde, W)
- 12: 80s - FDJ / communist party youth organisation (Mainstream, E)
- 13: 80s - Swords into ploughshares (Avantgarde, E)
- 14: 90s - digital media kids (Mainstream, E+W)
- 15: 90s - ecological awareness (Avantgarde, E+W)
In [64]:
# Let's have a look at it
df.PRAEGENDE_JUGENDJAHRE.sample(10)
Out[64]:
In [65]:
# Let's have a look at it
df.PRAEGENDE_JUGENDJAHRE.nunique()
Out[65]:
Fig7 : CountPlot of PRAEGENDE_JUGENDJAHRE categories¶
In [66]:
plt.style.use('ggplot')
fig, ax1 = plt.subplots(figsize=(14,8))
sns.countplot(data=df,
x='PRAEGENDE_JUGENDJAHRE',
color="darkslategray",
ax=ax1);
In [67]:
df.PRAEGENDE_JUGENDJAHRE.value_counts()
Out[67]:
In [68]:
df.PRAEGENDE_JUGENDJAHRE.isnull().sum()
Out[68]:
In [69]:
df.PRAEGENDE_JUGENDJAHRE.isnull().mean() * 100
Out[69]:
Observations
Level 14 (90s - digital media kids (Mainstream, E+W))leads with close to 172952 rows followed bylevel 8 (70s - family orientation (Mainstream, E+W)) with 134165 observationsLevels 7 (60s - opponents to the building of the Wall (Avantgarde, E)), 13 ((80s - Swords into ploughshares (Avantgarde, E))and2 (40s - reconstruction years (Avantgarde, E+W))have the least observations with 3888, 5300 and 7340 rows respectively- The feature has 16121
NaNs- 2.16 % of the total.
New Variables Data encoding:
ENG_DECADE:
Person’s decade of youth.
- 1: 40s
- 2: 50s
- 3: 60s
- 4: 70s
- 5: 80s
- 6: 90s
ENG_MOVEMENT:
Person’s movement alignment.
- 1: Mainstream
- 2: Avantgarde
In [70]:
decade_dict = {1: 1, 2: 1, 3: 2, 4: 2, 5: 3, 6: 3, 7: 3,
8: 4, 9: 4, 10: 5, 11: 5, 12: 5, 13: 5,
14: 6, 15: 6
}
movement_dict = {1: 1, 3: 1, 5: 1, 8: 1, 10: 1, 12: 1, 14: 1,
2: 2, 4: 2, 6: 2, 7: 2, 9: 2, 11: 2, 13: 2, 15: 2
}
# It appears that map is approximately 10x faster than replace.
# Ref: https://stackoverflow.com/questions/20250771/remap-values-in-pandas-column-with-a-dict
df['ENG_DECADE'] = df['PRAEGENDE_JUGENDJAHRE'].map(decade_dict)
df['ENG_MOVEMENT'] = df['PRAEGENDE_JUGENDJAHRE'].map(movement_dict)
Fig8 : CountPlot of ENG_DECADE categories¶
In [71]:
fig, ax1 = plt.subplots(figsize=(14,8))
sns.countplot(data=df,
x='ENG_DECADE',
palette="Greens",
ax=ax1);
In [72]:
df.ENG_DECADE.value_counts()
Out[72]:
In [73]:
# Sanity check
df.ENG_DECADE.isnull().sum()
Out[73]:
Observations
- In the new feature
Level 6 90sleads with close to 212488 rows followed bylevel 4 70swith 166658 observations consistently with the original feature. Levels 1: 40s, and2 50shave the least observations with 26923 and 71852 rows respectively.- The feature as expected still has 16121
NaNs- 2.16 % of the feature's total values.
Fig9 : CountPlot of ENG_MOVEMENT categories¶
In [74]:
fig, ax1 = plt.subplots(figsize=(14,8))
sns.countplot(data=df,
x='ENG_MOVEMENT',
palette="Greens_r",
ax=ax1);
In [75]:
df.ENG_MOVEMENT.value_counts()
Out[75]:
In [76]:
# Sanity check
df.ENG_MOVEMENT.isnull().sum()
Out[76]:
Observations
- In this new binary feature
1: Mainstreamleads with close to 564070 rows almost 3.5 times more than2: Avantgardewith 166918 observations. - The feature as expected still has 16121
NaNs- 2.16 % of the total.
Investigate "CAMEO_INTL_2015" and engineer two new variables.¶
Data Dictionary:
CAMEO_INTL_2015: Wealth / Life Stage Typology, mapped to international code
- -1: unknown
- 11: Wealthy Households - Pre-Family Couples & Singles
- 12: Wealthy Households - Young Couples With Children
- 13: Wealthy Households - Families With School Age Children
- 14: Wealthy Households - Older Families & Mature Couples
- 15: Wealthy Households - Elders In Retirement
- 21: Prosperous Households - Pre-Family Couples & Singles
- 22: Prosperous Households - Young Couples With Children
- 23: Prosperous Households - Families With School Age Children
- 24: Prosperous Households - Older Families & Mature Couples
- 25: Prosperous Households - Elders In Retirement
- 31: Comfortable Households - Pre-Family Couples & Singles
- 32: Comfortable Households - Young Couples With Children
- 33: Comfortable Households - Families With School Age Children
- 34: Comfortable Households - Older Families & Mature Couples
- 35: Comfortable Households - Elders In Retirement
- 41: Less Affluent Households - Pre-Family Couples & Singles
- 42: Less Affluent Households - Young Couples With Children
- 43: Less Affluent Households - Families With School Age Children
- 44: Less Affluent Households - Older Families & Mature Couples
- 45: Less Affluent Households - Elders In Retirement
- 51: Poorer Households - Pre-Family Couples & Singles
- 52: Poorer Households - Young Couples With Children
- 53: Poorer Households - Families With School Age Children
- 54: Poorer Households - Older Families & Mature Couples
- 55: Poorer Households - Elders In Retirement
- XX: unknown
In [77]:
# Let's have a look at it
df.CAMEO_INTL_2015.head(10)
Out[77]:
Fig10 : CountPlot of ENG_MOVEMENT categories¶
In [78]:
fig, ax1 = plt.subplots(figsize=(14,8))
levels_sort = [str(int(l)) for l in np.sort(df.CAMEO_INTL_2015.astype(float).unique())[:-1]]
color_map = dict(zip())
sns.countplot(data=df,
x='CAMEO_INTL_2015',
color = 'darkslategray',
order = levels_sort,
ax=ax1);
In [79]:
df.CAMEO_INTL_2015.value_counts()
Out[79]:
In [80]:
df.CAMEO_INTL_2015.isnull().sum()
Out[80]:
In [81]:
df.CAMEO_INTL_2015.isnull().mean() * 100
Out[81]:
Observations
- Level
51: Poorer Households - Pre-Family Couples & Singlesleads with 128033 rows followed by level41: Less Affluent Households - Pre-Family Couples & Singles) with 87902 observations, almost 1/3 less. - Levels
33: Comfortable Households - Families With School Age Children,32: Comfortable Households - Young Couples With Childrenand35 (40s - reconstruction years (Avantgarde, E+W))have the least observations with 9161, 9777 and 9882 rows respectively - The feature has 3231
NaNs- 4.32% of it's total.
"CAMEO_INTL_2015" combines information on two axes: wealth and life stage. Break up the two-digit codes by their 'tens'-place and 'ones'-place digits into two new ordinal variables (which, for the purposes of this project, is equivalent to just treating them as their raw numeric values)
New Variables Data encoding:
ENG_WEALTH:
Household's wealth
- 1: Wealthy Households
- 2: Prosperous Households
- 3: Comfortable Households
- 4: Less Affluent Households
- 5: Poorer Households
ENG_LIFE_STAGE:
Person's Life Stage
- 1: Pre-Family Couples & Singles
- 2: Young Couples With Children
- 3: Families With School Age Children
- 4: Older Families & Mature Couples
- 5: Elders In Retirement
In [82]:
def breakup_codes(x, digit):
''' Breaks up two-digit codes by their
'tens'-place and 'ones'-place digits
into two new ordinal variables.
Leaves NaNs unchanged'''
if not pd.isna(x):
if digit == 'first':
return int(str(x)[0])
elif digit == 'second':
return int(str(x)[1])
return x
In [83]:
df['ENG_WEALTH'] = df['CAMEO_INTL_2015'].apply(breakup_codes, digit='first')
df['ENG_LIFE_STAGE'] = df['CAMEO_INTL_2015'].apply(breakup_codes, digit='second')
Fig11 : CountPlot of ENG_WEALTH categories¶
In [84]:
fig, ax1 = plt.subplots(figsize=(14,8))
sns.countplot(data=df,
x='ENG_WEALTH',
color = 'teal',
ax=ax1);
In [85]:
df.ENG_WEALTH.value_counts()
Out[85]:
In [86]:
# Sanity check
df.ENG_WEALTH.isnull().sum()
Out[86]:
Observations
- In the new feature Level
5: Poorer Householdsleads with close to 213978 rows followed by level4: Less Affluent Householdswith 181187 observations consistently with the original feature. - Levels
3: Comfortable Households, and1: Wealthy Householdshave the least observations with 62642 and 111831 rows respectively.
Fig12 : CountPlot of ENG_LIFE_STAGE categories¶
In [87]:
fig, ax1 = plt.subplots(figsize=(14,8))
sns.countplot(data=df,
x='ENG_LIFE_STAGE',
color = 'c',
ax=ax1);
In [88]:
df.ENG_LIFE_STAGE.value_counts()
Out[88]:
In [89]:
# Sanity check
df.ENG_LIFE_STAGE.isnull().sum()
Out[89]:
In [90]:
mixed_feat_cols = feat_info[feat_info.type == 'mixed'].attribute
mixed_feat_cols
Out[90]:
In [91]:
# How many levels per feature
df[mixed_feat_cols].nunique()
Out[91]:
In [92]:
# Drop mixed dtype features
df = df.drop(mixed_feat_cols, axis=1)
df.head()
Out[92]:
Observations
- In the new feature Level
1: Pre-Family Couples & Singlesleads with 232220 rows followed by level4: Older Families & Mature Coupleswith 220077 observations. - Levels
2: Young Couples With Children, and3: Families With School Age Childrenhave the least observations with 71973 and 108061 rows respectively.
Summary¶
The dataset contains 6 mixed variables, 5 of which are binary and 13 are multilevel.
The mixed feature PRAEGENDE_JUGENDJAHRE: Dominating movement of person's youth (avantgarde vs. mainstream; east vs. west) (1: 40s - war years (Mainstream, E+W), 15: 90s - ecological awareness (Avantgarde, E+W)), with the following observations Fig7:
- Level
14 (90s - digital media kids (Mainstream, E+W))leads with close to 172952 rows followed by level8 (70s - family orientation (Mainstream, E+W)) with 134165 observations - Levels
7 (60s - opponents to the building of the Wall (Avantgarde, E)),3 ((80s - Swords into ploughshares (Avantgarde, E))and2 (40s - reconstruction years (Avantgarde, E+W))have the least observations with 3888, 5300 and 7340 rows respectively - The feature has 16121
NaNs- 2.16 % of it's total.
and two new variables have been engineered:
ENG_DECADE: Person’s decade of youth. (1: 40s , 6: 90s) , Fig8 observations:- In the new feature
Level 6 90sleads with close to 212488 rows followed bylevel 4 70swith 166658 observations consistently with the original feature. Levels 1: 40s, and2 50shave the least observations with 26923 and 71852 rows respectively.
- In the new feature
ENG_MOVEMENT: Person’s movement alignment (1: Mainstream, 2: Avantgarde ), Fig9 observations:- In this new binary feature
1: Mainstreamleads with close to 564070 rows almost 3.5 times more than2: Avantgardewith 166918 observations. - The feature as expected still has 16121
NaNs- 2.16 % of the total.
- In this new binary feature
The mixed feature CAMEO_INTL_2015: German CAMEO: Wealth / Life Stage Typology, mapped to international code (11: Wealthy Households - Pre-Family Couples & Singles, 55: Poorer Households - Elders In Retirement), has been explored, with the following observations Fig10:
- Level
51: Poorer Households - Pre-Family Couples & Singlesleads with 128033 rows followed by level41: Less Affluent Households - Pre-Family Couples & Singles) with 87902 observations, almost 1/3 less. - Levels
33: Comfortable Households - Families With School Age Children,32: Comfortable Households - Young Couples With Childrenand35 (40s - reconstruction years (Avantgarde, E+W))have the least observations with 9161, 9777 and 9882 rows respectively - The feature has 3231
NaNs- 4.32% of it's total.
and two new variables have been engineered:
ENG_WEALTH: Household's wealth (1: Wealthy Households, 5: Poorer Households) , Fig11 observations:- In the new feature Level
5: Poorer Householdsleads with close to 213978 rows followed by level4: Less Affluent Householdswith 181187 observations consistently with the original feature. - Levels
3: Comfortable Households, and1: Wealthy Householdshave the least observations with 62642 and 111831 rows respectively.
- In the new feature Level
ENG_MOVEMENT: Person's Life Stage (Pre-Family Couples & Singles, 5: Elders In Retirement), Fig12 observations:- In the new feature Level
1: Pre-Family Couples & Singlesleads with 232220 rows followed by level4: Older Families & Mature Coupleswith 220077 observations. - Levels
2: Young Couples With Children, and3: Families With School Age Childrenhave the least observations with 71973 and 108061 rows respectively.
- In the new feature Level
The rest of the mixed features are:
| Feature | # of Levels |
|---|---|
| LP_LEBENSPHASE_FEIN | 40 |
| LP_LEBENSPHASE_GROB | 12 |
| WOHNLAGE | 8 |
| PLZ8_BAUMAX | 5 |
which we chose to drop in this iteration of the analysis and keep things more straightforward and not increase the number of variables too much at this phase.
Complete Feature Selection¶
If there are other re-engineering tasks you need to perform, make sure you take care of them here. (Dealing with missing data will come in step 2.1.)¶
Do whatever you need to in order to ensure that the dataframe only contains the columns that should be passed to the algorithm functions.¶
In [93]:
df.columns
Out[93]:
In [94]:
df.shape
Out[94]:
In [95]:
# Check if all colums are of numeric dtype
from pandas.api.types import is_numeric_dtype
#should be equal to the no of cols if all True
assert(sum([is_numeric_dtype(df[col]) for col in df.columns]) == df.shape[1])
print('Tests OK')
Create a Cleaning Function¶
In [96]:
def clean_data(df):
"""
Perform feature trimming, re-encoding, and engineering for demographics
data
INPUT: Demographics DataFrame
OUTPUT: Trimmed and cleaned demographics DataFrame
"""
# Put in code here to execute all main cleaning steps:
# convert missing value codes into NaNs, ...
df = convert_missingdata_encodings_to_nans(df)
print('Convert missing value codes into NaN, OK')
# remove selected columns , ...
df = drop_selected_columns(df)
# remove selected rows and return both dfs
df, df_cust_highna = drop_selected_rows(df)
print('Remove selected columns and rows, OK')
# select, re-encode, and engineer column values.
df = reencode_categorical_feat(df)
df = engineer_mixed_feat(df)
print('Select, re-encode, and engineer column value, OK')
# Return the cleaned dataframe and the high NaN Df.
return df, df_cust_highna
def convert_missingdata_encodings_to_nans(df, debug=False):
'''Convert to nans using the missing_data_encode dictionary'''
for col in missing_data_encode.keys():
if debug:
print(col, 'nans before: ', df[col].isnull().sum())
df.loc[df[col].isin(missing_data_encode[col]), col] = np.nan
if debug:
print(col, 'nans after: ', df[col].isnull().sum(), '\n')
return df
def drop_selected_columns(df):
''''''
return df.drop(columns=cols_to_remove, axis=1)
def drop_selected_rows(df):
'''split the dataset into two parts:
One with a percntage of NaN per row over the threshhold
and one below. Return both.'''
threshhold = 0.1
df['row_nan_perc'] = df.isnull().mean(axis=1)
df_lowna = df.loc[df['row_nan_perc'] < threshhold, :].\
drop(columns='row_nan_perc', axis=1)
df_highna = df.loc[df['row_nan_perc'] > threshhold, :].\
drop(columns='row_nan_perc', axis=1)
return df_lowna, df_highna
def reencode_categorical_feat(df):
''''''
# Reencode the binary OST_WEST_KZ to numeric categories
df['OST_WEST_KZ'] = df['OST_WEST_KZ'].replace({'O': 0,
'W': 1})
# choose the categorical features with more than 2 levels
df_cat_cols = feat_info[feat_info.type == 'categorical'].attribute
cat_nunique = df[df_cat_cols].nunique()
cat_to_encode = cat_nunique[cat_nunique > 2]
# DROP THE MULTI-LEVEL CATEGORICAL FEATURES
df = df.drop(columns=cat_to_encode.index, axis=0)
return df
def engineer_mixed_feat(df):
'''Engineer new variables from the
mixed features. Drop he old features'''
# Engineer two new variables from 'PRAEGENDE_JUGENDJAHRE'
decade_dict = {1: 1, 2: 1, 3: 2, 4: 2, 5: 3, 6: 3, 7: 3,
8: 4, 9: 4, 10: 5, 11: 5, 12: 5, 13: 5,
14: 6, 15: 6
}
movement_dict = {1: 1, 3: 1, 5: 1, 8: 1, 10: 1, 12: 1, 14: 1,
2: 2, 4: 2, 6: 2, 7: 2, 9: 2, 11: 2, 13: 2, 15: 2
}
df['ENG_DECADE'] = df['PRAEGENDE_JUGENDJAHRE'].map(decade_dict)
df['ENG_MOVEMENT'] = df['PRAEGENDE_JUGENDJAHRE'].map(movement_dict)
df.drop
# Engineer two new variables from 'CAMEO_INTL_2015'
df['ENG_WEALTH'] = df['CAMEO_INTL_2015'].apply(
breakup_codes, digit='first')
df['ENG_LIFE_STAGE'] = df['CAMEO_INTL_2015'].apply(
breakup_codes, digit='second')
# Drop mixed dtype features
mixed_feat_cols = feat_info[feat_info.type == 'mixed'].attribute
df = df.drop(mixed_feat_cols, axis=1)
return df
In [97]:
def test_clean_data(df):
'''A simple test function to test that the cleaning function
is at least producing an identical Dataframe as the one we created
TODO: proper unit tests'''
df_test = pd.read_csv('Udacity_AZDIAS_Subset.csv', delimiter=';')
print('Load the dataset, OK')
pd.testing.assert_frame_equal(clean_data(df_test)[0], df)
print('Tests Passed!')
test_clean_data(df)
Feature Transformation¶
Apply Feature Scaling¶
Explore and clean the dataset of all NaN values¶
In [98]:
# Percentage of nans per column
(df.isnull().mean(axis=0) * 100).sort_values(ascending=False)[:15]
Out[98]:
W_KEIT_KIND_HH has the highest percentage of missing values (6%) with the next REGIOTYP and KKK at 5.4%.
In [99]:
# Percentage of nans per row
nan_perc = (df.isnull().mean(axis=1) * 100).sort_values(ascending=False)
#How many rows have nans?
nan_perc[nan_perc > 0].shape[0]
Out[99]:
In [100]:
# as a percentage of rows
print('Total percentage of row that have NaNs : ',
round(nan_perc[nan_perc > 0].shape[0]/ df.shape[0] * 100, 4),
'%')
In [101]:
df.isnull().sum().sum()
Out[101]:
In [102]:
np.product(df.shape)
Out[102]:
In [103]:
# How many nans as perc?
print('Total percentage of Nans : ',
round(df.isnull().sum().sum() / np.product(df.shape) * 100, 2),
'%')
In [104]:
# Fit the Scaler without NaNs
df_drop = df.dropna(axis=0)
scaler = StandardScaler()
scaler.fit(df_drop)
Out[104]:
In [105]:
# Imput missing values in the full df with the mean
imp = Imputer(missing_values='NaN', strategy='mean', axis=0)
X_impute = imp.fit_transform(df)
In [106]:
X_impute
Out[106]:
In [107]:
# Transform the df with the fitted scaler
X_scale = scaler.transform(X_impute)
In [108]:
# Check the distribution of the feature with the highest pec of NaNs
df_scale = pd.DataFrame(X_scale, columns=df.columns)
df_scale.W_KEIT_KIND_HH.describe()
Out[108]:
In [109]:
# check the distributions
df_scale.describe()
Out[109]:
In [110]:
# How many nans as perc after?
print('Total percentage of Nans after drop: ',
round(df_drop.isnull().sum().sum() / np.product(df_drop.shape) * 100, 2),
'%')
Summary¶
The total percentage of NaNs in the dataset is 0.67% while the percentage of rows that contain NanNs is 16,37%.
Droping them would discard a significant amount of data in order to handle much smaller amount of actually missing data.
Imputing them with the mean will result them having a 0 value after using StandardScaler(), possibly reducing the variance of the resulting variables.
An alternative which was used is to fit the StandardScaler to the data without the NaNs and then reintroduce them, imputing them with the mean and transforming them with the fitted StandardScaler object thus resulting in a distributions that retain a higher amount of information.
Dimensionality Reduction¶
Investigate the variance accounted for by each principal component.¶
In [111]:
X_scale.shape
Out[111]:
In [112]:
# With all the feaures
pca = PCA().fit(X_scale)
Fig13 : Explained Variance per Principal Component¶
In [113]:
num_components=len(pca.explained_variance_ratio_)
inds = np.arange(num_components)
vals = pca.explained_variance_ratio_
fig, ax = plt.subplots(figsize=(12,9))
plt.bar(inds, vals, color='b');
Fig14 : Explained Cumulative Variance per Principal Component¶
In [114]:
fig, ax = plt.subplots(figsize=(16,10))
plt.plot(np.cumsum(pca.explained_variance_ratio_), marker='o', color='b')
ax.set_xticks(inds+1, minor=False);
The numbers that stand out are 12, 18 and 25.
In [115]:
def scree_plot(pca, c, ann_f=6, figsize=(12, 9), pct_change=False):
'''
Creates a scree plot associated with the principal components
INPUT: pca - the result of instantian of PCA in scikit learn
OUTPUT:
None
'''
num_components = len(pca.explained_variance_ratio_)
ind = np.arange(num_components) + 1
if pct_change:
vals = -pd.Series(pca.explained_variance_ratio_).pct_change()
# the first item is NaN
vals = vals[1:].tolist()
ind = ind[1:]
else:
vals = pca.explained_variance_ratio_
plt.figure(figsize=figsize)
ax = plt. subplot(111)
cumvals = np.cumsum(vals)
ax.bar(ind, vals, color=c)
ax.plot(ind, cumvals, color=c)
for i in range(num_components):
try: # Bandage aid - need to fix
ax.annotate(r"%s%%" % (
(str(vals[i] * 100)[:4])),
(ind[i] + 0.2, vals[i]),
va="bottom", ha="center", fontsize=ann_f)
except IndexError:
pass
ax.xaxis.set_tick_params(width=0)
ax.set_xticks(ind, minor=False)
ax.yaxis.set_tick_params(width=2, length=12)
ax.set_xlabel("Principal Component")
ax.set_ylabel("Variance Explained (%)")
plt.title('Explained Variance Per Principal Component')
Fig15 : Combined Plot of cumulative Explained Variance and Cumulative explained Variance and per Principal Component (First 25)¶
In [116]:
# Re-apply PCA to the data while selecting for number of components to retain.
pca_25 = PCA(n_components=25)
X_pca_25 = pca_25.fit(X_scale)
scree_plot(pca_25, c='slategray', ann_f=8, figsize=(12,9))
Observations
There seems to be a stabilization of the variance explained and a decrease in the ratio of cumulative variance explained round 12 components.
In [117]:
np.cumsum(pca.explained_variance_ratio_[:12])[-1]
Out[117]:
The cumulative explained variance of the first 12 components is 66.21%.
Final dimensionality reduction¶
In [118]:
# Fit on X for 12 components
pca = PCA(n_components=12)
X_pca = pca.fit_transform(X_scale)
In [119]:
X_pca.shape
Out[119]:
In [120]:
df_pca = pd.DataFrame(X_pca)
df_pca.head()
Out[120]:
In [121]:
df_pca.describe()
Out[121]:
Summary¶
We applied the dimensionality reduction PCA() method to the scaled dataset and by observing the combined plot of explained variance and cumulative explained variance and per principal component cumulative explained variance plot (Fig15) we chose to keep 12 components as at this number the seems to start a decrease in the ratio of cumulative variance explained as new components are added.
The cumulative explained variance of the first 12 components is 66.21%.
Interpret Principal Components¶
In [122]:
def map_pca_weights_to_feats(pca, df, comp_no):
'''Map pca weights to individual features
and return two pd.Series on with the highest
positive weights and one with the lowest negative
weights'''
weights = pd.DataFrame(np.round(pca.components_, 4), columns=df.keys())
component = weights.iloc[comp_no - 1, :]
comp_pos = component[component > 0].sort_values(ascending=False)
comp_neg = component[component < 0].sort_values(ascending=True)
return comp_pos, comp_neg
In [123]:
def create_bar_table(comp_pos, comp_neg):
'''Create and display a conditionally styled pandas dataframe
for the interpertation of the positive and negative weights
of PCA and centroid distances for KMeans'''
head = """
<table>
"""
row = ""
for serie in [comp_pos[:15],comp_neg[:15]]:
s = serie.copy()
s.name=''
row += "<td>{}</td>".format(s.to_frame().style.bar(
align='mid',
color=['#d65f5f', '#5fba7d'],
width=100).render()
)
row += '</tr>'
head += row
head+= """
</table>"""
display(HTML(head))
Map weights for the first principal component to corresponding feature names.¶
In [124]:
# the weights df
weights = pd.DataFrame(np.round(pca.components_, 4), columns = df_drop.keys())
weights.head()
Out[124]:
In [125]:
pd.Series(X_pca[:, 1]).describe()
Out[125]:
Fig16 : Table-Bar Plot of the 1st Principal Component Weights¶
In [126]:
comp_pos1, comp_neg1 = map_pca_weights_to_feats(pca, df_drop, 1)
create_bar_table(comp_pos1, comp_neg1)
Most Positive
- PLZ8_ANTG3: Number of 6-10 family houses in the PLZ8 region (0: no 6-10 family homes, 3: high share of 6-10 family homes)
- Higher share of medium apartment buildings
- PLZ8_ANTG4 : Number of 10+ family houses in the PLZ8 region (0: no 10+ family home, 2: high share of 10+ family homes)
- Higher share of larger apartment buildings
- HH_EINKOMMEN_SCORE: Estimated household net income (1: highest income, 6: very low income)
- Lower incomes
- ENG_WEALTH: Household wealth - Engineered (1: Wealthy Households, 5: Poorer Households)
- Poorer Households
- ORTSGR_KLS9: Size of community (1: <= 2,000 inhabitants, 9: > 700,000 inhabitants)
- Bigger community size
- EWDICHTE : Density of households per square kilometer (1: less than 34 households per km^2, 6: more than 999 households per km^2)
- Denser populated areas
Most Negative
- MOBI_REGIO : Movement patterns (1: very high movement, 6: none)
- Lower movement patterns
- PLZ8_ANTG1: Number of 1-2 family houses in the PLZ8 region (0: no 1-2 family homes, 4: very high share of 1-2 family homes)
- Higher share of of 1-2 family homes
- KBA05_ANTG1 : Number of 1-2 family houses in the microcell (0: no 1-2 family homes, 4: very high share of 1-2 family homes)
- Higher share of of 1-2 family homes
- FINANZ_MINIMALIST : Financial typology - MINIMALIST: low financial interest (1: very high, 5: very low)
- Lower financial minimalism
The first component component seems to capture the population density, and financial affluence of the person
PLZ8_ANTG3, PLZ8_ANTG4 , ENG_WEALTH, HH_EINKOMMEN_SCORE, RTSGR_KLS9 and EWDICHTE have high positive weights (People that live in more densely populated areas and with lower financial affluence).
The assumption is supported by the negative weights of MOBI_REGIO , KBA05_ANTG1, PLZ8_ANTG1 and FINANZ_MINIMALIST (People that live in less densely populated areas and with higher financial affluence).
Map weights for the second principal component to corresponding feature names.¶
Fig17 : Table-Bar Plot of the 2nd Principal Component Weights¶
In [127]:
comp_pos2, comp_neg2 = map_pca_weights_to_feats(pca, df_drop, 2)
create_bar_table(comp_pos2, comp_neg2)
Most Positive
- ALTERSKATEGORIE_GROB: Estimated age based on given name analysis (1: < 30 years old, 4: > 60 years old, 9: uniformly distributed)
- Older or evenly distributed
- FINANZ_VORSORGER: Financial typology - Be prepared (1: Very high, 5: Very low)
- Lower financial preparation
- SEMIO_ERL: Personality typology - event-oriented ( 1: highest affinity, 7: lowest affinity)
- Lower event oriented affinity
- SEMIO_LUST: Personality typology - sensual-minded ( 1: highest affinity, 7: lowest affinity)
- Lower sensual-minded affinity
- RETOURTYP_BK_S: Return type (1: influenceable Crazy-Shopper, 5: determined Minimal-Returner)
- More conservative return type
- FINANZ_HAUSBAUER: Financial typology - home ownership (1: Very high, 5: Very low)
- Lower home ownership
Most Negative
- SEMIO_REL: Personality typology - religious ( 1: highest affinity, 7: lowest affinity)
- Lower religious affinity
- ENG_DECADE: Person’s decade of youth (1: 40s , 6: 90s)
- Younger generation
- FINANZ_SPARER: Financial typology - money-saver (1: Very high, 5: Very low)
- Lower money saving
- FINANZ_UNAUFFAELLIGER: Financial typology - inconspicuous (1: Very high, 5: Very low)
- Less inconspicuous
- SEMIO_TRADV: Personality typology - tradional-minded ( 1: highest affinity, 7: lowest affinity)
- Less traditional minded
The second component seems to capture the age, generation and culture (financial, buying and other) of the person.
RETOURTYP_BK_S, SEMIO_ERL and SEMIO_LUST ALTERSKATEGORIE_GROB and FINANZ_VORSORGER have high positive weights and (older people with lower financial preparation and home ownership and less sensual and event oriented).
The assumption is supported by the negative weights of SEMIO_REL, ENG_DECADE, FINANZ_SPARER and SEMIO_TRADV (younger generation, with lower religious affinity, lower money saving, less inconspicuous and less traditional minded).
Map weights for the third principal component to corresponding feature names.¶
Fig18 : Table-Bar Plot of the 3rd Principal Component Weights¶
In [128]:
comp_pos3, comp_neg3 = map_pca_weights_to_feats(pca, df_drop, 3)
create_bar_table(comp_pos3, comp_neg3)
Most Positive
- SEMIO_VERT: Personality typology - VERT: dreamful (1: highest affinity, 7: lowest affinity)
- Less dreamful
- SEMIO_SOZ: Personality typology - SOZ: socially-minded (1: highest affinity, 7: lowest affinity)
- Less socially-minded
- SEMIO_FAM: Personality typology - family-minded ( 1: highest affinity, 7: lowest affinity)
- Less family-minded
- SEMIO_KULT: Personality typology - cultural-minded ( 1: highest affinity, 7: lowest affinity)
- Less cultural-minded
- FINANZ_MINIMALIST : Financial typology - MINIMALIST: low financial interest (1: very high, 5: very low)
- Lower financial minimalism
- RETOURTYP_BK_S: Return type (1: influenceable Crazy-Shopper, 5: determined Minimal-Returner)
- More conservative return type
Most Negative
- ANREDE_KZ: Gender ( 1: male, 2: female)
- Female
- SEMIO_KAEM: Personality typology - KAEM: combative attitude ( 1: highest affinity, 7: lowest affinity)
- Less combative attitude
- SEMIO_DOM: Personality typology - DOM: dominant-minded ( 1: highest affinity, 7: lowest affinity)
- Less dominant-minded
- SEMIO_KRIT: Personality typology - KRIT: critical-minded ( 1: highest affinity, 7: lowest affinity)
- Less critical-minded
- SEMIO_RAT: Personality typology - RAT: rational ( 1: highest affinity, 7: lowest affinity)
- Less rational
- FINANZ_ANLEGER : Financial typology - ANLEGER: investor (1: very high, 5: very low)
- Less prone to investing
The third component seems to capture the gender and their main corresponding personality types (dreamfulness, dominance, rationality, social and family attitude) of the person.
SEMIO_VERT, SEMIO_SOZ and SEMIO_FAM, SEMIO_KULT and FINANZ_MINIMALIST and RETOURTYP_BK_S have high positive weights and in conjunction with the very strong negative weight of ANREDE_KZ can refer to (Men that are less dreamful, less socially-minded, less cultural-minded not financial minimalists and conservative return shopper types).
The assumption is supported by the negative weights of ANREDE_KZ, SEMIO_KAEM, SEMIO_DOM, SEMIO_KRIT, SEMIO_RAT and FINANZ_ANLEGER (Women that are less dominant-minded, less critical-minded, less rational have less combative attitude and are less prone to investing).
Summary¶
1.The first component component seems to capture the population density, and financial status of the person - Fig16.
PLZ8_ANTG3, PLZ8_ANTG4 , ENG_WEALTH, HH_EINKOMMEN_SCORE, RTSGR_KLS9 and EWDICHTE have high positive weights correlated with people that live in more densely populated areas and with lower financial affluence.
The assumption is supported by the negative weights of MOBI_REGIO , KBA05_ANTG1, PLZ8_ANTG1 and FINANZ_MINIMALIST correlated with People that live in less densely populated areas and with higher financial affluence.
2.The second component seems to capture the age, generation and culture (financial, buying and other) of the person - Fig17.
RETOURTYP_BK_S, SEMIO_ERL and SEMIO_LUST ALTERSKATEGORIE_GROB and FINANZ_VORSORGER have high positive weights correlated with older people with lower financial preparation and home ownership that are less sensual and event oriented.
The assumption is supported by the negative weights of SEMIO_REL, ENG_DECADE, FINANZ_SPARER and SEMIO_TRADV correlated with (younger generation, with lower religious affinity, lower money saving, less inconspicuous and less traditional minded).
3.The third component seems to capture the gender and their main corresponding personality types (dreamfulness, dominance, rationality, social and family attitude) of the person - Fig18.
SEMIO_VERT, SEMIO_SOZ and SEMIO_FAM, SEMIO_KULT and FINANZ_MINIMALIST and RETOURTYP_BK_S have high positive weights and in conjunction with the very strong negative weight of ANREDE_KZ can refer to (men that are less dreamful, less socially-minded, less cultural-minded not financial minimalists and conservative return shopper types).
The assumption is supported by the negative weights of ANREDE_KZ, SEMIO_KAEM, SEMIO_DOM, SEMIO_KRIT, SEMIO_RAT and FINANZ_ANLEGER correlated with (women that are less dominant-minded, less critical-minded, less rational, have less combative attitude and are less prone to investing).
Clustering¶
Apply Clustering to General Population¶
Investigate the change in within-cluster distance across number of clusters.¶
Fig19 : Total Distance to Cluster center per cluster center.¶
In [493]:
#THIS TAKES A LONG TIME - DON'T RUN EVERY TIME!
# compute the average within-cluster distances.
# Over a number of different cluster counts...
scores = {}
for k in range(2, 31):
# run k-means clustering on the data and...
scores[k] = np.abs(KMeans(n_clusters=k, n_jobs=-1).fit(X_pca).score(X_pca))
# 1print('Calculated: ', k)
# Plot relationship plot
fig, ax = plt.subplots(figsize=(16,10))
ax = pd.Series(scores).plot(marker='o', color='slategray')
ax.set_xticks(np.arange(2, 31), minor=False);
ax.set_xlabel("# of Clusters")
ax.set_ylabel("Total Distance to centroid");
Final general population clustering¶
In [129]:
# Chosen number of clusters
K=8
In [130]:
# Re-fit the k-means model with the selected number of clusters and obtain
# cluster predictions for the general population demographics data.
kmeans = KMeans(n_clusters=K, n_jobs=-1, random_state=0).fit(X_pca)
In [131]:
kmeans.labels_.shape
Out[131]:
In [132]:
kmeans.cluster_centers_.shape
Out[132]:
Summary¶
The curve in the Distance to Number of Clusters relationship plot (Fig19) is quite even and the the most obvious elbow in the curve where adding additional clusters less effective in reducing the within-cluster distance can be inferred to be 8 and this number has been used as the number of clusters to fit the the KMeans algorithm on the transformed demographics data and predict cluster labels.
Another elbow can be observed around the 12 clusters but the k = 8 was preferred in this iteration in order to keep the number of clusters lower and see if meaningful separation can be observed with this number.
Apply All Steps to the Customer Data¶
In [133]:
# Load in the customer demographics data.
customers = pd.read_csv('Customers.csv', delimiter=';')
In [134]:
customers.head()
Out[134]:
Data cleaning¶
In [135]:
df_cust, df_cust_highna = clean_data(customers)
In [136]:
df_cust.shape
Out[136]:
In [137]:
df_cust_highna.shape
Out[137]:
Feature transformation¶
In [138]:
# Impute
X_cust_imp = imp.transform(df_cust)
In [139]:
# Scale
X_cust_sc = scaler.transform(X_cust_imp)
Dimensionality reduction¶
In [140]:
X_cust_pca = pca.transform(X_cust_sc)
X_cust_pca.shape
Out[140]:
In [141]:
df_cust_pca = pd.DataFrame(X_cust_pca, columns=np.arange(1, 13))
df_cust_pca.describe()
Out[141]:
Clustering¶
In [142]:
X_cust_kmeans = kmeans.predict(X_cust_pca)
In [143]:
np.unique(X_cust_kmeans)
Out[143]:
In [144]:
X_cust_kmeans.shape
Out[144]:
Compare Customer Data to Demographics Data¶
Compare proportions¶
Compare the proportion of data in each cluster for the customer data to the proportion of data in each cluster for the general population.
In [146]:
# General population dataset kmean labels
gen_clust_labels = kmeans.labels_
# Add the high NaN dataset as a cluster with label -1
gen_clust_labels_na = np.append(gen_clust_labels, [-1] * df_highna.shape[0])
# Make proportion table
gen_clust_freq = itemfreq(gen_clust_labels_na)
gen_proportion_col = gen_clust_freq[:, 1] / gen_clust_freq[:, 1].sum()
gen_clust_freq = np.c_[gen_clust_freq, gen_proportion_col]
# customer dataset
cust_clust_labels = X_cust_kmeans
# Add the high NaN dataset as a cluster with label -1
cust_clust_labels_na = np.append(cust_clust_labels, [-1] * df_cust_highna.shape[0])
# Make proportion table
cust_clust_freq = itemfreq(cust_clust_labels_na)
cust_proportion_col = cust_clust_freq[:, 1] / cust_clust_freq[:, 1].sum()
cust_clust_freq = np.c_[cust_clust_freq, cust_proportion_col]
Fig20 : Barplots of Cluster proportions at the Population and Customer Datasets¶
In [147]:
# Plot
fig, axes = plt.subplots(1,2, figsize=(16, 10), sharey=True)
sns.barplot(x=gen_clust_freq[:, 0], y=gen_clust_freq[:, 2], color='b', ax=axes[0]);
sns.barplot(x=cust_clust_freq[:, 0], y=cust_clust_freq[:, 2], color='r', ax=axes[1]);
axes[0].set_title('General Population')
axes[1].set_title('Customer Data');
Fig21 : Barplot of Cluster proportions at the Population and Customer Datasets¶
In [148]:
# create plot
fig, ax = plt.subplots(figsize=(16, 10))
index = np.append([-1], np.arange(K))
bar_width = 0.35
opacity = 0.8
gen = plt.bar(index, gen_clust_freq[:, 2], bar_width,
alpha=opacity,
color='b',
label='General Population')
cust = plt.bar(index + bar_width, cust_clust_freq[:, 2], bar_width,
alpha=opacity,
color='g',
label='Customer Data')
plt.xlabel('# of Cluster')
plt.ylabel('Proportion of cluster')
plt.title('Proportions of Clusters for General and Customer Data')
plt.xticks(index + bar_width, index)
plt.legend()
plt.tight_layout()
sns.set(style="whitegrid")
sns.despine();
Fig22 : Barplots of difference of Cluster proportions at the Customer and Population Datasets¶
In [149]:
#Calculate differece of proportions
prop_diff = cust_clust_freq[:, 2] - gen_clust_freq[:, 2]
positive= prop_diff > 0
# create plot
fig, ax = plt.subplots(figsize=(14, 8))
index = np.append([-1], np.arange(K))
bar_width = 0.35
opacity = 0.8
diff = plt.bar(index, prop_diff, bar_width,
alpha=opacity,
color=np.vectorize({True: 'k', False: 'r'}.get)(positive)
)
plt.xlabel('# of Cluster')
plt.ylabel('Difference of proportions')
plt.title('Difference of Proportions of Clusters for General and Customer Data')
plt.xticks(index)
plt.tight_layout()
sns.set(style="whitegrid")
sns.despine();
Observations
- Cluster 3 is the most overrepresented and is actually dominating the customer data and cluster 6 the most underrepresented.
- The high-missing NaN cluster is quite overrepresented in the customer data which suggests a possibly methodical problem in the data acquisition of the customer data.
In [150]:
# Add cluster labels
df_cust_pca['cluster'] = X_cust_kmeans
Fig23: Pairplots of all the principal components colored by cluster¶
In [151]:
g = sns.pairplot(df_cust_pca,
vars = np.arange(1, 13),
diag_kind = 'kde',
hue="cluster",
palette="Paired");
Observations
- By looking at the density plots we see that there seems to be a much clearer separation of clusters in the first 3 components and a more uniform spread of clusters at the rest of the components.
- There is an interesting binary separation in most scatter plots that involve the third component which agrees with our assumption of it being a latent gender feature.
- We have not include the high-NAN in this plot as it refers to different dataset characteristics.
Let's have a closer look at the first three principal components.
In [152]:
# Rename components for better reading
df_cust_pca = df_cust_pca.rename(columns = {1: 'pop_density-finance_status',
2: 'age-generation',
3: 'gender'}
)
Fig24: Pairplots of the first 3 principal components colored by cluster.¶
In [153]:
g = sns.pairplot(df_cust_pca,
vars = ['pop_density-finance_status', 'age-generation', 'gender'],
diag_kind = 'kde',
hue="cluster",
palette="Paired");
Observations
- Cluster No 3 (dark green) seems to be strongly correlated to lower values of the 1st component (Population Density - Financial status) , values around 0 for the 2nd component (Age, Generation and Culture) and higher values of the 3rd component (Gender and their main corresponding personality types).
- Cluster No 6 (ligh orange) seems to be the opposite of Cluster no 3 and correlates with higher values of the first component, values around 0 for the 2nd component (Age, Generation and Culture) and lower values for the 3rd principal component(Gender and their main corresponding personality)
In [155]:
def map_kmeans_weights_to_feats(kmeans, df, clust_no):
'''Map pca weights to individual features
and return two pd.Series on with the highest
positive weights and one with the lowest negative
weights'''
weights = pd.DataFrame(np.round(kmeans.cluster_centers_, 4), columns=df.keys())
centroid = weights.iloc[clust_no, :]
cent_pos = centroid[centroid > 0].sort_values(ascending=False)
cent_neg = centroid[centroid < 0].sort_values(ascending=True)
return cent_pos, cent_neg
Summary¶
Fig20, Fig21 and Fig22 observations:
- Cluster 3 is the most overrepresented and is actually dominating the customer data and cluster 6 the most underrepresented.
- The high-missing NaN cluster is quite overrepresented in the customer data which suggests a possibly methodical problem in the data acquisition of the customer data.
Fig22 observations:
- By looking at the density plots we see that there seems to be a much clearer separation of clusters in the first 3 components and a more uniform spread of clusters at the rest of the components.
- There is an interesting binary separation in most scatter plots that involve the 3rd component which agrees with our assumption of it being a latent gender feature.
What kinds of people are part of a cluster that is overrepresented in the customer data compared to the general population?¶
In [156]:
# Overrepresented
over_pos, over_neg = map_kmeans_weights_to_feats(kmeans, df_cust_pca.iloc[:, :-1], 3)
create_bar_table(over_pos, over_neg)
- Looking at the centroid distances table with the highly negative value for
pop_density-finance_statusand the strongly positive value forgender, confirms the pairplot observations and strengthens the assumption that the overrepresented group are: Men that are less dreamful, less socially-minded, less cultural-minded, not financial minimalists, conservative return shopper types, that have a higher financial status and live in less densely populated areas. The negative weight for theage-generationprincipal component in conjunction with the attitude characteristics of these men suggests that they could be of relative younger age.
What kinds of people are part of a cluster that is underrepresented in the customer data compared to the general population?¶
In [157]:
# Underrepresented
under_pos, under_neg = map_kmeans_weights_to_feats(kmeans, df_cust_pca.iloc[:, :-1], 6)
create_bar_table(under_pos, under_neg)
- Looking at the centroid distances table with the highly positive value for
pop_density-finance_statusand the strongly negative value forgenderalong with negative weight for theage-generationprincipal component confirms the pairplot observations and strengthens the assumption that the underrepresented group are: Younger women that are less dominant-minded, less critical-minded, less rational, have a less combative attitude and are less prone to investing and live in densely populated areas with a lower financial status. They have lower religious affinity, lower money saving culture and are less traditional minded.
References¶
[1]:https://stackoverflow.com/questions/20250771/remap-values-in-pandas-column-with-a-dict
[2]:https://stackoverflow.com/questions/16992713/translate-every-element-in-numpy-array-
[3]:https://stackoverflow.com/questions/36653443/retaining-nan-values-after-get-dummies-in-pandas
[4]:http://scikit-learn.org/stable/modules/clustering.html#clustering-evaluation
[5]:http://greenteapress.com/thinkstats2/html/thinkstats2010.html#sec94
[6]:https://github.com/jadianes/data-science-your-way/tree/master/03-dimensionality-reduction-and-clustering
[7]:https://www.mailman.columbia.edu/research/population-health-methods/cluster-analysis-using-k-means