In this section we will be going over LINEAR REGRESSION. We'll be going over how to use the scikit-learn regression model, as well as how to train the regressor using the fit() method, and how to predict new labels using the predict() method. We'll be analyzing a data set consisting of house prices in Boston. We'll start off with a single variable linear regression using numpy and then move on to using scikit learn. We'll do an overview of the mathematics behind the method we're using, but mostly we'll dive deeper into pratical "hands-on" coding lessons.
In this section we will be working through linear regression with the following steps:
Step 1: Getting and setting up the data.
Step 2: Visualizing current data.
Step 3: The mathematics behind the Least Squares Method.
Step 4: Using Numpy for a Univariate Linear Regression.
Step 5: Getting the error.
Step 6: Using scikit learn to implement a multivariate regression.
Step 7: Using Training and Validation.
Step 8: Predicting Prices
Step 9 : Residual PlotsWe'll start by looking a an example of a dataset from scikit-learn. First we'll import our usual data analysis imports, then sklearn's built-in boston dataset.
In [9]:
import numpy as np
import pandas as pd
from pandas import Series,DataFrame
Imports for plotting
In [10]:
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
%matplotlib inline
Now import dataset from scikit learn as well as the linear_model module. Note: You may have to run a download, scikit learn will read an error and prompt you to if you don't have the datasets.
In [11]:
from sklearn.datasets import load_boston
Next we'll download the data set
In [12]:
# Load the housing dataset
boston = load_boston()
Let's see what the data set contains
In [13]:
print(boston.DESCR)
You should always try to do a quick visualization fo the data you have. Let's go ahead an make a histogram of the prices.
In [14]:
# Histogram of prices (this is the target of our dataset)
plt.hist(boston.target,bins=50)
#label
plt.xlabel('Price in $1000s')
plt.ylabel('Number of houses')
Out[14]:
Interesting, now let's see a scatter plot of one feature, versus the target. In this case we'll use the housing price versus the number of rooms in the dwelling.
In [15]:
# Plot the column at the 5 index (Labeled RM)
plt.scatter(boston.data[:,5],boston.target)
#label
plt.ylabel('Price in $1000s')
plt.xlabel('Number of rooms')
Out[15]:
Great! Now we can make out a slight trend that price increases along with the number of rooms in that house, which intuitively makes sense! Now let's use scikit learn to see if we can fit the data linearly.
Let's try to do the following:
1.) Use pandas to transform the boston dataset into a DataFrame:
2.) Then use seaborn to perform an lmplot on that DataFrame to reproduce the scatter plot with a linear fit line.
In [16]:
# reset data as pandas DataFrame
boston_df = DataFrame(boston.data)
# label columns
boston_df.columns = boston.feature_names
#show
boston_df.head()
Out[16]:
Now let's add the target of the boston data set, the price. We'll create a new column in our DataFrame.
In [17]:
# Set price column for target
boston_df['Price'] = boston.target
Now let's see the resultign DataFrame!
In [18]:
# Show result
boston_df.head()
Out[18]:
Now, you might be reminded of the seaborn lmplot function we used during the visualization lectures. You could use it here to do a linear fit automatically!
In [19]:
# Using seabron to create a linear fit
sns.lmplot('RM','Price',data = boston_df)
Out[19]:
However, we won't be able to do this when we move to more complicated regression models, so we'll stay focused on using the scikit learn library!
In this particular lecture we'll use the least squares method as the way to estimate the coefficients. Here's a quick breakdown of how this method works mathematically:
Take a quick look at the plot we created above using seaborn. Now consider each point, and know that they each have a coordinate in the form (X,Y). Now draw an imaginary line between each point and our current "best-fit" line. We'll call the distanace between each point and our current best-fit line, D. To get a quick image of what we're currently trying to visualize, take a look at the picture below:
In [20]:
# Quick display of image form wikipedia
from IPython.display import Image
url = 'http://upload.wikimedia.org/wikipedia/commons/thumb/b/b0/Linear_least_squares_example2.svg/220px-Linear_least_squares_example2.svg.png'
Image(url)
Out[20]:
Now as before, we're labeling each green line as having a distance D, and each red point as having a coordinate of (X,Y). Then we can define our best fit line as the line having the property were: $$ D_{1}^2 + D_{2}^2 + D_{3}^2 + D_{4}^2 + ....+ D_{N}^2$$
So how do we find this line? The least-square line approximating the set of points:
$$ (X,Y)_{1},(X,Y)_{2},(X,Y)_{3},(X,Y)_{4},(X,Y)_{5}, $$has the equation: $$ Y = a_{0} +a_{1}X $$ this is basically just a rewritten form of the standard equation for a line: $$Y=mx+b$$
We can solve for these constants a0 and a1 by simultaneously solving these equations: $$ \Sigma Y = a_{0}N + a_{1}\Sigma X $$ $$ \Sigma XY = a_{0}\Sigma X + a_{1}\Sigma X^2 $$
These are called the normal equations for the least squares line. There are further steps that can be taken in rearranging these equations to solve for y, but we'll let scikit-learn do the rest of the heavy lifting here. If you want further informatino on the mathematics of the above formulas, check out this great video.
For now, we'll use numpy to do a simple single variable linear regression. Afterwards we'll unleash the power of scikit learn to do a full multivariate linear regression.
Numpy has a built in Least Square Method in its linear algebra library. We'll use this first for our Univariate regression and then move on to scikit learn for out Multi variate regression.
We will start by setting up the X and Y arrays for numpy to take in. An important note for the X array: Numpy expects a two-dimensional array, the first dimension is the different example values, and the second dimension is the attribute number. In this case we have our value as the mean number of rooms per house, and this is a single attribute so the second dimension of the array is just 1. So we'll need to create a (506,1) shape array. There are a few ways to do this, but an easy way to do this is by using numpy's built-in vertical stack tool, vstack.
In [21]:
# Set up X as median room values
X = boston_df.RM
# Use v to make X two-dimensional
X = np.vstack(boston_df.RM)
# Set up Y as the target price of the houses.
Y = boston_df.Price
Now that we have our X and Y, let's go ahead and use numpy to create the single variable linear regression.
We know that a line has the equation: $$y=mx+b$$ which we can rewrite using matrices: $$y=Ap$$ where: $$A = \begin{bmatrix}x & 1\end{bmatrix}$$ and $$p= \begin{bmatrix}m \\b\end{bmatrix}$$
This is the same as the first equation if you carry out the linear algebra. So we'll start by creating the A matrix using numpy. We'll do this by creating a matrix in the form [X 1], so we'll call every value in our original X using a list comprehension and then set up an array in the form [X 1]
In [22]:
# Create the X array in the form [X 1]
X = np.array( [ [value,1] for value in X ] )
Great! Now we can get the best fit values!
In [23]:
# Now get out m and b values for our best fit line
m, b = np.linalg.lstsq(X, Y)[0]
Finally let's plot it all together! Note that we use the original format of the boston information. We only did our matrix transformations to utilize the numpy least square method.
In [24]:
# First the original points, Price vs Avg Number of Rooms
plt.plot(boston_df.RM,boston_df.Price,'o')
# Next the best fit line
x= boston_df.RM
plt.plot(x, m*x + b,'r',label='Best Fit Line')
Out[24]:
Great! We've just completed a single variable regression using the least squares method with Python! Let's see if we can find the error in our fitted line. Checking out the documentation here, we see that the resulting array has the total squared error. For each element, it checks the the difference between the line and the true value (our original D value), squares it, and returns the sum of all these. This was the summed D^2 value we discussed earlier.
It's probably easier to understand the root mean squared error, which is similar to the standard deviation. In this case, to find the root mean square error we divide by the number of elements and then take the square root. There is also an issue of bias and an unbiased regression, but we'll delve into those topics later.
For now let's see how we can get the root mean squared error of the line we just fitted.
In [25]:
# Get the resulting array
result = np.linalg.lstsq(X,Y)
# Get the total error
error_total = result[1]
# Get the root mean square error
rmse = np.sqrt(error_total/len(X) )
# Print
print("The root mean squared error was %.2f " %rmse)
Since the root mean square error (RMSE) corresponds approximately to the standard deviation we can now say that the price of a house won't vary more than 2 times the RMSE 95% of the time. Note: Review the Normal Distribution Appendix lecture if this doesn't make sense to you or check out this link.
Thus we can reasonably expect a house price to be within $13,200 of our line fit.
Now, we'll keep moving along with using scikit learn to do a multi variable regression. This will be a similar apporach to the above example, but sci kit learn will be able to take into account more than just a single data variable effecting the target!
We'll start by importing the linear regression library from the sklearn module.
The sklearn.linear_model.LinearRegression class is an estimator. Estimators predict a value based on the observed data. In scikit-learn, all estimators implement the fit() and predict() methods. The former method is used to learn the parameters of a model, and the latter method is used to predict the value of a response variable for an explanatory variable using the learned parameters. It is easy to experiment with different models using scikit-learn because all estimators implement the fit and predict methods.
In [26]:
# Import for Linear Regression
import sklearn
from sklearn.linear_model import LinearRegression
Next, we create a LinearRegression object, afterwards, type lm. then press tab to see the list of methods availble on this object.
In [27]:
# Create a LinearRegression Object
lreg = LinearRegression()
The functions we will be using are:
lreg.fit() which fits a linear model
lreg.predict() which is used to predict Y using the linear model with estimated coefficients
lreg.score() which returns the coefficient of determination (R^2). A measure of how well observed outcomes are replicated by the model, learn more about it here
We'll start the multi variable regression analysis by seperating our boston dataframe into the data columns and the target columns:
In [28]:
# Data Columns
X_multi = boston_df.drop('Price',1)
# Targets
Y_target = boston_df.Price
Finally, we're ready to pass the X and Y using the linear regression object.
In [29]:
# Implement Linear Regression
lreg.fit(X_multi,Y_target)
Out[29]:
Let's go ahead check the intercept and number of coefficients.
In [30]:
print(' The estimated intercept coefficient is %.2f ' %lreg.intercept_)
In [31]:
print(' The number of coefficients used was %d ' % len(lreg.coef_))
Great! So we have basically made an equation for a line, but instead of just oneo coefficient m and an intercept b, we now have 13 coefficients. To get an idea of what this looks like check out the documentation for this equation: $$ y(w,x) = w_0 + w_1 x_1 + ... + w_p x_p $$
Where $$w = (w_1, ...w_p)$$ as the coefficients and $$ w_0 $$ as the intercept
What we'll do next is set up a DataFrame showing all the Features and their estimated coefficients obtained form the linear regression.
In [32]:
# Set a DataFrame from the Features
coeff_df = DataFrame(boston_df.columns)
coeff_df.columns = ['Features']
# Set a new column lining up the coefficients from the linear regression
coeff_df["Coefficient Estimate"] = pd.Series(lreg.coef_)
# Show
coeff_df
Out[32]:
Just like we initially plotted out, it seems the highest correlation between a feature and a house price was the number of rooms.
Now let's move on to Predicting prices!
In a dataset a training set is implemented to build up a model, while a validation set is used to validate the model built. Data points in the training set are excluded from the validation set. The correct way to pick out samples from your dataset to be part either the training or validation (also called test) set is randomly.
Fortunately, scikit learn has a built in function specifically for this called train_test_split.
The parameters passed are your X and Y, then optionally test_size parameter, representing the proportion of the dataset to include in the test split. As well a train_size parameter. ou can learn more about these parameters here
In [33]:
# Grab the output and set as X and Y test and train data sets!
X_train, X_test, Y_train, Y_test = sklearn.cross_validation.train_test_split(X,boston_df.Price)
Let's go ahead and see what the output of the train_test_split was:
In [40]:
# Residual plot of all the dataset using seaborn
sns.residplot('RM', 'Price', data = boston_df)
Out[40]:
In [34]:
# Print shapes of the training and testing data sets
print(X_train.shape, X_test.shape, Y_train.shape, Y_test.shape)
Great! Now that we have our training and testing sets we can continue on to predicint gprices based on the multiple variables.
Now that we have our training and testing sets, let's go ahead and try to use them to predict house prices. We'll use our training set for the prediction and then use our testing set for validation.
In [35]:
# Create our regression object
lreg = LinearRegression()
# Once again do a linear regression, except only on the training sets this time
lreg.fit(X_train,Y_train)
Out[35]:
Now run a prediction on both the X training set and the testing set.
In [36]:
# Predictions on training and testing sets
pred_train = lreg.predict(X_train)
pred_test = lreg.predict(X_test)
Now we will get the mean square error
In [37]:
print("Fit a model X_train, and calculate MSE with Y_train: %.2f" % np.mean((Y_train - pred_train) ** 2))
print("Fit a model X_train, and calculate MSE with X_test and Y_test: %.2f" %np.mean((Y_test - pred_test) ** 2))
It looks like our mean square error between our training and testing was pretty close. But how do we actually visualize this?
In regression analysis, the difference between the observed value of the dependent variable (y) and the predicted value (ŷ) is called the residual (e). Each data point has one residual, so that:
$$Residual = Observed\:value - Predicted\:value $$You can think of these residuals in the same way as the D value we discussed earlier, in this case however, there were multiple data points considered.
A residual plot is a graph that shows the residuals on the vertical axis and the independent variable on the horizontal axis. If the points in a residual plot are randomly dispersed around the horizontal axis, a linear regression model is appropriate for the data; otherwise, a non-linear model is more appropriate.
Residual plots are a good way to visualize the errors in your data. If you have done a good job then your data should be randomly scattered around line zero. If there is some strucutre or pattern, that means your model is not capturing some thing. There could be an interaction between 2 variables that you're not considering, or may be you are measuring time dependent data. If this is the case go back to your model and check your data set closely.
So now let's go ahead and create the residual plot. For more info on the residual plots check out this great link.
In [39]:
# Scatter plot the training data
train = plt.scatter(pred_train,(Y_train-pred_train),c='b',alpha=0.5)
# Scatter plot the testing data
test = plt.scatter(pred_test,(Y_test-pred_test),c='r',alpha=0.5)
# Plot a horizontal axis line at 0
plt.hlines(y=0,xmin=-10,xmax=50)
#Labels
plt.legend((train,test),('Training','Test'),loc='lower left')
plt.title('Residual Plots')
Out[39]:
Great! Looks like there aren't any major patterns to be concerned about, it may be interesting to check out the line occuring towards the bottom right, but overall the majority of the residuals seem to be randomly allocated above and below the horizontal. We could also use seaborn to create these plots: