Module/Package import
In [1]:
import numpy as np # numpy module for linear algebra
import pandas as pd # pandas module for data manipulation
import matplotlib.pyplot as plt # module for plotting
import seaborn as sns # another module for plotting
In [2]:
import warnings # to handle warning messages
warnings.filterwarnings('ignore')
In [3]:
from sklearn.linear_model import LinearRegression # package for linear model
import statsmodels.api as sm # another package for linear model
import statsmodels.formula.api as smf
import scipy as sp
In [4]:
from sklearn.model_selection import train_test_split # split data into training and testing sets
Dataset import
The dataset that we will be using is the meuse dataset.
As described by the author of the data: "This data set gives locations and topsoil heavy metal concentrations, along with a number of soil and landscape variables at the observation locations, collected in a flood plain of the river Meuse, near the village of Stein (NL). Heavy metal concentrations are from composite samples of an area of approximately 15 m $\times$ 15 m."
In [5]:
soil = pd.read_csv("soil.csv") # data import
soil.head() # check if read in correctly
Out[5]:
In [6]:
soil.shape # rows x columns
Out[6]:
In [7]:
print(soil.describe())
In [8]:
index = pd.isnull(soil).any(axis = 1)
soil = soil[-index]
soil = soil.reset_index(drop = True)
In [9]:
soil.shape
Out[9]:
In [10]:
n = soil.shape[1]
Exploratory Data Analysis (EDA)¶
1. Correlation HeatmapIn [11]:
plt.figure(figsize = (10, 8))
sns.heatmap(soil.corr(), annot = True, square = True, linewidths = 0.1)
plt.ylim(n-1, 0)
plt.xlim(0, n-1)
plt.title("Pearson Correlation Heatmap between variables")
plt.show()
In [12]:
plt.figure(figsize = (10, 8))
corr = soil.corr()
mask = np.zeros_like(corr)
mask[np.triu_indices_from(mask)] = True
with sns.axes_style("white"):
sns.heatmap(corr, mask = mask, linewidths = 0.1, vmax = .3, square = True)
plt.ylim(n-1, 0)
plt.xlim(0, n-1)
plt.title("correlation heatmap without symmetric information")
plt.show()
In [13]:
fig, axs = plt.subplots(figsize = (8, 12))
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0.2)
plt.subplot(2, 1, 1)
# correlation heatmap without annotation
sns.heatmap(soil.corr(), linewidths = 0.1, square = True, cmap = "YlGnBu")
plt.ylim(n-1, 0)
plt.xlim(0, n-1)
plt.title("correlation heatmap without annotation")
plt.subplot(2, 1, 2)
# correlation heatmap with annotation
sns.heatmap(soil.corr(), linewidths = 0.1, square = True, annot = True, cmap = "YlGnBu")
plt.ylim(n-1, 0)
plt.xlim(0, n-1)
plt.title("correlation heatmap with annotation")
plt.show()
2. BoxplotsIn [14]:
plt.figure(figsize = (20, 10))
soil.iloc[:, 2:].boxplot()
# soil.iloc[:, 2:].plot(kind = "box") # 2nd method
plt.title("boxplot for each variable", fontsize = "xx-large")
plt.show()
In [15]:
fig, axs = plt.subplots(2, int(n/2), figsize = (20, 10))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'pink'] # to set color
for i, var in enumerate(soil.columns.values):
if var == "landuse":
axs[1, i-int(n/2)].set_xlabel(var, fontsize = 'large')
continue
if i < int(n/2):
axs[0, i].boxplot(soil[var])
axs[0, i].set_xlabel(var, fontsize = 'large')
axs[0, i].scatter(x = np.tile(1, soil.shape[0]), y = soil[var], color = colors[i])
else:
axs[1, i-int(n/2)].boxplot(soil[var])
axs[1, i-int(n/2)].set_xlabel(var, fontsize = 'large')
axs[1, i-int(n/2)].scatter(x = np.tile(1, soil.shape[0]), y = soil[var],
color = colors[i-int(n/2)])
plt.suptitle("boxplot of variables", fontsize = 'xx-large')
plt.show()
3. scatterplotsIn [16]:
fig, axs = plt.subplots(2, int(n/2), figsize = (20, 10))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k'] # to set color
for i, var in enumerate(soil.columns.values):
if var == "landuse":
axs[1, i-int(n/2)].set_xlabel(var, fontsize = 'large')
continue
if i < int(n/2):
axs[0, i].scatter(soil[var], soil["lead"], color = colors[i])
axs[0, i].set_xlabel(var, fontsize = 'large')
else:
axs[1, i-int(n/2)].scatter(soil[var], soil["lead"], color = colors[i-int(n/2)])
axs[1, i-int(n/2)].set_xlabel(var, fontsize = 'large')
plt.suptitle("scatterplot of lead vs. predictors", fontsize = 'xx-large')
plt.show()
4. histograms and density plotsIn [17]:
fig, axs = plt.subplots(2, int(n/2), figsize = (20, 10))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'pink'] # to set color
for i, var in enumerate(soil.columns.values):
if var == "landuse":
axs[1, i-int(n/2)].set_xlabel(var, fontsize = 'large')
continue
if i < int(n/2):
#axs[0, i].hist(soil[var], color = colors[i], alpha = 0.8)
#axs[0, i].set_xlabel(var, fontsize = 'large')
# add density plot
sns.distplot(soil[var], color = colors[i], ax = axs[0, i])
#kde = gaussian_kde(soil[var])
#xaxis = np.linspace(min(soil[var]), max(soil[var]), 20)
#axs[0, i].plot(xaxis, kde(xaxis) * 100)
else:
sns.distplot(soil[var], color = colors[i-int(n/2)], ax = axs[1, i-int(n/2)])
plt.suptitle("histogram and density plot of each variable", fontsize = 'xx-large')
plt.show()
Variable Selection and Modeling
Based on the scatterplot of lead vs. predictors:
cadmium,copper,zincandomseem to have strong positive correlation with leadelev,distanddist.mseem to have strong negative correlation with lead- categorical variable
limeseems to be a good indicator predictor variable
Therefore, we now build a linear model using 6 predictors: cadmium, copper, zinc, elev, dist and lime.
In [18]:
# variables = ["cadmium", "copper", "zinc", "elev", "dist", "dist.m", "lime"]
variables = ["cadmium", "copper", "zinc", "elev", "dist", "lime"]
In [19]:
x = soil[variables]
y = soil["lead"]
Modeling - building a linear model¶
split the dataset into training and testing set
In [20]:
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size = 0.33, random_state = 42)
build a linear model using sklearn package
In [21]:
model = LinearRegression().fit(X_train, y_train)
In [22]:
# R^2 on training set
R2_train = model.score(X_train, y_train)
R2_train
Out[22]:
In [23]:
# R^2 on testing set
R2_test = model.score(X_test, y_test)
R2_test
Out[23]:
In [24]:
# value of intercepts
for var, coef in zip(variables, model.coef_):
print(var, coef)
build a linear model using statsmodel package
In [25]:
df_train = pd.concat([X_train, y_train], axis = 1) # build a dataframe for training set
In [26]:
fullModel = smf.ols("lead ~ cadmium + copper + zinc + elev + lime", data = df_train)
In [27]:
fullModel = fullModel.fit()
print(fullModel.summary())
diagnostics plots
In [28]:
# diagnosticsPlots function definition
def diagnosticsPlots(predictions, truevalues):
'''
This function will plot the following two diagonistics plots:
1. residuals vs. fitted values plot (plot on the left)
2. qqplot of the residuals (plot on the right)
parameters required:
- predictions: predicted values using the linear regression model
- truevalues: true values of the response variable
required modules/packages:
- matplotlib.pyplot (plt)
- scipy (sp)
Author: Kaixin Wang
Created: October 2019
'''
residuals = truevalues - predictions
fig, axs = plt.subplots(figsize = (12, 2.5))
plt.subplot(1, 2, 1)
# residuals vs. fitted values
plt.scatter(predictions, residuals)
plt.plot(predictions, np.tile(np.mean(residuals), residuals.shape[0]), "r-")
plt.xlim([np.min(predictions) - 2.5, np.max(predictions) + 2.5])
plt.title("residuals vs. fitted values")
plt.ylabel("residuals")
plt.xlabel("fitted values")
plt.legend(["E[residuals]"])
print("Average of residuals: ", np.mean(residuals))
# qqplot
plot = plt.subplot(1, 2, 2)
sp.stats.probplot(residuals, plot = plot, fit = True)
plt.show()
In [29]:
pred_val = fullModel.fittedvalues.copy()
true_val = y_train.copy()
diagnosticsPlots(pred_val, true_val)
training and testing Root-Mean-Square Error (RMSE)
The root-mean-square error of a prediction is calcuated as the following:
$$\text{RMSE} = \frac{\sqrt{\sum_{i=1}^n (y_i - \hat {y_i}) ^ 2}}{n}$$where $y_i$ is the true value of the response, $\hat {y_i}$ is the fitted value of the response, n is the size of the input
In [30]:
def RMSETable(train_y, train_fitted, train_size, train_R2, test_y, test_fitted, test_size, test_R2):
'''
This function creates a function that returns a dataframe in the following format:
-------------------------------------------------
RMSE R-squared size
training set train_RMSE train_R2 n_training
testing set test_RMSE tesint_R2 n_testing
-------------------------------------------------
parameters required:
- train_y: true values of the response in the training set
- train_fitted: fitted values of the response in the training set
- train_size: size of the training set
- train_R2: R-squared on the training set
- test_y: true values of the response in the testing set
- test_fitted: fitted values of the response in the testing set
- test_size: size of the testing size
- test_R2: R-squared on the testing set
Author: Kaixin Wang
Created: October 2019
'''
train_RMSE = np.sqrt(sum(np.power(train_y - train_fitted, 2))) / train_size
test_RMSE = np.sqrt(sum(np.power(test_y - test_fitted, 2))) / test_size
train = [train_RMSE, train_size, train_R2]
test = [test_RMSE, test_size, test_R2]
return pd.DataFrame([train, test], index = ["training set", "testing set"], columns = ["RMSE", "R-squared", "size"])
In [31]:
table1 = RMSETable(y_train.values, model.predict(X_train), R2_train, len(y_train),
y_test.values, model.predict(X_test), R2_test, len(y_test))
print(table1)
Model selection¶
Since we are using 5 predictor variables for 152 entries, we can consider using fewer predictors to improve the $R_{adj}^2$, where is a version of $R^2$ that penalizes model complexity.
Selecting fewer predictor will also help with preventing overfitting.
Therefore, observing that zinc, copper and cadmium have relatively stronger correlation with lead, we now build a model using only these three predictors.
using sklearn linear regression module
In [32]:
variables = ["cadmium", "lime", "zinc"]
x = soil[variables]
y = soil["lead"]
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size = 0.33, random_state = 42)
In [33]:
model2 = LinearRegression().fit(X_train, y_train)
# R^2 on training set
R2_train = model2.score(X_train, y_train)
print("Training set R-squared", R2_train)
# R^2 on testing set
R2_test = model2.score(X_test, y_test)
print("Testing set R-squared", R2_test)
In [34]:
# value of intercepts
for var, coef in zip(variables, model2.coef_):
print(var, coef)
using statsmodel ordinary linear regrssion module
In [35]:
df_train = pd.concat([X_train, y_train], axis = 1) # build a dataframe for training set
reducedModel = smf.ols("lead ~ cadmium + C(lime) + zinc", data = df_train)
reducedModel = reducedModel.fit()
print(reducedModel.summary())
diagnostics plots
In [36]:
pred_val = reducedModel.fittedvalues.copy()
true_val = y_train.copy()
diagnosticsPlots(pred_val, true_val)
In [37]:
table2 = RMSETable(y_train.values, model2.predict(X_train), R2_train, len(y_train),
y_test.values, model2.predict(X_test), R2_test, len(y_test))
print(table2)
model evaluation and comparison¶
- model with 5 predictors (full model)
In [38]:
print(table1)
- model with 3 predictors (reduced model)
In [39]:
print(table2)
In [40]:
table = sm.stats.anova_lm(reducedModel, fullModel)
table
Out[40]:
Comparing the tables obtained above:
- The full model with 6 predictors has lower test RMSE and higher $R^2$ value than the reduced model with 3 predictors
- However, a higher $R^2$ doesn't mean that we favor the full model, since $R^2$ will always increase as the number of predictors used increases. Instead, we will take a look at $R^2_{adj}$, the adjusted $R^2$, which is calculated as
where $n$ is the sample size and $p$ is the number of predictors used in the model.
- Since $R^2_{adj}$ of the full model is 0.945, while the $R^2_{adj}$ of the reduced model is 0.941, we observe that the model doens't obtain very much improvement by adding three additional predictors.
Summary of the module¶
By building two linear regression models using 6 and 3 predictors respectively, we observe that predictor variables cadmium, copper and zinc have very strong positive correlation with the response variable lead. In addition, we also observed that by adding three more predictors, dist, elev and lime, the $R^2_{adj}$ didn't get much improvement comparing to the model with only 3 predictors.
Therefore, the final model that we decided to adopt is the reduced model with 3 predictors:
$$\text{Lead} \sim \text{Cadmium} + \text{lime} + \text{Zinc}$$Reference:
Dataset: meuse dataset
- P.A. Burrough, R.A. McDonnell, 1998. Principles of Geographical Information Systems. Oxford University Press.
- Stichting voor Bodemkartering (STIBOKA), 1970. Bodemkaart van Nederland : Blad 59 Peer, Blad 60 West en 60 Oost Sittard: schaal 1 : 50 000. Wageningen, STIBOKA.
To access this dataset in R:
install.packages("sp")
library(sp)
data(meuse)© Kaixin Wang, updated October 2019