import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from nose.tools import *
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.linear_model import ElasticNet
from sklearn.metrics import make_scorer
from scipy import stats
import seaborn as sns
Part 1. Loading packages
The following packages have been imported NymPy, Pandas, Matplotlib, Scipy, Seaborn. Sklearn is the most widely used package for the machine learning process. The fallowing subpackages have been used:
- train_test_split
- linear_model
- model_selection
- metrics
- tree
- ensemble
- preprocessing
life_data = pd.read_csv('data/lifeExpectancy.csv', sep = ',')
life_data = life_data.drop('Year', axis = 1)
Part 2. Reading the data
The data is saved as a csv file as LifeExpectancy.csv and it is read and stored in the life_data variable.The Year column is dropped as it will not be used in the analysis. Below the first 5 rows are shown. The data contains 21 columns and 2938 rows with the header row. The table contains data about:
- Country
- Status
- Life Expectancy
- Adult Mortality
- Alcohol
- percentage expenditure
- Hepatitis B
- Measles
- BMI
- under-five deaths
- Polio
- Total expenditure
- Diphtheria
- HIV/AIDS
- GDP
- Population
- thinness 1-19 years
- thinness 5-9 years
- Income composition of resources
- Schooling
With the exclution of Country name and Status(either developed or developing) all of the data is numeric. The values are either in years, precentages, millions or dollars in the case of Gross Domestic Product (GDP)
life_data.head()
status = pd.get_dummies(life_data.Status)
life_data = pd.concat([life_data, status], axis = 1)
life_data = life_data.drop(['Status'], axis=1)
life_data.rename(columns = {'Deloping' : '0', 'Developed' : 1})
life_data = life_data.groupby('Country').mean()
life_data.head()
Part 3. EDA
life_data.columns
plt.scatter(life_data[' HIV/AIDS'], life_data['Life expectancy '])
plt.xlabel('HIV/AIDS')
plt.ylabel('Life expectancy')
plt.scatter(life_data.GDP, life_data['Life expectancy '])
plt.xlabel('GDP')
plt.ylabel('Life expectancy')
plt.scatter(life_data[' BMI '], life_data['Life expectancy '])
plt.xlabel('BMI')
plt.ylabel('Life expectancy')
plt.scatter(life_data['under-five deaths '], life_data['Life expectancy '])
plt.xlabel('under-five deaths')
plt.ylabel('Life expectancy')
plt.scatter(life_data['Alcohol'], life_data['Life expectancy '])
plt.xlabel('Alcohol')
plt.ylabel('Life expectancy')
plt.scatter(life_data['Adult Mortality'], life_data['Life expectancy '])
plt.xlabel('Adult Mortality')
plt.ylabel('Life expectancy')
plt.scatter(life_data['Schooling'], life_data['Life expectancy '])
plt.xlabel('Schooling')
plt.ylabel('Life expectancy')
plt.scatter(life_data['percentage expenditure'], life_data['Life expectancy '])
plt.xlabel('Percentage Healhcare expenditure')
plt.ylabel('Life expectancy')
Using the scatter plot we plot the Life Expectancy against some other varyables to see if there is any correlation between them.
There seem to be a positive correlation between The Percentage of Healthcare Expenditure, Schooling, GDP and BMI and Life Expectancy, while there is a negative one between Adult Mortality, AIDS and Life Expectancy,there does not seem to have any correlation between Alcohol, under 5 years – old deaths and Life Expectancy.
plt.figure(figsize = (14, 10))
sns.heatmap(life_data.corr(), annot = True)
Now we will plot the correlation matrix visulizing it with a heatmap. The legend tells that the warmer colors show higher and positive correlation, while the colder low or negative.
There is a very high correlation between thinness of 5-9 year-old and that of 1-19 year-old. Also between population and infant deaths, under 5 deaths, another is between schooling and income composition of resources. On the other hand Life expectancy and Adult Mortality are very highly negatively correlated.
Part 4. Preprocessing the data
The raw data is not suitable for us to start building a model so some preprocessing will be done. First the Status of the country is turned into numberical with the get_dummies function, so we get 2 new columns. The original column is being dropped. Second the data is being grouped by the country and we find the mean values during the 2000 – 2015 year period.
Then the Life expectancy column is removed to form the life_labels variable or the output, and the rest is stored as the life_features variable. Now we consider that we have some null values in the table, the isnull function has been used to find the with the boolean True. Below that the number of null values are displayed in each separate column. It is mostly situated in the Population and GDP columns.
Now the missing values are filled with the mean of its respective column. This will create some distortions, but the other option in removing parts of the table will shrink the data so it will be avoided here because the number of rows is not that high. The final shape for the life_features is 193 rows to 20 columns. Finally considering the large differences in the values of the columns, there will be some scaling with the MinMaxScaler function.
Now we will split the data into a training part of 70% and a testing of 30%. Cross validation will be initialized with the creation of 5 fold split.
life_labels = life_data['Life expectancy ']
life_features = life_data.drop('Life expectancy ', axis = 1)
life_features.isnull().head()
life_features.isnull().sum()
life_labels.isnull().sum()
life_features.fillna(value = life_features.mean(), inplace = True)
life_labels.fillna(value = life_labels.mean(), inplace = True)
stats.describe(life_features[1:])
min_max_scaler = MinMaxScaler()
life_features = min_max_scaler.fit_transform(life_features)
life_features
life_features_train, life_features_test, life_labels_train, life_labels_test = train_test_split(
life_features, life_labels, train_size = 0.7, test_size = 0.3)
Part 5 Linear Regression and additions
linear_model = LinearRegression()
linear_model.fit(life_features_train, life_labels_train)
print('R_square score on the training: %.2f' % linear_model.score(life_features_train, life_labels_train))
linear_model_predict = linear_model.predict(life_features_test)
print('Coefficients: \n', linear_model.coef_)
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, linear_model_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, linear_model_predict))
print('R_square score: %.2f' % r2_score(life_labels_test, linear_model_predict))
Several algorithms will be tried out. First the classical linear regression.
The model is fitted first on the trained data the R square is 0.92 on the training data.Later its R square is checked on the testing data. The score is 87 % in the iteration of writing. We also calculate the MAE (the modulus between the predicted and the real value) at 2.32 and the MSE (the same only put to the power of 2) at 9.8. Now we will try to have some changes to the initial model. We will use both the properties of the Ridge regression and the Lasso and eventually the ElasticNet to see if the score can be improved.
scoring = make_scorer(r2_score)
grid_cv = GridSearchCV(Ridge(),
param_grid={'alpha': range(0, 10), 'max_iter' : [10, 100, 1000]},
scoring=scoring, cv=5, refit=True)
grid_cv.fit(life_features_train, life_labels_train)
print("Best Parameters: " + str(grid_cv.best_params_))
result = grid_cv.cv_results_
print("R^2 score on training data: %.2f" %grid_cv.score(life_features_train, life_labels_train))
print("R^2 score: %.2f"
% r2_score(life_labels_test, grid_cv.best_estimator_.predict(life_features_test)))
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, linear_model_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, linear_model_predict))
The Ridge Regression uses L2 regularization to minimize the weights of the coefficients, this is controlled with the hyperparameter alpha. Increasing the value of alpha decreases the weights of the coefficients.
We performed a grid search with cross-validation on the grid regression with alpha varying between 0 and 10 and having 3 max iterations of 10, 100 and 1000. Finally the best parameters here are alpha = 1, and max iterations = 10. The R square on the training data is 90% compared to 92 % on the standard linear model There is some improvement of the R square on the test data statistic 91% vs 87% the errors stayed the same. The MAE is 1.83 better and MSE is 6.05 again lower the classic linear regression
scoring = make_scorer(r2_score)
grid_cv = GridSearchCV(Lasso(),
param_grid={'alpha': range(0, 10), 'max_iter' : [10, 100, 1000]},
scoring=scoring, cv=5, refit=True)
grid_cv.fit(life_features_train, life_labels_train)
print("Best Parameters: " + str(grid_cv.best_params_))
result = grid_cv.cv_results_
print("R^2 score on training data: %.2f" % grid_cv.score(life_features_train, life_labels_train))
print("R^2 score: %.2f"
% r2_score(life_labels_test, grid_cv.best_estimator_.predict(life_features_test)))
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, linear_model_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, linear_model_predict))
The Lasso Regression uses the L1 regularization which reduces some coefficent weights to zero. Again the regularization is controlled with aplha, decreasing it reduces the weights, also at the same time the max iterations hyperparameter needs to be increased.
The best parameters for the Lasso Regression are alpha = 0, and max iterations set to 100. The R square is a little higher at 92 % on the training data and 92 % on the test data the errors continue to stay the same as the Ridge regression.
scoring = make_scorer(r2_score)
grid_cv = GridSearchCV(ElasticNet(),
param_grid={'alpha': range(0, 10), 'max_iter' : [10, 100, 1000], 'l1_ratio' : [0.1, 0.4, 0.8]},
scoring=scoring, cv=5, refit=True)
grid_cv.fit(life_features_train, life_labels_train)
print("Best Parameters: " + str(grid_cv.best_params_))
result = grid_cv.cv_results_
print("R^2 score on training data: %.2f" % grid_cv.score(life_features_train, life_labels_train))
print("R^2 score: %.2f"
% r2_score(life_labels_test, grid_cv.best_estimator_.predict(life_features_test)))
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, linear_model_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, linear_model_predict))
The ElasticNet combines both L1 and L2 regularization of the previous regression types.
Finally the grid search concerning ElasticNet shows that it performs as well as the Lasso regression so far with alpha = 0, max iterations = 100 and the L1 ration set at 0.1. The R square and the errors remain the same as the Lasso Regression
Part 5 Linear Regression with Polynomial Features
The Polynomial Features function has been used to get the interactions of the input variables only to the power of 2
quad_feature_transformer = PolynomialFeatures(2, interaction_only = True)
quad_feature_transformer.fit(life_features_train)
life_features_train_quad = quad_feature_transformer.transform(life_features_train)
life_features_test_quad = quad_feature_transformer.transform(life_features_test)
poly_model_quad = LinearRegression()
poly_model_quad.fit(life_features_train_quad, life_labels_train)
accuracy_score_quad = poly_model_quad.score(life_features_train_quad, life_labels_train)
print(accuracy_score_quad)
The Linear Regression is being tested on the training data with the new Polynomial Features, the accuracy is 1, meaning the model has learned all the data. An prediction for the life_features_test_quad variable has been made.
poly_model_quad_predict = poly_model_quad.predict(life_features_test_quad)
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, poly_model_quad_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, poly_model_quad_predict))
print('R_square score: %.2f' % r2_score(life_labels_test, poly_model_quad_predict))
All the errors are significantly higher then the privious models, and the R square is in this case negative (this happens only in sklearn). This is the worst performing model for now.
Part 5 Decision Tree Regression
decision_tree_model = DecisionTreeRegressor()
decision_tree_fit = decision_tree_model.fit(life_features_train, life_labels_train)
decision_tree_score = cross_val_score(decision_tree_fit, life_features_train, life_labels_train, cv = 5)
print("mean cross validation score: %.2f" % np.mean(decision_tree_score))
print("score without cv: %.2f" % decision_tree_fit.score(life_features_train, life_labels_train))
print("R^2 score on the test data %.2f"% r2_score(life_labels_test, decision_tree_fit.predict(life_features_test)))
Now we will try the Decision Tree Regression. Cross Validation has been performed
The R square on the training data is 1 meaning that the algorithm has learned the data by heart, with the cross validation the figure declines to 77% and using the test date we get 80%
Now we use the algorithm to predict the values of the life_features_test.
decision_tree_model_predict = decision_tree_model.predict(life_features_test)
scoring = make_scorer(r2_score)
grid_cv = GridSearchCV(DecisionTreeRegressor(),
param_grid={'min_samples_split': range(2, 10)},
scoring=scoring, cv=5, refit=True)
grid_cv.fit(life_features_train, life_labels_train)
grid_cv.best_params_
print("Best Parameters: " + str(grid_cv.best_params_))
result = grid_cv.cv_results_
print("R^2 score on training data: %.2f" % grid_cv.best_estimator_.score(life_features_train, life_labels_train))
print("R^2 score: %.2f"
% r2_score(life_labels_test, grid_cv.best_estimator_.predict(life_features_test)))
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, decision_tree_model_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, decision_tree_model_predict))
After performing grid search with minimun samples split in the range between 2 and 10 we get the best split of 3.
The R square on the training data is 98%, the algorithms has nearly leaned the data by hearth. On the test data we get R square of 81 %, the MAE is 2.71 and MSE is 16.59. This method is worse than the Elastic Net Regression
Part 6 Random Forest Regression
random_forest_model = RandomForestRegressor()
random_forest_fit = random_forest_model.fit(life_features_train, life_labels_train)
random_forest_score = cross_val_score(random_forest_fit, life_features_train, life_labels_train, cv = 5)
print("mean cross validation score: %.2f"
% np.mean(random_forest_score))
print("score without cv: %.2f"
% random_forest_fit.score(life_features_train, life_labels_train))
print("R^2 score on the test data %.2f"
%r2_score(life_labels_test, random_forest_fit.predict(life_features_test)))
Now we will use the Random Forest Regression
The alogorithm has learned 98% on the training data without cross validation and 88% with, the value is 92 % on the test data.
random_forest_model_predict = random_forest_model.predict(life_features_test)
scoring = make_scorer(r2_score)
grid_cv = GridSearchCV(RandomForestRegressor(),
param_grid={'min_samples_split': range(2, 10)},
scoring=scoring, cv=5, refit=True)
grid_cv.fit(life_features_train, life_labels_train)
grid_cv.best_params_
result = grid_cv.cv_results_
print("Best Parameters: " + str(grid_cv.best_params_))
result = grid_cv.cv_results_
print("R^2 score on training data: %.2f" % grid_cv.best_estimator_.score(life_features_train, life_labels_train))
print("R^2 score: %.2f"
% r2_score(life_labels_test, grid_cv.best_estimator_.predict(life_features_test)))
print("Mean squared error: %.2f"
% mean_squared_error(life_labels_test, random_forest_model_predict))
print("Mean absolute error: %.2f"
% mean_absolute_error(life_labels_test, random_forest_model_predict))
After performing grid search in the range 2 to 10 for the minimun samples split we get 2 for the split. The R square on the training is 98 % while on the test it is 91% and a MAE of 1.89 and MSE of 6.71.
Part 7 Conclution
After comparing all the algorithms we can conclude the Lasso and the Elastic Net Regression offer which are the same:
- Best Parameters: {‘alpha’: 0, ‘max_iter’: 10}
- R square on the test data of 92%
- MAE of 1.83
- MSE of 6.05
Part 8 Sources
The following sources have been used:
- https://www.kaggle.com/kumarajarshi/life-expectancy-who/data
- Introduction to Machine Learning with Python by Andreas C. Müller & Sarah Guido
- Labs of the couse
- stack overflow
- Lectures of the course