Following is the Python program I wrote to fulfill the last assignment of the Machine Learning for Data Analysis online course.
I decided to use Jupyter Notebook as it is a pretty way to write code and present results.
I decided to use the same research question than for the previous assignment on Lasso regression analysis.
Using the Gapminder database, I would like to see the variables that are the most influencing income per person (2010 GDP per capita in constant 2000 US$). It will therefore be my test variable.
The cluster analysis will be carried out on the following variables:
The countries for which data are missing will be discarded. As missing data in Gapminder database are replaced directly by NaN no special data treatment is needed.
In [1]:
# Magic command to insert the graph directly in the notebook
%matplotlib inline
# Load a useful Python libraries for handling data
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from IPython.display import Markdown, display
from sklearn.cross_validation import train_test_split
from sklearn import preprocessing
from sklearn.cluster import KMeans
In [2]:
# Read the data
data_filename = r'gapminder.csv'
data = pd.read_csv(data_filename)
data = data.set_index('country')
General information on the Gapminder data
In [3]:
display(Markdown("Number of countries: {}".format(len(data))))
In [4]:
explanatory_vars = ['employrate', 'urbanrate', 'polityscore', 'lifeexpectancy', 'internetuserate', 'relectricperperson']
test_var = 'incomeperperson'
constructor_dict = dict()
for var in explanatory_vars + [test_var, ]:
constructor_dict[var] = pd.to_numeric(data[var], errors='coerce')
numeric_data = pd.DataFrame(constructor_dict, index=data.index).dropna()
display(Markdown("Number of countries after discarding countries with missing data: {}".format(len(numeric_data))))
In [5]:
predictors = numeric_data[explanatory_vars]
target = numeric_data[test_var]
# Standardize predictors to have mean=0 and std=1
std_predictors = predictors.copy()
for var in std_predictors.columns:
std_predictors[var] = preprocessing.scale(std_predictors[var].astype('float64'))
# Check standardization
std_predictors.describe()
Out[5]:
The table above proves the explanatory variables have been standardized (i.e. mean=0 and std=1).
Next the data will be split in two sets; the training set (70% of the data) and the test set (the remaining 30%).
In [6]:
# split data into train and test sets
cluster_train, cluster_test = train_test_split(std_predictors,
test_size=.3,
random_state=123)
We will now perform a k-means cluster analysis for 1 to 10 clusters.
In [7]:
from scipy.spatial.distance import cdist
clusters = range(10)
meandist = list()
for k in clusters:
model = KMeans(n_clusters=k+1)
model.fit(cluster_train)
# clusterassign = model.predict(cluster_train)
meandist.append(sum(np.min(cdist(cluster_train,
model.cluster_centers_,
'euclidean'),
axis=1))
/ cluster_train.shape[0])
In [8]:
plt.plot(clusters, meandist)
plt.xlabel('Number of clusters')
plt.ylabel('Average distance')
plt.title('Selecting k with the Elbow Method');
The elbow seems to be at 2 clusters in the figure above. But we will try the 3 clusters solution as it is a better case to use the canonical variables visualization.
In [9]:
model3=KMeans(n_clusters=3)
model3.fit(cluster_train)
clusassign=model3.predict(cluster_train)
In [10]:
# plot clusters
color_map = {0 : 'r', 1 : 'b', 2 : 'g'}
def color(x):
return color_map[x]
colors = list(map(color, model3.labels_))
from sklearn.decomposition import PCA
pca_2 = PCA(2)
plot_columns = pca_2.fit_transform(cluster_train)
plt.scatter(x=plot_columns[:,0], y=plot_columns[:,1], c=colors)
plt.xlabel('Canonical variable 1')
plt.ylabel('Canonical variable 2')
plt.title('Scatterplot of Canonical Variables for 3 Clusters');
From the evolution of the distance with the number of clusters, it was seen that 2 or 3 clusters seems the best division. The visualization above of the three clusters (group 0 in red, 1 in blue and 2 in green) tends to favor a two clusters approach due to the overlapping of group 1 and 2.
In [11]:
new_clus = pd.Series(model3.labels_,
index=cluster_train.index,
name='cluster').to_frame()
new_clus.T
Out[11]:
In [12]:
merged_train = cluster_train.merge(new_clus, left_index=True, right_index=True)
merged_train.head()
Out[12]:
In [13]:
counts = merged_train.cluster.value_counts()
counts.name = '# countries'
counts.index.name = 'cluster'
counts.sort_index().to_frame()
Out[13]:
In [14]:
merged_train.groupby('cluster').mean()
Out[14]:
If we set employment rate aside, the group 0 gathers the better countries (high urban rate, democratic system, high life expectancy, high internet use rate and high electricity consumption) when the groups 1 and 2 gather the countries not having those characteristics.
The difference between the group 1 and 2 is mainly on the employment rate. The group 1 having a important employment rate when the group 2 has a low one.
The tree clusters will be tested against our test variable income per person. To confirm the accuracy of the clusters obtained.
In [15]:
target_train, target_test = train_test_split(target,
test_size=.3,
random_state=123)
merged_train_all=merged_train.merge(target_train.to_frame(),
left_index=True,
right_index=True)
sub1 = merged_train_all[[test_var, 'cluster']].dropna()
import statsmodels.formula.api as smf
import statsmodels.stats.multicomp as multi
income = smf.ols(formula='{} ~ C(cluster)'.format(test_var), data=sub1).fit()
income.summary()
Out[15]:
The ANOVA test confirms there is a significant variation of income per person with cluster (p-value = 1.46e-10 << 0.05).
In [16]:
display(Markdown('means for income by cluster :'))
sub1.groupby('cluster').mean()
Out[16]:
In [17]:
display(Markdown('standard deviations for income by cluster : '))
sub1.groupby('cluster').std()
Out[17]:
The group 0 has the highest income per person (even if we substract the standard deviation, it stays higher than the average of groups 1 and 2).
Group 1 has a lower average than group 2. But due to their respective standard deviation, they have a large overlapping confirming the idea of merging them for a better description.
In [18]:
mc1 = multi.MultiComparison(sub1[test_var], sub1['cluster'])
res1 = mc1.tukeyhsd()
res1.summary()
Out[18]:
The multiple comparison above confirms that the cluster analysis would be better keeping only two clusters as groups 1 and 2 have no significant difference.
In [19]:
meandist_test = list()
for k in clusters:
model = KMeans(n_clusters=k+1)
model.fit(cluster_train)
meandist_test.append(sum(np.min(cdist(cluster_test,
model.cluster_centers_,
'euclidean'),
axis=1))
/ cluster_test.shape[0])
plt.plot(clusters, meandist, label='Train data')
plt.plot(clusters, meandist_test, label='Test data')
plt.legend()
plt.xlabel('Number of clusters')
plt.ylabel('Average distance')
plt.title('Selecting k with the Elbow Method');
The test data will be plotted in the canonical variables space using lighter colors to quickly visualize their dispersion compare to the training data.
In [20]:
test_assign = model3.predict(cluster_test)
color_map = {0 : 'lightcoral', 1 : 'lightblue', 2 : 'lawngreen'}
def color(x):
return color_map[x]
color_test = list(map(color, test_assign))
plot_col_test = pca_2.transform(cluster_test)
plt.scatter(x=plot_columns[:,0], y=plot_columns[:,1], c=colors)
plt.scatter(x=plot_col_test[:,0], y=plot_col_test[:,1], c=color_test)
plt.xlabel('Canonical variable 1')
plt.ylabel('Canonical variable 2')
plt.title('Scatterplot of Canonical Variables for 3 Clusters');
From the two figures above, the clustering in the training dataset seems to be confirmed by the test dataset. This confirms the clustering obtained on the training set.
This assignment allows me to apply the unsupervised k-means cluster analysis on the Gapminder dataset. For the selected data, it seems that the countries are best described in a two clusters model; one have high quality of life (high urban rate, high life expectancy, high electricity consumption, high democratic system and high internet use rate) and the other one have lower quality of life. That cluster analysis was confirmed when testing for the income per person on the generated clusters.