Following is the Python program I wrote to fulfill the third 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.
Using the Gapminder database, I would like to see the variables that are the most influencing the income per person (2010 GDP per capita in constant 2000 US$) from the following list of variables:
For the question I'm interested in, 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.linear_model import LassoLarsCV
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']
response_var = 'incomeperperson'
constructor_dict = dict()
for var in explanatory_vars + [response_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))))
The number of countries has unfortunately been almost divided by two due to missing data. Nevertheless the k-fold cross-validation technique will be used.
In [5]:
predictors = numeric_data[explanatory_vars]
target = numeric_data[response_var]
# Standardize predictors to have mean=0 and std=1
from sklearn import preprocessing
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 (30% of the data). Then the lasso regression analysis will be carried out with 10 folds in the k-fold cross-validation process.
In [6]:
# split data into train and test sets
pred_train, pred_test, tar_train, tar_test = train_test_split(std_predictors, target,
test_size=.3, random_state=123)
# specify the lasso regression model
model=LassoLarsCV(cv=10, precompute=False).fit(pred_train,tar_train)
# print variable names and regression coefficients
pd.Series(model.coef_, index=explanatory_vars, name='Regression coefficients').to_frame()
Out[6]:
From the lasso regression, the employ rate and the urban rate have been discarded of the explanatory variables list. Then the most influential variables are in order :
In [7]:
# plot coefficient progression
plt.semilogx(model.alphas_, model.coef_path_.T)
plt.legend(explanatory_vars)
plt.axvline(model.alpha_, linestyle='--', color='k',
label='alpha CV')
plt.ylabel('Regression Coefficients')
plt.xlabel(r'$\alpha$')
plt.title('Regression Coefficients Progression for Lasso Paths');
The plot above showing the evolution of the regression coefficient proves the order of appearance of the four remaining explanatory variables; the internet use rate starting right from the beginning as it is the most infuential.
In [8]:
# plot mean square error for each fold
plt.semilogx(model.cv_alphas_, model.cv_mse_path_, ':')
plt.semilogx(model.cv_alphas_, model.cv_mse_path_.mean(axis=-1), 'k',
label='Average across the folds', linewidth=2)
plt.axvline(model.alpha_, linestyle='--', color='k',
label='alpha CV')
plt.legend()
plt.xlabel(r'$\alpha$')
plt.ylabel('Mean squared error')
plt.title('Mean squared error on each fold');
Unexpectedly one fold have a rising error with the lasso regression analysis progression. A country may have an unusual combination of parameters explaining it.
In [9]:
# MSE from training and test data
from sklearn.metrics import mean_squared_error
train_error = mean_squared_error(tar_train, model.predict(pred_train))
test_error = mean_squared_error(tar_test, model.predict(pred_test))
In [10]:
# R-square from training and test data
rsquared_train=model.score(pred_train,tar_train)
rsquared_test=model.score(pred_test,tar_test)
In [11]:
pd.DataFrame({'MSE' : [train_error, test_error], 'R-square' : [rsquared_train, rsquared_test]},
index=['Training data', 'Test data'])
Out[11]: